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');