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

[inform, 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 - {(‘solveCobraLP’), ‘cvx’, ‘nonconvex’} interface called to do the consistency check
  • method.solver - {(default solver),’gurobi’,’mosek’}
  • method.param - solver specific parameter structure

Outputs

• inform – 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
  • model.SConsistentRxnBool - n x 1 boolean vector non-exchange reaction involving a stoichiometrically consistent metabolite

findMassLeaksAndSiphons(model, metBool, rxnBool, modelBoundsFlag, params, 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, params, 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
• params – structure with fields:
  • params.epsilon - (1e-4)
  • params.eta - (feasTol*100), smallest nonzero mass leak/siphon
  • params.theta - (0.5) parameter of capped l1 approximation
  • params.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

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 – 1e-4, 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)

Inputs

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

Inputs

• 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 inputs

• massBalanceCheck – {(0), 1} where:

  • 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*100) min nonzero mass, 1/epsilon = max mass

OUTPUT: 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 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

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, params)

Maximises the cardinality of the conservation vector:

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

The l0 norm is approximated by capped l1 norm. The resulting problem is a DC program

Usage

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

Input

• S – m x n stoichiometric matrix

Optional input

• params – structure with:
  • params.nbMaxIteration - Stopping criteria - maximal number of iteration (Default value 1000)
  • params.epsilon - 1/epsilon is the largest molecular mass considered (Default value 1e-4)
  • params.zeta - Stopping criteria - threshold (Default value 1e-6)
  • params.theta - Parameter of capped l1 approximation (Default value 0.5)

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
  • 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, params, 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, params, printLevel)

Input

• S – m x n stoichiometric matrix

Optional inputs

• params – structure with:
  • params.epsilon - (1e-4) 1/epsilon is the largest flux expected
  • params.eta - (feasTol * 100), cutoff for mass leak/siphon
  • params.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 tru 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