Analysis#

MATLAB functions in RAVEN/src/analysis of the RAVEN toolbox.

Functions#

Function	Summary
`analyzeSampling`	Compare flux change significance with expression change.
`compareFluxes`	Find the reactions whose flux differs between two flux vectors.
`findGeneDeletions`	Delete genes and track the resulting fluxes.
`FSEOF`	Flux Scanning based on Enforced Objective Flux.
`getAllowedBounds`	Return the minimal and maximal fluxes through reactions.
`getAllSubGraphs`	Get all metabolic subgraphs in a model.
`getEssentialRxns`	Calculate the essential reactions for a solvable model.
`getFluxZ`	Calculate Z scores between two sets of random flux distributions.
`getMinimalMedium`	Find the minimal set of uptake reactions for growth.
`haveFlux`	Check which reactions can carry a flux.
`randomSampling`	Sample the flux solution space (entry point for all samplers).
`reporterMetabolites`	Identify metabolites around which transcriptional changes occur.
`runDynamicFBA`	Perform dynamic FBA using the static optimisation approach.
`runPhenotypePhasePlane`	Run phenotype phase plane analysis and plot the results.
`runProductionEnvelope`	Calculate the byproduct secretion envelope.
`runRobustnessAnalysis`	Perform robustness analysis for a reaction and objective.
`runSimpleOptKnock`	Find deletions that couple metabolite production to growth.
`sampleACHR`	Artificially Centered Hit-and-Run flux sampling.
`sampleChebyshevCenter`	Centre of the largest ball inscribed in a polytope.
`sampleCHRR`	Coordinate Hit-and-Run with Rounding flux sampling.
`sampleMaxVolEllipse`	Maximum-volume inscribed ellipsoid of a polytope.
`sampleWarmupPoints`	FVA warmup vertices for ACHR sampling.
`traceFluxPath`	Find the highest-flux-fraction path between two reactions.
`walkFluxes`	Interactively navigate a flux distribution reaction by reaction.

Reference#

analyzeSampling#

function scores=analyzeSampling(Tex, df, solutionsA, solutionsB, varargin)

Compare flux change significance with expression change.

Compares the significance of change in flux between two conditions with the significance of change in gene expression.

Parameters

Tex : double: a vector of t-scores for the change in gene expression for each reaction. The score is the Student t between the two conditions, or is calculated from a p-value by the inverse of the error function. Transcripts that increased in expression level have positive values and transcripts that decreased have negative values. A p-value does not give the direction of the change, so a score calculated from a p-value takes its sign from the direction of the expression change.
df : double: the degrees of freedom in the t-test.
solutionsA : double: random solutions for the reference condition (as generated by randomSampling).
solutionsB : double: random solutions for the test condition (as generated by randomSampling).

Name-Value Arguments

printResults : logical: true to print the 10 reactions with the highest probability in each category (default false).

Returns

scores : double

a Nx3 column matrix with the probabilities of a reaction:

changing both in flux and expression in the same direction
changing in expression but not in flux
changing in flux but not in expression or changing in opposed directions in flux and expression

Examples

scores = analyzeSampling(Tex, df, solutionsA, solutionsB, ...
                         printResults);

compareFluxes#

function result = compareFluxes(model, fluxes1, fluxes2, varargin)

Find the reactions whose flux differs between two flux vectors.

Compares two flux vectors from the same model and returns all reactions that changed, sorted by absolute flux difference. Classifies reactions as turned on, turned off, reversed direction, or changed magnitude.

Parameters

model : struct: a model structure.
fluxes1 : double: reference flux vector (condition A, such as wild-type).
fluxes2 : double: comparison flux vector (condition B, such as a mutant or knock-in).

Name-Value Arguments

cutoff : double: minimum |flux| for a reaction to be active, and minimum |flux2 - flux1| for a reaction to be changed (default 1e-8).
metaboliteList : cell: cell array of metabolite names. Only reactions involving at least one of these metabolites are compared; all other reactions are excluded from every field of the result. Matching is against model.metNames and is case insensitive. A name that matches no metabolite is reported as a warning (default all reactions).
nMax : double: maximum rows printed in the summary table (default 20).
verbose : logical: print the summary table (default true).

Returns

result : struct

structure with fields:

turnedOn : cell array with the IDs of reactions that are active in fluxes2 and not in fluxes1.
turnedOff : cell array with the IDs of reactions that are active in fluxes1 and not in fluxes2.
flipped : cell array with the IDs of reactions that are active in both conditions and change sign.
changed : struct with fields for all reactions with |flux2 - flux1| above cutoff, sorted by |flux2 - flux1| descending:
- rxn : cell array with the reaction IDs.
- flux1 : double array with the flux in condition A.
- flux2 : double array with the flux in condition B.
- absDelta : double array with |flux2 - flux1|.
- relChange : double array with (flux2 - flux1)/|flux1|. NaN when |flux1| is below cutoff.
- type : cell array with 'on', 'off', 'flip' or '' (changed magnitude) for each reaction.

