Added sim that runs without growth rate control with the new metaboli…

…sm, since it will conflict with future changes.

ecoli/composites/ecoli_configs/fba_redux_classic.json

@@ -0,0 +1,34 @@
+{
+    "experiment_id" : "fba_redux",
+    "total_time" : 10,
+    "swap_processes" : {
+        "ecoli-metabolism" : "ecoli-metabolism-redux-classic"
+    },
+    "exclude_processes": ["exchange_data"],
+    "flow": {
+        "ecoli-metabolism-redux": [["ecoli-chromosome-structure"]],
+        "ecoli-mass-listener": [["ecoli-metabolism-redux"]],
+        "RNA_counts_listener": [["ecoli-metabolism-redux"]],
+        "rna_synth_prob_listener": [["ecoli-metabolism-redux"]],
+        "monomer_counts_listener": [["ecoli-metabolism-redux"]],
+        "dna_supercoiling_listener": [["ecoli-metabolism-redux"]],
+        "replication_data_listener": [["ecoli-metabolism-redux"]],
+        "rnap_data_listener": [["ecoli-metabolism-redux"]],
+        "unique_molecule_counts": [["ecoli-metabolism-redux"]],
+        "ribosome_data_listener": [["ecoli-metabolism-redux"]]
+    },
+    "raw_output": false,
+    "operons": true,
+    "trna_charging": false,
+    "ppgpp_regulation": false,
+    "mass_distribution": true,
+    "trna_attenuation": false,
+    "variable_elongation_transcription": true,
+    "variable_elongation_translation": false,
+    "mechanistic_translation_supply": false,
+    "mechanistic_aa_transport": false,
+    "translation_supply": false,
+    "aa_supply_in_charging": false,
+    "adjust_timestep_for_charging": false,
+    "disable_ppgpp_elongation_inhibition": true
+}

ecoli/experiments/metabolism_redux_sim.py

@@ -156,6 +156,39 @@ def run_ecoli_with_metabolism_redux(
     sim.initial_state = get_state_from_file(path=f'data/{initial_state_file}.json')
     sim.raw_output = raw_output
     sim.save = save
+
+
+    sim.build_ecoli()
+    sim.run()
+
+    query = []
+    folder = f'out/fbagd/{name}_{total_time}_{datetime.date.today()}/'
+    # save_sim_output(folder, query, sim, save_model=True)
+
+
+def run_ecoli_with_metabolism_redux_classic(
+        filename='fba_redux_classic',
+        total_time=4,
+        divide=True,
+        initial_state_file='wcecoli_t0',
+        progress_bar=True,
+        log_updates=False,
+        emitter='timeseries',
+        name='fba-redux',
+        raw_output=False,
+        save=False,
+        # save_times=4,
+):
+    # filename = 'default'
+    sim = EcoliSim.from_file(CONFIG_DIR_PATH + filename + '.json')
+    sim.total_time = total_time
+    sim.divide = divide
+    sim.progress_bar = progress_bar
+    sim.log_updates = log_updates
+    sim.emitter = emitter
+    sim.initial_state = get_state_from_file(path=f'data/{initial_state_file}.json')
+    sim.raw_output = raw_output
+    sim.save = save
     # sim.save_times = [4]
 
     # # simplify working with uptake
@@ -172,8 +205,7 @@ def run_ecoli_with_metabolism_redux(
 
     query = []
     folder = f'out/fbagd/{name}_{total_time}_{datetime.date.today()}/'
-    save_sim_output(folder, query, sim, save_model=True)
-
+    # save_sim_output(folder, query, sim, save_model=True)
 
 
 @pytest.mark.slow
@@ -268,6 +300,7 @@ def run_ecoli_with_default_metabolism(
     '0': run_metabolism,
     '1': run_metabolism_composite,
     '2': run_ecoli_with_metabolism_redux,
+    '2a': run_ecoli_with_metabolism_redux_classic,
     '3': test_ecoli_with_metabolism_redux,
     '4': test_ecoli_with_metabolism_redux_div,
     '5': run_ecoli_with_default_metabolism,

ecoli/library/sim_data.py

@@ -384,6 +384,7 @@ def get_config_by_name(self, name, time_step=1, parallel=False):
             'ecoli-protein-degradation': self.get_protein_degradation_config,
             'ecoli-metabolism': self.get_metabolism_config,
             'ecoli-metabolism-redux': self.get_metabolism_redux_config,
+            'ecoli-metabolism-redux-classic': self.get_metabolism_redux_config,
             'ecoli-chromosome-replication': self.get_chromosome_replication_config,
             'ecoli-mass': self.get_mass_config,
             'ecoli-mass-listener': self.get_mass_listener_config,

ecoli/processes/__init__.py

@@ -15,6 +15,7 @@
 from ecoli.processes.protein_degradation import ProteinDegradation
 from ecoli.processes.metabolism import Metabolism
 from ecoli.processes.metabolism_redux import MetabolismRedux
+from ecoli.processes.metabolism_redux_classic import MetabolismReduxClassic
 from ecoli.processes.chromosome_replication import ChromosomeReplication
 from ecoli.processes.stubs.exchange_stub import Exchange
 from ecoli.processes.listeners.mass_listener import MassListener
@@ -69,6 +70,7 @@
 process_registry.register(ProteinDegradation.name, ProteinDegradation)
 process_registry.register(Metabolism.name, Metabolism)
 process_registry.register(MetabolismRedux.name, MetabolismRedux)
+process_registry.register(MetabolismReduxClassic.name, MetabolismReduxClassic)
 process_registry.register(ChromosomeReplication.name, ChromosomeReplication)
 process_registry.register(MassListener.name, MassListener)
 process_registry.register(DnaSupercoiling.name, DnaSupercoiling)

ecoli/processes/metabolism_redux.py

@@ -236,7 +236,8 @@ def __init__(self, parameters):
         self.network_flow_model = NetworkFlowModel(
             self.stoichiometry, self.metabolite_names,
             self.reaction_names, self.homeostatic_metabolites,
-            self.kinetic_constraint_reactions)
+            self.kinetic_constraint_reactions, self.parameters['get_mass'],
+            self.gam.asNumber(), self.active_constraints_mask)
 
         # important bulk molecule names
         self.catalyst_ids = self.parameters['catalyst_ids']
@@ -373,18 +374,8 @@ def ports_schema(self):
                     'target_kinetic_fluxes': [],
                     'target_kinetic_bounds': [],
                     'reaction_catalyst_counts': [],
-                    'maintenance_target': []
-                }),
-
-                'enzyme_kinetics': listener_schema({
-                    'metabolite_counts_init': 0,
-                    'metabolite_counts_final': 0,
-                    'enzyme_counts_init': 0,
-                    'counts_to_molar': 1.0,
-                    'actual_fluxes': [],
-                    'target_fluxes': [],
-                    'target_fluxes_upper': [],
-                    'target_fluxes_lower': []})
+                    'maintenance_target': 0
+                })
             },
 
             'global_time': {'_default': 0},
