Analysis¶

MATLAB functions in RAVEN/analysis of the RAVEN toolbox. Help text is collected from the source of the tracked branch.

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.
`followChanged`	Print reactions whose fluxes differ from a reference case.
`followFluxes`	Print reactions with fluxes in a specified interval.
`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.
`getMinNrFluxes`	Find the minimal set of fluxes that satisfy the model.
`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¶

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.

Input arguments:

Name	Type	Description	Default
`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).	required
`df`	`double`	the degrees of freedom in the t-test.	required
`solutionsA`	`double`	random solutions for the reference condition (as generated by randomSampling).	required
`solutionsB`	`double`	random solutions for the test condition (as generated by randomSampling).	required

Name-value arguments:

Name	Type	Description	Default
`printResults`	`logical`	prints the most significant reactions in each category (default false).

Output arguments:

Name	Type	Description
`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¶

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.

Input arguments:

Name	Type	Description	Default
`model`	`struct`	RAVEN model structure.	required
`fluxes1`	`double`	reference flux vector (condition A, e.g. wild-type).	required
`fluxes2`	`double`	comparison flux vector (condition B, e.g. mutant or knockin).	required

Name-value arguments:

Name	Type	Description
`cutoff`	`double`	minimum \|flux\| to consider a reaction active (default 1e-8).
`nMax`	`double`	maximum rows printed in the summary table (default 20).
`verbose`	`logical`	print the summary table (default true).

Output arguments:

Name	Type	Description
`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)

See also walkFluxes, traceFluxPath, modelSummary

findGeneDeletions¶

Delete genes and track the resulting fluxes.

Deletes genes, optimizes the model, and keeps track of the resulting fluxes. This is used for identifying gene deletion targets.

Input arguments:

Name	Type	Description	Default
`model`	`struct`	a model structure.	required

Name-value arguments:

Name	Type	Description
`testType`	`char`	single/double gene deletions/over expressions. Over expression is only available if using MOMA (default "sgd"): "sgd" : single gene deletion "dgd" : double gene deletion "sgo" : single gene over expression "dgo" : double gene over expression
`analysisType`	`char`	determines whether to use FBA ("fba") or MOMA ("moma") in the optimization (default "fba").
`refModel`	`struct`	MOMA works by fitting the flux distributions of two models to be as similar as possible. The most common application is where there is a reference model with some fluxes constrained from experimental data. This model is required when using MOMA.
`oeFactor`	`double`	a factor by which the fluxes should be increased if a gene is overexpressed (default 10).

Output arguments:

Name	Type	Description
`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/overexpressed 2 : proved lethal in sgd (single gene deletion) 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,analysisType,refModel,oeFactor);

followChanged¶

Print reactions whose fluxes differ from a reference case.

Prints fluxes and reactions for each of the reactions that result in different fluxes compared to the reference case.

Input arguments:

Name	Type	Description	Default
`model`	`struct`	a model structure.	required
`fluxesA`	`double`	flux vector for the test case.	required
`fluxesB`	`double`	flux vector for the reference test.	required

Name-value arguments:

Name	Type	Description
`cutOffChange`	`double`	reactions where the fluxes differ by less than this many percent won't be printed (default 10^-8).
`cutOffFlux`	`double`	reactions where the absolute value of both fluxes are below this value won't be printed (default 10^-8).
`cutOffDiff`	`double`	reactions where the fluxes differ by less than cutOffDiff won't be printed (default 10^-8).
`metaboliteList`	`cell`	cell array of metabolite names. Only reactions involving any of these metabolites will be printed.

Examples:

followChanged(model,fluxesA,fluxesB,cutOffChange,cutOffFlux,...
    cutOffDiff,metaboliteList);

followFluxes¶

Print reactions with fluxes in a specified interval.

Prints fluxes and reactions for each of the reactions that result in fluxes within the specified interval.

Input arguments:

Name	Type	Description	Default
`model`	`struct`	a model structure.	required
`fluxesA`	`double`	flux vector for the test case.	required
`lowerFlux`	`double`	only reactions with fluxes above this cutoff value are displayed.	required

Name-value arguments:

Name	Type	Description	Default
`upperFlux`	`double`	only reactions with fluxes below this cutoff value are displayed (default Inf).
`fluxesB`	`double`	flux vector for the reference case.

Output arguments:

Name	Type	Description
`errorFlag`	`double`	set to 1 if upperFlux is not larger than lowerFlux, otherwise empty.

Examples:

errorFlag=followFluxes(model,fluxesA,lowerFlux,upperFlux,fluxesB);

FSEOF¶

FSEOF 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.

Input arguments:

Name	Type	Description	Default
`model`	`struct`	a model structure.	required
`biomassRxn`	`char`	reaction ID of the biomass formation or growth reaction.	required
`targetRxn`	`char`	reaction ID of the target reaction.	required

Name-value arguments:

