feat: check charge balance in getElementalBalance

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)

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
 

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

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

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