@@ -409,10 +400,14 @@ def next_update(self, timestep, states):
         # Initialize indices
         if self.homeostatic_metabolite_idx is None:
             bulk_ids = states['bulk']['id']
-            self.homeostatic_metabolite_idx = bulk_name_to_idx(self.homeostatic_metabolites, bulk_ids)
+            self.homeostatic_metabolite_idx = bulk_name_to_idx(
+                self.homeostatic_metabolites, bulk_ids)
             self.catalyst_idx = bulk_name_to_idx(self.catalyst_ids, bulk_ids)
-            self.kinetics_enzymes_idx = bulk_name_to_idx(self.kinetic_constraint_enzymes, bulk_ids)
-            self.kinetics_substrates_idx = bulk_name_to_idx(self.kinetic_constraint_substrates, bulk_ids)
+            self.kinetics_enzymes_idx = bulk_name_to_idx(
+                self.kinetic_constraint_enzymes, bulk_ids)
+            self.kinetics_substrates_idx = bulk_name_to_idx(
+                self.kinetic_constraint_substrates, bulk_ids)
+            self.aa_idx = bulk_name_to_idx(self.aa_names, bulk_ids)
 
         unconstrained = states['environment']['exchange_data']['unconstrained']
         constrained = states['environment']['exchange_data']['constrained']
@@ -463,16 +458,13 @@ def next_update(self, timestep, states):
         ngam_target = total_ngam.asNumber()
 
         # binary kinetic targets - sum up enzyme counts for each reaction
-        reaction_catalyst_counts = np.array([sum([current_catalyst_counts[enzyme_idx]
-                                                  for enzyme_idx in enzymes_idx])
-                                             if len(enzymes_idx) > 0 else -1
-                                             for enzymes_idx in self.catalyzed_rxn_enzymes_idx])
+        reaction_catalyst_counts = self.catalysis_matrix.dot(current_catalyst_counts)
         # Get reaction indices whose fluxes should be set to zero
         # because there are no enzymes to catalyze the rxn
         binary_kinetic_idx = np.where(reaction_catalyst_counts == 0)[0]
-        
+
         # TODO: Figure out how to handle changing media ID