Examples

solA = solveLP(model);
modelAna = setParam(model, 'ub', 'O2exchange', 0);
solB = solveLP(modelAna);
result = compareFluxes(model, solA.x, solB.x)
result = compareFluxes(model, solA.x, solB.x, 'nMax', 50)
result = compareFluxes(model, solA.x, solB.x, 'metaboliteList', {'ATP','NADH'})

See also

walkFluxes, traceFluxPath, modelSummary

findGeneDeletions#

function [genes, fluxes, originalGenes, details, grRatioMuts]=findGeneDeletions(model,varargin)

Delete genes and track the resulting fluxes.

Deletes single genes or gene pairs, solves the model by FBA after each deletion, and returns the resulting fluxes. The output identifies gene deletion targets.

Parameters

model : struct: a model structure.

Name-Value Arguments

testType : char

type of gene deletion (default "sgd"):

"sgd" : single gene deletion
"dgd" : double gene deletion

Returns

genes : double

a matrix with the genes that were deleted in each simulation (the gene indexes in originalGenes). Each row corresponds to a column in fluxes.

fluxes : double

a matrix with the resulting fluxes. Double deletions that result in an unsolvable problem have all zero flux. Single deletions that result in an unsolvable problem are indicated in details instead.

originalGenes : cell

the genes in the input model.

details : double

a vector with one element per gene in originalGenes that states whether the gene was deleted and why. Genes that only belong to dead-end reactions, for example, are not deleted:

1 : deleted, and the problem was solved (a lethal deletion with zero growth also gets 1)
2 : the single-deletion problem had no optimal solution (infeasible or failed); fluxes and grRatioMuts are 0
3 : redundant (not assigned)
4 : involved in dead-end reaction

grRatioMuts : double

growth rate ratio between mutated strain and wild type, matching the originalGenes(genes) mutants. The elements do not map directly to model.genes, unlike the output of COBRA getEssentialGenes. The mapping to model.genes is obtained with:

grRatio=zeros(1,numel(model.genes));
grRatio(genes)=grRatioMuts;

Examples

[genes, fluxes, originalGenes, details, grRatioMuts]=...
    findGeneDeletions(model,testType);

FSEOF#

function targets=FSEOF(model,biomassRxn,targetRxn,varargin)

Flux Scanning based on Enforced Objective Flux.

Implements the Flux Scanning based on Enforced Objective Flux algorithm. This function writes a tab-delimited file or prints to the command window. If an output is specified (targets), it also returns a structure indicating for each model reaction whether FSEOF identifies it as a target and the slope of the reaction when switching from biomass formation to product formation.

Parameters

model : struct: a model structure.
biomassRxn : char: reaction ID of the biomass formation or growth reaction.
targetRxn : char: reaction ID of the target reaction.

Name-Value Arguments

iterations : double: number of enforced target flux levels, equally spaced up to coefficient times the maximum target flux (default 10).
coefficient : double: ratio of the maximum target reaction flux that is enforced at the last iteration, less than 1 (default 0.9).
outputFile : char: output filename (default prints to command window).
corrThreshold : double: minimum absolute Pearson correlation coefficient between the absolute flux of a reaction and the enforced product flux for the reaction to be selected as a target. Values are in [0, 1]. The filter keeps reactions whose flux is close to linear in the enforced product flux, and a higher value keeps fewer reactions (default 0.9).

Returns

targets : struct

structure with information for identified targets, with fields:

logical : logical array, true for amplification targets (identical to amplify).
amplify : logical array, reactions whose absolute flux increases with enforced product flux.
knockdown : logical array, reactions whose absolute flux decreases but stays above zero at the highest enforced product flux.
knockout : logical array, reactions whose absolute flux decreases to zero at the highest enforced product flux.
slope : numeric array with per-reaction regression slopes of absolute flux versus enforced product flux (for all reactions).

Examples

targets = FSEOF(model, biomassRxn, targetRxn, iterations, ...
                coefficient, outputFile);

getAllowedBounds#

function [minFluxes, maxFluxes, exitFlags]=getAllowedBounds(model,varargin)

Return the minimal and maximal fluxes through reactions.

Parameters

model : struct: a model structure.

Name-Value Arguments

rxns : cell or logical or double: either a cell array of reaction IDs, a logical vector with the same number of elements as reactions in the model, or a vector of reaction indexes (default model.rxns).
runParallel : logical: run the calculations in parallel (default true). Parallel processing requires the MATLAB Parallel Computing Toolbox; without it, the calculations run serially regardless of this setting. For a few reactions, the overhead of parallel processing exceeds the time it saves.

Returns

