Skip to content

sigmaFitter, getStandardKcat and getSubsetEcModel: plateau, subsystem case, empty subsystems, light ecModels #15

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 30, 2026
Merged

sigmaFitter, getStandardKcat and getSubsetEcModel: plateau, subsystem case, empty subsystems, light ecModels #15

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
7087741
sigmaFitter: return the start of the growth plateau on every solver
edkerk Oct 28, 2026
793593c
getStandardKcat, getSubsetEcModel: handle subSystem case and empty su…
edkerk Oct 29, 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
14 changes: 12 additions & 2 deletions src/change_model/getSubsetEcModel.m
View file
Open in desktop
Original file line number Diff line number Diff line change
Expand Up @@ -50,7 +50,11 @@
if ismember('standard', genesToRemove)
genesToRemove(strcmp(genesToRemove, 'standard')) = [];
end
smallEcModel = removeGenes(bigEcModel,genesToRemove,true,true,false);
if isempty(genesToRemove)
smallEcModel = bigEcModel;
else
smallEcModel = removeGenes(bigEcModel,genesToRemove,true,true,false);
end

% Remove genes from ec-structure
enzToRemove = ismember(smallEcModel.ec.genes,genesToRemove);
Expand All @@ -68,7 +72,13 @@
smallEcModel = removeReactions(smallEcModel,~keepRxns, true, true, true);

% Remove reactions from ec-structure
ecRxnsToRemove = ~ismember(smallEcModel.ec.rxns,smallEcModel.rxns);
% The reaction identifiers of a light ecModel carry a 4 character prefix
if smallEcModel.ec.geckoLight
ecModelRxns = extractAfter(smallEcModel.ec.rxns,4);
else
ecModelRxns = smallEcModel.ec.rxns;
end
ecRxnsToRemove = ~ismember(ecModelRxns,smallEcModel.rxns);

smallEcModel.ec.rxns(ecRxnsToRemove) = [];
smallEcModel.ec.kcat(ecRxnsToRemove) = [];
Expand Down
59 changes: 18 additions & 41 deletions src/gather_kcats/getStandardKcat.m
View file
Open in desktop
Original file line number Diff line number Diff line change
Expand Up @@ -133,13 +133,14 @@
end
% Map ec-rxns to model.rxns
[~,rxnIdx] = ismember(modRxns,model.rxns);
% Choose first subSystem
enzSubSystems = flattenCell(model.subSystems(rxnIdx));
enzSubSystems = enzSubSystems(:,1);
% Choose first subSystem ('' when a reaction has none)
enzSubSystems = cellfun(@firstSubSystem, model.subSystems(rxnIdx), 'UniformOutput', false);
if ~all(cellfun(@isempty, enzSubSystems))

% Make list of unique subsystems, and which rxns are linked to them
[enzSubSystem_names, ~, rxnToSub] = unique(enzSubSystems);
% Make list of unique subsystems, and which rxns are linked to them.
% Names that differ only in letter case are one subSystem.
[~, firstOccurrence, rxnToSub] = unique(lower(enzSubSystems));
enzSubSystem_names = enzSubSystems(firstOccurrence);
% Calculate the median kcat per subSystem, counting only reactions
% with a real (nonzero, non-NaN) kcat toward both the median itself
% and the threshold that decides whether a subSystem-specific value
Expand Down Expand Up @@ -289,8 +290,9 @@
% all(), is the right test -- all() only ever succeeds when the model
% has exactly one subSystem in total, silently falling back to the
% model-wide standardKcat for every reaction otherwise.
kcatSubSystemIdx = strcmpi(enzSubSystem_names, model.subSystems{rxnIdx}(1));
if any(kcatSubSystemIdx)
rxnSubSystem = firstSubSystem(model.subSystems{rxnIdx});
kcatSubSystemIdx = strcmpi(enzSubSystem_names, rxnSubSystem);
if ~isempty(rxnSubSystem) && any(kcatSubSystemIdx)
model.ec.kcat(end+1) = kcatSubSystem(kcatSubSystemIdx);
else
model.ec.kcat(end+1) = standardKcat;
Expand Down Expand Up @@ -319,39 +321,14 @@
end
end

