From 334cbf0f314de7bd3ca8087d56c31e5a33852d7d Mon Sep 17 00:00:00 2001 From: Eduard Kerkhoven Date: Fri, 17 Jul 2026 11:21:09 +0200 Subject: [PATCH] Fix inverted allowExcretion constraint sense in ftINITInternalAlg allowExcretion set csense 'L' on the metabolite rows, i.e. S*v <= 0. That allows a metabolite to be consumed without ever being produced -- free uptake, the exact opposite of excretion, and enough for the MILP to switch a reaction on by conjuring its substrates out of nothing. Since ftINIT locks step-1 picks in as essentialRxns, such a reaction is then kept for the rest of the run. The intended constraint is S*v >= 0 ('G'): a metabolite may be produced in excess of what is consumed, and the surplus is implicitly excreted. This is what the function's own docstring describes, what the parameter name says, and what runINIT does for the same flag -- it sets the metabolite row's upper bound to inf (metUB=inf(nMets,1)) in the slack formulation, which is S*v >= 0. Verified against optimizeProb directly: with a single row S*v and one producing reaction, csense 'L' gives max v = 0 (production forbidden) and 'G' gives max v = 10 (production allowed). Also give the else branch one csense character per row. A scalar '=' only works because gurobi broadcasts it; optimizeProb maps csense elementwise for glpk ('E'->'S') and cobra, so the scalar does not survive dispatch. Honest scope note: the synthetic model in tINIT does not discriminate between the two senses -- it has exchange reactions, which mask the effect, and '1+1', '2+1' and '2+0' all return the same reactions before and after this change. The fix rests on the constraint semantics above, not on an observed output change, and no test in the suite regresses. Tests: only '1+1' and 'full' were ever exercised, so the paper's 2-step series had no coverage. Adds ftINITSeriesVariantsRun pinning '2+1' and '2+0'; '2+0' keeps the GPR-less transport R2, as its docstring says it should. --- INIT/ftINITInternalAlg.m | 11 +++++++++-- testing/function_tests/tINIT.m | 26 ++++++++++++++++++++++++++ 2 files changed, 35 insertions(+), 2 deletions(-) diff --git a/INIT/ftINITInternalAlg.m b/INIT/ftINITInternalAlg.m index 19c0ddc4..0b990920 100644 --- a/INIT/ftINITInternalAlg.m +++ b/INIT/ftINITInternalAlg.m @@ -392,10 +392,17 @@ metVarInd = (1:nMetabolMets) + (length(prob.vartype) - nMetVars); if allowExcretion - prob.csense = [repmat('L', 1, length(milpModel.mets)), ... + %S*v >= 0: a metabolite may be produced in excess of what is consumed, + %and the surplus is implicitly excreted. 'L' would say S*v <= 0, which + %lets a metabolite be consumed without ever being produced -- free + %uptake, the opposite of excretion, and enough for the MILP to switch a + %reaction on by conjuring its substrates out of nothing. + prob.csense = [repmat('G', 1, length(milpModel.mets)), ... repmat('E', 1, length(prob.b) - length(milpModel.mets))]; else - prob.csense = '='; + %One character per row: optimizeProb maps csense elementwise for glpk and + %cobra, so a scalar '=' only happens to work on gurobi. + prob.csense = repmat('E', 1, length(prob.b)); end params.intTol = 10^-7; %This value is very important. If set too low diff --git a/testing/function_tests/tINIT.m b/testing/function_tests/tINIT.m index 5c77724c..71c5e1ac 100644 --- a/testing/function_tests/tINIT.m +++ b/testing/function_tests/tINIT.m @@ -259,6 +259,32 @@ function ftINITFullVsThreeStepRuns(testCase) testCase.verifyTrue(all(contains(mres2.rxns, expResult))); end + function ftINITSeriesVariantsRun(testCase) + % Only '1+1' and 'full' were ever exercised, so the 2-step series + % from the paper -- and the allowExcretion constraint they lean on + % -- had no coverage at all. + testCase.assumeMILPSolver(); + testModel = getTstModel(); + testParams = struct(); + evalc('prepData = prepINITModel(testModel, {}, {}, false, {}, ''s'');'); + arrayData.genes = testModel.genes; + arrayData.tissues = {'a'}; + arrayData.levels = getExprForRxnScore(getTstModelRxnScores()); + arrayData.threshold = 1; + + steps = getINITSteps([], '2+1'); + evalc(['resModel = ftINIT(prepData,arrayData.tissues{1},[],[],' ... + 'arrayData,[],steps,true,true,testParams,false);']); + % Same answer as the '1+1' series in ftINITPipelineRuns. + testCase.verifyEqual(resModel.rxns, {'R1';'R4';'R6';'R8';'R9';'R10'}); + + steps = getINITSteps([], '2+0'); + evalc(['resModel = ftINIT(prepData,arrayData.tissues{1},[],[],' ... + 'arrayData,[],steps,true,true,testParams,false);']); + % '2+0' skips step 3, so the GPR-less transport R2 survives. + testCase.verifyEqual(resModel.rxns, {'R1';'R2';'R4';'R6';'R8';'R9';'R10'}); + end + end end