From 80ec339771f07ea86553342617010b82ec45ede4 Mon Sep 17 00:00:00 2001 From: Arthur Goldberg Date: Mon, 20 Nov 2017 12:37:17 -0500 Subject: [PATCH] find minimal reaction network, using dead end species and reaction gaps clean up Sphinx docstrings, partly --- tests/test_prepare.py | 220 ++++++++++++++++++++++++++++++++---------- wc_lang/prepare.py | 185 ++++++++++++++++++++++------------- 2 files changed, 286 insertions(+), 119 deletions(-) diff --git a/tests/test_prepare.py b/tests/test_prepare.py index 09e4a4c..69579a7 100644 --- a/tests/test_prepare.py +++ b/tests/test_prepare.py @@ -12,7 +12,7 @@ from wc_lang.core import (Model, Submodel, ObjectiveFunction, Reaction, SpeciesType, Species, Compartment, ReactionParticipant, RateLaw, RateLawEquation, RateLawDirection, SubmodelAlgorithm, - Concentration, BiomassComponent, BiomassReaction) + Concentration, BiomassComponent, BiomassReaction, SpeciesTypeType) from wc_lang.io import Reader from wc_lang.prepare import PrepareModel, CheckModel @@ -33,21 +33,22 @@ def setUp(self): BiomassReaction.objects.reset() # read and initialize a model self.model = Reader().run(self.MODEL_FILENAME) + self.dfba_submodel = Submodel.objects.get_one(id='submodel_1') self.prepare_model = PrepareModel(self.model) + self.id_idx = 0 def test_create_dfba_exchange_rxns(self): - dfba_submodel = Submodel.objects.get_one(id='submodel_1') EXTRACELLULAR_COMPARTMENT_ID = config_wc_lang['EXTRACELLULAR_COMPARTMENT_ID'] self.assertEqual( - self.prepare_model.create_dfba_exchange_rxns(dfba_submodel, EXTRACELLULAR_COMPARTMENT_ID), 2) + self.prepare_model.create_dfba_exchange_rxns(self.dfba_submodel, EXTRACELLULAR_COMPARTMENT_ID), 2) # should add these exchange reactions: # -> specie_1[e] # -> specie_2[e] EXCHANGE_RXN_ID_PREFIX = config_wc_lang['EXCHANGE_RXN_ID_PREFIX'] species_found = set() - for rxn in dfba_submodel.reactions: + for rxn in self.dfba_submodel.reactions: if EXCHANGE_RXN_ID_PREFIX in rxn.id: self.assertEqual(-float('inf'), rxn.min_flux) self.assertEqual(float('inf'), rxn.max_flux) @@ -57,33 +58,32 @@ def test_create_dfba_exchange_rxns(self): species_found.add(participant.species) self.assertEqual(species_found, - set(Species.get(['specie_1[e]', 'specie_2[e]'], dfba_submodel.get_species()))) + set(Species.get(['specie_1[e]', 'specie_2[e]'], self.dfba_submodel.get_species()))) def test_confirm_dfba_submodel_obj_func(self): - dfba_submodel = Submodel.objects.get_one(id='submodel_1') confirm_dfba_submodel_obj_func = self.prepare_model.confirm_dfba_submodel_obj_func - self.assertEqual(confirm_dfba_submodel_obj_func(dfba_submodel), None) + self.assertEqual(confirm_dfba_submodel_obj_func(self.dfba_submodel), None) - dfba_submodel.algorithm = None + self.dfba_submodel.algorithm = None with self.assertRaises(ValueError) as context: - confirm_dfba_submodel_obj_func(dfba_submodel) + confirm_dfba_submodel_obj_func(self.dfba_submodel) self.assertIn("not a dfba submodel", str(context.exception)) - dfba_submodel.algorithm = SubmodelAlgorithm.dfba + self.dfba_submodel.algorithm = SubmodelAlgorithm.dfba - dfba_submodel.objective_function = None - self.assertEqual(confirm_dfba_submodel_obj_func(dfba_submodel), None) - # dfba_submodel should be using its biomass reaction as its objective function - self.assertEqual(dfba_submodel.objective_function.expression, dfba_submodel.biomass_reaction.id) - self.assertEqual(dfba_submodel.objective_function.reactions, []) - self.assertEqual(dfba_submodel.objective_function.biomass_reaction_coefficients, [1.0]) + self.dfba_submodel.objective_function = None + self.assertEqual(confirm_dfba_submodel_obj_func(self.dfba_submodel), None) + # self.dfba_submodel should be using its biomass reaction as its objective function + self.assertEqual(self.dfba_submodel.objective_function.expression, + self.dfba_submodel.biomass_reaction.id) + self.assertEqual(self.dfba_submodel.objective_function.reactions, []) + self.assertEqual(self.dfba_submodel.objective_function.biomass_reaction_coefficients, [1.0]) def test_parse_dfba_submodel_obj_func(self): - dfba_submodel = Submodel.objects.get_one(id='submodel_1') parse_dfba_submodel_obj_func = self.prepare_model.parse_dfba_submodel_obj_func - of = dfba_submodel.objective_function + of = self.dfba_submodel.objective_function # list of tests and expected results # (test, [(coef, reaction_id), ... ] [(coeff, biomass_reaction_id), ... ]) @@ -100,7 +100,7 @@ def test_parse_dfba_submodel_obj_func(self): test, reaction_results, biomass_results = test_and_result of.expression = test - (reactions, biomass_reactions) = parse_dfba_submodel_obj_func(dfba_submodel) + (reactions, biomass_reactions) = parse_dfba_submodel_obj_func(self.dfba_submodel) for coeff,reaction in reaction_results: self.assertIn((coeff,reaction), reactions) self.assertEqual(len(reactions), len(reaction_results)) @@ -119,15 +119,14 @@ def test_parse_dfba_submodel_obj_func(self): for error_input,msg in error_inputs: of.expression = error_input with self.assertRaises(ValueError) as context: - parse_dfba_submodel_obj_func(dfba_submodel) + parse_dfba_submodel_obj_func(self.dfba_submodel) self.assertIn(msg, str(context.exception)) def test_assign_linear_objective_fn(self): - dfba_submodel = Submodel.objects.get_one(id='submodel_1') - of = dfba_submodel.objective_function + of = self.dfba_submodel.objective_function of.expression = 'Metabolism_biomass + reaction_1 + reaction_2*2.0' - (reactions, biomass_reactions) = self.prepare_model.parse_dfba_submodel_obj_func(dfba_submodel) - PrepareModel.assign_linear_objective_fn(dfba_submodel, reactions, biomass_reactions) + (reactions, biomass_reactions) = self.prepare_model.parse_dfba_submodel_obj_func(self.dfba_submodel) + PrepareModel.assign_linear_objective_fn(self.dfba_submodel, reactions, biomass_reactions) self.assertEqual(of.biomass_reactions[0].id, 'Metabolism_biomass') self.assertEqual(of.biomass_reaction_coefficients[0], 1.0) coeffs_n_ids = zip(of.reaction_coefficients, [r.id for r in of.reactions]) @@ -135,28 +134,26 @@ def test_assign_linear_objective_fn(self): self.assertIn((2.0, 'reaction_2'), coeffs_n_ids) def test_apply_default_dfba_submodel_flux_bounds(self): - dfba_submodel = Submodel.objects.get_one(id='submodel_1') self.assertEqual( - self.prepare_model.apply_default_dfba_submodel_flux_bounds(dfba_submodel), (1,1)) - test_non_rev = dfba_submodel.reactions.create( + self.prepare_model.apply_default_dfba_submodel_flux_bounds(self.dfba_submodel), (1,1)) + test_non_rev = self.dfba_submodel.reactions.create( id='__test_1', reversible=False ) - test_rev = dfba_submodel.reactions.create( + test_rev = self.dfba_submodel.reactions.create( id='__test_2', reversible=True ) - self.assertEqual(self.prepare_model.apply_default_dfba_submodel_flux_bounds(dfba_submodel), + self.assertEqual(self.prepare_model.apply_default_dfba_submodel_flux_bounds(self.dfba_submodel), (2,2)) - self.prepare_model.apply_default_dfba_submodel_flux_bounds(dfba_submodel) + self.prepare_model.apply_default_dfba_submodel_flux_bounds(self.dfba_submodel) self.assertEqual(test_non_rev.max_flux, test_rev.max_flux) self.assertEqual(-test_rev.min_flux, test_rev.max_flux) - self.assertEqual(self.prepare_model.apply_default_dfba_submodel_flux_bounds(dfba_submodel), + self.assertEqual(self.prepare_model.apply_default_dfba_submodel_flux_bounds(self.dfba_submodel), (0,0)) def test_run(self): - dfba_submodel = Submodel.objects.get_one(id='submodel_1') - of = dfba_submodel.objective_function + of = self.dfba_submodel.objective_function of.expression = 'Metabolism_biomass + reaction_1 + reaction_2*2.0' self.prepare_model.run() self.assertTrue(of.linear) @@ -165,6 +162,130 @@ def test_run(self): self.assertFalse(of.linear) +class TestGapFinding(unittest.TestCase): + + def setUp(self): + # make model + self.model = Model(id='model') + comp = self.model.compartments.create(id='comp') + self.species = [] + self.num_species = 20 + for i in range(1, self.num_species+1): + spec_type = self.model.species_types.create(id='spec_type_{}'.format(i), + type=SpeciesTypeType.metabolite) + self.species.append(Species(species_type=spec_type, compartment=comp)) + self.dfba_submodel = self.model.submodels.create( + id='metabolism', algorithm=SubmodelAlgorithm.dfba) + + self.id_idx = 0 + self.prepare_model = PrepareModel(self.model) + + def next_id(self): + self.id_idx += 1 + return "rxn_{}".format(self.id_idx) + + def make_reaction(self, submodel, reactant, product, reversible=True): + rxn = submodel.reactions.create(id=self.next_id(), reversible=reversible) + rxn.participants.create(species=reactant, coefficient=-1) + rxn.participants.create(species=product, coefficient=1) + + def create_reaction_network(self, submodel, network_type, **kwargs): + # make networks of reactions with 1 reactant and 1 product + # first delete all Reactions + submodel.reactions = [] + if network_type == 'ring': + # kwargs options: size, reversible + species = self.species + if len(species) < kwargs['size']: + self.fail("not enough species, len(species) < kwargs['size']") + for r_idx in range(kwargs['size']): + product_idx = (r_idx+1) % kwargs['size'] + self.make_reaction(submodel, species[r_idx], species[product_idx], kwargs['reversible']) + else: + self.Fail("Unknown network type: {}".format(network_type)) + + def test_get_inactive_reactions(self): + # make ring of 3 irreversible reactions + self.create_reaction_network(self.dfba_submodel, 'ring', **{'size':3, 'reversible':False}) + + # no dead end species -> no inactive reactions + self.assertEqual(self.prepare_model.get_inactive_reactions(self.dfba_submodel, (set(), set())), []) + + # one dead end species -> 2 inactive reactions + first_specie = self.species[0] + dead_end_species = set([first_specie]) + inactive_reactions = self.prepare_model.get_inactive_reactions(self.dfba_submodel, + (set(), dead_end_species)) + self.assertEqual(len(inactive_reactions), 2) + self.assertIn(self.dfba_submodel.reactions[0], inactive_reactions) + self.assertIn(self.dfba_submodel.reactions[-1], inactive_reactions) + + def test_find_dead_end_species(self): + prep_mdl = self.prepare_model + + # make ring of 4 irreversible reactions + self.create_reaction_network(self.dfba_submodel, 'ring', **{'size':4, 'reversible':False}) + + # ring with no inactive reactions -> no dead end species + species_not_consumed, species_not_produced = prep_mdl.find_dead_end_species(self.dfba_submodel, set()) + self.assertFalse(species_not_consumed) + self.assertFalse(species_not_produced) + + # ring with first reaction missing -> + # species_not_consumed = first reaction's reactant + # species_not_produced = first reaction's product + for part in self.dfba_submodel.reactions[0].participants: + if part.coefficient == -1: + reactant = part.species + if part.coefficient == 1: + product = part.species + del self.dfba_submodel.reactions[0] + species_not_consumed, species_not_produced = prep_mdl.find_dead_end_species(self.dfba_submodel, set()) + self.assertEqual(species_not_consumed.pop(), reactant) + self.assertEqual(species_not_produced.pop(), product) + self.assertFalse(species_not_consumed) + self.assertFalse(species_not_produced) + + # make ring of 4 irreversible reactions + self.create_reaction_network(self.dfba_submodel, 'ring', **{'size':4, 'reversible':False}) + # ring with first reaction inactive -> + # species_not_consumed = first reaction's reactant + # species_not_produced = first reaction's product + species_not_consumed, species_not_produced = prep_mdl.find_dead_end_species(self.dfba_submodel, + set([self.dfba_submodel.reactions[0]])) + self.assertEqual(species_not_consumed.pop(), reactant) + self.assertEqual(species_not_produced.pop(), product) + self.assertFalse(species_not_consumed) + self.assertFalse(species_not_produced) + + # make ring of reversible reactions + self.create_reaction_network(self.dfba_submodel, 'ring', **{'size':3, 'reversible':True}) + # ring with first reaction missing -> all species produced and consumed + del self.dfba_submodel.reactions[0] + species_not_consumed, species_not_produced = prep_mdl.find_dead_end_species(self.dfba_submodel, set()) + self.assertFalse(species_not_consumed) + self.assertFalse(species_not_produced) + + def test_identify_dfba_submodel_rxn_gaps(self): + prep_mdl = self.prepare_model + size = 4 + kwargs = {'size':size, 'reversible':False} + # ring of 4 irreversible reactions -> no dead end species or inactive reactions + self.create_reaction_network(self.dfba_submodel, 'ring', **kwargs) + (not_consumed, not_produced), inactive_rxns = prep_mdl.identify_dfba_submodel_rxn_gaps(self.dfba_submodel) + self.assertFalse(not_consumed) + self.assertFalse(not_produced) + self.assertFalse(inactive_rxns) + + # ring of 4 irreversible reactions with one missing -> all species dead end and all reactions inactive + del self.dfba_submodel.reactions[0] + (not_consumed, not_produced), inactive_rxns = prep_mdl.identify_dfba_submodel_rxn_gaps(self.dfba_submodel) + species_in_ring = set(self.species[0:size]) + self.assertEqual(not_consumed, species_in_ring) + self.assertEqual(not_produced, species_in_ring) + self.assertEqual(sorted(inactive_rxns, key=lambda x: x.id), + sorted(self.dfba_submodel.reactions, key=lambda x: x.id)) + class TestCheckModel(unittest.TestCase): MODEL_FILENAME = os.path.join(os.path.dirname(__file__), 'fixtures', 'test_check_model_model.xlsx') @@ -174,20 +295,19 @@ def setUp(self): model.objects.reset() # read a wc model self.model = Reader().run(self.MODEL_FILENAME) + self.dfba_submodel = Submodel.objects.get_one(id='dfba_submodel') self.check_model = CheckModel(self.model) def test_check_dfba_submodel_1(self): - dfba_submodel = Submodel.objects.get_one(id='dfba_submodel') - self.assertEqual(self.check_model.check_dfba_submodel(dfba_submodel), []) + self.assertEqual(self.check_model.check_dfba_submodel(self.dfba_submodel), []) # delete a reaction's min flux reaction1 = Reaction.objects.get_one(id='reaction_1') reaction1.min_flux = np.nan - errors = self.check_model.check_dfba_submodel(dfba_submodel) + errors = self.check_model.check_dfba_submodel(self.dfba_submodel) self.assertIn("Error: no min_flux for reaction 'reaction_name_1' in submodel", errors[0]) def test_check_dfba_submodel_2(self): - dfba_submodel = Submodel.objects.get_one(id='dfba_submodel') # violate reaction.min_flux <= reaction.max_flux reaction1 = Reaction.objects.get_one(id='reaction_1') @@ -197,39 +317,37 @@ def test_check_dfba_submodel_2(self): reaction2 = Reaction.objects.get_one(id='reaction_2') reaction2.min_flux = 1 - errors = self.check_model.check_dfba_submodel(dfba_submodel) + errors = self.check_model.check_dfba_submodel(self.dfba_submodel) self.assertIn("Error: max_flux < min_flux ({} < {}) for reaction '{}' in submodel".format( reaction1.max_flux, reaction1.min_flux, reaction1.name), errors[0]) self.assertIn("Error: 0 < min_flux ({}) for reversible reaction '{}' in submodel".format( reaction2.min_flux, reaction2.name), errors[1]) def test_check_dfba_submodel_3(self): - dfba_submodel = Submodel.objects.get_one(id='dfba_submodel') # remove all BiomassComponents from the BiomassReaction - dfba_submodel.biomass_reaction.biomass_components = [] - errors = self.check_model.check_dfba_submodel(dfba_submodel) - self.assertIn("Error: submodel '{}' uses dfba but lacks a biomass reaction".format(dfba_submodel.name), + self.dfba_submodel.biomass_reaction.biomass_components = [] + errors = self.check_model.check_dfba_submodel(self.dfba_submodel) + self.assertIn("Error: submodel '{}' uses dfba but lacks a biomass reaction".format(self.dfba_submodel.name), errors[0]) # remove the BiomassReaction - dfba_submodel.biomass_reaction = None - errors = self.check_model.check_dfba_submodel(dfba_submodel) - self.assertIn("Error: submodel '{}' uses dfba but lacks a biomass reaction".format(dfba_submodel.name), + self.dfba_submodel.biomass_reaction = None + errors = self.check_model.check_dfba_submodel(self.dfba_submodel) + self.assertIn("Error: submodel '{}' uses dfba but lacks a biomass reaction".format(self.dfba_submodel.name), errors[0]) # remove the objective function - dfba_submodel.objective_function = None - errors = self.check_model.check_dfba_submodel(dfba_submodel) - self.assertIn("Error: submodel '{}' uses dfba but lacks an objective function".format(dfba_submodel.name), + self.dfba_submodel.objective_function = None + errors = self.check_model.check_dfba_submodel(self.dfba_submodel) + self.assertIn("Error: submodel '{}' uses dfba but lacks an objective function".format(self.dfba_submodel.name), errors[0]) def test_check_dfba_submodel_4(self): - dfba_submodel = Submodel.objects.get_one(id='dfba_submodel') # remove a reaction to test that all species used in biomass reactions are defined - del dfba_submodel.reactions[-1] - errors = self.check_model.check_dfba_submodel(dfba_submodel) + del self.dfba_submodel.reactions[-1] + errors = self.check_model.check_dfba_submodel(self.dfba_submodel) self.assertEquals(len(errors), 1) six.assertRegex(self, errors[0], "Error: undefined species '.*' in biomass reaction '.*' used by submodel") diff --git a/wc_lang/prepare.py b/wc_lang/prepare.py index 1461f71..3e60be0 100644 --- a/wc_lang/prepare.py +++ b/wc_lang/prepare.py @@ -1,4 +1,4 @@ -""" Prepare a WC model for further processing, such as export or simulation +""" Prepare a WC model for further processing, such as export or simulation. :Author: Arthur Goldberg :Date: 2017-10-22 @@ -25,6 +25,8 @@ EXTRACELLULAR_COMPARTMENT_ID = config_wc_lang['EXTRACELLULAR_COMPARTMENT_ID'] +# TODO: distinguish between preparing and analyzing a model: analyze includes gap finding + class PrepareModel(object): '''Statically prepare a model @@ -35,9 +37,9 @@ class PrepareModel(object): * Missing concentrations * Create implicit exchange reactions for dFBA submodels - * Identify gaps in dFBA submodel reaction networks * Ensure that dFBA submodels have objective functions * Apply default flux bounds to the reactions in dFBA submodels + * Identify dead end species and reaction network gaps in dFBA submodels ''' def __init__(self, model): self.model = model @@ -51,14 +53,6 @@ def run(self): EXTRACELLULAR_COMPARTMENT_ID) warn("{} exchange reactions created for submodel '{}'.".format(reactions_created, submodel.name)) - ''' - reaction_gaps = self.identify_dfba_submodel_reaction_gaps(submodel) - if reaction_gaps: - warn("{} gaps the reaction network for submodel '{}'.".format(len(reaction_gaps), - submodel.name)) - for gap in reaction_gaps: - warn("{}".format(gap)) - ''' self.confirm_dfba_submodel_obj_func(submodel) (min_bounds_set, max_bounds_set) = self.apply_default_dfba_submodel_flux_bounds(submodel) warn("{} minimum and {} maximum default flux bounds set for submodel '{}'.".format( @@ -77,17 +71,17 @@ def run(self): def create_dfba_exchange_rxns(self, submodel, extracellular_compartment_id): '''Create exchange reactions for a dFBA submodel's reaction network. - * To represent FBA's mathematical assumption that it models a closed system, create - 'implicit' forward exchange reactions that synthesize all extracellular metabolites. + To represent FBA's mathematical assumption that it models a closed system, create + 'implicit' forward exchange reactions that synthesize all extracellular metabolites. - # TODO: - * To model how other pathways consume metabolites generated by metabolism, create 'implicit' + # TODO: To model how other pathways consume metabolites generated by metabolism, create 'implicit' reactions which exchange these metabolites between a dFBA metabolism submodel and the other pathway(s)/submodel(s). Algorithm to synthesize extracellular metabolites: - E = the set of all extracellular metabolites used by the submodel - generate a "-> e" reaction for each e in E in the submodel + + | E = the set of all extracellular metabolites used by the submodel + | generate a "-> e" reaction for each e in E in the submodel Args: submodel (`Submodel`): a DFBA submodel @@ -97,7 +91,7 @@ def create_dfba_exchange_rxns(self, submodel, extracellular_compartment_id): ValueError: if `submodel` is not a dFBA submodel Returns: - (:obj:`int`): the number of reactions created + :obj:`int`: the number of reactions created ''' if submodel.algorithm != SubmodelAlgorithm.dfba: raise ValueError("submodel '{}' not a dfba submodel".format(submodel.name)) @@ -121,17 +115,25 @@ def create_dfba_exchange_rxns(self, submodel, extracellular_compartment_id): return reaction_number-1 - def identify_dfba_submodel_reaction_gaps(self, submodel): - '''Create reactions to fill gaps in a dFBA submodel's reaction network. + def identify_dfba_submodel_rxn_gaps(self, submodel): + '''Identify gaps in a dFBA submodel's reaction network + + Species that are not consumed or not produced indicate gaps in the reaction network. + These can be found by a static analysis of the model. Reactions that use species that + are not produced or produce species that are not consumed must eventually have zero flux. + A reaction network can be reduced to a minimal network of reactions that can all + have positive fluxes. - Considering only reactions used by this submodel, generate gap filling reactions - that consume species that are not consumed and produce species that are not produced. These - reactions enable an FBA solver to obtain a steady state solution. + Algorithm:: - Algorithm: - S = the set of all species used by submodel - generate a "-> s" reaction for each s in S that is not produced by a reaction in submodel - generate an "s ->" reaction for each s in S that is not consumed by a reaction in submodel + all_gap_species = get_gap_species([]) + delta_gap_species = all_gap_species + while delta_gap_species: + all_gap_reactions = get_gap_reactions(all_gap_species) + tmp_gap_species = all_gap_species + all_gap_species = get_gap_species(all_gap_reactions) + delta_gap_species = all_gap_species - tmp_gap_species + return (all_gap_species, all_gap_reactions) Args: submodel (`Submodel`): a DFBA submodel @@ -140,52 +142,92 @@ def identify_dfba_submodel_reaction_gaps(self, submodel): ValueError: if `submodel` is not a dFBA submodel Returns: - (:obj:`int`): the number of reactions created + :obj:`tuple`: + + * :obj:`set` of :obj:`Species`: `Species` not in the minimal reaction network + * :obj:`set` of :obj:`Reaction`: `Reaction`s not in the minimal reaction network ''' if submodel.algorithm != SubmodelAlgorithm.dfba: raise ValueError("submodel '{}' not a dfba submodel".format(submodel.name)) - reaction_number = 1 + all_dead_end_species = PrepareModel.find_dead_end_species(submodel, set()) + delta_dead_end_species = all_dead_end_species + inactive_reactions = set() + while any(delta_dead_end_species): + inactive_reactions = PrepareModel.get_inactive_reactions(submodel, all_dead_end_species) + tmp_not_consumed, tmp_not_produced = all_dead_end_species + all_dead_end_species = PrepareModel.find_dead_end_species(submodel, inactive_reactions) + all_not_consumed, all_not_produced = all_dead_end_species + delta_dead_end_species = (all_not_consumed-tmp_not_consumed, all_not_produced-tmp_not_produced) + return(all_dead_end_species, inactive_reactions) + @staticmethod + def find_dead_end_species(submodel, inactive_reactions): + '''Find the dead end species in a reaction network + + Given a set of inactive reactions in submodel, determine species that are not consumed by + any reaction, or are not produced by any reaction. Costs :math:`O(n*p)`, where :math:`n` is + the number of reactions in `submodel` and :math:`p` is the maximum number of participants in + a reaction. + + Args: + submodel (`Submodel`): a DFBA submodel + inactive_reactions (`set` of `Reaction`): the inactive reactions in `submodel` + + Returns: + :obj:`tuple`: + + * :obj:`set` of :obj:`Species`: the species that are not consumed + * :obj:`set` of :obj:`Species`: the species that are not produced + ''' species = submodel.get_species() - species_not_produced = set(species) species_not_consumed = set(species) + species_not_produced = set(species) for rxn in submodel.reactions: + if rxn in inactive_reactions: + continue if rxn.reversible: for part in rxn.participants: - species_not_produced.discard(part.species) species_not_consumed.discard(part.species) + species_not_produced.discard(part.species) else: for part in rxn.participants: if part.coefficient<0: species_not_consumed.discard(part.species) elif 0 specie" reaction - new_rxn = submodel.reactions.create( - id = "{}_{}".format(EXCHANGE_RXN_ID_PREFIX, reaction_number), - name = "{}_{}".format(EXCHANGE_RXN_NAME_PREFIX, reaction_number)) - reaction_number += 1 - new_rxn.participants.create(species=specie, coefficient = 1) - - for specie in species_not_consumed: - # generate a "specie ->" reaction - new_rxn = submodel.reactions.create( - id = "{}_{}".format(EXCHANGE_RXN_ID_PREFIX, reaction_number), - name = "{}_{}".format(EXCHANGE_RXN_NAME_PREFIX, reaction_number)) - reaction_number += 1 - new_rxn.participants.create(species=specie, coefficient = -1) + @staticmethod + def get_inactive_reactions(submodel, dead_end_species): + '''Find the inactive reactions in a reaction network - return reaction_number-1 + Given the dead end species in a reaction network, find the reactions that must eventually + become inactive. Reactions that consume species which are not produced must become inactive. + And reactions that produce species which are not consumed must become inactive to prevent + the copy numbers of those species from growing without bound. + Costs :math:`O(n*p)`, where :math:`n` is the number of reactions in `submodel` and :math:`p` + is the maximum number of participants in a reaction. + + Args: + submodel (:obj:`Submodel`): a DFBA submodel + dead_end_species (:obj:`tuple`): + + * :obj:`set` of :obj:`Species`: the `Species` that are not consumed by any `Reaction` in `submodel` + * :obj:`set` of :obj:`Species`: the `Species` that are not produced by any `Reaction` in `submodel` + + Returns: + :obj:`set` of :obj:`Reaction`: the inactive reactions in `submodel`'s reaction network + ''' + species_not_consumed, species_not_produced = dead_end_species + inactive_reactions = [] + for rxn in submodel.reactions: + for part in rxn.participants: + if (part.species in species_not_consumed or + part.species in species_not_produced): + inactive_reactions.append(rxn) + break + return inactive_reactions def confirm_dfba_submodel_obj_func(self, submodel): '''Ensure that a dFBA submodel has an objective function @@ -234,9 +276,10 @@ def parse_dfba_submodel_obj_func(self, submodel): submodel (`Submodel`): a dFBA submodel Returns: - :obj:`(`list`, `list`)`: a pair of lists representing the objective's linear form; - (`list` of (coeff, id) pairs for reactions, `list` of (coeff, id) pairs for biomass - reactions) + :obj:`tuple`: a pair of lists representing the objective's linear form + + * obj:`list`: (coeff, id) pairs for reactions + * obj:`list`: (coeff, id) pairs for biomass reactions Raises: ValueError: if `submodel` is not a dFBA submodel @@ -402,13 +445,16 @@ def apply_default_dfba_submodel_flux_bounds(self, submodel): * max_flux = default_max_flux_bound Args: - submodel (`Submodel`): a dFBA submodel + submodel (:obj:`Submodel`): a dFBA submodel Raises: ValueError: if `submodel` is not a dFBA submodel Returns: - :obj:`tuple` of (`int`, `int`): counts of min and max flux bounds set to the default + :obj:`tuple`: + + * obj:`int`: number of min flux bounds set to the default + * obj:`int`: number of max flux bounds set to the default ''' if submodel.algorithm != SubmodelAlgorithm.dfba: raise ValueError("submodel '{}' not a dfba submodel".format(submodel.name)) @@ -465,7 +511,7 @@ class CheckModel(object): * Rate laws transcode and evaluate without error * All reactants in each submodel's reactions are in the submodel's compartment - Other properties to check: + Other properties to be checked: * The model does not contain dead-end species which are only consumed or produced * Reactions are balanced @@ -476,6 +522,7 @@ class CheckModel(object): * Ensure that Reaction and BiomassReaction ids don't overlap; can then simplify ObjectiveFunction.deserialize() # TODO: implement these, and expand the list of properties + # TODO: fix doc string formatting ''' def __init__(self, model): @@ -496,7 +543,8 @@ def check_dfba_submodel(self, submodel): '''Check the inputs to a DFBA submodel Ensure that: - * All regular DFBA reactions have min flux and max flux with appropriate values + + * All DFBA reactions have min flux and max flux with appropriate values * The DFBA submodel contains a biomass reaction and an objective function * All species used in biomass reactions are defined @@ -504,8 +552,8 @@ def check_dfba_submodel(self, submodel): submodel (`Submodel`): a DFBA submodel Returns: - :obj:`list` of `str`: if no errors, returns an empty `list`, otherwise a `list` of - error messages + :obj:`list` of :obj:`str`: if no errors, returns an empty `list`; otherwise a `list` of + error messages ''' errors = [] for reaction in submodel.reactions: @@ -545,14 +593,15 @@ def check_dynamic_submodel(self, submodel): '''Check the inputs to a dynamic submodel Ensure that: + * All reactions have rate laws for the appropriate directions Args: submodel (`Submodel`): a dynamic (SSA or ODE) submodel Returns: - :obj:`list` of :obj:`str` if no errors, returns an empty `list`, otherwise a `list` of - error messages + :obj:`list` of :obj:`str`: if no errors, returns an empty `list`; otherwise a `list` of + error messages ''' errors = [] for reaction in submodel.reactions: @@ -579,8 +628,8 @@ def transcode_and_check_rate_law_equations(self): Ensure that all rate law equations can be transcoded and evaluated. Returns: - :obj:`list` of `str`: if no errors, returns an empty `list`, otherwise a `list` of - error messages + :obj:`list` of `str`: if no errors, returns an empty `list`; otherwise a `list` of + error messages ''' errors = [] species = self.model.get_species() @@ -605,8 +654,8 @@ def verify_reactant_compartments(self): '''Verify that all reactants in each submodel's reactions are in the submodel's compartment Returns: - :obj:`list` of `str`: if no errors, returns an empty `list`, otherwise a `list` of - error messages + :obj:`list` of `str`: if no errors, returns an empty `list`; otherwise a `list` of + error messages ''' errors = [] for lang_submodel in self.model.get_submodels():