diff --git a/src/enzyme_usage/getEnzymeBottlenecks.m b/src/enzyme_usage/getEnzymeBottlenecks.m index f345476f..83f6dbe5 100644 --- a/src/enzyme_usage/getEnzymeBottlenecks.m +++ b/src/enzyme_usage/getEnzymeBottlenecks.m @@ -1,9 +1,20 @@ function bottlenecks = getEnzymeBottlenecks(model, varargin) % getEnzymeBottlenecks Top-N enzymes limiting the current objective. % -% Solves the model and ranks every enzyme by the absolute shadow price of -% its prot_ mass-balance constraint --- the enzymes whose extra -% availability would help the objective the most. +% Solves the model and ranks every enzyme by the change of the objective +% after the upper bound of its usage_prot_ reaction is raised by +% foldChange and the model is solved again. The objective is maximised, so +% the gain is the rise of the maximised objective. Only enzymes whose usage runs at +% its upper bound are solved again; raising a bound that does not bind +% cannot change the optimum, so every other enzyme has an objectiveGain of 0 +% without a solve. +% +% The shadow price of the prot_ mass-balance constraint is reported for +% reference and does not enter the order: at a degenerate optimum, which is +% common in enzyme-constrained models, several constraints are active and +% the solver assigns the dual to one of them depending on the basis. +% Enzymes with equal objectiveGain (rounded to 9 decimals) keep their order +% in model.ec.enzymes. % % Parameters % ---------- @@ -16,31 +27,54 @@ % -------------------- % top : double % number of enzymes to return (default 10). +% foldChange : double +% a value how much to increase the usage upper bound of an enzyme that +% runs at its upper bound; must be greater than 1 (default 2). % % Returns % ------- % bottlenecks : table -% one row per returned enzyme, sorted by descending absolute shadow -% price, with columns uniprot, gene, shadowPrice, flux, capUsage -% (abs(flux)/upperBound, NaN if upperBound is 0) and upperBound. +% one row per returned enzyme, sorted by descending objectiveGain, with +% columns uniprot, gene, objectiveGain (objective after raising the +% bound minus the objective of the unchanged model; 0 for enzymes not at +% their upper bound, for gains below 1e-10 and for solves that are not +% optimal), shadowPrice, flux, capUsage (abs(flux)/upperBound, NaN if +% upperBound is 0), upperBound and multipleActive. An enzyme is at its +% upper bound when upperBound is above 0 and capUsage is at least +% 1 - 1e-6. multipleActive is true when at least one other constraint +% is active: another enzyme at its upper bound, or prot_pool_exchange at +% its upper bound (within 1e-6 relative). Where it is true, the dual is +% not informative about a single enzyme. The table property +% UserData.nResolved holds the number of enzymes that were solved again +% and UserData.activeConstraints the number of active constraints. +% When no enzyme is at its usage bound a warning (getEnzymeBottlenecks:noEnzymeAtBound) +% says that no single enzyme limits the objective; every objectiveGain is then 0 and +% the rows are in model order. % % Examples % -------- % bottlenecks = getEnzymeBottlenecks(ecModel); -% bottlenecks = getEnzymeBottlenecks(ecModel, 'top', 5); +% bottlenecks = getEnzymeBottlenecks(ecModel, 'top', 5, 'foldChange', 3); -p = parseGECKOargs(varargin, {'top', []}); +p = parseGECKOargs(varargin, {'top', []; 'foldChange', []}); top = p.top; +foldChange = p.foldChange; if isempty(top) top = 10; end +if isempty(foldChange) + foldChange = 2; +end +if ~isnumeric(foldChange) || ~isscalar(foldChange) || ~(foldChange > 1) + error('getEnzymeBottlenecks:foldChange', 'foldChange must be greater than 1') +end if isfield(model.ec, 'geckoLight') && model.ec.geckoLight error(['getEnzymeBottlenecks requires a full ecModel, not gecko-light: ' ... 'light models have no per-enzyme prot_ / usage_prot_ machinery.']) end -sol = solveLP(model); +[sol, hs] = solveLP(model); if sol.stat ~= 1 error('getEnzymeBottlenecks: solver status %d; cannot rank bottlenecks.', sol.stat) end @@ -65,8 +99,46 @@ capUsage = abs(flux) ./ upperBound; capUsage(upperBound == 0) = NaN; -bottlenecks = table(uniprot, gene, shadowPrice, flux, capUsage, upperBound); -[~, order] = sort(abs(bottlenecks.shadowPrice), 'descend'); +nEnz = numel(uniprot); +atCap = upperBound > 0 & capUsage >= 1 - 1e-6; +base = sol.f; +objectiveGain = zeros(nEnz, 1); +for i = find(atCap)' + tempModel = model; + tempModel.ub(rxnIdx(i)) = upperBound(i) * foldChange; + [tempSol, hs] = solveLP(tempModel, 0, [], hs); + if tempSol.stat == 1 + gain = tempSol.f - base; + if gain >= 1e-10 + objectiveGain(i) = gain; + end + end +end + +activeConstraints = nnz(atCap); +poolActive = false; +poolIdx = find(strcmp(model.rxns, 'prot_pool_exchange'), 1); +if ~isempty(poolIdx) && isfinite(model.ub(poolIdx)) && model.ub(poolIdx) > 0 ... + && abs(sol.x(poolIdx)) >= model.ub(poolIdx) * (1 - 1e-6) + activeConstraints = activeConstraints + 1; + poolActive = true; +end +if ~any(atCap) + if poolActive + state = 'the protein pool exchange is at its bound'; + else + state = 'no capacity constraint is active'; + end + warning('getEnzymeBottlenecks:noEnzymeAtBound', ... + ['getEnzymeBottlenecks: no enzyme is at its usage bound, so no single enzyme limits ' ... + 'the objective; %s. Every objectiveGain is 0 and the rows are in model order.'], state) +end +multipleActive = (activeConstraints - double(atCap)) >= 1; + +bottlenecks = table(uniprot, gene, objectiveGain, shadowPrice, flux, capUsage, upperBound, multipleActive); +[~, order] = sortrows([-round(objectiveGain, 9), (1:nEnz)']); bottlenecks = bottlenecks(order, :); bottlenecks = bottlenecks(1:min(top, height(bottlenecks)), :); +bottlenecks.Properties.UserData.nResolved = nnz(atCap); +bottlenecks.Properties.UserData.activeConstraints = activeConstraints; end diff --git a/test/unit_tests/geckoCoreFunctionTests.m b/test/unit_tests/geckoCoreFunctionTests.m index 308b830f..4b92695c 100644 --- a/test/unit_tests/geckoCoreFunctionTests.m +++ b/test/unit_tests/geckoCoreFunctionTests.m @@ -946,26 +946,19 @@ function testWriteOpenKineticsPredictorInputReactionSubsetIndex_tc0023(testCase) end -function testGetEnzymeBottlenecksRanksByShadowPrice_tc0024(testCase) - % Verifies getEnzymeBottlenecks solves the model, ranks enzymes by the - % absolute shadow price of their prot_ mass-balance constraint, and - % returns the top N rows. Checks the primal-derived columns (flux, - % capUsage, upperBound) exactly, that rows are sorted by non-increasing - % |shadowPrice|, and that 'top' limits the row count -- but not the exact - % shadow-price value for an enzyme carrying zero flux at the optimum, - % since that value is LP-dual-degenerate and can vary with the solver's - % platform/presolve path. - geckoPath = findGECKOroot; - adapter = ModelAdapterManager.getAdapter(fullfile(geckoPath,'test','unit_tests','ecTestGEM', 'TestGEMAdapter.m')); - model = getGeckoTestModel(); - model.lb(strcmp(model.rxns,'R2')) = 0; model.ub(strcmp(model.rxns,'R2')) = 0; - model.lb(strcmp(model.rxns,'R4')) = 0; model.ub(strcmp(model.rxns,'R4')) = 0; - ecModel = makeEcModel(model, false, adapter); - ecModel = getECfromGEM(ecModel); - ecModel.ec.kcat(:) = 10; - ecModel.ec.source(:) = {'manual'}; - ecModel = applyKcatConstraints(ecModel); - ecModel = setProtPoolSize(ecModel, 0.5, 0.5, 0.5, adapter); +function testGetEnzymeBottlenecksRanksByObjectiveGain_tc0024(testCase) + % Verifies getEnzymeBottlenecks solves the model and returns the top N + % rows with the columns objectiveGain, shadowPrice, flux, capUsage, + % upperBound and multipleActive. In this model no usage reaction runs at + % its upper bound and the protein pool is the only active constraint, so + % no enzyme is solved again, every objectiveGain is 0, the order is the + % order of model.ec.enzymes, and every row is marked multipleActive. + % Checks the primal-derived columns (flux, capUsage, upperBound) exactly, + % and that 'top' limits the row count -- but not the shadow-price value + % for an enzyme carrying zero flux at the optimum, since that value is + % LP-dual-degenerate and can vary with the solver's platform/presolve + % path. + ecModel = buildBottleneckTestModel(0.5); sol = solveLP(ecModel); verifyEqual(testCase, sol.stat, 1) @@ -973,21 +966,196 @@ function testGetEnzymeBottlenecksRanksByShadowPrice_tc0024(testCase) bottlenecks = getEnzymeBottlenecks(ecModel); verifyEqual(testCase, height(bottlenecks), 5) + verifyEqual(testCase, bottlenecks.Properties.VariableNames, ... + {'uniprot','gene','objectiveGain','shadowPrice','flux','capUsage','upperBound','multipleActive'}) + + verifyEqual(testCase, bottlenecks.uniprot, {'P1';'P2';'P3';'P4';'P5'}) + verifyEqual(testCase, bottlenecks.objectiveGain, zeros(5, 1)) + verifyEqual(testCase, bottlenecks.flux, [0;0;0;55.555555555556;69.444444444444], 'AbsTol', 1e-9) + verifyEqual(testCase, bottlenecks.capUsage, [0;0;0;0.055555555555556;0.069444444444444], 'AbsTol', 1e-9) + verifyEqual(testCase, bottlenecks.upperBound, repmat(1000, 5, 1), 'AbsTol', 1e-9) + verifyTrue(testCase, all(bottlenecks.multipleActive)) + verifyEqual(testCase, bottlenecks.Properties.UserData.nResolved, 0) + verifyEqual(testCase, bottlenecks.Properties.UserData.activeConstraints, 1) - % Primal-derived values: unique regardless of dual degeneracy. - [~, order] = ismember({'P1','P2','P3','P4','P5'}, bottlenecks.uniprot); - verifyEqual(testCase, bottlenecks.flux(order), [0;0;0;55.555555555556;69.444444444444], 'AbsTol', 1e-9) - verifyEqual(testCase, bottlenecks.capUsage(order), [0;0;0;0.055555555555556;0.069444444444444], 'AbsTol', 1e-9) - verifyEqual(testCase, bottlenecks.upperBound(order), repmat(1000, 5, 1), 'AbsTol', 1e-9) - - % Sorted by non-increasing |shadow price| --- true whichever degenerate - % vertex the solver lands on. - verifyTrue(testCase, all(diff(abs(bottlenecks.shadowPrice)) <= 1e-9)) - - % top limits the row count to exactly what was asked for, whatever the - % tie-break among degenerate rows put there. top2 = getEnzymeBottlenecks(ecModel, 'top', 2); verifyEqual(testCase, height(top2), 2) + verifyEqual(testCase, top2.uniprot, {'P1';'P2'}) +end + + +function testGetEnzymeBottlenecksOrdersByResolvedGain_tc0277(testCase) + % Verifies that of two enzymes running at their upper bound, the one + % whose bound limits growth is ranked first with the growth gain of a + % solve with its bound doubled (positive for a maximised growth + % objective), and the one whose raised bound does not change growth (P3, + % whose usage is fixed by equal lower and upper bounds) has a gain of 0. + % Both enzymes are solved again, so nResolved is 2 and, with two active + % constraints, both rows are multipleActive. Only P5 and P3 are asserted + % as active; which of P1, P2 and P4 carry flux in the parallel routes is + % not unique. + ecModel = buildBottleneckTestModel(100, false); + ecModel.ub(strcmp(ecModel.rxns, 'usage_prot_P5')) = 20; + p3 = strcmp(ecModel.rxns, 'usage_prot_P3'); + ecModel.lb(p3) = 1; + ecModel.ub(p3) = 1; + + base = solveLP(ecModel); + verifyEqual(testCase, base.stat, 1) + raised = ecModel; + raised.ub(strcmp(raised.rxns, 'usage_prot_P5')) = 40; + expectedGain = solveLP(raised).f - base.f; + verifyGreaterThan(testCase, expectedGain, 1) + + bottlenecks = getEnzymeBottlenecks(ecModel); + verifyEqual(testCase, bottlenecks.uniprot{1}, 'P5') + verifyEqual(testCase, bottlenecks.objectiveGain(1), expectedGain, 'AbsTol', 1e-9) + verifyEqual(testCase, bottlenecks.objectiveGain(2:end), zeros(4, 1)) + row = strcmp(bottlenecks.uniprot, 'P3'); + verifyEqual(testCase, bottlenecks.capUsage(row), 1, 'AbsTol', 1e-9) + verifyTrue(testCase, all(bottlenecks.multipleActive(1) & bottlenecks.multipleActive(row))) + verifyEqual(testCase, bottlenecks.Properties.UserData.nResolved, 2) + verifyEqual(testCase, bottlenecks.Properties.UserData.activeConstraints, 2) + % rows with equal gain keep the order of model.ec.enzymes + verifyEqual(testCase, bottlenecks.uniprot(2:end), {'P1';'P2';'P3';'P4'}) +end + + +function testGetEnzymeBottlenecksRespectsFoldChange_tc0278(testCase) + % Verifies the foldChange argument sets the factor by which the upper + % bound of an enzyme at its bound is raised: the gain of P4 equals the + % growth change of a solve with its bound multiplied by 1.5, by 2 + % (default) and by 3, and the gain differs between factors. + ecModel = buildBottleneckTestModel(100); + ecModel.ub(strcmp(ecModel.rxns, 'usage_prot_P4')) = 16; + base = solveLP(ecModel).f; + + gains = zeros(1, 3); + factors = [1.5 2 3]; + for i = 1:3 + raised = ecModel; + raised.ub(strcmp(raised.rxns, 'usage_prot_P4')) = 16 * factors(i); + expected = solveLP(raised).f - base; + b = getEnzymeBottlenecks(ecModel, 'foldChange', factors(i)); + verifyEqual(testCase, b.uniprot{1}, 'P4') + verifyEqual(testCase, b.objectiveGain(1), expected, 'AbsTol', 1e-9) + gains(i) = b.objectiveGain(1); + end + verifyTrue(testCase, all(diff(gains) > 0)) + default = getEnzymeBottlenecks(ecModel); + verifyEqual(testCase, default.objectiveGain(1), gains(2), 'AbsTol', 1e-9) +end + + +function testGetEnzymeBottlenecksNonBindingEnzymeHasZeroGain_tc0279(testCase) + % Verifies an enzyme carrying flux below its upper bound is not solved + % again and has objectiveGain 0, while the one enzyme at its bound has a + % positive gain. nResolved is 1 and the protein pool is not active, so + % the enzyme at its bound is not multipleActive, while every other row + % is, because one constraint besides its own is active. + ecModel = buildBottleneckTestModel(100); + ecModel.ub(strcmp(ecModel.rxns, 'usage_prot_P4')) = 16; + ecModel.ub(strcmp(ecModel.rxns, 'usage_prot_P5')) = 40; + + bottlenecks = getEnzymeBottlenecks(ecModel); + verifyEqual(testCase, bottlenecks.uniprot{1}, 'P4') + verifyGreaterThan(testCase, bottlenecks.objectiveGain(1), 1) + verifyFalse(testCase, bottlenecks.multipleActive(1)) + row = strcmp(bottlenecks.uniprot, 'P5'); + verifyEqual(testCase, bottlenecks.flux(row), 20, 'AbsTol', 1e-9) + verifyEqual(testCase, bottlenecks.capUsage(row), 0.5, 'AbsTol', 1e-9) + verifyEqual(testCase, bottlenecks.objectiveGain(row), 0) + verifyTrue(testCase, bottlenecks.multipleActive(row)) + verifyEqual(testCase, bottlenecks.Properties.UserData.nResolved, 1) + verifyEqual(testCase, bottlenecks.Properties.UserData.activeConstraints, 1) +end + + +function testGetEnzymeBottlenecksErrorsWhenNotOptimal_tc0280(testCase) + % Verifies getEnzymeBottlenecks raises an error when the model cannot be + % solved to optimality (here: the growth reaction is forced above what + % the protein pool allows). + ecModel = buildBottleneckTestModel(0.5); + ecModel.lb(ecModel.c == 1) = 1000; + verifyNotEqual(testCase, solveLP(ecModel).stat, 1) + try + getEnzymeBottlenecks(ecModel); + message = ''; + catch err + message = err.message; + end + verifyTrue(testCase, startsWith(message, 'getEnzymeBottlenecks: solver status')) +end + + +function testGetEnzymeBottlenecksRejectsInvalidFoldChange_tc0281(testCase) + % Verifies foldChange values of 1, below 1, NaN and non-scalar raise + % getEnzymeBottlenecks:foldChange. + ecModel = buildBottleneckTestModel(100); + for value = {1, 0.5, 0, NaN, [2 3]} + verifyError(testCase, @() getEnzymeBottlenecks(ecModel, 'foldChange', value{1}), ... + 'getEnzymeBottlenecks:foldChange') + end +end + + +function testGetEnzymeBottlenecksWarnsWhenOnlyThePoolBinds_tc0282(testCase) + % Verifies that with every usage upper bound out of reach and a binding + % protein pool no enzyme is at its bound: the warning names the pool, + % nResolved is 0, every objectiveGain is 0 and the rows are in model order. + ecModel = buildBottleneckTestModel(100); + usageIdx = startsWith(ecModel.rxns, 'usage_prot_'); + ecModel.ub(usageIdx) = 1e6; + + lastwarn(''); + bottlenecks = verifyWarning(testCase, @() getEnzymeBottlenecks(ecModel), ... + 'getEnzymeBottlenecks:noEnzymeAtBound'); + verifyEqual(testCase, bottlenecks.Properties.UserData.nResolved, 0) + verifyEqual(testCase, bottlenecks.Properties.UserData.activeConstraints, 1) + verifyEqual(testCase, bottlenecks.objectiveGain, zeros(height(bottlenecks), 1)) + verifyEqual(testCase, bottlenecks.uniprot, ecModel.ec.enzymes(1:height(bottlenecks))) + verifyTrue(testCase, all(bottlenecks.multipleActive)) + verifyNotEmpty(testCase, strfind(lastwarn, 'the protein pool exchange is at its bound')) +end + + +function testGetEnzymeBottlenecksWarnsWhenNoConstraintIsActive_tc0283(testCase) + % Verifies that with the growth reaction R5 capped at 10 and the pool and + % usage bounds out of reach the warning says that no capacity constraint + % is active. + ecModel = buildBottleneckTestModel(1e6); + usageIdx = startsWith(ecModel.rxns, 'usage_prot_'); + ecModel.ub(usageIdx) = 1e6; + ecModel.ub(strcmp(ecModel.rxns, 'R5')) = 10; + + lastwarn(''); + bottlenecks = verifyWarning(testCase, @() getEnzymeBottlenecks(ecModel), ... + 'getEnzymeBottlenecks:noEnzymeAtBound'); + verifyEqual(testCase, bottlenecks.Properties.UserData.activeConstraints, 0) + verifyFalse(testCase, any(bottlenecks.multipleActive)) + verifyNotEmpty(testCase, strfind(lastwarn, 'no capacity constraint is active')) +end + + +function ecModel = buildBottleneckTestModel(poolSize, blockR2) + % Test model with kcat 10 for all enzymes, R4 blocked, R2 blocked unless + % blockR2 is false, and a protein pool of poolSize. + if nargin < 2 + blockR2 = true; + end + geckoPath = findGECKOroot; + adapter = ModelAdapterManager.getAdapter(fullfile(geckoPath,'test','unit_tests','ecTestGEM', 'TestGEMAdapter.m')); + model = getGeckoTestModel(); + if blockR2 + model.lb(strcmp(model.rxns,'R2')) = 0; model.ub(strcmp(model.rxns,'R2')) = 0; + end + model.lb(strcmp(model.rxns,'R4')) = 0; model.ub(strcmp(model.rxns,'R4')) = 0; + ecModel = makeEcModel(model, false, adapter); + ecModel = getECfromGEM(ecModel); + ecModel.ec.kcat(:) = 10; + ecModel.ec.source(:) = {'manual'}; + ecModel = applyKcatConstraints(ecModel); + ecModel = setProtPoolSize(ecModel, poolSize, 0.5, 0.5, adapter); end