diff --git a/ecoli/composites/ecoli_configs/fba_redux_classic.json b/ecoli/composites/ecoli_configs/fba_redux_classic.json new file mode 100644 index 000000000..853723d17 --- /dev/null +++ b/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 +} diff --git a/ecoli/experiments/metabolism_redux_sim.py b/ecoli/experiments/metabolism_redux_sim.py index 7ff09e9bd..52c2e3467 100644 --- a/ecoli/experiments/metabolism_redux_sim.py +++ b/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, diff --git a/ecoli/library/sim_data.py b/ecoli/library/sim_data.py index ec33dc0b2..0ec7bb247 100644 --- a/ecoli/library/sim_data.py +++ b/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, diff --git a/ecoli/processes/__init__.py b/ecoli/processes/__init__.py index 2b37fa746..de54ed894 100644 --- a/ecoli/processes/__init__.py +++ b/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) diff --git a/ecoli/processes/metabolism_redux.py b/ecoli/processes/metabolism_redux.py index 36e41b4f8..e6fe681d9 100644 --- a/ecoli/processes/metabolism_redux.py +++ b/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." diff --git a/ecoli/processes/metabolism_redux_classic.py b/ecoli/processes/metabolism_redux_classic.py index 36983f410..0b009a23d 100644 --- a/ecoli/processes/metabolism_redux_classic.py +++ b/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."