Skip to content

getEnzymeBottlenecks ranks by the re-solved change of the objective #20

New issue

Have a question about this project? Sign up for a free GitHub account to open an issue and contact its maintainers and the community.

Sign up for GitHub

By clicking “Sign up for GitHub”, you agree to our terms of service and privacy statement. We’ll occasionally send you account related emails.

Already on GitHub? Sign in to your account

Merged
edkerk merged 2 commits into from
Oct 31, 2026
Merged

getEnzymeBottlenecks ranks by the re-solved change of the objective #20

Show file tree
Hide file tree
Changes from all commits
Commits
Show all changes
2 commits
Select commit Hold shift + click to select a range
e2245cf
Rank getEnzymeBottlenecks by the re-solved change of the objective
edkerk Oct 31, 2026
84fe193
getEnzymeBottlenecks warns when no enzyme is at its usage bound
edkerk Oct 31, 2026
File filter

Filter by extension

Filter by extension
Conversations
Failed to load comments.
Loading
Jump to
Jump to file
Failed to load files.
Loading
Diff view
Diff view
94 changes: 83 additions & 11 deletions src/enzyme_usage/getEnzymeBottlenecks.m
View file
Open in desktop
Original file line number Diff line number Diff line change
@@ -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_<id> 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_<id> 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_<id> 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
% ----------
Expand All @@ -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_<id> / usage_prot_<id> machinery.'])
end

sol = solveLP(model);
[sol, hs] = solveLP(model);
if sol.stat ~= 1
error('getEnzymeBottlenecks: solver status %d; cannot rank bottlenecks.', sol.stat)
end
Expand All @@ -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
232 changes: 200 additions & 32 deletions test/unit_tests/geckoCoreFunctionTests.m
View file
Open in desktop
Original file line number Diff line number Diff line change
Expand Up @@ -946,48 +946,216 @@ function testWriteOpenKineticsPredictorInputReactionSubsetIndex_tc0023(testCase)
end


function testGetEnzymeBottlenecksRanksByShadowPrice_tc0024(testCase)
% Verifies getEnzymeBottlenecks solves the model, ranks enzymes by the
% absolute shadow price of their prot_<id> 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)
verifyEqual(testCase, sol.f, 50, 'AbsTol', 1e-9)

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


Expand Down
Loading