minFluxes : double: minimal allowed fluxes.
maxFluxes : double: maximal allowed fluxes.
exitFlags : double: a matrix with one row per reaction. Column 1 holds the solver status of the minimisation and column 2 the status of the maximisation (1 if a flux could be calculated).

Notes

In cases where no solution can be calculated, NaN is returned.

Examples

[minFluxes, maxFluxes, exitFlags] = getAllowedBounds(model, rxns, runParallel);

getAllSubGraphs#

function subGraphs=getAllSubGraphs(model)

Get all metabolic subgraphs in a model.

Two metabolites are connected if they share a reaction.

Parameters

model : struct: a model structure.

Returns

subGraphs : logical: a logical matrix with one row per metabolite and one column per subgraph. An element is true if the metabolite belongs to the subgraph. The columns are ordered by decreasing subgraph size.

Examples

subGraphs = getAllSubGraphs(model);

getEssentialRxns#

function [essentialRxns, essentialRxnsIndexes]=getEssentialRxns(model,varargin)

Calculate the essential reactions for a solvable model.

Parameters

model : struct: a model structure.

Name-Value Arguments

ignoreRxns : cell: cell array of reaction IDs which should not be checked (default {}).

Returns

essentialRxns : cell: cell array with the IDs of the essential reactions.
essentialRxnsIndexes : double: vector with the indexes of the essential reactions.

Notes

Essential reactions are those which, when constrained to 0, result in an infeasible problem.

Examples

[essentialRxns, essentialRxnsIndexes] = getEssentialRxns(model, ignoreRxns);

getFluxZ#

function Z=getFluxZ(solutionsA, solutionsB)

Calculate Z scores between two sets of random flux distributions.

Parameters

solutionsA : double: random solutions for the reference condition (as generated by randomSampling).
solutionsB : double: random solutions for the test condition (as generated by randomSampling).

Returns

Z : double: a vector with one Z-score per reaction: the difference between the mean fluxes of solutionsB and solutionsA, divided by the square root of the sum of their variances. A positive sign indicates an increase of the flux in the second condition with respect to the first, a negative sign a decrease. Values are limited to the range -100 to 100.

Examples

Z = getFluxZ(solutionsA, solutionsB);

See also

randomSampling

getMinimalMedium#

function [medium, mediumIdx] = getMinimalMedium(model, varargin)

Find the minimal set of uptake reactions for growth.

Uses MILP to find the smallest set of exchange reactions that can carry uptake flux (lb < 0) while the model achieves a target growth rate. The MILP formulation gives a medium with the globally minimal number of uptake reactions.

Parameters

model : struct: a model structure.

Name-Value Arguments

minGrowth : double: minimum objective value (growth rate) the medium must support. If empty, 10 % of the optimum objective value under the current bounds, which requires one extra FBA call. A value given explicitly avoids that call (for example 'minGrowth', 0.1) (default empty).
verbose : logical: print the resulting minimal medium (default true).

Returns

medium : cell: cell array of reaction IDs that form the minimal medium. Empty if no feasible medium can support the requested growth.
mediumIdx : double: column vector of reaction indices corresponding to medium.

Examples

[med, idx] = getMinimalMedium(model)
[med, idx] = getMinimalMedium(model, 'minGrowth', 0.05)
[med, idx] = getMinimalMedium(model, 'verbose', false)

Notes

Only exchange reactions with lb < 0 are candidate uptake reactions. Exchange reactions with lb >= 0 are not included.

The MILP formulation: minimise sum(y_i) s.t. S v = 0 c' v >= minGrowth v_i >= lb_i * y_i for each candidate uptake reaction i lb_j <= v_j <= ub_j for all reactions j y_i in {0, 1}

A y_i = 0 forces reaction i to carry non-negative flux (no net uptake); y_i = 1 allows it to carry flux as low as lb_i.

The number of uptake reactions is minimal, but the medium is not necessarily unique: when several media of the same size support the growth rate, the solver returns one of them, and another solver or solver version can return a different one.

The MILP is solved with the solver set by setRavenSolver. RAVEN's bundled GLPK solves it for small models. The RAVEN:glpkMILP warning issued with GLPK states that GLPK works for small problems and that Gurobi is much faster and more reliable for larger ones (setRavenSolver('gurobi')).

haveFlux#

function I=haveFlux(model,varargin)

Check which reactions can carry a flux.

Checks which reactions can carry a (positive or negative) flux under the model's current bounds; other media or bounds are not explored. Unlike getAllowedBounds, which minimises and maximises every reaction, it skips a reaction once an earlier solution has given it a flux, and minimises only reversible reactions that still have none.

Parameters

model : struct: a model structure.

Name-Value Arguments

cutOff : double: the absolute flux value that a reaction has to exceed to be identified as carrying flux (default 10^-6).
rxns : cell or logical or double: either a cell array of IDs, a logical vector with the same number of elements as reactions in the model, or a vector of indexes (default model.rxns).
seed : double: seed for the random number generator that shuffles the order in which reactions are tested. The order is then reproducible across runs (default: the current rng state is left untouched).