Name	Type	Description
`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).

Output arguments:

Name	Type	Description
`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¶

Return the minimal and maximal fluxes through reactions.

Input arguments:

Name	Type	Description	Default
`model`	`struct`	a model structure.	required

Name-value arguments:

Name	Type	Description	Default
`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).

Output arguments:

Name	Type	Description
`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¶

Get all metabolic subgraphs in a model.

Two metabolites are connected if they share a reaction.

Input arguments:

Name	Type	Description	Default
`model`	`struct`	a model structure.	required

Output arguments:

Name	Type	Description
`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¶

Calculate the essential reactions for a solvable model.

Input arguments:

Name	Type	Description	Default
`model`	`struct`	a model structure.	required

Name-value arguments:

Name	Type	Description	Default
`ignoreRxns`	`cell`	cell array of reaction IDs which should not be checked (default {}).

Output arguments:

Name	Type	Description
`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¶

Calculate Z scores between two sets of random flux distributions.

Input arguments:

Name	Type	Description	Default
`solutionsA`	`double`	random solutions for the reference condition (as generated by randomSampling).	required
`solutionsB`	`double`	random solutions for the test condition (as generated by randomSampling).	required

Output arguments:

Name	Type	Description
`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¶

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.

Input arguments:

Name	Type	Description	Default
`model`	`struct`	RAVEN model structure.	required

Name-value arguments:

Name	Type	Description	Default
`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).

Output arguments:

Name	Type	Description
`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.

getMinNrFluxes¶

Find the minimal set of fluxes that satisfy the model.

Uses mixed integer linear programming to find the minimal set of fluxes that satisfy the model.

Input arguments:

Name	Type	Description	Default
`model`	`struct`	a model structure.	required

Name-value arguments:

Name	Type	Description
`toMinimize`	`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 indexes for the reactions that should be minimized (default model.rxns).
`params`	`struct`	obsolete option.
`scores`	`double`	vector of weights for the reactions. Negative scores should not have flux. Positive scores are not possible in this implementation, and they are changed to max(scores(scores<0)). Must have the same dimension as toMinimize (find(toMinimize) if it is a logical vector) (default -1 for all reactions).

Output arguments:

Name	Type	Description
`x`	`double`	the corresponding fluxes for the full model.
`I`	`double`	the indexes of the reactions in toMinimize that were used in the solution.
`exitFlag`	`double`	exit status: 1 : optimal solution found -1 : no feasible solution found -2 : optimization time out

Examples:

[x, I, exitFlag] = getMinNrFluxes(model, toMinimize, params, scores);

Notes

Uses 1000 mmol/gDW/h as an arbitary large flux. Could possibly cause problems if the fluxes in the model are larger than that.

haveFlux¶

Check which reactions can carry a flux.

Checks which reactions can carry a (positive or negative) flux. Is used as a faster version of getAllowedBounds if it is only interesting whether the reactions can carry a flux or not.

Input arguments:

Name	Type	Description	Default
`model`	`struct`	a model structure.	required

Name-value arguments:

Name	Type	Description	Default
`cutOff`	`double`	the flux value that a reaction has to carry to be identified as positive (default 10^-8).
`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).

Output arguments:

Name	Type	Description
`I`	`logical`	logical array with true if the corresponding reaction can carry a flux.

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.

randomSampling¶

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. This was randomSampling's historical behaviour; it is no longer the default.

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.

Input arguments:

Name	Type	Description	Default
`model`	`struct`	a model structure.	required

Name-value arguments:

Name	Type	Description
`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 in order 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).

Output arguments:

Name	Type	Description
`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.

Notes

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.

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¶

Identify metabolites around which transcriptional changes occur.

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

Input arguments:

Name	Type	Description	Default
`model`	`struct`	a model structure.	required
`genes`	`cell`	a cell array of gene names (should match with model.genes).	required
`genePValues`	`double`	P-values for differential expression of the genes.	required

Name-value arguments:

Name	Type	Description
`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).

Output arguments:

Name	Type	Description
`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 : Z-scores for differential expression around each metabolite in "mets". metPValues : P-values for differential expression around each metabolite in "mets". 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.

runDynamicFBA¶

Perform dynamic FBA using the static optimization approach.

Input arguments:

Name	Type	Description	Default
`model`	`struct`	a model structure.	required
`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).	required
`initConcentrations`	`double`	initial concentrations of substrates (matching substrateRxns).	required
`initBiomass`	`double`	initial biomass (must be non-zero).	required
`timeStep`	`double`	time step size.	required
`nSteps`	`double`	maximum number of time steps.	required
`plotRxns`	`cell`	cell array with exchange reaction identifiers for substrates whose concentration should be plotted.	required
`exclUptakeRxns`	`cell`	cell array with exchange reaction identifiers for substrates whose concentration does not change (e.g. co2, o2, h2o, h).	required

Output arguments:

