Skip to content

test: port invariants asserted by the raven-toolbox suite #683

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 1 commit into from
Jul 17, 2026
Merged

test: port invariants asserted by the raven-toolbox suite #683

Show file tree
Hide file tree
Changes from all commits
Commits
Show all changes
1 commit
Select commit
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
76 changes: 76 additions & 0 deletions testing/function_tests/tConditions.m
View file
Open in desktop
Original file line number Diff line number Diff line change
Expand Up @@ -12,5 +12,81 @@ function applyConditionSetsBounds(testCase)
testCase.verifyEqual(m2.lb(idx), 0);
end

function applyConditionResetExchangesZeroesUptake(testCase)
% The prelude reopens exchanges to (0, 1000), i.e. no uptake.
m = testCase.model;
[~, exchangeRxns] = getExchangeRxns(m);
m.lb(exchangeRxns) = -37;
m.ub(exchangeRxns) = 42;
condition.prelude.reset_exchanges = 'both';
evalc('m2 = applyCondition(m, condition);');
testCase.verifyEqual(unique(m2.lb(exchangeRxns)), 0);
testCase.verifyEqual(unique(m2.ub(exchangeRxns)), 1000);
end

function applyConditionRemoveMetZeroesCoefficient(testCase)
% remove_mets drops a metabolite from the cofactor pseudoreaction.
m = tConditions.cofactorModel();
condition.cofactor_pseudoreaction = struct('rxn_id', 'R_cofactor', ...
'remove_mets', {{struct('met', 'nadh_c')}});
evalc('m2 = applyCondition(m, condition);');
testCase.verifyEqual(full(m2.S(strcmp(m2.mets,'nadh_c'), 1)), 0);
end

function applyConditionRecomputesChargeBalance(testCase)
% charge_balance_met takes whatever coefficient makes the
% pseudoreaction charge neutral, recomputed after the removals.
m = tConditions.cofactorModel();
condition.cofactor_pseudoreaction = struct('rxn_id', 'R_cofactor', ...
'remove_mets', {{struct('met', 'nadh_c')}}, ...
'charge_balance_met', 'h_c');
evalc('m2 = applyCondition(m, condition);');
% Left after removing NADH: atp_c (-1, charge -4) and amp_c
% (+1, charge -2), a residual of +2, so H+ must enter at -2.
testCase.verifyEqual(full(m2.S(strcmp(m2.mets,'h_c'), 1)), -2);
% which is to say the reaction ends up charge neutral
idx = find(m2.S(:,1));
testCase.verifyEqual(sum(full(m2.S(idx,1)) .* m2.metCharges(idx)), 0);
end

function applyConditionBiomassDeltaAddsToCoefficient(testCase)
% A stoichiometry delta adds to the existing coefficient rather
% than replacing it.
m = tConditions.cofactorModel();
m.rxns{2} = 'R_biomass'; m.rxnNames{2} = 'R_biomass';
before = full(m.S(strcmp(m.mets,'atp_c'), 2));
condition.biomass_stoichiometry_delta = struct('rxn_id', 'R_biomass', ...
'add', {{struct('met', 'atp_c', 'coef', -0.5)}});
evalc('m2 = applyCondition(m, condition);');
testCase.verifyEqual(full(m2.S(strcmp(m2.mets,'atp_c'), 2)), before - 0.5);
end

end

methods (Static, Access = private)

function m = cofactorModel()
% R_cofactor consumes ATP and NADH; R_other consumes ATP.
m = struct();
m.id = 'cofac';
m.rxns = {'R_cofactor';'R_other'};
m.rxnNames = {'R_cofactor';'R_other'};
m.mets = {'atp_c';'nadh_c';'amp_c';'h_c'};
m.metNames = {'ATP';'NADH';'AMP';'H+'};
m.metComps = [1;1;1;1];
m.metCharges= [-4;-2;-2;1];
m.comps = {'c'};
m.compNames = {'cytosol'};
% R_cofactor R_other
m.S = sparse([ -1 -1; % atp_c
-1 0; % nadh_c
1 0; % amp_c
0 0]); % h_c
m.lb = [0;0]; m.ub = [1000;1000]; m.rev = [0;0]; m.c = [0;0];
m.b = zeros(4,1);
m.genes = {}; m.grRules = {'';''}; m.rxnGeneMat = sparse(2,0);
m.metFormulas = {'C10';'C21';'C10';'H'};
end

end
end
118 changes: 118 additions & 0 deletions testing/function_tests/tManipulation.m
View file
Open in desktop
Original file line number Diff line number Diff line change
Expand Up @@ -113,6 +113,86 @@ function convertToIrrevAllIrreversible(testCase)
testCase.verifyGreaterThanOrEqual(numel(m2.rxns), numel(testCase.model.rxns));
end

function convertToIrrevSplitsBoundsAndStoichiometry(testCase)
% A reversible reaction with bounds (-500,1000) keeps (0,1000) and
% its stoichiometry, while the _REV copy gets (0,500) and the
% negated stoichiometry.
m = tManipulation.twoMetModel();
m.rev = 1; m.lb = -500; m.ub = 1000;
m2 = convertToIrrev(m);