Returns

I : logical: logical array with true if the corresponding reaction can carry a flux under the tested bounds.

Examples

I = haveFlux(model, cutOff, rxns);

Notes

Bounds of +/- Inf are replaced with +/- 10000 before solving.

Dead-end reactions are removed (simplifyModel) before solving. A dead-end reaction with a forced flux (lb > 0 or ub < 0) makes the model infeasible, but it is removed like any other dead end, so the remaining reactions are tested as if it were absent. A forced flux elsewhere in an infeasible model gives false for all reactions.

randomSampling#

function [solutions, goodRxns, info]=randomSampling(model,varargin)

Sample the flux solution space (entry point for all samplers).

Dispatches to one of three sampling methods via the 'method' argument: 'achr' (default) Artificially Centered Hit-and-Run; near-uniform MCMC sampling of the polytope interior (see sampleACHR). 'chrr' Coordinate Hit-and-Run with Rounding; near-uniform MCMC with maximum-volume-ellipsoid rounding, which improves mixing on thin or ill-conditioned polytopes such as enzyme-constrained models (see sampleCHRR). 'randomObjective' the random-objective vertex method of Bordel et al. (2010) PLoS Comput Biol (doi:10.1371/journal.pcbi.1000859): each sample optimises a random objective of nObjectives reactions with random weights, returning a polytope vertex.

The 'achr' and 'chrr' methods draw from the (near-)uniform distribution over the interior; the 'randomObjective' method draws vertices. Arguments below marked [randomObjective] apply only to that method; 'thinning' and 'nBurnin' apply only to the MCMC methods.

Parameters

model : struct: a model structure.

Name-Value Arguments

nSamples : double: the number of solutions to return (default 1000).
method : char: sampling method: 'achr' (default), 'chrr', or 'randomObjective'.
thinning : double: MCMC steps between recorded samples for 'achr' and 'chrr' (default 100).
nBurnin : double: burn-in steps discarded before the first sample for 'chrr' (default 1000).
replaceBoundsWithInf : logical: [randomObjective] replace the largest upper bounds with Inf and the smallest lower bounds with -Inf (the latter only if the smallest lower bound is negative). This is required to get solutions without loops if the model has arbitrarily large bounds, for example 1000 and -1000. A model with only biologically relevant bounds can use false (default true).
supressErrors : logical: if false, randomSampling raises an error after 99 consecutive failed attempts to find a non-zero solution that is not involved in loops. Failures can result from constraints that are too relaxed (such as unlimited glucose uptake) or too strict (such as many narrow constraints). If true, the attempts continue without an error (default false).
runParallel : logical: run the calculations in parallel (default true). Requires the MATLAB Parallel Computing Toolbox; without it, the calculations run serially regardless of this setting.
goodRxns : double: vector of indexes of those reactions that are not involved in loops and can be used as random objective functions, as generated by a previous run of randomSampling on the same model (default empty).
minFlux : logical: determines if a second optimisation is performed for each random sample, which minimises the sum of absolute fluxes to prevent loops (default false).
nObjectives : double: number of reactions included in the random objective function per sample iteration (default 2).
seed : double: random seed for reproducibility. If empty, the current MATLAB random state is used (default empty).

Returns

solutions : double: matrix with the solutions.
goodRxns : double: [randomObjective] vector of indexes of those reactions that are not involved in loops or always carry zero flux and can be used as random objective functions. Empty ([]) for the 'achr' and 'chrr' methods.
info : struct: [chrr] the diagnostics from sampleCHRR, including mveConverged (whether the maximum-volume ellipsoid rounding reached tolerance). Empty ([]) for the 'achr' and 'randomObjective' methods.

Notes

With 'randomObjective', the solutions are generated by optimising (with random weights) a random set of nObjectives reactions (default 2). The sign of each weight is chosen at random, which selects maximisation or minimisation of the reaction; reactions with an upper bound below 0 always receive a positive weight. The default 'achr' and the 'chrr' methods use no random objectives; they take hit-and-run steps through the interior of the solution space.

If the model is a GECKO v3+ ecModel, then usage_prot reactions are not selected as random objectives; the objectives are drawn from the metabolic reactions.

Examples

% Near-uniform interior sampling (default, ACHR):
solutions = randomSampling(model, 1000);
% Coordinate hit-and-run with rounding:
solutions = randomSampling(model, 1000, 'method', 'chrr', 'seed', 1);
% Random-objective vertex method:
solutions = randomSampling(model, 1000, 'method', 'randomObjective');

See also

sampleACHR, sampleCHRR, sampleMaxVolEllipse

reporterMetabolites#

function repMets=reporterMetabolites(model,genes,genePValues,varargin)

Identify metabolites around which transcriptional changes occur.

