From d6ab2f4bc421ccebc1930ae0f1ca5e0f88614c5a Mon Sep 17 00:00:00 2001 From: Eduard Kerkhoven Date: Fri, 17 Jul 2026 13:10:58 +0200 Subject: [PATCH] feat: check charge balance in getElementalBalance RAVEN checked charge nowhere: getElementalBalance is elemental only and parseFormulas strips +/- from formulas. getElementalBalance now also reports chargeStatus (1 balanced, 0 unbalanced, -1 unknown) and chargeResidual, as sibling fields rather than folded into balanceStatus, which removeBadRxns and printModelStats read as an elemental verdict. Both are masked to the metabolites that participate in each reaction. S is sparse and 0*NaN is NaN, so an unmasked sum would let a single unset charge anywhere in the model poison every reaction's residual (measured in R2024b). Meanwhile RAVEN had two functions that rebalance charge, and both summed with 'omitnan': model.S(Hc, rxnPos) = -sum(model.S(:, rxnPos) .* model.metCharges, 'omitnan'); addMets defaults metCharges to NaN, so 'omitnan' treated an unknown charge as neutral and silently wrote the wrong coefficient onto the balancing metabolite. Both now refuse, naming the metabolites whose charge is unset: an unknown charge balance is not a zero one. --- biomass/scaleBiomassPseudoreaction.m | 24 ++++++++++- conditions/applyCondition.m | 18 +++++++- queries/getElementalBalance.m | 38 ++++++++++++++++- testing/function_tests/tBiomass.m | 47 +++++++++++++++++++++ testing/function_tests/tQueries.m | 61 ++++++++++++++++++++++++++++ 5 files changed, 185 insertions(+), 3 deletions(-) diff --git a/biomass/scaleBiomassPseudoreaction.m b/biomass/scaleBiomassPseudoreaction.m index f565800b..4209bc89 100644 --- a/biomass/scaleBiomassPseudoreaction.m +++ b/biomass/scaleBiomassPseudoreaction.m @@ -57,7 +57,29 @@ % Rebalance H+ to keep charge neutrality. Hc = find(strcmp(model.mets, biomassConfig.proton_met)); model.S(Hc, rxnPos) = 0; -model.S(Hc, rxnPos) = -sum(model.S(:, rxnPos) .* model.metCharges, 'omitnan'); +model.S(Hc, rxnPos) = -chargeResidual(model, rxnPos); +end + +function residual = chargeResidual(model, rxnPos) +% Sum of charges over the metabolites that take part in a reaction. Only the +% participants may be summed: S is sparse and 0*NaN is NaN, so any unset +% charge in the model would otherwise poison the result. Summing with +% 'omitnan' instead would treat an unset charge as neutral and silently +% rebalance the reaction to the wrong coefficient. +metIdx = find(model.S(:, rxnPos)); +if ~isfield(model, 'metCharges') + error('scaleBiomassPseudoreaction:noCharges', ... + ['Cannot rebalance charge: the model has no metCharges field. ' ... + 'Provide metabolite charges, or rescale without charge balancing.']); +end +unknown = metIdx(isnan(model.metCharges(metIdx))); +if ~isempty(unknown) + error('scaleBiomassPseudoreaction:unknownCharge', ... + ['Cannot rebalance charge of %s: the following metabolites have no ' ... + 'charge, so the charge balance is unknown rather than zero:\n\t%s'], ... + model.rxns{rxnPos}, strjoin(model.mets(unknown), ', ')); +end +residual = sum(full(model.S(metIdx, rxnPos)) .* model.metCharges(metIdx)); end function comp = findComponent(cfg, name) diff --git a/conditions/applyCondition.m b/conditions/applyCondition.m index 940c5d2f..c3f1f58b 100644 --- a/conditions/applyCondition.m +++ b/conditions/applyCondition.m @@ -88,8 +88,24 @@ if isfield(cp, 'charge_balance_met') balanceIdx = find(strcmp(model.mets, cp.charge_balance_met)); model.S(balanceIdx, cofacIdx) = 0; + % Only the participating metabolites may be summed: S is sparse and + % 0*NaN is NaN, while 'omitnan' would treat an unset charge as + % neutral and silently write the wrong balancing coefficient. + metIdx = find(model.S(:, cofacIdx)); + if ~isfield(model, 'metCharges') + error('applyCondition:noCharges', ... + ['Cannot charge balance %s: the model has no metCharges ' ... + 'field.'], cp.rxn_id); + end + unknown = metIdx(isnan(model.metCharges(metIdx))); + if ~isempty(unknown) + error('applyCondition:unknownCharge', ... + ['Cannot charge balance %s: the following metabolites have ' ... + 'no charge, so the charge balance is unknown rather than ' ... + 'zero:\n\t%s'], cp.rxn_id, strjoin(model.mets(unknown), ', ')); + end model.S(balanceIdx, cofacIdx) = ... - -sum(model.S(:, cofacIdx) .* model.metCharges, 'omitnan'); + -sum(full(model.S(metIdx, cofacIdx)) .* model.metCharges(metIdx)); end end diff --git a/queries/getElementalBalance.m b/queries/getElementalBalance.m index c3a7b4a7..f7b88f8b 100755 --- a/queries/getElementalBalance.m +++ b/queries/getElementalBalance.m @@ -28,7 +28,14 @@ % % - balanceStatus : 1 if the reaction is balanced, 0 if it is % unbalanced, -1 if it could not be balanced due to missing -% information, -2 if it could not be balanced due to an error +% information, -2 if it could not be balanced due to an error. +% Elemental only; charge is reported separately below +% - chargeStatus : 1 if the reaction is charge balanced, 0 if it is +% not, -1 if any participating metabolite has no charge (or the +% model has no metCharges field), in which case the charge balance +% is unknown rather than zero +% - chargeResidual : the sum of charges over the reaction, NaN where +% chargeStatus is -1 % - elements : struct with fields abbrevs (cell array with % abbreviations for all used elements) and names (cell array with % the names for all used elements) @@ -115,6 +122,33 @@ %The remaining ones are all balanced balanceStructure.balanceStatus(isnan(balanceStructure.balanceStatus))=1; +%Charge balance. This is reported separately from balanceStatus, which +%callers such as removeBadRxns and printModelStats read as a purely +%elemental verdict. +balanceStructure.chargeStatus=zeros(numel(model.rxns),1); +balanceStructure.chargeResidual=nan(numel(model.rxns),1); +if ~isfield(model,'metCharges') + balanceStructure.chargeStatus(:)=-1; +else + for j=1:numel(model.rxns) + %Only the participating metabolites may be touched: S is sparse and + %0*NaN is NaN, so a single unset charge anywhere in the model would + %otherwise poison every reaction. Summing with 'omitnan' instead + %would be worse, silently reporting an unknown residual as 0. + idx=find(model.S(:,j)); + if isempty(idx) || any(isnan(model.metCharges(idx))) + balanceStructure.chargeStatus(j)=-1; + else + balanceStructure.chargeResidual(j)=sum(full(model.S(idx,j)).*model.metCharges(idx)); + if abs(balanceStructure.chargeResidual(j))>10^-8 %Roundoff error + balanceStructure.chargeStatus(j)=0; + else + balanceStructure.chargeStatus(j)=1; + end + end + end +end + %Print warnings toPrint=[]; if printUnbalanced==true @@ -150,6 +184,8 @@ rxns = getIndexes(model,rxns,'rxns'); [~,i] = sort(rxns); balanceStructure.balanceStatus(i) = balanceStructure.balanceStatus; + balanceStructure.chargeStatus(i) = balanceStructure.chargeStatus; + balanceStructure.chargeResidual(i) = balanceStructure.chargeResidual; balanceStructure.leftComp(i,:) = balanceStructure.leftComp; balanceStructure.rightComp(i,:) = balanceStructure.rightComp; end diff --git a/testing/function_tests/tBiomass.m b/testing/function_tests/tBiomass.m index b13ae3f2..535de476 100644 --- a/testing/function_tests/tBiomass.m +++ b/testing/function_tests/tBiomass.m @@ -33,6 +33,26 @@ function scaleBiomassPseudoreactionErrorsOnMissingComponent(testCase) 'protein', 0.9), ?MException); end + function scaleBiomassPseudoreactionRebalancesProton(testCase) + % Rescaling the substrates rebalances H+ so the pseudoreaction + % stays charge neutral. + [m, cfg] = tBiomass.protonToyModel(); + out = scaleBiomassPseudoreaction(m, cfg, 'protein', 0.5); + % x_c: -0.5 * charge -2 = +1, so H+ must come in at -1 + testCase.verifyEqual(full(out.S(strcmp(out.mets,'x_c'), 1)), -0.5, 'AbsTol', 1e-9); + testCase.verifyEqual(full(out.S(strcmp(out.mets,'h_c'), 1)), -1, 'AbsTol', 1e-9); + end + + function scaleBiomassPseudoreactionRefusesUnknownCharge(testCase) + % An unset charge on a participating metabolite leaves the charge + % balance unknown. Treating it as neutral would silently write the + % wrong proton coefficient, so the rescale must refuse instead. + [m, cfg] = tBiomass.protonToyModel(); + m.metCharges(strcmp(m.mets,'x_c')) = NaN; + testCase.verifyError(@() scaleBiomassPseudoreaction(m, cfg, 'protein', 0.5), ... + 'scaleBiomassPseudoreaction:unknownCharge'); + end + function fitParametersRunsWhenQuadprogAvailable(testCase) testCase.assumeDependency(exist('quadprog','file')==2, ... 'Optimization Toolbox (quadprog)'); @@ -46,6 +66,33 @@ function fitParametersRunsWhenQuadprogAvailable(testCase) end + methods (Static, Access = private) + function [m, cfg] = protonToyModel() + % One pseudoreaction: x_c -> protein, with h_c available to carry + % the charge balance. + m = struct(); + m.id = 'toy'; + m.rxns = {'R1'}; + m.rxnNames = {'protein pseudoreaction'}; + m.mets = {'x_c';'h_c';'prot_c'}; + m.metNames = {'x';'h';'protein'}; + m.metComps = [1;1;1]; + m.metCharges= [-2;1;0]; + m.comps = {'c'}; + m.compNames = {'cytosol'}; + m.S = sparse([-1;0;1]); + m.lb = 0; m.ub = 1000; m.rev = 0; m.c = 1; m.b = zeros(3,1); + m.grRules = {''}; m.genes = {}; m.rxnGeneMat = sparse(1,0); + + cfg = struct(); + cfg.biomass_rxn = 'R1'; + cfg.proton_met = 'h_c'; + cfg.components = {struct('name','protein', ... + 'pseudoreaction_name','protein pseudoreaction', ... + 'mass_strategy','mw')}; + end + end + methods (Access = private) function cfg = biomassConfig(testCase) biomassRxn = testCase.model.rxns{find(testCase.model.c == 1, 1)}; diff --git a/testing/function_tests/tQueries.m b/testing/function_tests/tQueries.m index d9800601..72ded278 100644 --- a/testing/function_tests/tQueries.m +++ b/testing/function_tests/tQueries.m @@ -79,6 +79,67 @@ function getElementalBalanceEmptyRxnIsUnbalanced(testCase) testCase.verifyLessThanOrEqual(bs.balanceStatus(end), -1); end + function getElementalBalanceReportsChargeBalance(testCase) + % A charge-imbalanced reaction is reported as such, while the + % elemental verdict is left alone. + m = testCase.model; + m.mets = {'a';'b'}; + m.metNames = {'a';'b'}; + m.metComps = [1;1]; + m.metFormulas = {'H';'H'}; + m.metCharges = [0;1]; + m.comps = {'c'}; + m.compNames = {'cytosol'}; + m.rxns = {'R1'}; + m.rxnNames = {'R1'}; + m.S = sparse([-1;1]); + m.lb = 0; m.ub = 1000; m.rev = 0; m.c = 0; m.b = zeros(2,1); + m.grRules = {''}; m.genes = {}; m.rxnGeneMat = sparse(1,0); + + bs = getElementalBalance(m); + testCase.verifyEqual(bs.chargeStatus(1), 0); + testCase.verifyEqual(bs.chargeResidual(1), 1, 'AbsTol', 1e-9); + % a -> b is elementally balanced (H on both sides) + testCase.verifyEqual(bs.balanceStatus(1), 1); + + % Balancing the charge flips chargeStatus, not balanceStatus + m.metCharges = [1;1]; + bs = getElementalBalance(m); + testCase.verifyEqual(bs.chargeStatus(1), 1); + testCase.verifyEqual(bs.chargeResidual(1), 0, 'AbsTol', 1e-9); + end + + function getElementalBalanceChargeUnknownIsNotZero(testCase) + % An unset charge on a participating metabolite makes the charge + % balance unknown, and must not be reported as balanced. An unset + % charge on a metabolite that does not participate must not leak + % into the reaction's residual. + m = testCase.model; + m.mets = {'a';'b';'spectator'}; + m.metNames = {'a';'b';'spectator'}; + m.metComps = [1;1;1]; + m.metFormulas = {'H';'H';'H'}; + m.metCharges = [1;1;NaN]; + m.comps = {'c'}; + m.compNames = {'cytosol'}; + m.rxns = {'R1'}; + m.rxnNames = {'R1'}; + m.S = sparse([-1;1;0]); + m.lb = 0; m.ub = 1000; m.rev = 0; m.c = 0; m.b = zeros(3,1); + m.grRules = {''}; m.genes = {}; m.rxnGeneMat = sparse(1,0); + + % The NaN belongs to a metabolite outside the reaction + bs = getElementalBalance(m); + testCase.verifyEqual(bs.chargeStatus(1), 1); + testCase.verifyEqual(bs.chargeResidual(1), 0, 'AbsTol', 1e-9); + + % Now the NaN is on a participant: unknown, not balanced + m.metCharges = [1;NaN;0]; + bs = getElementalBalance(m); + testCase.verifyEqual(bs.chargeStatus(1), -1); + testCase.verifyTrue(isnan(bs.chargeResidual(1))); + end + function getExchangeRxnsConsistent(testCase) [exch, idx] = getExchangeRxns(testCase.model); testCase.verifyClass(exch, 'cell');