Source code for ecoli.processes.metabolism

"""
==========
Metabolism
==========

Encodes molecular simulation of microbial metabolism using flux-balance analysis.

This process demonstrates how metabolites are taken up from the environment
and converted into other metabolites for use in other processes.

NOTE:
- In wcEcoli, metabolism only runs after all other processes have completed
and internal states have been updated (deriver-like, no partitioning necessary)
"""

from typing import Any, Optional
import warnings

import numpy as np
import numpy.typing as npt
from scipy.sparse import csr_matrix
from unum import Unum
from vivarium.core.process import Step
from vivarium.library.units import units as vivunits

from ecoli.processes.registries import topology_registry
from ecoli.library.schema import numpy_schema, bulk_name_to_idx, counts, listener_schema
from wholecell.utils import units
from wholecell.utils.random import stochasticRound
from wholecell.utils.modular_fba import FluxBalanceAnalysis
from reconstruction.ecoli.dataclasses.process.metabolism import REVERSE_TAG


# Register default topology for this process, associating it with process name
NAME = "ecoli-metabolism"
TOPOLOGY = {
    "bulk": ("bulk",),
    # Non-partitioned counts
    "bulk_total": ("bulk",),
    "listeners": ("listeners",),
    "environment": {
        "_path": ("environment",),
        "exchange": ("exchange",),
    },
    "boundary": ("boundary",),
    "polypeptide_elongation": ("process_state", "polypeptide_elongation"),
    "global_time": ("global_time",),
    "timestep": ("timestep",),
    "next_update_time": ("next_update_time", "metabolism"),
}
topology_registry.register(NAME, TOPOLOGY)

COUNTS_UNITS = units.mmol
VOLUME_UNITS = units.L
MASS_UNITS = units.g
TIME_UNITS = units.s
CONC_UNITS = COUNTS_UNITS / VOLUME_UNITS
CONVERSION_UNITS = MASS_UNITS * TIME_UNITS / VOLUME_UNITS
GDCW_BASIS = units.mmol / units.g / units.h

USE_KINETICS = True