Applies the Reporter Metabolites algorithm, which scores each metabolite by the differential expression of the genes around it.

Parameters

model : struct: a model structure.
genes : cell: a cell array of gene names (matching model.genes).
genePValues : double: P-values for differential expression of the genes. Values must lie in [0, 1]; NaN values are ignored, any other value raises an error.

Name-Value Arguments

printResults : logical: true to print the top 20 Reporter Metabolites to the screen (default false).
outputFile : char: the results are printed to this file (default none).
geneFoldChanges : double: log-fold changes for the genes. If supplied, Reporter Metabolites are also calculated for the upregulated genes (fold change >= 0) alone and for the downregulated genes (fold change < 0) alone, in addition to the full test (default none).

Returns

repMets : struct

an array of structures with the following fields:

test : a string that describes the genes used to calculate the Reporter Metabolites ("all", "only up", or "only down"). The tests "only up" and "only down" are only calculated if geneFoldChanges is supplied.
mets : a cell array of metabolite IDs for the metabolites for which a score could be calculated.
metNames : a cell array with the names of the metabolites in "mets".
metZScores : background-corrected Z-scores for differential expression around each metabolite in "mets". The entries of all fields are sorted by this score, highest first.
metPValues : one-sided (upper-tail) P-values of "metZScores".
metNGenes : number of neighbouring genes for each metabolite in "mets".
meanZ : average Z-scores for the genes around each metabolite in "mets".
stdZ : standard deviations of the Z-scores around each metabolite in "mets".

Examples

repMets = reporterMetabolites(model, genes, genePValues, ...
    printResults, outputFile, geneFoldChanges);

Notes

For details about the algorithm, see Patil KR, Nielsen J, Uncovering transcriptional regulation of metabolism by using metabolic network topology. Proc. Natl Acad. Sci. USA 2005;102:2685-2689.

The background correction uses the closed-form mean and standard deviation of the aggregate Z-score of n genes drawn at random from all scored genes, so the results are deterministic.

runDynamicFBA#

function [concentrationMatrix, excRxnNames, timeVec, biomassVec] = runDynamicFBA(model, substrateRxns, initConcentrations, initBiomass, timeStep, nSteps, plotRxns, exclUptakeRxns)

Perform dynamic FBA using the static optimisation approach.

Parameters

model : struct: a model structure.
substrateRxns : cell: cell array with exchange reaction identifiers for substrates that are initially in the media, whose concentration may change (for example not h2o or co2).
initConcentrations : double: initial concentrations of substrates (matching substrateRxns).
initBiomass : double: initial biomass (must be non-zero).
timeStep : double: time step size.
nSteps : double: maximum number of time steps.
plotRxns : cell: cell array with exchange reaction identifiers for substrates whose concentration should be plotted.
exclUptakeRxns : cell: cell array with exchange reaction identifiers for substrates whose concentration does not change (for example co2, o2, h2o, h).

Returns

concentrationMatrix : double: matrix with extracellular metabolite concentrations.
excRxnNames : cell: cell array with exchange reaction identifiers that match the metabolites included in the concentrationMatrix.
timeVec : double: vector of time points.
biomassVec : double: vector with biomass concentrations.

Examples

[concentrationMatrix, excRxnNames, timeVec, biomassVec] = ...
    runDynamicFBA(model, substrateRxns, initConcentrations, ...
    initBiomass, timeStep, nSteps, plotRxns, exclUptakeRxns);

Notes

If no initial concentration is given for a substrate that has an open uptake in the model (that is, model.lb < 0), the initial concentration is set to 1000. If the uptake bound calculated from the concentration exceeds the maximum uptake rate specified in the model for that nutrient, and the maximum uptake rate specified is > 0, the maximum uptake rate specified in the model is used as the uptake bound.

Modified from COBRA Toolbox dynamicFBA.m.

runPhenotypePhasePlane#

function [growthRates, shadowPrices1, shadowPrices2] = runPhenotypePhasePlane(model, controlRxn1, controlRxn2, varargin)

Run phenotype phase plane analysis and plot the results.

Runs phenotype phase plane analysis and plots the results. The first plot is a 3D surface plot showing the phenotype phase plane, the other two plots show the shadow prices of the metabolites from the two control reactions, which define the phases.

Parameters

model : struct: a model structure.
controlRxn1 : char: reaction identifier of the first reaction to be plotted.
controlRxn2 : char: reaction identifier of the second reaction to be plotted.

Name-Value Arguments

nPts : double: the number of points to plot in each dimension (default 50).
range1 : double: the range [from 0 to range1] of reaction 1 to plot (default 20).
range2 : double: the range [from 0 to range2] of reaction 2 to plot (default 20).

Returns

growthRates : double: a matrix of maximum growth rates.
shadowPrices1 : double: a matrix with shadow prices for reaction 1.
shadowPrices2 : double: a matrix with shadow prices for reaction 2.

