fix: include new metabolites in gapfilling

cobra/flux_analysis/gapfilling.py

@@ -118,12 +118,23 @@ def extend_model(self, exchange_reactions=False, demand_reactions=True):
         """
         for rxn in self.universal.reactions:
             rxn.gapfilling_type = 'universal'
+        new_metabolites = self.universal.metabolites.query(
+            lambda metabolite: metabolite not in self.model.metabolites
+                                                           )
+        self.model.add_metabolites(new_metabolites)
+        existing_exchanges = []
+        for rxn in self.universal.boundary:
+            existing_exchanges = existing_exchanges + \
+                [met.id for met in list(rxn.metabolites)]
+
         for met in self.model.metabolites:
             if exchange_reactions:
-                rxn = self.universal.add_boundary(
-                    met, type='exchange_smiley', lb=-1000, ub=0,
-                    reaction_id='EX_{}'.format(met.id))
-                rxn.gapfilling_type = 'exchange'
+                # check for exchange reaction in model already
+                if met.id not in existing_exchanges:
+                    rxn = self.universal.add_boundary(
+                        met, type='exchange_smiley', lb=-1000, ub=0,
+                        reaction_id='EX_{}'.format(met.id))
+                    rxn.gapfilling_type = 'exchange'
             if demand_reactions:
                 rxn = self.universal.add_boundary(
                     met, type='demand_smiley', lb=0, ub=1000,
@@ -158,7 +169,7 @@ def add_switches_and_objective(self):
                     for r in self.model.reactions)
         prob = self.model.problem
         for rxn in self.model.reactions:
-            if not hasattr(rxn, 'gapfilling_type') or rxn.id.startswith('DM_'):
+            if not hasattr(rxn, 'gapfilling_type'):
                 continue
             indicator = prob.Variable(
                 name='indicator_{}'.format(rxn.id), lb=0, ub=1, type='binary')
@@ -233,6 +244,9 @@ def fill(self, iterations=1):
 
     def validate(self, reactions):
         with self.original_model as model:
+            mets = [x.metabolites for x in reactions]
+            all_keys = set().union(*(d.keys() for d in mets))
+            model.add_metabolites(all_keys)
             model.add_reactions(reactions)
             model.slim_optimize()
             return (model.solver.status == OPTIMAL and

cobra/test/test_flux_analysis/test_gapfilling.py

@@ -51,6 +51,34 @@ def test_gapfilling(salmonella):
     assert len(result[1]) == 1
     assert {i[0].id for i in result} == {"EX_b", "EX_c"}
 
+    # # Gapfilling solution adds metabolites not present in original model
+    # test for when demand = T
+    # a demand reaction must be added to clear new metabolite
+    universal_noDM = Model()
+    a2b = Reaction("a2b")
+    universal_noDM.add_reactions([a2b])
+    a2b.build_reaction_from_string("a --> b + d", verbose=False)
+    result = gapfill(m, universal_noDM,
+                     exchange_reactions=False,
+                     demand_reactions=True)[0]
+    # add reaction a2b and demand reaction to clear met d
+    assert len(result) == 2
+    assert "a2b" in [x.id for x in result]
+
+    # test for when demand = False
+    # test for when metabolites are added to the model and
+    # must be cleared by other reactions in universal model
+    # (i.e. not necessarily a demand reaction)
+    universal_withDM = universal_noDM.copy()
+    d_dm = Reaction("d_dm")
+    universal_withDM.add_reactions([d_dm])
+    d_dm.build_reaction_from_string("d -->", verbose=False)
+    result = gapfill(m, universal_withDM,
+                     exchange_reactions=False,
+                     demand_reactions=False)[0]
+    assert len(result) == 2
+    assert "a2b" in [x.id for x in result]
+
     # somewhat bigger model
     universal = Model("universal_reactions")
     with salmonella as model: