fix: predicted pathways containing exchange reactions for native compounds (#234)

A change of `model.exchanges` in upstream cobrapy has led to the
incorporation of DM_ reactions for native compounds (H+, H2O, CO2, etc.)
into the `PathwayPredictor.model` (host model + universal model).

I have replaced `model.exchanges` with `model.boundary` (which has
the same behavior as `model.exchanges` had in the past).

Author: phantomas1234

cameo/core/utils.py

@@ -48,7 +48,7 @@ def medium(model):
     reaction_names = []
     lower_bounds = []
     upper_bounds = []
-    for ex in model.exchanges:
+    for ex in model.boundary:
         metabolite = list(ex.metabolites.keys())[0]
         coeff = ex.metabolites[metabolite]
         if coeff * ex.lower_bound > 0:
@@ -107,7 +107,7 @@ def load_medium(model, medium_def, copy=False, delimiter="\t"):
 
 def _load_medium_from_dict(model, medium_def):
     assert isinstance(medium_def, dict)
-    for ex_reaction in model.exchanges:
+    for ex_reaction in model.boundary:
         ex_reaction.lower_bound = medium_def.get(ex_reaction.id, 0)
 
 
@@ -126,7 +126,7 @@ def _load_medium_from_file(model, file_path, delimiter="\t"):
 
 def _load_medium_from_dataframe(model, medium_df):
     assert isinstance(medium_df, DataFrame)
-    for ex_reaction in model.exchanges:
+    for ex_reaction in model.boundary:
         if ex_reaction.id in medium_df.reaction_id.values:
             medium_row = medium_df[medium_df.reaction_id == ex_reaction.id]
             ex_reaction.lower_bound = medium_row.lower_bound.values[0]

cameo/flux_analysis/analysis.py

@@ -169,7 +169,7 @@ def find_blocked_reactions(model):
 
     """
     with model:
-        for exchange in model.exchanges:
+        for exchange in model.boundary:
             exchange.bounds = (-9999, 9999)
         fva_solution = flux_variability_analysis(model)
     return frozenset(

cameo/flux_analysis/structural.py

@@ -302,7 +302,7 @@ def __init__(self, model, reactions=None, c=1e-5, copy=True, change_bounds=True)
             self._elementary_mode_generator = self.__generate_elementary_modes()
 
     def __set_exchange_bounds(self):
-        exchanges = self.model.exchanges
+        exchanges = self.model.boundary
         min_bound = min(exchange.lower_bound for exchange in exchanges)
         max_bound = max(exchange.upper_bound for exchange in exchanges)
         for exchange in exchanges:

cameo/flux_analysis/util.py

@@ -49,7 +49,7 @@ def remove_infeasible_cycles(model, fluxes, fix=()):
     """
     with model:
         # make sure the original object is restored
-        exchange_reactions = model.exchanges
+        exchange_reactions = model.boundary
         exchange_ids = [exchange.id for exchange in exchange_reactions]
         internal_reactions = [reaction for reaction in model.reactions if reaction.id not in exchange_ids]
         for exchange in exchange_reactions:

cameo/strain_design/deterministic/flux_variability_based.py

@@ -197,8 +197,8 @@ def __init__(self, design_space_model, objective, variables=None, reference_mode
         reference_blocked_reactions = find_blocked_reactions_nullspace(self.reference_model, self.reference_nullspace)
         self.exclude += [reaction.id for reaction in reference_blocked_reactions]
 
-        self.exclude += [reaction.id for reaction in self.design_space_model.exchanges]
-        self.exclude += [reaction.id for reaction in self.reference_model.exchanges]
+        self.exclude += [reaction.id for reaction in self.design_space_model.boundary]
+        self.exclude += [reaction.id for reaction in self.reference_model.boundary]
 
         self.exclude += [reaction.id for reaction in self.design_space_model.reactions
                          if _BIOMASS_RE_.match(reaction.id)]
@@ -856,7 +856,7 @@ def run(self, target=None, max_enforced_flux=0.9, number_of_results=10, exclude=
         ndecimals = config.ndecimals
 
         # Exclude list
-        exclude = list(exclude) + model.exchanges
+        exclude = list(exclude) + model.boundary
         exclude_ids = [target.id]
         for reaction in exclude:
             if isinstance(reaction, Reaction):

cameo/strain_design/deterministic/linear_programming.py

@@ -157,7 +157,7 @@ def _reduce_to_nullspace(self, reactions):
     def _build_problem(self, exclude_reactions, use_nullspace_simplification):
         logger.debug("Starting to formulate OptKnock problem")
 
-        self.essential_reactions = find_essential_reactions(self._model, processes=1).union(self._model.exchanges)
+        self.essential_reactions = find_essential_reactions(self._model, processes=1).union(self._model.boundary)
         if exclude_reactions:
             self.exclude_reactions = set.union(
                 self.essential_reactions,

cameo/strain_design/heuristic/evolutionary/optimization.py

@@ -678,7 +678,7 @@ def __init__(self, reactions=None, essential_reactions=None, use_nullspace_simpl
         if use_nullspace_simplification:
             ns = nullspace(create_stoichiometric_array(self.model))
             dead_ends = set(find_blocked_reactions_nullspace(self.model, ns=ns))
-            exchanges = set(self.model.exchanges)
+            exchanges = set(self.model.boundary)
             reactions = [
                 r for r in self.model.reactions
                 if (r not in exchanges) and (
@@ -694,7 +694,7 @@ def __init__(self, reactions=None, essential_reactions=None, use_nullspace_simpl
         else:
             groups = None
             to_keep = set(r.id for r in self.model.reactions)
-            to_keep.difference_update(r.id for r in self.model.exchanges)
+            to_keep.difference_update(r.id for r in self.model.boundary)
             to_keep.difference_update(self.essential_reactions)
             to_keep = list(to_keep)
 

cameo/strain_design/pathway_prediction/pathway_predictor.py

@@ -246,7 +246,7 @@ def __init__(self, model, universal_model=None, mapping=None, compartment_regexp
         except SolverNotFound:
             logger.info('cplex not available for pathway predictions.')
 
-        self.new_reactions = self._extend_model(model.exchanges)
+        self.new_reactions = self._extend_model(model.boundary)
 
         logger.debug("Adding adapter reactions to connect model with universal model.")
         self.adpater_reactions = util.create_adapter_reactions(model.metabolites, self.universal_model,
@@ -421,7 +421,7 @@ def _add_switches(self, reactions):
         self._y_vars_ids = [var.name for var in y_vars]
 
     def _extend_model(self, original_exchanges):
-        for exchange in self.model.exchanges:
+        for exchange in self.model.boundary:
             if len(exchange.reactants) > 0 >= exchange.lower_bound:
                 exchange.upper_bound = 999999.
 
@@ -431,7 +431,7 @@ def _extend_model(self, original_exchanges):
                                       r in original_exchanges for m, coeff in six.iteritems(r.metabolites)
                                       if len(r.metabolites) == 1 and coeff < 0 < r.upper_bound]
 
-        universal_exchanges = self.universal_model.exchanges
+        universal_exchanges = self.universal_model.boundary
         for reaction in self.universal_model.reactions:
             if reaction in self.model.reactions:
                 continue

tests/conftest.py

@@ -42,7 +42,7 @@ def iaf1260(data_directory):
 @pytest.fixture(scope="session")
 def universal_model(data_directory):
     universal = load_model(join(data_directory, 'iJO1366.xml'), sanitize=False)
-    universal.remove_reactions(universal.exchanges)
+    universal.remove_reactions(universal.boundary)
     return universal
 
 

tests/test_api.py

@@ -27,7 +27,7 @@
 MODELS = os.path.dirname(models.__file__)
 
 UNIVERSALMODEL = load_model(os.path.join(MODELS, 'json/iJO1366.json'))
-UNIVERSALMODEL.remove_reactions(UNIVERSALMODEL.exchanges)
+UNIVERSALMODEL.remove_reactions(UNIVERSALMODEL.boundary)
 
 
 def test_api():