Examples

[growthRates, shadowPrices1, shadowPrices2] = ...
    runPhenotypePhasePlane(model, controlRxn1, controlRxn2, nPts, ...
    range1, range2);

Notes

Modified from COBRA Toolbox phenotypePhasePlane.m.

runProductionEnvelope#

function [biomassValues, targetValues] = runProductionEnvelope(model, targetRxn, biomassRxn, varargin)

Calculate the byproduct secretion envelope.

Parameters

model : struct: a model structure.
targetRxn : char: identifier of target metabolite production reaction.
biomassRxn : char: identifier of biomass reaction.

Name-Value Arguments

nPts : double: number of biomass values at which the target bounds are calculated (default 20).

Returns

biomassValues : double: nPts biomass values, evenly spaced from the minimum to the maximum biomass flux.
targetValues : double: a matrix with two rows: the lower bounds (row 1) and the upper bounds (row 2) of the target reaction flux at each biomass value. NaN where the problem has no solution.

Examples

[biomassValues, targetValues] = runProductionEnvelope(model, ...
    targetRxn, biomassRxn, nPts);

Notes

Modified from COBRA Toolbox productionEnvelope.m.

runRobustnessAnalysis#

function [controlFlux, objFlux] = runRobustnessAnalysis(model, controlRxn, varargin)

Perform robustness analysis for a reaction and objective.

Varies the flux of a reaction of interest from its minimum to its maximum and calculates the optimal value of an objective of interest at each flux.

Parameters

model : struct: a model structure.
controlRxn : char: reaction of interest whose value is to be controlled.

Name-Value Arguments

nPoints : double: number of control flux values (default 20).
objRxn : char: reaction identifier of the objective to be maximised (default: the objective defined in the model).
plotRedCost : logical: whether the reduced cost of the control reaction is also plotted on the right axis (default false).

Returns

controlFlux : double: flux values of the reaction of interest, ranging from its minimum to its maximum value.
objFlux : double: optimal values of objective reaction at each control reaction flux value.

Examples

[controlFlux, objFlux] = runRobustnessAnalysis(model, controlRxn, ...
    nPoints, objRxn);

Notes

Modified from COBRA Toolbox robustnessAnalysis.m.

runSimpleOptKnock#

function out = runSimpleOptKnock(model, targetRxn, biomassRxn, varargin)

Find deletions that couple metabolite production to growth.

Checks all gene or reaction deletions for growth-coupled metabolite production, by testing combinations of up to maxNumKO deletions. It tries a deeper deletion only after a deletion whose growth rate is above minGrowth. A combination is reported if the growth rate exceeds minGrowth and the product of growth rate and production rate exceeds the value without the last added deletion (the wild type for the first deletion) by more than 1 %. The method is not defined as a MILP and is therefore slow.

Parameters

model : struct: a model structure.
targetRxn : char: identifier of target reaction.
biomassRxn : char: identifier of biomass reaction.

Name-Value Arguments

deletions : cell: cell array with gene or reaction identifiers that should be considered for knockout (default model.rxns).
genesOrRxns : char: string indicating whether deletions parameter is given with 'genes' or 'rxns' identifiers (default 'rxns').
maxNumKO : double: maximum number of simultaneous knockouts (default 1).
minGrowth : double: minimum growth rate (default 0.05).

Returns

out : struct

structure with deletion strategies that result in growth-coupled production, with fields:

KO : cell array with gene(s) or reaction(s) to be deleted.
growthRate : vector with growth rates after deletion.
prodRate : vector with production rates after deletion.
score : vector with the product of growth rate and production rate.

Examples

out = runSimpleOptKnock(model, targetRxn, biomassRxn, deletions, ...
    genesOrRxns, maxNumKO, minGrowth);

sampleACHR#

function solutions = sampleACHR(model, nSamples, thinning, warmup, seed)

Artificially Centered Hit-and-Run flux sampling.

Draws nSamples flux vectors approximately uniformly from the flux polytope {v : Sv = 0, lb <= v <= ub}, using the ACHR algorithm (Kaufman & Smith 1998). No rounding is performed, so on strongly elongated polytopes (for example enzyme-constrained models) mixing is slower than with sampleCHRR, which rounds the polytope. Reactions whose bounds or warmup values coincide never move. The remaining reactions are corrected towards Sv = 0 every few steps and before each recorded sample, with a sparse Cholesky factor of the row-normalised SS' (no dense null-space basis), so the rounding error in Sv does not grow with the number of steps and samples stay inside the bounds. The guarantee S*v = 0 holds for a well-conditioned S and warmup points that are feasible for the model. When the largest residual of a recorded sample exceeds 1e-9 times max(1, max(abs(v))) (nearly dependent rows of S, a noisy or infeasible warmup), sampleACHR issues the warning sampleACHR:steadyStateNotRestored.

Parameters