fwd = strcmp(m2.rxns, 'R1');
rev = strcmp(m2.rxns, 'R1_REV');
testCase.verifyTrue(any(rev));
testCase.verifyEqual(full(m2.lb(fwd)), 0);
testCase.verifyEqual(full(m2.ub(fwd)), 1000);
testCase.verifyEqual(full(m2.lb(rev)), 0);
testCase.verifyEqual(full(m2.ub(rev)), 500);
testCase.verifyEqual(full(m2.S(:,fwd)), [-1;1]);
testCase.verifyEqual(full(m2.S(:,rev)), [1;-1]);
end

function convertToIrrevReverseCopyInheritsGrRule(testCase)
% The _REV copy is catalysed by the same genes as the forward one.
m = tManipulation.twoMetModel();
m.rev = 1; m.lb = -500; m.ub = 1000;
m.genes = {'g1'}; m.grRules = {'g1'}; m.rxnGeneMat = sparse(1,1,1);
m2 = convertToIrrev(m);
testCase.verifyEqual(m2.grRules{strcmp(m2.rxns,'R1_REV')}, 'g1');
end

function findDuplicateRxnsIgnoreDirection(testCase)
% a -> b and b -> a are the same reaction run backwards, so they
% group by default and stay separate when direction matters.
m = tManipulation.twoMetModel();
m.rxns = {'R1';'R2'};
m.rxnNames = {'R1';'R2'};
m.S = sparse([-1 1; 1 -1]);
m.lb = [0;0]; m.ub = [1000;1000]; m.rev = [0;0]; m.c = [0;0];
m.grRules = {'';''}; m.rxnGeneMat = sparse(2,0);

pairs = findDuplicateRxns(m);
testCase.verifyEqual(pairs, [1 2]);

pairs = findDuplicateRxns(m, 'ignoreDirection', false);
testCase.verifyEmpty(pairs);
end

function changeGrRulesAppendsToExistingRule(testCase)
% replace=false must OR the new rule onto the existing one rather
% than overwrite it, and add the new gene to the model.
m2 = changeGrRules(testCase.model, 'ACKr', 'b9999', false);
rule = m2.grRules{strcmp(m2.rxns, 'ACKr')};
testCase.verifySubstring(rule, 'b9999');
testCase.verifySubstring(rule, ' or ');
% the gene the reaction already had must still be in the rule
oldRule = testCase.model.grRules{strcmp(testCase.model.rxns, 'ACKr')};
oldGene = regexp(oldRule, 'b\d+', 'match', 'once');
testCase.verifySubstring(rule, oldGene);
testCase.verifyTrue(ismember('b9999', m2.genes));
end

function copyToCompsDeleteOriginalIsAMove(testCase)
% deleteOriginal turns the copy into a move: the reaction count is
% unchanged and the original compartment's copy is gone.
evalc('copied = copyToComps(testCase.model, {''p''}, ''rxns'', ''ACKr'');');
evalc(['moved = copyToComps(testCase.model, {''p''}, ''rxns'', ''ACKr'', ' ...
'''deleteOriginal'', true);']);
testCase.verifyEqual(numel(copied.rxns), numel(testCase.model.rxns) + 1);
testCase.verifyEqual(numel(moved.rxns), numel(testCase.model.rxns));
end

function mergeModelsMetParamDecidesUnification(testCase)
% The same metabolite under two ids unifies when matching on names
% and stays distinct when matching on ids.
a = tManipulation.namedMetModel('glc_c', 'A');
b = tManipulation.namedMetModel('glucose_c', 'B');
evalc('byName = mergeModels({a; b});');
evalc('byId = mergeModels({a; b}, ''metParam'', ''mets'');');
testCase.verifyEqual(nnz(strcmp(byName.metNames, 'Glucose')), 1);
testCase.verifyEqual(nnz(strcmp(byId.metNames, 'Glucose')), 2);
end

function copyToCompsAddsCompartment(testCase)
evalc('m2 = copyToComps(testCase.model, {''p''}, ''ACKr'');');
testCase.verifyTrue(ismember('p', m2.comps));
Expand Down Expand Up @@ -252,4 +332,42 @@ function standardizeGrRulesReturnsRules(testCase)
end

end

methods (Static, Access = private)

function m = twoMetModel()
% Single reaction a -> b in one compartment.
m = struct();
m.id = 'toy';
m.rxns = {'R1'};
m.rxnNames = {'R1'};
m.mets = {'a';'b'};
m.metNames = {'a';'b'};
m.metComps = [1;1];
m.comps = {'c'};
m.compNames = {'cytosol'};
m.S = sparse([-1;1]);
m.lb = 0; m.ub = 1000; m.rev = 0; m.c = 0; m.b = zeros(2,1);
m.genes = {}; m.grRules = {''}; m.rxnGeneMat = sparse(1,0);
m.metFormulas = {'C';'C'};
end

function m = namedMetModel(glucoseId, modelId)
% Glucose[c] under a caller-chosen id, consumed by one reaction.
m = struct();
m.id = modelId;
m.rxns = {['R_' modelId]};
m.rxnNames = m.rxns;
m.mets = {glucoseId; ['product_' modelId]};
m.metNames = {'Glucose'; ['Product' modelId]};
m.metComps = [1;1];
m.comps = {'c'};
m.compNames = {'cytosol'};
m.S = sparse([-1;1]);
m.lb = 0; m.ub = 1000; m.rev = 0; m.c = 0; m.b = zeros(2,1);
m.genes = {}; m.grRules = {''}; m.rxnGeneMat = sparse(1,0);
m.metFormulas = {'C6H12O6';'C'};
end

end
end