Stoichconsistency¶
 checkStoichiometricConsistency(model, printLevel, method)[source]¶
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)[source]¶
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)[source]¶
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)[source]¶
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)[source]¶
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)[source]¶
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)[source]¶
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)[source]¶
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)[source]¶
 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)[source]¶
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)[source]¶
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)[source]¶
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