fix: gapFillTopological seed detection and zero-substrate reactions

gapfilling/gapFillTopological.m

@@ -31,7 +31,8 @@
 % --------------------------------
 % 'seeds' (default [])
 %     Metabolite IDs (cell array) available from the medium. Default:
-%     metabolites involved in uptake exchange reactions (lb < 0).
+%     metabolites of every exchange reaction whose bounds let it supply the
+%     metabolite, i.e. getExchangeRxns 'uptake' plus 'reverse'.
 % 'targets' (default [])
 %     Metabolite IDs (cell array) that should be produced. Default:
 %     substrates of the objective (biomass) reaction.
@@ -67,9 +68,22 @@
 
 % ---- Identify seed metabolites ----
 if isempty(seeds)
-    % Default: metabolites in exchange reactions that allow uptake (lb < 0)
-    [~, exchIdx] = getExchangeRxns(model);
-    uptakeIdx = exchIdx(model.lb(exchIdx) < 0);
+    % Default: metabolites that an exchange reaction can supply. Ask
+    % getExchangeRxns rather than testing lb < 0 here: RAVEN writes an uptake
+    % exchange either as "=> met" (addExchangeRxns 'in', so lb = 0 and a
+    % positive flux supplies the metabolite) or as "met =>" with lb < 0. Only
+    % the second has lb < 0, so a hand-rolled lb < 0 test silently misses
+    % every 'in' exchange.
+    %
+    % 'uptake' alone is not enough either: it means the bounds allow *only*
+    % uptake, so a reversible exchange ("met <=>") is classified 'reverse'
+    % even though it supplies the metabolite just as well. The union of the
+    % two is exactly "forward flux produces it, or reverse flux does", which
+    % is what analyse_topology's lower_bound < 0 test means under cobra's
+    % convention, where every exchange is written "met <=>".
+    [~, onlyUptakeIdx] = getExchangeRxns(model, 'uptake');
+    [~, reversibleIdx] = getExchangeRxns(model, 'reverse');
+    uptakeIdx = union(onlyUptakeIdx, reversibleIdx);
     if isempty(uptakeIdx)
         warning(['gapFillTopological: no uptake exchange reactions found. ' ...
             'Provide seeds manually via the ''seeds'' option.']);
@@ -150,6 +164,20 @@
 qHead  = 1;
 reachable(seedMetIdx) = true;
 
+% Reactions with no substrates have nothing to wait for and fire straight
+% away. They are in no subOf list, so the countdown below never reaches them
+% and their products would otherwise stay unreachable forever. This is where
+% "=> met" uptake exchanges enter the scope.
+for j = find(subCountFwd == 0)'
+    firedFwd(j) = true;
+    for p = rxnProds{j}
+        if ~reachable(p)
+            reachable(p) = true;
+            queue(end+1) = p; %#ok<AGROW>
+        end
+    end
+end
+
 while qHead <= numel(queue)
     m = queue(qHead);
     qHead = qHead + 1;

testing/function_tests/tGapfilling.m

@@ -140,6 +140,37 @@ function gapFillTopologicalIdentifiesGaps(testCase)
             testCase.verifyNumElements(result.reachableMets, numel(gapModel.mets));
         end
 
+        function gapFillTopologicalReachesThroughInExchange(testCase)
+            % A -> B -> C fed by an 'in' exchange on A. addExchangeRxns writes
+            % 'in' as "=> A", which has lb = 0, so a seed test of lb < 0 finds
+            % no uptake at all and every metabolite comes back blocked.
+            m = testCase.chainModel();
+            m = addExchangeRxns(m, 'in', {'A'});
+            evalc(['result = gapFillTopological(m, m, ''targets'', ' ...
+                '{''A'',''B'',''C''}, ''verbose'', false);']);
+            testCase.verifyTrue(all(result.reachableMets));
+            testCase.verifyEmpty(result.blockedMets);
+        end
+
+        function gapFillTopologicalReachesThroughReversibleExchange(testCase)
+            % "A <=>" supplies A too, but is 'reverse' rather than 'uptake'.
+            m = testCase.chainModel();
+            m = addExchangeRxns(m, 'both', {'A'});
+            evalc(['result = gapFillTopological(m, m, ''targets'', ' ...
+                '{''A'',''B'',''C''}, ''verbose'', false);']);
+            testCase.verifyTrue(all(result.reachableMets));
+        end
+
+        function gapFillTopologicalBlocksWithoutExchange(testCase)
+            % Without any exchange nothing is producible: the positive control
+            % for the two tests above.
+            m = testCase.chainModel();
+            evalc(['result = gapFillTopological(m, m, ''seeds'', {}, ''targets'', ' ...
+                '{''A'',''B'',''C''}, ''verbose'', false);']);
+            testCase.verifyFalse(any(result.reachableMets));
+            testCase.verifyNumElements(result.blockedMets, 3);
+        end
+
         function gapFillMILPRepairsGrowth(testCase)
             testCase.assumeMILPSolver();
             modelDB = testCase.model; modelDB.id = 'DB';
@@ -190,4 +221,18 @@ function fillGapsDispatchesGapfillMILP(testCase)
         end
 
     end
+
+    methods (Access = private)
+        function m = chainModel(~)
+            % A -> B -> C, no exchanges; the caller adds the ones it needs.
+            m = struct();
+            m.rxns = {'R1';'R2'}; m.rxnNames = {'R1';'R2'};
+            m.mets = {'A';'B';'C'}; m.metNames = {'A';'B';'C'}; m.metComps = [1;1;1];
+            m.comps = {'c'}; m.compNames = {'c'};
+            m.S = sparse([-1 0; 1 -1; 0 1]);
+            m.lb = [0;0]; m.ub = [1000;1000]; m.rev = [0;0]; m.c = [0;1];
+            m.b = [0;0;0];
+            m.genes = {}; m.grRules = {'';''}; m.rxnGeneMat = sparse(2,0);
+        end
+    end
 end