[docs] class Metabolism(Step): """Metabolism Process""" name = NAME topology = TOPOLOGY defaults = { "get_import_constraints": lambda u, c, p: (u, c, []), "nutrientToDoublingTime": {}, "use_trna_charging": False, "include_ppgpp": False, "mechanistic_aa_transport": False, "aa_names": [], "aa_targets_not_updated": set(), "import_constraint_threshold": 0, "exchange_molecules": [], "current_timeline": None, "media_id": "minimal", "imports": {}, "metabolism": {}, "ngam": 8.39 * units.mmol / (units.g * units.h), "avogadro": 6.02214076e23 / units.mol, "cell_density": 1100 * units.g / units.L, "dark_atp": 33.565052868380675 * units.mmol / units.g, "cell_dry_mass_fraction": 0.3, "get_biomass_as_concentrations": lambda doubling_time: {}, "ppgpp_id": "ppgpp", "get_ppGpp_conc": lambda media: 0.0, "exchange_data_from_media": lambda media: [], "get_masses": lambda exchanges: [], "doubling_time": 44.0 * units.min, "amino_acid_ids": {}, "linked_metabolites": None, "aa_exchange_names": [], "removed_aa_uptake": [], "seed": 0, # TODO: For testing, remove later (perhaps after modifying sim data) "reduce_murein_objective": False, "base_reaction_ids": [], "fba_reaction_ids_to_base_reaction_ids": [], "time_step": 1, } def __init__(self, parameters=None): super().__init__(parameters) # Use information from the environment and sim self.get_import_constraints = self.parameters["get_import_constraints"] self.nutrientToDoublingTime = self.parameters["nutrientToDoublingTime"] self.use_trna_charging = self.parameters["use_trna_charging"] self.include_ppgpp = self.parameters["include_ppgpp"] self.mechanistic_aa_transport = self.parameters["mechanistic_aa_transport"] self.current_timeline = self.parameters["current_timeline"] self.media_id = self.parameters["media_id"] self.exchange_molecules = self.parameters["exchange_molecules"] self.environment_molecules = sorted( [mol[:-3] for mol in self.exchange_molecules] ) # Create model to use to solve metabolism updates self.model = FluxBalanceAnalysisModel( self.parameters, timeline=self.current_timeline, include_ppgpp=self.include_ppgpp, ) # Save constants self.nAvogadro = self.parameters["avogadro"] self.cellDensity = self.parameters["cell_density"] # Track updated AA concentration targets with tRNA charging self.aa_targets = {} self.aa_targets_not_updated = self.parameters["aa_targets_not_updated"] self.aa_names = self.parameters["aa_names"] # Comparing two values with units is faster than converting units # and comparing magnitudes self.import_constraint_threshold = ( self.parameters["import_constraint_threshold"] * vivunits.mM ) # Molecules with concentration updates for listener self.linked_metabolites = self.parameters["linked_metabolites"] doubling_time = self.nutrientToDoublingTime.get( self.media_id, self.nutrientToDoublingTime["minimal"] ) update_molecules = list( self.model.getBiomassAsConcentrations(doubling_time).keys() ) if self.use_trna_charging: update_molecules += [ aa for aa in self.aa_names if aa not in self.aa_targets_not_updated ] update_molecules += list(self.linked_metabolites.keys()) if self.include_ppgpp: update_molecules += [self.model.ppgpp_id] self.conc_update_molecules = sorted(update_molecules) self.aa_exchange_names = self.parameters["aa_exchange_names"] self.removed_aa_uptake = self.parameters["removed_aa_uptake"] self.aa_environment_names = [aa[:-3] for aa in self.aa_exchange_names] self.seed = self.parameters["seed"] self.random_state = np.random.RandomState(seed=self.seed) # TODO: For testing, remove later (perhaps after modifying sim data) self.reduce_murein_objective = self.parameters["reduce_murein_objective"] # Helper indices for Numpy indexing self.metabolite_idx = None # Get conversion matrix to compile individual fluxes in the FBA # solution to the fluxes of base reactions self.fba_reaction_ids = self.model.fba.getReactionIDs() self.base_reaction_ids = self.parameters["base_reaction_ids"] fba_reaction_ids_to_base_reaction_ids = self.parameters[ "fba_reaction_ids_to_base_reaction_ids" ] self.externalMoleculeIDs = self.model.fba.getExternalMoleculeIDs() self.outputMoleculeIDs = self.model.fba.getOutputMoleculeIDs() self.kineticTargetFluxNames = self.model.fba.getKineticTargetFluxNames() self.homeostaticTargetMolecules = self.model.fba.getHomeostaticTargetMolecules() fba_reaction_id_to_index = { rxn_id: i for (i, rxn_id) in enumerate(self.fba_reaction_ids) } base_reaction_id_to_index = { rxn_id: i for (i, rxn_id) in enumerate(self.base_reaction_ids) } base_rxn_indexes = [] fba_rxn_indexes = [] v = [] for fba_rxn_id in self.fba_reaction_ids: base_rxn_id = fba_reaction_ids_to_base_reaction_ids[fba_rxn_id] base_rxn_indexes.append(base_reaction_id_to_index[base_rxn_id]) fba_rxn_indexes.append(fba_reaction_id_to_index[fba_rxn_id]) if fba_rxn_id.endswith(REVERSE_TAG): v.append(-1) else: v.append(1) base_rxn_indexes = np.array(base_rxn_indexes) fba_rxn_indexes = np.array(fba_rxn_indexes) v = np.array(v) shape = (len(self.base_reaction_ids), len(self.fba_reaction_ids)) self.reaction_mapping_matrix = csr_matrix( (v, (base_rxn_indexes, fba_rxn_indexes)), shape=shape ) def __getstate__(self): return self.parameters def __setstate__(self, state): self.__init__(state)
[docs] def ports_schema(self): ports = { "bulk": numpy_schema("bulk"), "bulk_total": numpy_schema("bulk"), "environment": { "media_id": {"_default": "", "_updater": "set"}, "exchange": { str(element): {"_default": 0} for element in self.environment_molecules }, "exchange_data": { "unconstrained": {"_default": {}}, "constrained": {"_default": set()}, }, }, "boundary": {"external": {"*": {"_default": 0 * vivunits.mM}}}, "listeners": { "mass": listener_schema( { "cell_mass": 0.0, "dry_mass": 0.0, "rna_mass": 0.0, "protein_mass": 0.0, } ), "fba_results": listener_schema( { "media_id": "", "conc_updates": ( [0.0] * len(self.conc_update_molecules), self.conc_update_molecules, ), "catalyst_counts": ( [0] * len(self.model.catalyst_ids), self.model.catalyst_ids, ), "translation_gtp": 0.0, "coefficient": 0.0, "unconstrained_molecules": ( [False] * len(self.exchange_molecules), self.exchange_molecules, ), "constrained_molecules": ( [False] * len(self.exchange_molecules), self.exchange_molecules, ), "uptake_constraints": ( [-1.0] * len(self.exchange_molecules), self.exchange_molecules, ), "delta_metabolites": ( [0] * len(self.model.metaboliteNamesFromNutrients), self.model.metaboliteNamesFromNutrients, ), "reaction_fluxes": ( [0.0] * len(self.fba_reaction_ids), self.fba_reaction_ids, ), "external_exchange_fluxes": ( [0.0] * len(self.externalMoleculeIDs), self.externalMoleculeIDs, ), "objective_value": 0.0, "shadow_prices": ( [0.0] * len(self.outputMoleculeIDs), self.outputMoleculeIDs, ), "reduced_costs": ( [0.0] * len(self.fba_reaction_ids), self.fba_reaction_ids, ), "target_concentrations": ( [0.0] * len(self.homeostaticTargetMolecules), self.homeostaticTargetMolecules, ), "homeostatic_objective_values": ( [0.0] * len(self.homeostaticTargetMolecules), self.homeostaticTargetMolecules, ), "kinetic_objective_values": ( [0.0] * len(self.kineticTargetFluxNames), self.kineticTargetFluxNames, ), "base_reaction_fluxes": ( [0.0] * len(self.base_reaction_ids), self.base_reaction_ids, ), # 'estimated_fluxes': {}, # 'estimated_homeostatic_dmdt': {}, # 'target_homeostatic_dmdt': {}, # 'estimated_exchange_dmdt': {}, # 'target_kinetic_fluxes': {}, # 'target_kinetic_bounds': {}, # 'target_maintenance_flux': 0 } ), "enzyme_kinetics": listener_schema( { "metabolite_counts_init": ( [0] * len(self.model.metaboliteNamesFromNutrients), self.model.metaboliteNamesFromNutrients, ), "metabolite_counts_final": ( [0] * len(self.model.metaboliteNamesFromNutrients), self.model.metaboliteNamesFromNutrients, ), "enzyme_counts_init": ( [0] * len(self.model.kinetic_constraint_enzymes), self.model.kinetic_constraint_enzymes, ), "counts_to_molar": 1.0, "actual_fluxes": ( [0.0] * len(self.model.kinetics_constrained_reactions), self.model.kinetics_constrained_reactions, ), "target_fluxes": ( [0.0] * len(self.model.kinetics_constrained_reactions), self.model.kinetics_constrained_reactions, ), "target_fluxes_upper": ( [0.0] * len(self.model.kinetics_constrained_reactions), self.model.kinetics_constrained_reactions, ), "target_fluxes_lower": ( [0.0] * len(self.model.kinetics_constrained_reactions), self.model.kinetics_constrained_reactions, ), "target_aa_conc": ([0.0] * len(self.aa_names), self.aa_names), } ), }, "polypeptide_elongation": { "aa_count_diff": { "_default": {}, "_emit": True, "_updater": "set", "_divider": "empty_dict", }, "gtp_to_hydrolyze": { "_default": 0.0, "_emit": True, "_updater": "set", "_divider": "zero", }, "aa_exchange_rates": { "_default": 0.0, "_emit": True, "_updater": "set", "_divider": "zero", }, }, "global_time": {"_default": 0.0}, "timestep": {"_default": self.parameters["time_step"]}, "next_update_time": { "_default": self.parameters["time_step"], "_updater": "set", "_divider": "set", }, } return ports
[docs] def update_condition(self, timestep, states): """ See :py:meth:`~ecoli.processes.partition.Requester.update_condition`. """ if states["next_update_time"] <= states["global_time"]: if states["next_update_time"] < states["global_time"]: warnings.warn( f"{self.name} updated at t=" f"{states['global_time']} instead of t=" f"{states['next_update_time']}. Decrease the " "timestep for the global clock process for more " "accurate timekeeping." ) return True return False
[docs] def next_update(self, timestep, states): # At t=0, convert all strings to indices if self.metabolite_idx is None: self.metabolite_idx = bulk_name_to_idx( self.model.metaboliteNamesFromNutrients, states["bulk"]["id"] ) self.catalyst_idx = bulk_name_to_idx( self.model.catalyst_ids, states["bulk"]["id"] ) self.kinetics_enzymes_idx = bulk_name_to_idx( self.model.kinetic_constraint_enzymes, states["bulk"]["id"] ) self.kinetics_substrates_idx = bulk_name_to_idx( self.model.kinetic_constraint_substrates, states["bulk"]["id"] ) self.aa_idx = bulk_name_to_idx(self.aa_names, states["bulk"]["id"]) timestep = states["timestep"] # Load current state of the sim # Get internal state variables metabolite_counts_init = counts(states["bulk"], self.metabolite_idx) catalyst_counts = counts(states["bulk"], self.catalyst_idx) kinetic_enzyme_counts = counts(states["bulk"], self.kinetics_enzymes_idx) kinetic_substrate_counts = counts(states["bulk"], self.kinetics_substrates_idx) translation_gtp = states["polypeptide_elongation"]["gtp_to_hydrolyze"] cell_mass = states["listeners"]["mass"]["cell_mass"] * units.fg dry_mass = states["listeners"]["mass"]["dry_mass"] * units.fg # Calculate state values cellVolume = cell_mass / self.cellDensity counts_to_molar = (1 / (self.nAvogadro * cellVolume)).asUnit(CONC_UNITS) # Coefficient to convert between flux (mol/g DCW/hr) basis and # concentration (M) basis coefficient = dry_mass / cell_mass * self.cellDensity * timestep * units.s # Get exchange constraints unconstrained = set(states["environment"]["exchange_data"]["unconstrained"]) constrained = states["environment"]["exchange_data"]["constrained"] # Determine updates to concentrations depending on the current state current_media_id = states["environment"]["media_id"] doubling_time = self.nutrientToDoublingTime.get( current_media_id, self.nutrientToDoublingTime[self.media_id] ) if self.include_ppgpp: # Sim does not include ppGpp regulation so do not update biomass # by RNA/protein ratio and include ppGpp concentration as a # homeostatic target conc_updates = self.model.getBiomassAsConcentrations(doubling_time) conc_updates[self.model.ppgpp_id] = self.model.getppGppConc( doubling_time ).asUnit(CONC_UNITS) else: # Sim includes ppGpp regulation so update biomass based on # RNA/protein ratio which is controlled by ppGpp and other growth # regulation rp_ratio = ( states["listeners"]["mass"]["rna_mass"] / states["listeners"]["mass"]["protein_mass"] ) conc_updates = self.model.getBiomassAsConcentrations( doubling_time, rp_ratio=rp_ratio ) if self.use_trna_charging: conc_updates.update( self.update_amino_acid_targets( counts_to_molar, states["polypeptide_elongation"]["aa_count_diff"], dict(zip(self.aa_names, counts(states["bulk_total"], self.aa_idx))), ) ) # Converted from units to make reproduction from listener data # accurate to model results (otherwise can have floating point diffs) conc_updates = { met: conc.asNumber(CONC_UNITS) for met, conc in conc_updates.items() } aa_uptake_package = None if self.mechanistic_aa_transport: aa_in_media = np.array( [ states["boundary"]["external"][aa_name] > self.import_constraint_threshold for aa_name in self.aa_environment_names ] ) aa_in_media[self.removed_aa_uptake] = False exchange_rates = ( states["polypeptide_elongation"]["aa_exchange_rates"] * timestep ).asNumber(CONC_UNITS / TIME_UNITS) aa_uptake_package = ( exchange_rates[aa_in_media], self.aa_exchange_names[aa_in_media], True, ) if self.parameters["reduce_murein_objective"]: conc_updates["CPD-12261[p]"] /= 2.27 # Update FBA problem based on current state # Set molecule availability (internal and external) self.model.set_molecule_levels( metabolite_counts_init, counts_to_molar, coefficient, current_media_id, unconstrained, constrained, conc_updates, aa_uptake_package, ) # Set reaction limits for maintenance and catalysts present self.model.set_reaction_bounds( catalyst_counts, counts_to_molar, coefficient, translation_gtp ) # Constrain reactions based on targets targets, upper_targets, lower_targets = self.model.set_reaction_targets( kinetic_enzyme_counts, kinetic_substrate_counts, counts_to_molar, timestep * units.s, ) # Solve FBA problem and update states n_retries = 3 fba = self.model.fba fba.solve(n_retries) # Internal molecule changes delta_metabolites = (1 / counts_to_molar) * ( CONC_UNITS * fba.getOutputMoleculeLevelsChange() ) metabolite_counts_final = np.fmax( stochasticRound( self.random_state, metabolite_counts_init + delta_metabolites.asNumber() ), 0, ).astype(np.int64) delta_metabolites_final = metabolite_counts_final - metabolite_counts_init # Environmental changes exchange_fluxes = CONC_UNITS * fba.getExternalExchangeFluxes() converted_exchange_fluxes = (exchange_fluxes / coefficient).asNumber(GDCW_BASIS) delta_nutrients = ( ((1 / counts_to_molar) * exchange_fluxes).asNumber().astype(int) ) # Write outputs to listeners unconstrained, constrained, uptake_constraints = self.get_import_constraints( unconstrained, constrained, GDCW_BASIS ) # below is used for comparing target and estimate between GD-FBA # and LP-FBA, no effect on model # maintenance_ngam = ((self.ngam * coefficient) / # (counts_to_molar*timestep)).asNumber() # # TODO (Cyrus) Add change in mass when implementing, # # currently counts/mass. # maintenance_gam = (self.gam).asNumber() # maintenance_gam_active = translation_gtp/timestep # maintenance_target = maintenance_ngam + maintenance_gam \ # + maintenance_gam_active # objective_counts = {str(key): int((self.model.homeostatic_objective[ # key] / counts_to_molar).asNumber()) - int(states['bulk']['count'][ # states['bulk']['id'] == key]) # for key in fba.getHomeostaticTargetMolecules()} # denom = counts_to_molar*timestep # kinetic_targets = {str(self.model.kinetics_constrained_reactions[i]): # int((targets[i] / denom).asNumber()) # for i in range(len(targets))} # target_kinetic_bounds = { # str(self.model.kinetics_constrained_reactions[i]): # (int((lower_targets[i] / denom).asNumber()), # int((upper_targets[i] / denom).asNumber())) # for i in range(len(targets))} # fluxes = fba.getReactionFluxes() / timestep # names = fba.getReactionIDs() # flux_dict = {str(names[i]): int((fluxes[i] / denom).asNumber()) # for i in range(len(names))} reaction_fluxes = fba.getReactionFluxes() / timestep update = { "bulk": [(self.metabolite_idx, delta_metabolites_final)], "environment": { "exchange": { str(molecule[:-3]): delta_nutrients[index] for index, molecule in enumerate(self.externalMoleculeIDs) } }, "listeners": { "fba_results": { "media_id": current_media_id, "conc_updates": [ conc_updates.get(m, 0) for m in self.conc_update_molecules ], "catalyst_counts": catalyst_counts, "translation_gtp": translation_gtp, "coefficient": coefficient.asNumber(CONVERSION_UNITS), "unconstrained_molecules": unconstrained, "constrained_molecules": constrained, "uptake_constraints": uptake_constraints, "delta_metabolites": delta_metabolites_final, "reaction_fluxes": reaction_fluxes, "external_exchange_fluxes": converted_exchange_fluxes, "objective_value": fba.getObjectiveValue(), "shadow_prices": fba.getShadowPrices( self.model.metaboliteNamesFromNutrients ), "reduced_costs": fba.getReducedCosts(fba.getReactionIDs()), "target_concentrations": [ self.model.homeostatic_objective[mol] for mol in fba.getHomeostaticTargetMolecules() ], "homeostatic_objective_values": fba.getHomeostaticObjectiveValues(), "kinetic_objective_values": fba.getKineticObjectiveValues(), "base_reaction_fluxes": self.reaction_mapping_matrix.dot( reaction_fluxes ), # Quite large, comment out to reduce emit size # 'estimated_fluxes': flux_dict , # 'estimated_homeostatic_dmdt': { # metabolite: delta_metabolites_final[index] # for index, metabolite in enumerate( # self.model.metaboliteNamesFromNutrients)}, # 'target_homeostatic_dmdt': objective_counts, # 'estimated_exchange_dmdt': { # molecule: delta_nutrients[index] # for index, molecule in enumerate( # fba.getExternalMoleculeIDs())}, # 'target_kinetic_fluxes': kinetic_targets, # 'target_kinetic_bounds': target_kinetic_bounds, # 'target_maintenance_flux': maintenance_target, }, "enzyme_kinetics": { "metabolite_counts_init": metabolite_counts_init, "metabolite_counts_final": metabolite_counts_final, "enzyme_counts_init": kinetic_enzyme_counts, "counts_to_molar": counts_to_molar.asNumber(CONC_UNITS), "actual_fluxes": fba.getReactionFluxes( self.model.kinetics_constrained_reactions ) / timestep, "target_fluxes": targets / timestep, "target_fluxes_upper": upper_targets / timestep, "target_fluxes_lower": lower_targets / timestep, "target_aa_conc": [ self.aa_targets.get(id_, 0.0) for id_ in self.aa_names ], }, }, "next_update_time": states["global_time"] + states["timestep"], } return update
[docs] def update_amino_acid_targets( self, counts_to_molar: Unum, count_diff: dict[str, float], amino_acid_counts: dict[str, float], ) -> dict[str, Unum]: """ Finds new amino acid concentration targets based on difference in supply and number of amino acids used in polypeptide_elongation. Skips updates to molecules defined in self.aa_targets_not_updated: - L-SELENOCYSTEINE: rare AA that led to high variability when updated Args: counts_to_molar: conversion from counts to molar Returns: ``{AA name (str): new target AA conc (float with mol/volume units)}`` """ if len(self.aa_targets): for aa, diff in count_diff.items(): if aa in self.aa_targets_not_updated: continue self.aa_targets[aa] += diff # TODO (Santiago): Improve targets update if self.aa_targets[aa] < 0: print( "Warning: updated amino acid target for " f"{aa} was negative - adjusted to be positive." ) self.aa_targets[aa] = 1.0 # First time step of a simulation so set target to current counts to # prevent concentration jumps between generations else: for aa, counts in amino_acid_counts.items(): if aa in self.aa_targets_not_updated: continue self.aa_targets[aa] = float(counts) conc_updates = { aa: counts * counts_to_molar for aa, counts in self.aa_targets.items() } # Update linked metabolites that will follow an amino acid for met, link in self.linked_metabolites.items(): conc_updates[met] = ( conc_updates.get(link["lead"], 0 * counts_to_molar) * link["ratio"] ) return conc_updates
[docs] class FluxBalanceAnalysisModel(object): """ Metabolism model that solves an FBA problem with modular_fba. """ def __init__( self, parameters: dict[str, Any], timeline: tuple[tuple[int, str], ...], include_ppgpp: bool = True, ): """ Args: parameters: parameters from simulation data timeline: timeline for nutrient changes during simulation (time of change, media ID), by default [(0.0, 'minimal')] include_ppgpp: if True, ppGpp is included as a concentration target """ nutrients = timeline[0][1] # Local sim_data references metabolism = parameters["metabolism"] self.stoichiometry = metabolism.reaction_stoich self.maintenance_reaction = metabolism.maintenance_reaction # Load constants self.ngam = parameters["ngam"] gam = parameters["dark_atp"] * parameters["cell_dry_mass_fraction"] self.exchange_constraints = metabolism.exchange_constraints self._biomass_concentrations = {} # type: dict self.getBiomassAsConcentrations = parameters["get_biomass_as_concentrations"] # Include ppGpp concentration target in objective if not handled # kinetically in other processes self.ppgpp_id = parameters["ppgpp_id"] self.getppGppConc = parameters["get_ppGpp_conc"] # go through all media in the timeline and add to metaboliteNames metaboliteNamesFromNutrients = set() conc_from_nutrients = ( metabolism.concentration_updates.concentrations_based_on_nutrients ) if include_ppgpp: metaboliteNamesFromNutrients.add(self.ppgpp_id) for time, media_id in timeline: exchanges = parameters["exchange_data_from_media"](media_id) metaboliteNamesFromNutrients.update( conc_from_nutrients(imports=exchanges["importExchangeMolecules"]) ) self.metaboliteNamesFromNutrients = list(sorted(metaboliteNamesFromNutrients)) exchange_molecules = sorted(parameters["exchange_molecules"]) molecule_masses = dict( zip( exchange_molecules, parameters["get_masses"](exchange_molecules).asNumber( MASS_UNITS / COUNTS_UNITS ), ) ) # Setup homeostatic objective concentration targets # Determine concentrations based on starting environment conc_dict = conc_from_nutrients( media_id=nutrients, imports=parameters["imports"] ) doubling_time = parameters["doubling_time"] conc_dict.update(self.getBiomassAsConcentrations(doubling_time)) if include_ppgpp: conc_dict[self.ppgpp_id] = self.getppGppConc(doubling_time) self.homeostatic_objective = dict( (key, conc_dict[key].asNumber(CONC_UNITS)) for key in conc_dict ) # TODO: For testing, remove later (perhaps after modifying sim data) if parameters["reduce_murein_objective"]: self.homeostatic_objective["CPD-12261[p]"] /= 2.27 # Include all concentrations that will be present in a sim for constant # length listeners for met in self.metaboliteNamesFromNutrients: if met not in self.homeostatic_objective: self.homeostatic_objective[met] = 0.0 # Data structures to compute reaction bounds based on enzyme # presence/absence self.catalyst_ids = metabolism.catalyst_ids self.reactions_with_catalyst = metabolism.reactions_with_catalyst i = metabolism.catalysis_matrix_I j = metabolism.catalysis_matrix_J v = metabolism.catalysis_matrix_V shape = (i.max() + 1, j.max() + 1) self.catalysis_matrix = csr_matrix((v, (i, j)), shape=shape) # Function to compute reaction targets based on kinetic parameters and # molecule concentrations self.get_kinetic_constraints = metabolism.get_kinetic_constraints # Remove disabled reactions so they don't get included in the FBA # problem setup kinetic_constraint_reactions = metabolism.kinetic_constraint_reactions constraintsToDisable = metabolism.constraints_to_disable self.active_constraints_mask = np.array( [(rxn not in constraintsToDisable) for rxn in kinetic_constraint_reactions] ) self.kinetics_constrained_reactions = list( np.array(kinetic_constraint_reactions)[self.active_constraints_mask] ) self.kinetic_constraint_enzymes = metabolism.kinetic_constraint_enzymes self.kinetic_constraint_substrates = metabolism.kinetic_constraint_substrates # Set solver and kinetic objective weight (lambda) solver = metabolism.solver kinetic_objective_weight = metabolism.kinetic_objective_weight kinetic_objective_weight_in_range = metabolism.kinetic_objective_weight_in_range # Disable kinetics completely if weight is 0 or specified in file above if not USE_KINETICS or kinetic_objective_weight == 0: objective_type = "homeostatic" self.use_kinetics = False kinetic_objective_weight = 0 else: objective_type = "homeostatic_kinetics_mixed" self.use_kinetics = True # Set up FBA solver # reactionRateTargets value is just for initialization, it gets reset # each timestep during evolveState fba_options = { "reactionStoich": metabolism.reaction_stoich, "externalExchangedMolecules": exchange_molecules, "objective": self.homeostatic_objective, "objectiveType": objective_type, "objectiveParameters": { "kineticObjectiveWeight": kinetic_objective_weight, "kinetic_objective_weight_in_range": kinetic_objective_weight_in_range, "reactionRateTargets": { reaction: 1 for reaction in self.kinetics_constrained_reactions }, "oneSidedReactionTargets": [], }, "moleculeMasses": molecule_masses, # The "inconvenient constant"--limit secretion (e.g., of CO2) "secretionPenaltyCoeff": metabolism.secretion_penalty_coeff, "solver": solver, "maintenanceCostGAM": gam.asNumber(COUNTS_UNITS / MASS_UNITS), "maintenanceReaction": metabolism.maintenance_reaction, } self.fba = FluxBalanceAnalysis(**fba_options) self.metabolite_names = { met: i for i, met in enumerate(self.fba.getOutputMoleculeIDs()) } self.aa_names_no_location = [x[:-3] for x in parameters["amino_acid_ids"]]
[docs] def update_external_molecule_levels( self, objective: dict[str, Unum], metabolite_concentrations: Unum, external_molecule_levels: npt.NDArray[np.float64], ) -> npt.NDArray[np.float64]: """ Limit amino acid uptake to what is needed to meet concentration objective to prevent use as carbon source, otherwise could be used as an infinite nutrient source. Args: objective: homeostatic objective for internal molecules (molecule ID: concentration in counts/volume units) metabolite_concentrations: concentration for each molecule in metabolite_names external_molecule_levels: current limits on external molecule availability Returns: Updated limits on external molecule availability TODO(wcEcoli): determine rate of uptake so that some amino acid uptake can be used as a carbon/nitrogen source """ external_exchange_molecule_ids = self.fba.getExternalMoleculeIDs() for aa in self.aa_names_no_location: if aa + "[p]" in external_exchange_molecule_ids: idx = external_exchange_molecule_ids.index(aa + "[p]") elif aa + "[c]" in external_exchange_molecule_ids: idx = external_exchange_molecule_ids.index(aa + "[c]") else: continue conc_diff = objective[aa + "[c]"] - metabolite_concentrations[ self.metabolite_names[aa + "[c]"] ].asNumber(CONC_UNITS) if conc_diff < 0: conc_diff = 0 if external_molecule_levels[idx] > conc_diff: external_molecule_levels[idx] = conc_diff return external_molecule_levels
[docs] def set_molecule_levels( self, metabolite_counts: npt.NDArray[np.int64], counts_to_molar: Unum, coefficient: Unum, current_media_id: str, unconstrained: set[str], constrained: set[str], conc_updates: dict[str, Unum], aa_uptake_package: Optional[ tuple[npt.NDArray[np.float64], npt.NDArray[np.str_], bool] ] = None, ): """ Set internal and external molecule levels available to the FBA solver. Args: metabolite_counts: counts for each metabolite with a concentration target counts_to_molar: conversion from counts to molar (counts/volume) coefficient: coefficient to convert from mmol/g DCW/hr to mM basis (mass*time/volume) current_media_id: ID of current media unconstrained: molecules that have unconstrained import constrained: molecules (keys) and their limited max uptake rates (mol / mass / time) conc_updates: updates to concentrations targets for molecules (mmol/L) aa_uptake_package: (uptake rates, amino acid names, force levels), determines whether to set hard uptake rates """ # Update objective from media exchanges external_molecule_levels, objective = self.exchange_constraints( self.fba.getExternalMoleculeIDs(), coefficient, CONC_UNITS, current_media_id, unconstrained, constrained, conc_updates, ) self.fba.update_homeostatic_targets(objective) self.homeostatic_objective = {**self.homeostatic_objective, **objective} # Internal concentrations metabolite_conc = counts_to_molar * metabolite_counts self.fba.setInternalMoleculeLevels(metabolite_conc.asNumber(CONC_UNITS)) # External concentrations external_molecule_levels = self.update_external_molecule_levels( objective, metabolite_conc, external_molecule_levels ) self.fba.setExternalMoleculeLevels(external_molecule_levels) if aa_uptake_package: levels, molecules, force = aa_uptake_package self.fba.setExternalMoleculeLevels( levels, molecules=molecules, force=force, allow_export=True )
[docs] def set_reaction_bounds( self, catalyst_counts: npt.NDArray[np.int64], counts_to_molar: Unum, coefficient: Unum, gtp_to_hydrolyze: float, ): """ Set reaction bounds for constrained reactions in the FBA object. Args: catalyst_counts: counts of enzyme catalysts counts_to_molar: conversion from counts to molar (counts/volume) coefficient: coefficient to convert from mmol/g DCW/hr to mM basis (mass*time/volume) gtp_to_hydrolyze: number of GTP molecules to hydrolyze to account for consumption in translation """ # Maintenance reactions # Calculate new NGAM flux = (self.ngam * coefficient).asNumber(CONC_UNITS) self.fba.setReactionFluxBounds( self.fba._reactionID_NGAM, lowerBounds=flux, upperBounds=flux, ) # Calculate GTP usage based on how much was needed in polypeptide # elongation in previous step. flux = (counts_to_molar * gtp_to_hydrolyze).asNumber(CONC_UNITS) self.fba.setReactionFluxBounds( self.fba._reactionID_polypeptideElongationEnergy, lowerBounds=flux, upperBounds=flux, ) # Set hard upper bounds constraints based on enzyme presence # (infinite upper bound) or absence (upper bound of zero) reaction_bounds = np.inf * np.ones(len(self.reactions_with_catalyst)) no_rxn_mask = self.catalysis_matrix.dot(catalyst_counts) == 0 reaction_bounds[no_rxn_mask] = 0 self.fba.setReactionFluxBounds( self.reactions_with_catalyst, upperBounds=reaction_bounds, raiseForReversible=False, )
[docs] def set_reaction_targets( self, kinetic_enzyme_counts: npt.NDArray[np.int64], kinetic_substrate_counts: npt.NDArray[np.int64], counts_to_molar: Unum, time_step: Unum, ) -> tuple[ npt.NDArray[np.float64], npt.NDArray[np.float64], npt.NDArray[np.float64] ]: """ Set reaction targets for constrained reactions in the FBA object. Args: kinetic_enzyme_counts: counts of enzymes used in kinetic constraints kinetic_substrate_counts: counts of substrates used in kinetic constraints counts_to_molar: conversion from counts to molar (counts/volume) time_step: current time step (time) Returns: 3-element tuple containing - **mean_targets**: mean target for each constrained reaction - **upper_targets**: upper target limit for each constrained reaction - **lower_targets**: lower target limit for each constrained reaction """ if self.use_kinetics: enzyme_conc = counts_to_molar * kinetic_enzyme_counts substrate_conc = counts_to_molar * kinetic_substrate_counts # Set target fluxes for reactions based on their most relaxed # constraint reaction_targets = self.get_kinetic_constraints(enzyme_conc, substrate_conc) # Calculate reaction flux target for current time step targets = (time_step * reaction_targets).asNumber(CONC_UNITS)[ self.active_constraints_mask, : ] lower_targets = targets[:, 0] mean_targets = targets[:, 1] upper_targets = targets[:, 2] # Set kinetic targets only if kinetics is enabled self.fba.set_scaled_kinetic_objective(time_step.asNumber(units.s)) self.fba.setKineticTarget( self.kinetics_constrained_reactions, mean_targets, lower_targets=lower_targets, upper_targets=upper_targets, ) else: lower_targets = np.zeros(len(self.kinetics_constrained_reactions)) mean_targets = np.zeros(len(self.kinetics_constrained_reactions)) upper_targets = np.zeros(len(self.kinetics_constrained_reactions)) return mean_targets, upper_targets, lower_targets
def test_metabolism_listener(): from ecoli.experiments.ecoli_master_sim import EcoliSim sim = EcoliSim.from_file() sim.total_time = 2 sim.raw_output = False sim.build_ecoli() sim.run() data = sim.query() reaction_fluxes = data["agents"]["0"]["listeners"]["fba_results"]["reaction_fluxes"] assert isinstance(reaction_fluxes[0], list) assert isinstance(reaction_fluxes[1], list) if __name__ == "__main__": test_metabolism_listener()