Name	Type	Description
`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¶

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.

Input arguments:

Name	Type	Description	Default
`model`	`struct`	a model structure.	required
`controlRxn1`	`char`	reaction identifier of the first reaction to be plotted.	required
`controlRxn2`	`char`	reaction identifier of the second reaction to be plotted.	required

Name-value arguments:

Name	Type	Description
`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).

Output arguments:

Name	Type	Description
`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¶

Calculate the byproduct secretion envelope.

Input arguments:

Name	Type	Description	Default
`model`	`struct`	a model structure.	required
`targetRxn`	`char`	identifier of target metabolite production reaction.	required
`biomassRxn`	`char`	identifier of biomass reaction.	required

Name-value arguments:

Name	Type	Description	Default
`nPts`	`double`	number of points in the plot (default 20).

Output arguments:

Name	Type	Description
`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¶

Perform robustness analysis for a reaction and objective.

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

Input arguments:

Name	Type	Description	Default
`model`	`struct`	a model structure.	required
`controlRxn`	`char`	reaction of interest whose value is to be controlled.	required

Name-value arguments:

Name	Type	Description
`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).

Output arguments:

Name	Type	Description
`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¶

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

Input arguments:

Name	Type	Description	Default
`model`	`struct`	a model structure.	required
`targetRxn`	`char`	identifier of target reaction.	required
`biomassRxn`	`char`	identifier of biomass reaction.	required

Name-value arguments:

Name	Type	Description
`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).

Output arguments:

Name	Type	Description
`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¶

Artificially Centered Hit-and-Run flux sampling.

Draws nSamples flux vectors from the flux polytope {v : S*v = 0, lb <= v <= ub} using ACHR (Kaufman & Smith 1998): from the current point, step towards a randomly chosen warmup vertex through the running centre of all points so far, sampling uniformly along the feasible chord. No rounding is performed, so on strongly elongated polytopes (e.g. enzyme-constrained models) mixing is slower than sampleCHRR; prefer CHRR there.

Because every search direction is a difference of feasible points (which all satisfy S*v = 0), each step stays exactly on the steady-state manifold.

Input arguments:

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

Output arguments:

Name	Type	Description
`solutions`	`double`	nRxns-by-nSamples matrix of sampled flux distributions.

See also

randomSampling, sampleCHRR, sampleWarmupPoints

sampleChebyshevCenter¶

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 : A*x <= b}. Used to initialise the interior-point MVE solver (sampleMaxVolEllipse) inside sampleCHRR.

Input arguments:

Name	Type	Description	Default
`A`	`double`	m-by-n constraint matrix.	required
`b`	`double`	m-by-1 right-hand side.	required

Output arguments:

Name	Type	Description
`x0`	`double`	n-by-1 Chebyshev centre (strictly interior point).

See also

sampleMaxVolEllipse, sampleCHRR

sampleCHRR¶

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

Reduce the polytope to a full-dimensional body via the nullspace of S (v = v0 + N*x); reactions with ~zero flux range are folded into the equality system so the reduced polytope is full-dimensional.
Round it with the maximum-volume inscribed ellipsoid (sampleMaxVolEllipse). Rounding makes mixing independent of how elongated the polytope is, which is the regime of enzyme-constrained (ecModel + proteomics) and flux-measured models.
Walk the rounded polytope with coordinate hit-and-run and map back.

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

Input arguments:

Name	Type	Description	Default
`model`	`struct`	a RAVEN model structure.	required
`nSamples`	`double`	number of flux vectors to return (default 1000).	required
`thinning`	`double`	Markov-chain steps between recorded samples (default 100).	required
`nBurnin`	`double`	burn-in steps discarded before the first recorded sample (default 1000).	required
`seed`	`double`	random seed for reproducibility (default: current RNG state).	required
`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).	required

Output arguments:

Name	Type	Description
`solutions`	`double`	nRxns-by-nSamples matrix of sampled flux distributions.
`info`	`struct with fields:`	.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¶

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.

Input arguments:

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

Output arguments:

Name	Type	Description
`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¶

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.

Input arguments:

Name	Type	Description	Default
`model`	`struct`	a RAVEN model structure.	required

Output arguments:

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

See also

sampleACHR, randomSampling

traceFluxPath¶

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.

Input arguments:

Name	Type	Description	Default
`model`	`struct`	RAVEN model structure.	required
`fluxes`	`double`	flux vector (length = numel(model.rxns)).	required
`fromRxn`	`char`	source reaction ID.	required
`toRxn`	`char`	target reaction ID.	required

Name-value arguments:

Name	Type	Description
`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 {}).

Output arguments:

Name	Type	Description
`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¶

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'.

Input arguments:

Name	Type	Description	Default
`model`	`struct`	RAVEN model structure.	required
`fluxes`	`double`	flux vector (length = numel(model.rxns)).	required
`startRxn`	`char or double`	starting reaction ID (string) or column index.	required

Name-value arguments:

Name	Type	Description	Default
`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