Subspaces¶
- checkScaling(model, estLevel, printLevel, matrixA)[source]¶
checks the scaling of the stoichiometric matrix and provides a recommendation on the precision of the solver
- USAGE
[precisionEstimate, scalingProperties] = checkScaling (model, estLevel, printLevel)
- INPUT
model – (the following fields are required - others can be supplied)
S - m x n Stoichiometric matrix
lb - n x 1 Lower bounds on net flux
ub - n x 1 Upper bounds on net flux
- OPTIONAL INPUTS
model –
c - n x 1 Linear objective coefficients
b - m x 1 change in concentration with time
csense - m x 1 character array with entries in {L,E,G} (The code is backward compatible with an m + k x 1 csense vector, where k is the number of coupling constraints)
C - k x n Left hand side of C*v <= d
d - k x 1 Right hand side of C*v <= d
dsense - k x 1 character array with entries in {L,E,G}
estLevel – level of estimation: crude, medium, fine (default)
printLevel – verbose level (default: 1). Level 0 is quiet.
- OUTPUTS
precisionEstimate – estimation of precision (string, double or quad, default: double)
solverRecommendation – cell array with recommendation of solver(s) to be used (default: {‘glpk’})
scalingProperties – structure with properties of scaling
.estLevel: crude, medium, fine (default)
.scltol: value between 0 and 1 (column or row ratio as large as possible)
.matrixA: name of matrix
.nMets: number of metabolites
.nRxns: number of reactions
.minS: minimum of all stoichiometric coefficients
.maxS: maximum of all stoichiometric coefficients
.rmin: minimum of all row scaling coefficients
.imin: index of the minimum of all row scaling coefficients
.rmax: maximum of all row scaling coefficients
.imax: index of the maximum of all row scaling coefficients
.cmin: minimum of all column scaling coefficients
.jmin: index of the minimum of all column scaling coefficients
.cmax: maximum of all column scaling coefficients
.jmax: index of the maximum of all column scaling coefficients
.minLB: minimum of the lower bound vector
.maxLB: maximum of the lower bound vector
.minUB: minimum of the upper bound vector
.maxUB: maximum of the upper bound vector
.ratioS: ratio of the maximum and minimum stoichiometric coefficients
.ratioS_orderOfMag: order of magnitude of the ratio of the maximum and minimum stoichiometric coefficients
.ratioL: ratio of the maximum and minimum values of the lower bound vector
.ratioL_orderOfMag: order of magnitude of the ratio of the maximum and minimum values of the lower bound vector
.ratioU: ratio of the maximum and minimum values of the upper bound vector
.ratioU_orderOfMag: order of magnitude of the ratio of the maximum and minimum values of the upper bound vector
.ratioR: ratio of the maximum and minimum row scaling coefficients
.ratioR_orderOfMag: order of magnitude of the ratio of the maximum and minimum row scaling coefficients
.ratioC: ratio of the maximum and minimum column scaling coefficients
.ratioC_orderOfMag: order of magnitude of the ratio of the maximum and minimum column scaling coefficients
- conservationAnalysis(model, massBalanced, printLevel, tol)[source]¶
Returns the left and right nullspaces of a given stoichiometric matrix in echelon form.
i.e. \(L S = 0\) where \([-Lzero I] Pl^T S = 0\); and \(S N = 0\) where \(S Pn [-Nzero; I] = S, [-Nzero^T I^T] Pn^T = 0\);
- USAGE
[L, N, Lzero, Nzero, Pl, Pn,iR, dR, iC, dC, rankS] = conservationAnalysis (model, massBalanced, printLevel, tol)
- INPUT
model – structure
model.S - m x n stoichiometric matrix
- OPTIONAL INPUTS
massBalanced – {(0), 1, -1}
0 = conservation analysis of entire model.S
1 = conservation analysis of mass balanced reactions only
-1 = conservation analysis of all except biomass reaction
model.SIntRxnBool – boolean vector indicating mass balanced reactions only. Optional if conservation analysis of massBalanced reactions only as otherwise findSExRxnInd will be used to find the mass imbalanced reactions
model.biomassRxnAbbr – String with abbreviation of biomass reaction
model.c – n x 1 linear objective function
printLevel – {(0), 1}
0 = Silent
1 = Print out conservation relations using metabolite and reaction abbreviations
model.mets – m x 1 cell array of metabolite abbreviations
model.rxns – n x 1 cell array of reaction abbreviations
tol – upper bound on tolerance of linear independence,default no greater than 1e-12
- OUTPUTS
L – Echelon form Left nullspace of S
N – Echelon form Right nullspace of S
iR – Boolean index of Independent rows
dR – Boolean index of Dependent rows
iC – Boolean index of Independent columns
dC – Boolean index of Dependent columns
See: Conservation analysis of large biochemical networks Ravishankar Rao Vallabhajosyula , Vijay Chickarmane and Herbert M. Sauro
- detMinSpan(model, params, vectors)[source]¶
Calculates the MinSpan vectors for a COBRA model. The algorithm determines a set of linearly independent basis vectors that span the nullspace that meet the criteria of the flux bounds of the network and minimizes the number of reactions used. The algorithm operates in an interative manner and checkes for convergence after each iteration. Parameters may be provided to skip convergence and terminate the problem when models are large. See Bordbar et al. Mol Syst Biol 2014 for more details. This algorithm has only been tested and requires Gurobi as MILP solver.
- INPUTS
model – COBRA model structure (requires S, lb, ub) Note: MinSpan calculations are done w/o biomass reaction and all reactions must be able to carry a flux = 0 for the trivial solution to be feasible. model is auto corrected by bounds but biomass must be removed manually
- OPTIONAL INPUTS
params – Optional parameters to calculate MinSpan. Determining the MinSpan is a NP hard calculation. For large models, the algorithm may not converge and parameters must be provided to stop the algorithm to provide an approximate solution.
.coverage - Number of iterations to run the algorithm, if not
converge (Default = 10) * .timeLimit - Time to spend on each MinSpan calculation (sec), (Default = 30) * .saveIntV - Save intermediate vectors in order to restart from latest iteration (Default = 0) * .cores - Number of cores to use (Default = 1)
vectors – Set of intermediate MinSpan vectors that may have not yet reached convergence, allowing to pickup calculation from last spot.
- OUTPUTS
finalVectors – MinSpan vectors for COBRA model
- determineSignMatrix(S, sorted)[source]¶
Determine the binaryform of the stoichiometric matrix S and the connectivity vectors
- INPUT
S – Stoichiometric matrix of size m x n (m rows, n columns)
sorted – Boolean flag to sort the connectivity vectors (default: false)
- OUTPUTS
Shat – Binary form of the stoichiometric matrix
Shatabs – Absolute value of Shat
mconnect – Compound connectivity (connectivity number; sum of the number of non-zero entries in a row)
nconnect – Reaction participation number (sum of the number of non-zero entries in a column)
mconnectin – Produced compound connectivity
mconnectout – Consumed compound connectivity
- fourFundamentalSubpaces(S, printLevel)[source]¶
Calculate linearly independent bases the four fundamental subspaces of S using LUSOL or SVD
- USAGE
[N, R, L, C, p, q, rankS] = fourFundamentalSubpaces (S, printLevel)
- INPUT
S – m x n stoichiometric matrix defined to be of rank r.
- OPTIONAL INPUT
printLevel – default = 1
- OUTPUTS
N – n x (n-r) null space basis The null space of S contains all the steady-state flux distributions allowable in the network. The steady-state is of much interest since most homeostatic states are close to being steady-states.
R – n x r row space basis: The row space of S contains all the dynamic flux distributions of a network, and thus the thermodynamic driving forces that change the rate of reaction activity.
L – m x (m-r) left null space basis: The left null space of S contains all the conservation relationships, or time-invariants, that a network contains. The sum of conserved metabolites or conserved metabolic pools do not change with time and are combinations of concentration variables.
C – m x r column space basis: The column space of S contains all the possible time derivatives of the concentration vector, and thus how the thermodynamic driving forces move the concentration state of the network.
rankS – r, rank of the stoichiometric matrix
p – p(1:rankS) gives indices to the linearly independent columns of S
p(rankS + 1:end) gives indices to the linearly dependent columns of S
q – q(1:rankS) gives indices to the linearly independent rows of S
q(rankS + 1:end) gives indices to the linearly dependent rows of S
- gmscale(A, iprint, scltol)[source]¶
Geometric-Mean Scaling finds the scale values for the m x n sparse matrix A.
- USAGE
[cscale, rscale] = gmscale (A, iprint, scltol)
- INPUTS
A (i, j) – contains entries of A.
iprint – > 0 requests messages to the screen (0 means no output).
scltol – should be in the range (0.0, 1.0). Typically scltol = 0.9. A bigger value like 0.99 asks gmscale to work a little harder (more passes).
- OUTPUTS
cscale, rscale – column vectors of column and row scales such that R (inverse) A C (inverse) should have entries near 1.0, where R= diag(rscale), C = diag(cscale).
An iterative procedure based on geometric means is used, following a routine written by Robert Fourer, 1979. Several passes are made through the columns and rows of A. The main steps are:
Compute \(aratio = max_j (max_i Aij / min_i Aij)\).
Divide each row i by \(sqrt( max_j Aij * min_j Aij)\).
Divide each column j by \(sqrt( max_i Aij * min_i Aij)\).
Compute sratio as in Step 1.
If \(sratio < aratio * scltol\), set \(aratio = sratio\) and repeat from Step 2.
To dampen the effect of very small elements, on each pass, a new row or column scale will not be smaller than sqrt(damp) times the largest (scaled) element in that row or column.
Use of the scales: To apply the scales to a linear program, \(min c^T x\) st \(A x = b\), \(l \leq x \leq u\), we need to define “barred” quantities by the following relations: A = R Abar C, b = R bbar, C cbar = c, C l = lbar, C u = ubar, C x = xbar.
This gives the scaled problem \(min\ cbar^T xbar\) st \(Abar\ xbar = bbar\), \(lbar \leq xbar \leq ubar\).