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 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)[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 constraint-based 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

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)[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