-        
+
         ## Determine updates to concentrations depending on the current state
         doubling_time = self.nutrient_to_doubling_time.get(
             states['environment']['media_id'],
@@ -533,7 +525,7 @@ def next_update(self, timestep, states):
             aa_uptake_package = (exchange_rates[aa_in_media],
                                  self.aa_exchange_names[aa_in_media], True)
 
-        # kinetic constraints 
+        # kinetic constraints
         # TODO (Cyrus) eventually collect isozymes in single reactions, map
         # enzymes to reacts via stoich instead of kinetic_constraint_reactions
         kinetic_enzyme_conc = self.counts_to_molar * kinetic_enzyme_counts
@@ -555,7 +547,7 @@ def next_update(self, timestep, states):
             kinetic_targets=enzyme_kinetic_boundaries,
             binary_kinetic_idx=binary_kinetic_idx,
             objective_weights=objective_weights,
-            solver=cp.GLOP)
+            aa_uptake_package=aa_uptake_package)
 
         self.reaction_fluxes = solution.velocities
         self.metabolite_dmdt = solution.dm_dt
@@ -704,7 +696,10 @@ def __init__(self,
         metabolites: Iterable[list],
         reactions: Iterable[list],
         homeostatic_metabolites: Iterable[str],
-        kinetic_reactions: Iterable[str]
+        kinetic_reactions: Iterable[str],
+        get_mass: Callable[[str], Unum],
+        gam: float = 0,
+        active_constraints_mask: np.ndarray = None,
     ):
         self.S_orig = csr_matrix(stoich_arr.astype(np.int64))
         self.S_exch = None
@@ -735,7 +730,7 @@ def set_up_exchanges(self,
         the system. Uptakes allow certain metabolites to also have flow into the system."""
         all_exchanges = exchanges.copy()
         all_exchanges.update(uptakes)
-        
+
         # All exchanges can secrete but only uptakes go in both directions
         self.S_exch = np.zeros((self.n_mets, len(exchanges) + len(uptakes)))
         self.exchanges = []
@@ -879,9 +874,11 @@ def test_network_flow_model():
 
     model.set_up_exchanges(exchanges=exchanges, uptakes=uptakes)
 
-    solution: FlowResult = model.solve(homeostatic_targets=list(homeostatic_metabolites.values()),
-                                       objective_weights={'secretion': 0.01, 'efficiency': 0.0001},
-                                       upper_flux_bound=100, solver=cp.GLOP)
+    solution: FlowResult = model.solve(
+        homeostatic_concs=list(homeostatic_metabolites.values()),
+        homeostatic_dm_targets=list(homeostatic_metabolites.values()),
+        objective_weights={'secretion': 0.01, 'efficiency': 0.0001},
+        upper_flux_bound=100)
 
     assert np.isclose(solution.velocities, np.array([1, 1, 0])).all() == True, "Network flow toy model did not converge to correct solution."
 

ecoli/processes/metabolism_redux_classic.py

@@ -7,6 +7,7 @@
 from scipy.sparse import csr_matrix
 
 from vivarium.core.process import Step
+from vivarium.library.units import units as vivunits
 
 from ecoli.library.schema import (numpy_schema, bulk_name_to_idx,
     listener_schema, counts)
@@ -27,7 +28,7 @@
 GDCW_BASIS = units.mmol / units.g / units.h
 
 
-NAME = 'ecoli-metabolism-redux'
+NAME = 'ecoli-metabolism-redux-classic'
 TOPOLOGY = topology_registry.access('ecoli-metabolism')
 # TODO (Cyrus) - Re-add when kinetics are added.
 # TOPOLOGY['kinetic_flux_targets'] = ('rates', 'fluxes')
@@ -40,7 +41,7 @@
 
 FREE_RXNS = ["TRANS-RXN-145", "TRANS-RXN0-545", "TRANS-RXN0-474"]
 
-class MetabolismRedux(Step):
+class MetabolismReduxClassic(Step):
     name = NAME
     topology = TOPOLOGY
 
@@ -216,6 +217,16 @@ def ports_schema(self):
                     'target_fluxes_lower': []})
             },
 
+            # these three were added in parca update, may be able to remove
+            'boundary': {
+                'external': {
+                    '*': {'_default': 0 * vivunits.mM}
+                }
+            },
+
+            'global_time': {'_default': 0},
+            'timestep': {'_default': self.parameters['time_step']},
+
             'first_update': {
                 '_default': True,
                 '_updater': 'set',
@@ -324,7 +335,7 @@ def next_update(self, timestep, states):
             kinetic_targets=target_kinetic_values,
             binary_kinetic_idx=binary_kinetic_idx,
             objective_weights=objective_weights,
-            solver=cp.MOSEK)
+            solver=cp.GLOP)
 
         self.reaction_fluxes = solution.velocities
         self.metabolite_dmdt = solution.dm_dt
@@ -541,7 +552,7 @@ def test_network_flow_model():
 
     solution: FlowResult = model.solve(homeostatic_targets=list(homeostatic_metabolites.values()),
                                        objective_weights={'secretion': 0.01, 'efficiency': 0.0001},
-                                       upper_flux_bound=100, solver=cp.MOSEK)
+                                       upper_flux_bound=100, solver=cp.GLOP)
 
     assert np.isclose(solution.velocities, np.array([1, 1, 0])).all() == True, "Network flow toy model did not converge to correct solution."