stoichConsistency¶

checkStoichiometricConsistency(model, printLevel, method)¶

Verification of stoichiometric consistency by checking for at least one strictly positive basis in the left nullspace of S. If S is not stoichiometrically consistent detect conserved and unconserved metabolites, by returning a maximal conservation vector, which is a non-negative basis with as many strictly positive entries as possible. This omits rows of S that are entirely zero, when any exchange reactions are removed. The strictly positive and zero entries in m correspond to the conserved and unconserved metabolites respectively.

Verification of stoichiometric consistency is as initially described in: A. Gevorgyan, M. G. Poolman, and D. A. Fell Detection of stoichiometric inconsistencies in biomolecular models. Bioinformatics, 24(19):2245–2251, 2008.

Detection of conserved and unconserved metabolites based on a new implementation by Nikos Vlassis & Ronan Fleming.

USAGE:

[isConsistent, m, model] = checkStoichiometricConsistency(model, printLevel, method)

INPUT:

model: structure with fields:

model.S - m x n stoichiometric matrix

model.mets - if exists, but SIntRxnBool does not exist, then `findSExRxnBool will attempt to identify exchange reactions. (optional)

model.SIntMetBool - m x 1 boolean vector indicating the metabolites that are thought to be exclusively involved in non-exchange reactions (optional)

model.SIntRxnBool n x 1 boolean vector indicating the reactions that are thought to be non-exchange reactions. If this is not present, then we will attempt to identify such reactions, model.rxns exists, otherwise, we will assume all reactions are supposed to be non-exchange reactions. (optional)

OPTIONAL INPUTS:

printLevel: {(0), 1} method: structure with fields:

method.interface - {(‘SDCCO’),’LP’, ‘MILP’, ‘DCCO’} interface called to do the consistency check

method.solver - {(default solver as specified by CBT_LP_SOLVER), or any other CBT compatible LP solver}

method.param - solver specific parameter structure

OUTPUTS:

isConsistent: Solver status in standardized form:

1 - Optimal solution (Stoichiometrically consistent)

2 - Unbounded solution (Should never happen)

0 - Infeasible (Stoichiometrically INconsistent)

-1 - No solution reported (timelimit, numerical problem etc)

m: m x 1 strictly positive vector in left nullspace: (empty if it does not exist)

model: structure with fields:

.SConsistentMetBool - m x 1 boolean vector indicating metabolites involved in the maximal consistent vector

.SConsistentRxnBool - n x 1 boolean vector non-exchange reaction involving a stoichiometrically consistent metabolite

findMassLeaksAndSiphons(model, metBool, rxnBool, modelBoundsFlag, param, printLevel)¶

Finds the metabolites in a network that either leak mass or act as a siphon for mass, with (default) or without the bounds on a model. The approach is to solve the problem:

\[\begin{split}max ~& ||y||_0 \\ s.t. ~& Sv - y = 0\end{split}\]

with either \(l \leq v \leq u\) or \(-\infty \leq v \leq \infty\)

and with either \(0 \leq y \leq \infty\) (semipositive net stoichiometry = leak) or \(-\infty \leq y \leq 0\) (seminegative net stoichiometry = siphon)

If there are any zero rows of S, then the corresponding entry in y is set to zero.

USAGE:

[leakMetBool, leakRxnBool, siphonMetBool, siphonRxnBool, leakY, siphonY, statp, statn] = findMassLeaksAndSiphons(model, metBool, rxnBool, modelBoundsFlag, param, printLevel)

INPUT:

model: structure with fields (only .S is mandatory)

.S - m x n stoichiometric matrix

.lb - Lower bounds

.ub - Upper bounds

.SConsistentMetBool - m x 1 boolean vector indicating consistent mets

.SConsistentRxnBool - m x 1 boolean vector indicating consistent rxns

OPTIONAL INPUTS:

metBool: m x 1 boolean vector of metabolites to test for leakage rxnBool: n x 1 boolean vector of reactions to test for leakage modelBoundsFlag: {0, (1)}

0 = set all reaction bounds to -inf, inf

1 = use reaction bounds provided by model.lb and .ub

param: structure with fields:

param.epsilon - (getCobraSolverParams(‘LP’, ‘feasTol’)*100)

param.eta - (feasTol*100), smallest nonzero mass leak/siphon

param.theta - (0.5) parameter of capped l1 approximation

param.method - {(‘quasiConcave’), ‘dc’} method of approximation

printLevel: {(0), 1, 2 = debug}

OUTPUTS:

leakRxnBool: m x 1 boolean of metabolites in a positive leakage mode leakRxnBool: n x 1 boolean of reactions exclusively involved in a positive leakage mode siphonMetBool: m x 1 boolean of metabolites in a negative leakage mode siphonRxnBool: n x 1 boolean of reactions exclusively involved in a negative leakage mode leakY: m x 1 boolean of metabolites in a positive leakage mode siphonY: m x 1 boolean of metabolites in a negative siphon mode statp: status (positive leakage modes)

1 = Solution found

2 = Unbounded

0 = Infeasible

-1 = Invalid input

statn: status (negative leakage modes)

1 = Solution found

2 = Unbounded

0 = Infeasible

-1 = Invalid input

Fleming RMT, Haraldsdottir HS, Le HM, Vuong PT, Hankemeier T, Thiele I. Cardinality optimisation in constraint-based modelling: Application to human metabolism, 2022 (submitted).

findMinimalLeakageMode(model, metBool, modelBoundsFlag, epsilon, printLevel)¶

Solve the problem

\[\begin{split}min ~& ||v||_0 + ||y||_0 \\ s.t. ~& Sv - y = 0, \\ ~& l \leq v \leq u\end{split}\]

with either \(0 \leq y\) (semipositive net stoichiometry) or \(y \leq 0\) (seminegative net stoichiometry)

USAGE:

[Vp, Yp, statp, Vn, Yn, statn] = findMinimalLeakageMode(model, metBool, modelBoundsFlag, epsilon, printLevel)

INPUT:

model: (the following fields are required: .S, .lb, .ub)

.S - m x n stoichiometric matrix

.lb - Lower bounds

.ub - Upper bounds

.SConsistentMetBool - m x 1 boolean vector indicating consistent mets

.SConsistentRxnBool - m x 1 boolean vector indicating consistent rxns

OPTIONAL INPUTS:

metBool: m x 1 boolean vector of metabolites to test for leakage modelBoundsFlag: {0, (1)} where:

0 = set all reaction bounds to -inf, inf

1 = use reaction bounds provided by model.lb and .ub

epsilon: getCobraSolverParams(‘LP’, ‘feasTol’)*100, smallest nonzero reaction flux in leakage mode printLevel: {(0), 1}

OUTPUTS

Vp: n x 1 vector (positive leakage modes) Yp: m x 1 vector (positive leakage modes) statp: status (positive leakage modes)

1 = Solution found

2 = Unbounded

0 = Infeasible

-1 = Invalid input

Vn: n x 1 vector (negative leakage modes) Yn: m x 1 vector (negative leakage modes) statn: status (negative leakage modes)

1 = Solution found

2 = Unbounded

0 = Infeasible

-1 = Invalid input

findMinimalLeakageModeMet(model, metBool, rxnBool, modelBoundsFlag, params, printLevel)¶

Finds a minimal set of leak (or siphon) metabolites and the corresponding minimal set of reactions involved. Test metabolites in metBool.

Solve the problem

min ~& ||v||_0 + ||y||_0 \ s.t. ~& Sv - y = 0, \

& l leq v leq u

with either \(l(rxnBool) > 0\) or \(u(rxnBool) < 0\) and either \(0 \leq y\) (semipositive net stoichiometry) or \(y \leq 0\) (seminegative net stoichiometry)

and \(1 \leq y(metBool)\) (semipositive net stoichiometry) or \(y(metBool) \leq 1\) (seminegative net stoichiometry)

USAGE:

[minLeakMetBool, minLeakRxnBool, minSiphonMetBool, minSiphonRxnBool, leakY, siphonY, statp, statn] = findMinimalLeakageModeMet(model, metBool, rxnBool, modelBoundsFlag, params, printLevel)

INPUT:

model: structure with fields (bools are not mandatory)

.S - m x n stoichiometric matrix

.lb - Lower bounds

.ub - Upper bounds

.SConsistentMetBool - m x 1 boolean vector indicating consistent mets

.SConsistentRxnBool - m x 1 boolean vector indicating consistent rxns

metBool: m x 1 boolean vector of metabolites to test for leakage

OPTIONAL INPUTS:

rxnBool: n x 1 boolean vector of reactions to test for leakage modelBoundsFlag: {0,(1)}

0 = set all reaction bounds to -inf, inf

1 = use reaction bounds provided by model.lb and .ub

params: structure with fields:

params.epsilon - (feasTol*100), smallest nonzero mass leak/siphon

params.monoMetMode - {(0), 1} boolean to test for leakage of only one metabolite

printLevel: {(0), 1}

OUTPUTS:

minleakRxnBool: m x 1 boolean of metabolites in a positive leakage mode minleakRxnBool: n x 1 boolean of reactions exclusively involved in a positive leakage mode minsiphonMetBool: m x 1 boolean of metabolites in a negative leakage mode imnsiphonRxnBool: n x 1 boolean of reactions exclusively involved in a negative leakage mode leakY: m x 1 boolean of metabolites in a positive leakage mode siphonY: m x 1 boolean of metabolites in a negative siphon mode statp: status (positive leakage modes)

1 = Solution found

2 = Unbounded

0 = Infeasible

-1 = Invalid input

statn: status (negative leakage modes)

1 = Solution found

2 = Unbounded

0 = Infeasible

-1 = Invalid input

findMinimalLeakageModeRxn(model, rxnBool, metBool, modelBoundsFlag, params, printLevel)¶

Finds a minimal set of leak (or siphon) metabolites and the corresponding minimal set of exclusive reactions involved. Test reactions in rxnBool.

Solve the problem: min ~& ||v||_0 + ||y||_0 \ s.t. ~& Sv - y = 0, \

& l leq v leq u
with either: 0 <= v(rxnBool)
or: v(rxnBool) <= 0
and: 1 <= y(metBool) (semipositive net stoichiometry)
or: y(metBool) <= -1 (seminegative net stoichiometry)

USAGE:

[minLeakMetBool, minLeakRxnBool, minSiphonMetBool, minSiphonRxnBool, leakY, siphonY, statp, statn] = findMinimalLeakageModeRxn(model, rxnBool, metBool, modelBoundsFlag, params, printLevel)

INPUT:

model: structure with fields (bools are not mandatory)

.S - m x n stoichiometric matrix

.lb - Lower bounds

.ub - Upper bounds

.SConsistentMetBool - m x 1 boolean vector indicating consistent mets

.SConsistentRxnBool - m x 1 boolean vector indicating consistent rxns

rxnBool: n x 1 boolean vector of reactions to test for leakage

OPTIONAL INPUTS:

metBool: m x 1 boolean vector of metabolites to test for leakage modelBoundsFlag: {0, (1)}

0 = set all reaction bounds to -inf, inf

1 = use reaction bounds provided by model.lb and .ub

params: structure with fields:

params.eta - (feasTol*100), smallest nonzero mass leak/siphon

params.monoRxnMode - {(1), 0} adds one stoichiometrically inconsistent reaction at a time

printLevel: {(0), 1}

OUTPUTS:

minLeakMetBool m x 1 boolean of metabolites in a positive leakage mode minLeakRxnBool: n x 1 boolean of reactions exclusively involved in a positive leakage mode minSiphonMetBool: m x 1 boolean of metabolites in a negative leakage mode minSiphonRxnBool: n x 1 boolean of reactions exclusively involved in a negative leakage mode leakY: m x 1 boolean of metabolites in a positive leakage mode siphonY: m x 1 boolean of metabolites in a negative siphon mode statp: status (positive leakage modes)

1 = Solution found

2 = Unbounded

0 = Infeasible

-1 = Invalid input

statn: status (negative leakage modes)

1 = Solution found

2 = Unbounded

0 = Infeasible

-1 = Invalid input

findStoichConsistentSubset(model, massBalanceCheck, printLevel, fileName, epsilon)¶

Finds the subset of S that is stoichiometrically consistent using an iterative cardinality optimisation approach

USAGE:

[SConsistentMetBool, SConsistentRxnBool, SInConsistentMetBool, SInConsistentRxnBool, unknownSConsistencyMetBool, unknownSConsistencyRxnBool, model] = findStoichConsistentSubset(model, massBalanceCheck, printLevel, fileName, epsilon)

INPUT:

model: structure with:

.S - m x n stoichiometric matrix

OPTIONAL INPUT:

model.SIntRxnBool n x 1 boolean of reactions heuristically though to be mass balanced massBalanceCheck: {0, (1)} mass and charge balance can be checked by looking at formulas

0 = heuristic detection of exchange reactions (using

findSExRxnInd) will be use to warm start algorithmic part * 1 = reactions that seem mass imbalanced will be used to warm start the algorithmic steps to find the stoichiometrically consistent part. model.metFormulas must exist for mass balance check

printLevel: verbose level fileName: char, used when writing inconsistent metabolites and

reactions to a file

epsilon: (feasTol*10) min nonzero mass, 1/epsilon = max mass

OUTPUT: SConsistentMetBool m x 1 boolean vector indicating stoichiometrically consistent mets SConsistentRxnBool n x 1 boolean vector indicating stoichiometrically consistent rxns SInConsistentMetBool m x 1 boolean vector indicating stoichiometrically inconsistent mets SInConsistentRxnBool n x 1 boolean vector indicating stoichiometrically inconsistent rxns unknownSConsistencyMetBool m x 1 boolean vector indicating unknown consistent mets (all zeros when algorithm converged perfectly!) unknownSConsistencyRxnBool n x 1 boolean vector indicating unknown consistent rxns (all zeros when algorithm converged perfectly!) model

.SConsistentMetBool m x 1 boolean vector indicating consistent mets .SConsistentRxnBool n x 1 boolean vector indicating consistent rxns .SInConsistentMetBool m x 1 boolean vector indicating inconsistent mets .SInConsistentRxnBool n x 1 boolean vector indicating inconsistent rxns .metUnknownInconsistentRemoveBool m x 1 boolean vector indicating removed mets .rxnUnknownInconsistentRemoveBool n x 1 boolean vector indicating removed rxns .unknownSConsistencyMetBool m x 1 boolean vector indicating unknown consistent mets (all zeros when algorithm converged perfectly!) .unknownSConsistencyRxnBool n x 1 boolean vector indicating unknown consistent rxns (all zeros when algorithm converged perfectly!) .SIntMetBool m x 1 boolean of metabolites heuristically though to be involved in mass balanced reactions .SIntRxnBool n x 1 boolean of reactions heuristically though to be mass balanced .metRemoveBool m x 1 boolean vector of metabolites removed to form stoichConsistModel .rxnRemoveBool n x 1 boolean vector of reactions removed to form stoichConsistModel

stoichConsistModel model with stoichiometrically inconsistent heuristically internal reactions removed and any stoichiometrically inconsistent metabolites removed.

leakTest(model, params, printLevel)¶

Solves the problem

\[\begin{split}min ~& \lambda ||v||_0 - \delta ||s||_0 \\ s.t. ~& Sv + s = 0, \\ ~& l \leq v \leq u.\end{split}\]

USAGE:

solution = leakTest(model, params, printLevel)

INPUTS:

model: structure with fields (bools are not mandatory)

.S - m x n stoichiometric matrix

.lb - Lower bounds

.ub - Upper bounds

params: structure with fields printLevel: verbose level

OUTPUT:

solution: Structure containing the following fields:

v - n x 1 vector

s - m x 1 vector

stat - status:

1 = Solution found

2 = Unbounded

0 = Infeasible

-1= Invalid input

maxCardinalityConservationVector(S, param)¶

Maximises the cardinality of the conservation vector:

\[\begin{split}max ~& ||l||_0 \\ st. ~& S^T l = 0 \\ ~& 0 \leq l\end{split}\]

When param.method = ‘optimizeCardinality’; then approximately solve the problem: max ~& ||l||_0 \ st. ~& S^T l = 0 \

~& 0 leq l leq 1/epsilon
When param.method = ‘quasiConcave’; then solve the linear problem: max 1’*z s.t S’*l = 0

z <= l 0 <= l <= 1/epsilon 0 <= z <= epsilon

When param.method = ‘optimizeCardinality’; then solve the difference of convex function problem were the l0 norm is approximated by the capped l1 norm. The resulting problem is solved with a DC program.

When param.method = ‘dc’; then solve the difference of convex function problem were the l0 norm is approximated by the capped l1 norm. The resulting problem is solved with a DC program. Should give the same answer as optimizeCardinality.

USAGE:

[maxConservationMetBool, maxConservationRxnBool, solution] = maxCardinalityConservationVector(S, param)

INPUT:

S: m x n stoichiometric matrix

OPTIONAL INPUTS:

param: structure with:

.nbMaxIteration - Stopping criteria - maximal number of iteration (Default value 1000)

.eta - Smallest value considered non-zero (Default value feasTol*10)

.epsilon - 1/epsilon is the largest molecular mass considered (Default value getCobraSolverParams(‘LP’, ‘feasTol’)*10)

.zeta - Stopping criteria - threshold (Default value 1e-6)

.theta - Parameter of capped l1 approximation (Default value 0.5)

.method - {‘quasiConcave’, (‘optimizeCardinality’)}

OUTPUTS:

maxConservationMetBool: m x 1 boolean for consistent metabolites maxConservationRxnBool: n x 1 boolean for reactions exclusively involving consistent metabolites solution: Structure containing the following fields:

l - m x 1 molecular mass vector

k - n x 1 reaction complex mass vector

stat - status:

1 = Solution found

2 = Unbounded

0 = Infeasible

-1= Invalid input

maxEntConsVector(SInt, printLevel)¶

USAGE:

[m, bool] = maxEntConsVector(SInt, printLevel)

INPUTS:: SInt: m x n matrix printLevel: verbose level
OUTPUTS:: m: primal solution bool: is 1 if m > 0

minCardinalityConservationRelaxationVector(S, param, printLevel)¶

DC programming for solving the cardinality optimization problem

\[\begin{split}min ~& \lambda ||x||_0 \\ s.t. ~& x + S^T z = 0 \\ ~& -\infty \leq x \leq \infty, \\ ~& 1 \leq z \leq 1 / \epsilon\end{split}\]

USAGE:

[relaxRxnBool, solutionRelax] = minCardinalityConservationRelaxationVector(S, param, printLevel)

INPUT:

S: m x n stoichiometric matrix

OPTIONAL INPUTS:

param: structure with:

param.epsilon - (getCobraSolverParams(‘LP’, ‘feasTol’)*100) 1/epsilon is the largest flux expected

param.eta - (feasTol * 100), cutoff for mass leak/siphon

param.nonRelaxBool - (false(n, 1)), n x 1 boolean vector for reactions not to relax

printLevel: verbose level

OUTPUTS:

relaxRxnBool: n x 1 boolean vector where true correspond to relaxation solutionRelax: structure with:

solutionRelax.stat - solution status

solutionRelax.x - n x 1 vector where nonzeros>eta correspond to relaxations

solutionRelax.z - m x 1 vector where positives correspond to molecular mass

optimalConservationVectors(S, lambda, delta)¶

DC programming for solving the cardinality optimization problem

\[\begin{split}min ~& \lambda ||x||_0 - \delta ||y||_0 \\ s.t. ~& x + S^T y = 0 \\ ~& 0 \leq y \leq 1e4\end{split}\]

USAGE:

solution = optimalConservationVectors(S, lambda, delta)

INPUT:: S: m x n stoichiometric matrix
OPTIONAL INPUTS:: lambda: default = 1 delta: default = 1
OUTPUT:: solution: result of the solved cardinality optimization problem

printMinimalLeakageMode(model, minMetBool, minRxnBool, y, printLevel, fileName)¶

Prints out the data on each leakage mode

USAGE:

printMinimalLeakageMode(model, minMetBool, minRxnBool, y, printLevel, fileName)

INPUTS:: model: model structure minMetBool: boolean of metabolites in a positive leakage mode minRxnBool: boolean of reactions exclusively involved in a positive leakage mode y: contains minMetBool printLevel: verbose level fileName: name of the file