model : struct: a model structure.
nSamples : double: number of flux vectors to return (default 1000).
thinning : double: Markov-chain steps between recorded samples (default 100).
warmup : double: nPoints-by-nRxns warmup vertices (default: generated by sampleWarmupPoints).
seed : double: random seed for reproducibility (default: current RNG state).

Returns

solutions : double: nRxns-by-nSamples matrix of sampled flux distributions.

See also

randomSampling, sampleCHRR, sampleWarmupPoints

sampleChebyshevCenter#

function x0 = sampleChebyshevCenter(A, b)

Centre of the largest ball inscribed in a polytope.

Solves the LP max r s.t. a_i'*x + ||a_i||r <= b_i, r >= 0, giving a strictly interior point of {x : Ax <= b}. Used to initialise the interior-point MVE solver (sampleMaxVolEllipse) inside sampleCHRR.

Parameters

A : double: m-by-n constraint matrix.
b : double: m-by-1 right-hand side.

Returns

x0 : double: n-by-1 Chebyshev centre (strictly interior point).

See also

sampleMaxVolEllipse, sampleCHRR

sampleCHRR#

function [solutions, info] = sampleCHRR(model, nSamples, thinning, nBurnin, seed, fixedWidthTol)

Coordinate Hit-and-Run with Rounding flux sampling.

Draws nSamples flux vectors approximately uniformly from the flux polytope {v : S*v = 0, lb <= v <= ub}, using the CHRR algorithm (Haraldsdottir et al. 2017, Bioinformatics 33:1741). The polytope is rounded before sampling, which improves mixing on strongly elongated polytopes (for example enzyme-constrained models) in comparison with sampleACHR, which does not round.

Constraints to condition on (for example a biomass lower bound or measured exchange fluxes) are set on the model before the call.

Parameters

model : struct: a model structure.
nSamples : double: number of flux vectors to return (default 1000).
thinning : double: Markov-chain steps between recorded samples (default 100).
nBurnin : double: burn-in steps discarded before the first recorded sample (default 1000).
seed : double: random seed for reproducibility (default: current RNG state).
fixedWidthTol : double: a reaction whose allowed flux range is narrower than this is treated as stoichiometrically fixed and added to the equality system at the middle of its range, which keeps the reduced polytope full-dimensional (default 1e-7).

Returns

solutions : double: nRxns-by-nSamples matrix of sampled flux distributions.
info : struct: .nDimensions : dimension of the sampled flux polytope .mveConverged : whether MVE rounding reached tolerance .fixedRxns : IDs of reactions added to the equality system as fixed

See also

randomSampling, sampleACHR, sampleMaxVolEllipse

sampleMaxVolEllipse#

function [x, E, converged] = sampleMaxVolEllipse(A, b, x0, maxiter, tol, reg)

Maximum-volume inscribed ellipsoid of a polytope.

