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 nonnegative 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 nonexchange reactions (optional)
model.SIntRxnBool n x 1 boolean vector indicating the reactions that are thought to be nonexchange 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 nonexchange 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 nonexchange 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 constraintbased 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

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 nonzero (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 1e6)
.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