From 3f00aaacaa7e066eca9790706a1a4cb1183dc5f1 Mon Sep 17 00:00:00 2001 From: Eduard Kerkhoven Date: Mon, 2 Nov 2026 23:42:17 +0100 Subject: [PATCH] getEnzymeUsage: demand ratio of the usage over the measured concentrations --- src/enzyme_usage/getEnzymeUsage.m | 31 +++++++++ src/enzyme_usage/reportEnzymeUsage.m | 14 +++++ test/VerificationMatrix.tsv | 4 +- test/unit_tests/geckoCoreFunctionTests.m | 80 ++++++++++++++++++++++++ 4 files changed, 127 insertions(+), 2 deletions(-) diff --git a/src/enzyme_usage/getEnzymeUsage.m b/src/enzyme_usage/getEnzymeUsage.m index ac09a085..21353e41 100644 --- a/src/enzyme_usage/getEnzymeUsage.m +++ b/src/enzyme_usage/getEnzymeUsage.m @@ -43,6 +43,25 @@ % usage_prot_ reaction in the model is left out, with the warning % getEnzymeUsage:missingUsageReaction. % - fluxes : vector of fluxes, copy of input fluxes. +% - measuredProtein : sum of the measured concentrations in model.ec.concs +% (mg/gDCW), over the enzymes with a non-NaN entry. +% - demandedProtein : sum of absUsage over the same enzymes (mg/gDCW). +% - demandRatio : demandedProtein / measuredProtein. NaN when no +% concentration is measured; Inf when the measured concentrations sum to 0 +% and the demand is positive. A ratio above 1 means that the solution uses +% more of the measured enzymes than was measured, which happens when the +% upper bounds of the usage reactions are above the measurements (for +% instance after flexibilizeEnzConcs). A ratio at or below 1 does not +% exclude that one enzyme exceeds its measurement, since the sum is over +% all measured enzymes; capUsage gives the ratio per enzyme. Enzymes +% without a measured concentration draw from the protein pool and are not +% counted. +% - measuredProteinActive, demandRatioActive : the same two values restricted +% to the enzymes with an absUsage above 1e-7 mg/gDCW. Measured enzymes +% that are inactive in the solution raise measuredProtein and can hide a +% conflict among the active ones. +% The five fields are computed over all enzymes with a usage reaction, +% whatever the value of zero. % % Examples % -------- @@ -61,6 +80,9 @@ error('This function does not work on GECKO light models.') end usageData.protID = ecModel.ec.enzymes; +concs = nan(numel(ecModel.ec.enzymes),1); +nConcs = min(numel(ecModel.ec.concs),numel(concs)); +concs(1:nConcs) = ecModel.ec.concs(1:nConcs); [~,rxnIdx] = ismember(strcat('usage_prot_',ecModel.ec.enzymes),ecModel.rxns); noUsageRxn = rxnIdx == 0; if any(noUsageRxn) @@ -68,6 +90,7 @@ 'The model has no usage reaction for the enzymes %s; they are left out.', ... strjoin(usageData.protID(noUsageRxn), ', ')) usageData.protID = usageData.protID(~noUsageRxn); + concs = concs(~noUsageRxn); rxnIdx = rxnIdx(~noUsageRxn); end @@ -76,6 +99,14 @@ usageData.capUsage = abs(usageData.absUsage./usageData.UB); usageData.fluxes = fluxes; +measured = ~isnan(concs); +active = usageData.absUsage > 1e-7; +usageData.measuredProtein = sum(concs(measured)); +usageData.measuredProteinActive = sum(concs(measured & active)); +usageData.demandedProtein = sum(usageData.absUsage(measured)); +usageData.demandRatio = usageData.demandedProtein / usageData.measuredProtein; +usageData.demandRatioActive = sum(usageData.absUsage(measured & active)) / usageData.measuredProteinActive; + if ~zero nonzero = usageData.absUsage>0; usageData.absUsage = usageData.absUsage(nonzero); diff --git a/src/enzyme_usage/reportEnzymeUsage.m b/src/enzyme_usage/reportEnzymeUsage.m index ff313280..4336f377 100644 --- a/src/enzyme_usage/reportEnzymeUsage.m +++ b/src/enzyme_usage/reportEnzymeUsage.m @@ -31,6 +31,14 @@ % the topAbsUsage enzymes with the highest absolute usage. % totalUsageFlux : double % total protein pool used, i.e. the upper bound of prot_pool_exchange. +% demandRatio : double +% demandRatio of usageData: the protein used of the enzymes with a +% measured concentration divided by their measured concentrations. +% Above 1, the solution needs more of those enzymes than was +% measured. NaN when usageData has no demandRatio field. +% demandRatioActive : double +% demandRatioActive of usageData: the same ratio over the enzymes +% that carry flux in the solution. % % Examples % -------- @@ -197,4 +205,10 @@ end usageReport.topAbsUsage = struct2table(topUsage); usageReport.totalUsageFlux = protPool; +usageReport.demandRatio = NaN; +usageReport.demandRatioActive = NaN; +if isfield(usageData, 'demandRatio') + usageReport.demandRatio = usageData.demandRatio; + usageReport.demandRatioActive = usageData.demandRatioActive; +end end diff --git a/test/VerificationMatrix.tsv b/test/VerificationMatrix.tsv index 6457127e..8f18ca58 100644 --- a/test/VerificationMatrix.tsv +++ b/test/VerificationMatrix.tsv @@ -11,8 +11,8 @@ 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 -enzyme_usage/getEnzymeUsage.m tc0507, tc0508 -enzyme_usage/reportEnzymeUsage.m tc0508, tc0509 +enzyme_usage/getEnzymeUsage.m tc0286, tc0507, tc0508 +enzyme_usage/reportEnzymeUsage.m tc0286, tc0508, tc0509 gather_kcats/fuzzyKcatMatching.m tc0010, tc0223, tc0224, tc0225, tc0226, tc0231, tc0239 gather_kcats/getStandardKcat.m tc0256, tc0257, tc0309, tc0310 gather_kcats/mergeDlkcatAndFuzzyKcats.m tc0011; tc0241 diff --git a/test/unit_tests/geckoCoreFunctionTests.m b/test/unit_tests/geckoCoreFunctionTests.m index e133e721..4a158f03 100644 --- a/test/unit_tests/geckoCoreFunctionTests.m +++ b/test/unit_tests/geckoCoreFunctionTests.m @@ -1291,6 +1291,86 @@ function testPfbaEnzymesFractionOfNegativeOptimum_tc0285(testCase) verifyGreaterThanOrEqual(testCase, solution.objectiveValue, -11) end +function testGetEnzymeUsageDemandRatio_tc0286(testCase) + % Verifies the demand ratio of getEnzymeUsage: the usage of the enzymes + % with a measured concentration divided by their measured + % concentrations, over all of them and over those that carry flux, NaN + % without a measured concentration, and the same whatever 'zero' is. + geckoPath = findGECKOroot; + adapter = ModelAdapterManager.getAdapter(fullfile(geckoPath,'test','unit_tests','ecTestGEM', 'TestGEMAdapter.m')); + model = getGeckoTestModel(); + ecModel = makeEcModel(model, false, adapter); + ecModel = getECfromGEM(ecModel); + ecModel.ec.kcat(:) = 100; + ecModel.ec.source(:) = {'manual'}; + ecModel = applyKcatConstraints(ecModel); + ecModel = setProtPoolSize(ecModel, [], [], [], adapter); + + [~, idx] = ismember(strcat('usage_prot_', ecModel.ec.enzymes), ecModel.rxns); + ecModel.ec.concs = nan(numel(ecModel.ec.enzymes), 1); + ecModel.ec.concs(1:3) = [2; 4; 1]; + fluxes = zeros(numel(ecModel.rxns), 1); + fluxes(idx(1)) = 3; + fluxes(idx(3)) = 1.5; + fluxes(idx(4)) = 7; + + usageData = getEnzymeUsage(ecModel, fluxes); + verifyEqual(testCase, usageData.measuredProtein, 7, 'AbsTol', 1e-12) + verifyEqual(testCase, usageData.demandedProtein, 4.5, 'AbsTol', 1e-12) + verifyEqual(testCase, usageData.demandRatio, 4.5/7, 'AbsTol', 1e-12) + verifyEqual(testCase, usageData.measuredProteinActive, 3, 'AbsTol', 1e-12) + verifyEqual(testCase, usageData.demandRatioActive, 1.5, 'AbsTol', 1e-12) + + withoutZero = getEnzymeUsage(ecModel, fluxes, 'zero', false); + verifyEqual(testCase, withoutZero.demandRatio, usageData.demandRatio, 'AbsTol', 1e-12) + verifyEqual(testCase, withoutZero.demandRatioActive, usageData.demandRatioActive, 'AbsTol', 1e-12) + + report = reportEnzymeUsage(ecModel, usageData); + verifyEqual(testCase, report.demandRatio, usageData.demandRatio, 'AbsTol', 1e-12) + verifyEqual(testCase, report.demandRatioActive, usageData.demandRatioActive, 'AbsTol', 1e-12) + + % A usage below the threshold counts as demand but not as active. + lowFluxes = fluxes; + lowFluxes(idx(2)) = 1e-8; + usageData = getEnzymeUsage(ecModel, lowFluxes); + verifyEqual(testCase, usageData.demandedProtein, 4.5 + 1e-8, 'AbsTol', 1e-12) + verifyEqual(testCase, usageData.measuredProteinActive, 3, 'AbsTol', 1e-12) + verifyEqual(testCase, usageData.demandRatioActive, 1.5, 'AbsTol', 1e-12) + + % An enzyme without a usage reaction is left out with its concentration. + reduced = removeReactions(ecModel, ecModel.rxns(idx(1))); + reducedFluxes = fluxes; + reducedFluxes(idx(1)) = []; + warning('off', 'getEnzymeUsage:missingUsageReaction') + restore = onCleanup(@() warning('on', 'getEnzymeUsage:missingUsageReaction')); + usageData = getEnzymeUsage(reduced, reducedFluxes); + verifyEqual(testCase, usageData.measuredProtein, 5, 'AbsTol', 1e-12) + verifyEqual(testCase, usageData.demandedProtein, 1.5, 'AbsTol', 1e-12) + verifyEqual(testCase, usageData.demandRatioActive, 1.5, 'AbsTol', 1e-12) + + % A usageData without the demand fields gives NaN in the report. + bare = rmfield(getEnzymeUsage(ecModel, fluxes), ... + {'measuredProtein','measuredProteinActive','demandedProtein','demandRatio','demandRatioActive'}); + report = reportEnzymeUsage(ecModel, bare); + verifyTrue(testCase, isnan(report.demandRatio)) + verifyTrue(testCase, isnan(report.demandRatioActive)) + + % A concentration vector shorter than the enzyme list counts as unmeasured. + short = ecModel; + short.ec.concs = [2; 4]; + usageData = getEnzymeUsage(short, fluxes); + verifyEqual(testCase, usageData.measuredProtein, 6, 'AbsTol', 1e-12) + short.ec.concs = []; + usageData = getEnzymeUsage(short, fluxes); + verifyTrue(testCase, isnan(usageData.demandRatio)) + + ecModel.ec.concs(:) = NaN; + usageData = getEnzymeUsage(ecModel, fluxes); + verifyEqual(testCase, usageData.measuredProtein, 0) + verifyTrue(testCase, isnan(usageData.demandRatio)) + verifyTrue(testCase, isnan(usageData.demandRatioActive)) +end + function ecModel = fixtureForRelaxProteomicsGreedy() % Builds the shared fixture used by the relaxProteomicsGreedy tests below: