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 most significant flux changes between two conditions.
`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 optimization 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`	Simple OptKnock for growth-coupled production.
`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. This score could be the Student t between the two conditions, or you can calculate it from a p-value (by computing the inverse of the so called error function). If you choose the second alternative you should be aware that the transcripts that increased in expression level should have positive values and those who decreased in expression level should have negative values (the p-values only tell you if the fluxes changed or not but not in which direction).
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: prints the most significant reactions 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 most significant flux changes between two conditions.

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: RAVEN model structure.
fluxes1 : double: reference flux vector (condition A, e.g. wild-type).
fluxes2 : double: comparison flux vector (condition B, e.g. mutant or knockin).

Name-Value Arguments

cutoff : double: minimum |flux| to consider a reaction active (default 1e-8).
metaboliteList : cell: cell array of metabolite names. Only reactions involving at least one of these metabolites are compared; every other reaction is left out of the result entirely. Matching is against model.metNames and is case insensitive. A name that matches no metabolite is reported as a warning rather than an error (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: .turnedOn — cell, rxn IDs newly active in fluxes2 .turnedOff — cell, rxn IDs silenced in fluxes2 .flipped — cell, rxn IDs that reversed sign between conditions .changed — struct with fields (all reactions with |dflux|>cutoff, sorted by |dflux| descending): .rxn — cell of reaction IDs .flux1 — double, flux in condition A .flux2 — double, flux in condition B .absDelta — double, |flux2 - flux1| .relChange — double, (flux2-flux1)/|flux1|; NaN when flux1 is near-zero (reaction emerged) .type — cell, 'on'|'off'|'flip'|'' per 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 genes, optimizes the model by FBA, and keeps track of the resulting fluxes. This is used for identifying gene deletion targets.

Parameters

model : struct: a model structure.

Name-Value Arguments

testType : char

single/double 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 optimization (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

simply the genes in the input model. Included for simple presentation of the output.

details : double

not all genes will be deleted in all analyses. It is for example not necessary to delete genes for dead end reactions. This is a vector with details about each gene in originalGenes and why or why not it was deleted:

1 : was 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, no longer used
4 : involved in dead-end reaction

grRatioMuts : double

growth rate ratio between mutated strain and wild type, matching the originalGenes(genes) mutants. Note that this does not directly map to model.genes, as is the case for COBRA getEssentialGenes. However, this can be obtained by afterwards running:

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 has been specified (targets), it will also generate a structure indicating for each model reaction whether it is identified by FSEOF 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 iterations (default 10).
coefficient : double: ratio of optimal target reaction flux, must be less than 1 (default 0.9).
outputFile : char: output filename (default prints to command window).
corrThreshold : double: minimum absolute Pearson correlation coefficient between a reaction absolute flux and the enforced product flux for the reaction to be selected as a target. Values in [0, 1]; higher values require a more linear response and suppress LP alternative- optima noise (default 0.9).

Returns

targets : struct

structure with information for identified targets, with fields:

logical : logical array, true for amplification targets (backward compatible alias for amplify).
amplify : logical array, reactions whose absolute flux increases with enforced product flux.
knockdown : logical array, reactions driven toward but not to zero.
knockout : logical array, reactions driven to zero at max product.
slope : numeric array with per-reaction regression slopes of absolute flux versus enforced product flux (all reactions, not just targets).

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: speed up calculations by parallel processing. This is not beneficial if allowed bounds are calculated for only a few reactions, as the overhead of parallel processing will take longer. It requires MATLAB Parallel Computing Toolbox. If this is not installed, the calculations will not be parallelized, regardless of what is indicated as runParallel (default true).

Returns

minFluxes : double: minimal allowed fluxes.
maxFluxes : double: maximal allowed fluxes.
exitFlags : double: exit flags for min/max for each of the reactions. True if it was possible to calculate a flux.

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 boolean matrix where the rows correspond to the metabolites and the columns to which subgraph they are assigned to. The columns are ordered so that larger subgraphs come first.

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 Z-scores that tells you for each reaction how likely it is for its flux to have increased (positive sign) or decreased (negative sign) in the second condition with respect to the first.

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 still allowing the model to achieve a target growth rate. Unlike iterative heuristics, the MILP formulation finds the globally minimal medium.

Parameters

model : struct: RAVEN model structure.

Name-Value Arguments

minGrowth : double: minimum objective value (growth rate) the medium must support. Defaults to 10 % of the unconstrained optimum computed automatically. Set explicitly to avoid an extra FBA call (e.g. 'minGrowth', 0.1).
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 considered as candidate uptake reactions. Exchange reactions that are always secretion-only (lb >= 0) or that are completely blocked 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 minimal number of uptakes is guaranteed, but the medium itself need not be 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 (with a RAVEN:glpkMILP warning); for larger models Gurobi is much faster and more reliable (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. Is used as a faster version of getAllowedBounds if it is only interesting whether the reactions can carry a flux or not.

Parameters

model : struct: a model structure.

Name-Value Arguments

cutOff : double: the flux value that a reaction has to carry to be identified as positive (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 metabolites in the model, or a vector of indexes (default model.rxns).
seed : double: seed for the random number generator used to shuffle the order in which reactions are tested, so that order is reproducible across runs (default: 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

If a model has +/- Inf bounds then those are replaced with an arbitary large value of +/- 10000 prior to 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, for better mixing on thin/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 maximises a small random objective, returning a polytope vertex.

The 'achr'/'chrr' methods draw the (near-)uniform interior distribution; the 'randomObjective' method draws diverse vertices. Arguments below marked [randomObjective] apply only to that method; 'thinning'/'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'/'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. This is needed to get solutions without loops if your model has, for example, 1000/-1000 as arbitrarily large bounds. If your model only has "biologically relevant" bounds, then set this to false (default true).
supressErrors : logical: the program will halt if it has problems finding non-zero solutions which are not involved in loops. This could be because the constraints on the model are too relaxed (such as unlimited glucose uptake) or too strict (such as too many and too narrow constraints) (default false).
runParallel : logical: speed up calculations by parallel processing. Requires the MATLAB Parallel Computing Toolbox. If this is not installed, the calculations will not be parallelized, regardless of what is indicated as runParallel (default true).
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 optimization should be performed for each random sample, to minimize the number of fluxes and thereby prevent loops. Typically, loops are averaged out when a large number of samples are taken, but this is not always the case (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 maximizing (with random weights) for a random set of nObjectives reactions (default 2). For reversible reactions it randomly chooses between maximizing and minimizing. 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 for sampling, instead focusing on sampling from the metabolic aspects that form the solution space.

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);
% Historical 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.

The Reporter Metabolites algorithm for identifying metabolites around which transcriptional changes occur.

Parameters

model : struct: a model structure.
genes : cell: a cell array of gene names (should match with model.genes).
genePValues : double: P-values for differential expression of the genes.

Name-Value Arguments

printResults : logical: true if the top 20 Reporter Metabolites should be printed 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, then Reporter Metabolites are calculated for only up/down-regulated genes 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 two latter structures are only calculated if geneFoldChanges are supplied.
mets : a cell array of metabolite IDs for the metabolites for which a score could be calculated.
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 optimization 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 (e.g. 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 (e.g. 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 (i.e. model.lb < 0) the concentration is assumed to be high enough to not be limiting. If the uptake rate for a nutrient is calculated to exceed the maximum uptake rate for that nutrient specified in the model and the max uptake rate specified is > 0, the maximum uptake rate specified in the model is used instead of the calculated uptake rate.

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 points in the plot (default 20).

Returns

biomassValues : double: biomass values for plotting.
targetValues : double: target upper and lower bounds for plotting.

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.

Performs robustness analysis for a reaction of interest and an objective of interest.

Parameters

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

Name-Value Arguments

nPoints : double: number of points to show on plot (default 20).
objRxn : char: reaction identifier of objective to be maximized (default uses the objective defined in the model).
plotRedCost : logical: whether reduced cost should also be plotted (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)

Simple OptKnock for growth-coupled production.

Simple OptKnock algorithm that checks all gene or reaction deletions for growth-coupled metabolite production, by testing all possible combinations. This is not defined as MILP, and is therefore slow (but simple).

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.

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 : S*v = 0, lb <= v <= ub}, using the ACHR algorithm (Kaufman & Smith 1998). No rounding is performed, so on strongly elongated polytopes (e.g. enzyme-constrained models) mixing is slower than sampleCHRR; prefer CHRR there.

Parameters

model : struct: a RAVEN 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, so mixing stays fast even on strongly elongated polytopes (e.g. enzyme-constrained models), unlike sampleACHR.

Set any constraints you want to condition on (e.g. a biomass lower bound or measured exchange fluxes) on the model before calling.

Parameters

model : struct: a RAVEN 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 folded into the equality system, keeping 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 folded in 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 mixes rapidly.

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 keeping the Newton solves well-conditioned (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. A false value is not fatal — the last iterate is still a valid (if looser) rounding.

Notes

Validated against analytic cases: 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 seed the hit-and-run directions.

Parameters

model : struct: a RAVEN model structure.

Name-Value Arguments

runParallel : logical: speed up calculations by parallel processing, as in getAllowedBounds, which solves the same per-reaction min/max 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.

By default the search excludes currency metabolites (ATP, NAD, H2O, etc.) so the returned path represents material (carbon-skeleton) flow rather than cofactor shortcuts. Set traceMaterial=false to disable this.

Parameters

model : struct: RAVEN 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 |flux| for a reaction to be considered flux-carrying (default 1e-8).
maxHops : double: maximum number of reactions in the path; longer paths are pruned (default 30).
verbose : logical: print the path to the command window (default true).
traceMaterial : logical: exclude common currency/cofactor metabolites (ATP, ADP, NAD, NADH, H2O, CoA, Pi, CO2, …) so the path follows material (carbon-skeleton) flux rather than energy-carrier shortcuts (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. Silently ignored when model.metFormulas is absent (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 (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 routes 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/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 (e.g. they are in parallel branches), an empty path is returned. A best-first (greedy) search ensures the highest-fraction path is found efficiently; a visited set prevents revisiting reactions.

walkFluxes#

function walkFluxes(model, fluxes, startRxn, varargin)

Interactively navigate a flux distribution reaction by reaction.

Starting from a chosen reaction, displays all flux-carrying reactions connected through shared metabolites, grouped by metabolite and labelled with each neighbour's role (produces / consumes). Select a neighbour by number to step to it; build a history and rewind with 'b'.

Parameters

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

Name-Value Arguments

cutoff : double: minimum |flux| to display a reaction as a neighbour (default 1e-8).
maxPerMet : double: maximum neighbour reactions displayed per metabolite (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 errors immediately under -batch, in deployed code, or otherwise without a real terminal, rather than blocking forever waiting for input that will never arrive.