Solves max log det E over centre x and SPD matrix E such that the ellipsoid {x + Es : ||s||_2 <= 1} is contained in the polytope {z : Az <= b}, using the Zhang & Gao (2003) primal-dual interior-point method (the regularised variant used in COBRA's CHRR sampler). This is the rounding step of sampleCHRR: the returned transform E maps the polytope to a well-conditioned body in which coordinate hit-and-run samples.

This function is pure linear algebra (no LP solver) except when x0 is empty, in which case a Chebyshev centre is obtained via sampleChebyshevCenter.

Parameters

A : double: m-by-n constraint matrix of the polytope {z : A*z <= b}. Requires m >= n+1 and a non-empty interior.
b : double: m-by-1 right-hand side.
x0 : double (optional): a strictly interior point (A*x0 < b). If empty, computed by sampleChebyshevCenter(A, b).
maxiter : double (optional): maximum interior-point iterations (default 150).
tol : double (optional): convergence tolerance (default 1e-6).
reg : double (optional): diagonal (Levenberg) regularisation of the Newton solves (default 1e-8).

Returns

x : double: n-by-1 ellipsoid centre.
E : double: n-by-n lower-triangular rounding transform with EE' = E2 (the SPD ellipsoid matrix). The rounding substitution is z = x + Ey.
converged : logical: whether tol was reached within maxiter. If false, the last iterate is returned as the centre and transform, which is a valid but looser rounding.

Notes

Analytic cases, tested in tSampling: a box maps to the unit ball (E = I); a sheared box {z : -1 <= Mz <= 1} gives EE' = inv(M)*inv(M)'; a triangle gives its Steiner inellipse.

Based on Y. Zhang & L. Gao (2003), SIAM J. Optim. 14:53-76; as used in Haraldsdottir et al. (2017) Bioinformatics 33:1741 (CHRR).

See also

sampleCHRR, sampleChebyshevCenter, randomSampling

sampleWarmupPoints#

function warmup = sampleWarmupPoints(model,varargin)

FVA warmup vertices for ACHR sampling.

Generates warmup points for sampleACHR by maximising and minimising each reaction in turn and storing the full flux distribution at each optimum. These vertices span the flux polytope and define the hit-and-run directions.

Parameters

model : struct: a model structure.

Name-Value Arguments

runParallel : logical: run the calculations in parallel, as in getAllowedBounds, which solves the same per-reaction minimisation and maximisation LPs (default true).

Returns

warmup : double: nPoints-by-nRxns matrix; each row is a vertex flux distribution. Duplicate rows are removed.

See also

sampleACHR, randomSampling, getAllowedBounds

traceFluxPath#

function [pathRxns, pathMets, cumFrac] = traceFluxPath(model, fluxes, fromRxn, toRxn, varargin)

Find the highest-flux-fraction path between two reactions.

Traces how flux flows from one reaction to another through the metabolic network, quantifying what fraction of the source reaction's flux reaches the target. At each metabolite junction the flux is split proportionally among all consuming reactions; the path with the highest cumulative fraction is returned, and printed if verbose is true.

By default the search excludes currency metabolites (for example ATP, NAD, H2O), which makes the returned path follow material flow. The setting traceMaterial=false disables this.

Parameters

model : struct: a model structure.
fluxes : double: flux vector (length = numel(model.rxns)).
fromRxn : char: source reaction ID.
toRxn : char: target reaction ID.

Name-Value Arguments

cutoff : double: minimum net flux through a metabolite (|S_ij * v_j|) for a consuming reaction to be considered (default 1e-8).
maxHops : double: paths of more than maxHops reactions are not extended further, so a returned path has at most maxHops + 1 reactions (default 30).
verbose : logical: print the path to the command window (default true).
traceMaterial : logical: exclude common currency and cofactor metabolites (ATP, ADP, NAD, NADH, H2O, CoA, Pi, CO2 and others); the path then follows material flux (default true).
carbonOnly : logical: additionally require each intermediate metabolite to contain at least one carbon atom, inferred from model.metFormulas. Drops H+, H2O, Pi, O2, NH3, CO2 that are not already caught by the currency list. Ignored with a warning when model.metFormulas is absent or empty (default false).
excludeMets : cell of char: extra metabolite names or IDs (after stripping compartment suffix) to exclude, in addition to the built-in currency list (which applies when traceMaterial is true) (default {}).

Returns

pathRxns : cell: cell array of reaction IDs on the best path, including fromRxn and toRxn. Empty if no path is found.
pathMets : cell: metabolite IDs at each junction (numel = numel(pathRxns)-1).
cumFrac : double: fraction of fromRxn's flux that passes through this exact path (0 = no path found, 1 = all flux reaches toRxn directly).

Examples

[p, m, f] = traceFluxPath(model, sol.x, 'PFK', 'ATPS4rpp')
% disable material filter to see all connections including ATP and ADP:
[p, m, f] = traceFluxPath(model, sol.x, 'PFK', 'ATPS4rpp', 'traceMaterial', false)
% add model-specific cofactors to the exclusion list:
[p, m, f] = traceFluxPath(model, sol.x, 'PFK', 'CS', 'excludeMets', {'acetyl-CoA'})

Notes

The algorithm follows metabolites in the forward direction only: it traces the metabolites that fromRxn produces and finds which consuming reactions eventually reach toRxn. If fromRxn and toRxn are not connected in the forward flux direction (for example they are in parallel branches), an empty path is returned. A best-first search, which expands the entry with the highest cumulative fraction first, finds the highest-fraction path; a path does not revisit a reaction.

walkFluxes#

function walkFluxes(model, fluxes, startRxn, varargin)

Interactively navigate a flux distribution reaction by reaction.

Starting from a chosen reaction, displays the flux-carrying reactions connected through shared metabolites, grouped by metabolite and labelled with each neighbour's role (produces or consumes). Entering the number of a neighbour steps to it. The visited reactions are kept in a history, and 'b' steps back through it.

Parameters

model : struct: a model structure.
fluxes : double: flux vector (length = numel(model.rxns)).
startRxn : char or double: starting reaction ID (string) or reaction index.

Name-Value Arguments

cutoff : double: minimum net flux through a metabolite (|S_ij * v_j|) for the metabolite to be listed for the current reaction and for a reaction to be displayed as a neighbour (default 1e-8).
maxPerMet : double: maximum neighbour reactions displayed per metabolite, selected by decreasing absolute flux (default 8).

Examples

sol = solveLP(model);
walkFluxes(model, sol.x, 'BIOMASS_Ec_iJO1366_core_53p95M')
walkFluxes(model, sol.x, 'PFK', 'cutoff', 1e-6)

Notes

Navigation commands: : step to that neighbour reaction b : go back to the previous reaction q : quit

This function is interactive and reads from the terminal at every step. It raises an error under -batch and in deployed code.