function Cflat = flattenCell(C,strFlag)
%FLATTENCELL Flatten a nested column cell array into a matrix cell array.
%
% CFLAT = FLATTENCELL(C) takes a column cell array in which one or more
% entries is a nested cell array, and flattens it into a 2D matrix cell
% array, where the nested entries are spread across new columns.
%
% CFLAT = FLATTENCELL(C,STRFLAG) if STRFLAG is TRUE, empty entries in the
% resulting CFLAT will be replaced with empty strings {''}. Default = FALSE
if nargin < 2
strFlag = false;
end

% determine which entries are cells
cells = cellfun(@iscell,C);

% determine number of elements in each nested cell
cellsizes = cellfun(@numel,C);
cellsizes(~cells) = 1; % ignore non-cell entries

% flatten single-entry cells
Cflat = C;
Cflat(cells & (cellsizes == 1)) = cellfun(@(x) x{1},Cflat(cells & (cellsizes == 1)),'UniformOutput',false);

% iterate through multi-entry cells
multiCells = find(cellsizes > 1);
for i = 1:length(multiCells)
cellContents = Cflat{multiCells(i)};
Cflat(multiCells(i),1:length(cellContents)) = cellContents;
end

% change empty elements to strings, if specified
if ( strFlag )
Cflat(cellfun(@isempty,Cflat)) = {''};
function name = firstSubSystem(subSystems)
% First subSystem of a reaction as a character vector; '' when the reaction
% has none.
if ischar(subSystems)
name = subSystems;
elseif iscell(subSystems) && ~isempty(subSystems) && ischar(subSystems{1})
name = subSystems{1};
else
name = '';
end
end
9 changes: 8 additions & 1 deletion src/kcat_tuning/sigmaFitter.m
View file
Open in desktop
Original file line number Diff line number Diff line change
Expand Up @@ -95,7 +95,14 @@
error = num2str(((growthRate-solution.x(objPos))/growthRate)*100);
sigParam = [sigParam; sigma];
end
[~, minIndx] = min(errors);
% Growth that is flat from some sigma on (a target above what the model can
% reach) is returned by the solver with rounding noise. The smallest sigma
% whose growth is within 1e-9 (relative) of the best growth is chosen, so the
% result is the start of the plateau on every solver. The tolerance is
% expressed on the percent error.
[minError, minIndx] = min(errors);
errTol = 1e-9*100*abs(objValues(minIndx))/abs(growthRate);
minIndx = find(errors <= minError + errTol, 1);
sigma = sigParam(minIndx);
% The loop above leaves the protein pool sized for the *last* sigma tried
% (i=100, i.e. sigma=1), not the best-fitting one found by the search --
Expand Down
5 changes: 3 additions & 2 deletions test/VerificationMatrix.tsv
View file
Open in desktop
Original file line number Diff line number Diff line change
Expand Up @@ -7,9 +7,10 @@ change_model/getComplexData.m tc0003, tc0004 - File download not covered by test
change_model/makeEcModel.m tc0001, tc0002, tc0210, tc0212, tc0233, tc0234
change_model/setProtPoolSize.m tc0005
change_model/getKcatAcrossIsoenzymes.m tc0011, tc0252
change_model/getSubsetEcModel.m tc0258, tc0259
change_model/applyCustomKcats.m tc0227, tc0228, tc0229, tc0230
gather_kcats/fuzzyKcatMatching.m tc0010, tc0223, tc0224, tc0225, tc0226, tc0231, tc0239
gather_kcats/getStandardKcat.m
gather_kcats/getStandardKcat.m tc0256, tc0257
gather_kcats/mergeDlkcatAndFuzzyKcats.m tc0011; tc0241
gather_kcats/readDLKcatOutput.m tcm0001, tcm0002 - The test doesn't really check the function, only that some output is read
gather_kcats/runDLKcat.m tcm0001, tcm0002 - The test doesn't really check the function, only that we get something
Expand Down Expand Up @@ -54,7 +55,7 @@ kcat_tuning/evotune/* Has only been tested by running it manually, not tested by
kcat_tuning/fillCarbonNum.m tc0055 (deprecated-alias check only, not a real usage test)
kcat_tuning/findTopLimitations.m Not explicitly tested - The code has existed for a long time
kcat_tuning/sensitivityTuning.m tcm0001, tc0251, tc0253, tc0254
kcat_tuning/sigmaFitter.m Not explicitly tested - The code has existed for a long time
kcat_tuning/sigmaFitter.m tc0020, tc0255
kcat_tuning/topUsedEnzymes.m Not explicitly tested - The code has existed for a long time Not sure this works anymore - written for the old model structure?
kcat_tuning/truncateValues.m Not explicitly tested - The code has existed for a long time
limit_proteins/constrainProtConcs.m tc0013
Expand Down
122 changes: 122 additions & 0 deletions test/unit_tests/parityRound3Tests.m
View file
Open in desktop
Original file line number Diff line number Diff line change
Expand Up @@ -120,6 +120,128 @@ function testGetKcatAcrossIsozymesUsesPositiveSiblingsOnly_tc0252(testCase)
verifyTrue(testCase, isnan(out.ec.kcat(strcmp(out.ec.rxns, 'R2_EXP_2'))))
end

function testSigmaFitterReturnsStartOfGrowthPlateau_tc0255(testCase)
% The growth target is above what the uptake-limited model reaches, so growth
% is flat from sigma 0.06 on. GLPK returns the plateau growth with rounding
% noise (49.999999999999986 at sigma 0.06, 50 at sigma 1). The fitted sigma is
% the start of the plateau, 0.06, whatever the noise, and the same on Gurobi.
previous = getpref('RAVEN', 'solver', 'glpk');
restore = onCleanup(@() setRavenSolver(previous)); %#ok<NASGU>
setRavenSolver('glpk');
sigmaGlpk = plateauSigma();
verifyEqual(testCase, sigmaGlpk, 0.06, 'AbsTol', 1e-12)
hasGurobi = false;
if exist('gurobi', 'file')
try
setRavenSolver('gurobi');
hasGurobi = true;
catch
end
end
assumeTrue(testCase, hasGurobi, 'Gurobi is not available')
verifyEqual(testCase, plateauSigma(), sigmaGlpk, 'AbsTol', 1e-12)
end

function testGetStandardKcatSubSystemNamesIgnoreCase_tc0256(testCase)
% R3 is in 'Alpha', R5 in 'alpha', the entries of R2 in 'Beta'. The names
% 'Alpha' and 'alpha' are one subSystem, so R6 ('ALPHA') gets the median of
% R3 and R5 (5) and not the model-wide median (100).
[ecModel, adapter] = standardKcatModel(struct('R5', 'alpha', 'R2', 'Beta', 'R6', 'ALPHA'));
ecModel = getStandardKcat(ecModel, 'modelAdapter', adapter, 'threshold', 2);
verifyEqual(testCase, ecModel.ec.kcat(strcmp(ecModel.ec.rxns, 'R6')), 5)
end

function testGetStandardKcatReactionWithoutSubSystem_tc0257(testCase)
% R7 has an empty subSystem while the other reactions have one: it gets the
% model-wide median kcat (100), R6 the one of its subSystem 'Alpha' (5).
[ecModel, adapter] = standardKcatModel(struct('R2', 'Beta', 'R7', ''));
ecModel = getStandardKcat(ecModel, 'modelAdapter', adapter, 'threshold', 2);
verifyEqual(testCase, ecModel.ec.kcat(strcmp(ecModel.ec.rxns, 'R7')), 100)
verifyEqual(testCase, ecModel.ec.kcat(strcmp(ecModel.ec.rxns, 'R6')), 5)
end

function testGetSubsetEcModelKeepsEverythingForTheFullModel_tc0258(testCase)
% The context model is the conventional model, so no gene is removed and the
% ecModel is returned as it is, for the full and the light ecModel. The
% reaction identifiers of the light ecModel have a prefix; they are matched
% by the reaction they stand for, so they are not all removed.
adapter = getAdapter();
model = loadConventionalGEM('modelAdapter', adapter);
for light = [false true]
bigEcModel = makeEcModel(model, light, adapter);
smallEcModel = getSubsetEcModel(bigEcModel, model);
verifyEqual(testCase, smallEcModel, bigEcModel)
verifyNotEmpty(testCase, smallEcModel.ec.rxns)
end
end

function testGetSubsetEcModelLightDropsEcEntriesOfRemovedRxns_tc0259(testCase)
% The context model lacks R2: the light ecModel loses the entry 001_R2 of
% its ec structure and keeps the others.
adapter = getAdapter();
model = loadConventionalGEM('modelAdapter', adapter);
small = removeReactions(model, {'R2'}, 'removeUnusedMets', true, 'removeUnusedGenes', true);
bigEcModel = makeEcModel(model, true, adapter);
smallEcModel = getSubsetEcModel(bigEcModel, small);
verifyFalse(testCase, ismember('001_R2', smallEcModel.ec.rxns))
verifyTrue(testCase, all(ismember({'001_R3', '001_R5'}, smallEcModel.ec.rxns)))
verifyEqual(testCase, size(smallEcModel.ec.rxnEnzMat, 1), numel(smallEcModel.ec.rxns))
end

function [ecModel, adapter] = standardKcatModel(subSystems)
% ecTestGEM with R6 and R7 added (no genes) and the subSystems given as a
% struct of reaction id to name (an empty name is no subSystem); reactions
% not listed are in 'Alpha'. The kcats
% are R3 2, R5 8, the entries of R2 100.
adapter = getAdapter();
model = loadConventionalGEM('modelAdapter', adapter);
rxnsToAdd.rxns = {'R6'; 'R7'};
rxnsToAdd.equations = {'m1[c] => m2[c]'; 'm1[c] => m2[c]'};
rxnsToAdd.grRules = {''; ''};
model = addRxns(model, rxnsToAdd, 3);
model.subSystems = cell(numel(model.rxns), 1);
for i = 1:numel(model.rxns)
name = 'Alpha';
if isfield(subSystems, model.rxns{i})
name = subSystems.(model.rxns{i});
end
if isempty(name)
model.subSystems{i} = {};
else
model.subSystems{i} = {name};
end
end
ecModel = makeEcModel(model, false, adapter);
ecModel = getECfromGEM(ecModel);
ecModel.ec.kcat(:) = 100;
ecModel.ec.kcat(strcmp(ecModel.ec.rxns, 'R3')) = 2;
ecModel.ec.kcat(strcmp(ecModel.ec.rxns, 'R5')) = 8;
ecModel.ec.source(:) = {'manual'};
end

function sigma = plateauSigma()
adapter = getAdapter();
model = loadConventionalGEM('modelAdapter', adapter);
model.lb(strcmp(model.rxns, 'S1')) = -50;
for id = {'R2', 'R4'}
idx = strcmp(model.rxns, id{1});
model.lb(idx) = 0; model.ub(idx) = 0;
end
ecModel = makeEcModel(model, false, adapter);
ecModel = getECfromGEM(ecModel);
kcats = struct('R3', 5, 'R5', 200);
for id = fieldnames(kcats)'
j = strcmp(ecModel.ec.rxns, id{1});
ecModel.ec.kcat(j) = kcats.(id{1});
ecModel.ec.source{j} = 'manual';
end
ecModel.c = double(strcmp(ecModel.rxns, 'R5'));
ecModel = applyKcatConstraints(ecModel);
ecModel = setProtPoolSize(ecModel, [], [], [], adapter);
[~, sigma] = sigmaFitter(ecModel, 'growthRate', 100, 'Ptot', 0.5, 'f', 4, ...
'makePlot', false, 'modelAdapter', adapter);
end

function [ecModel, tuneAdapter] = tuningModel()
% R3 limits growth, its enzyme P4 is ignored by the tests, so tuning raises R5.
adapter = getAdapter();
Expand Down
Loading