Steadycom

SteadyCom(modelCom, options, varargin)[source]

Find the maximum community growth rate at community steady-state using the SteadyCom algorithm

USAGE

[sol, result, LP, LPminNorm, indLP] = SteadyCom (modelCom, options, parameter, ‘param1’, value1, ‘param2’, value2, …)

INPUT

modelCom – A community COBRA model structure with the following fields (created using createMultipleSpeciesModel) (the first 5 fields are required, at least one of the last two is needed. Can be obtained using getMultiSpecisModelId):

  • S - Stoichiometric matrix

  • b - Right hand side

  • c - Objective coefficients

  • lb - Lower bounds

  • ub - Upper bounds

  • infoCom - structure containing community reaction info

  • indCom - the index structure corresponding to infoCom

OPTIONAL INPUTS

  • options – struct with the following possible fields: (for constraining individual growth rates and biomass amounts, default []):

    • GRfx - Fixed growth rate for organisms apart from the community (\(N_{organisms} * 1\) vector, NaN for unfixed growth rate, or [#organisms | value]) e.g. to fix organisms 2, 3 at growth rate 0.1, GRfx = [2, 0.1; 3, 0.1];

    • BMcon - Biomass constraint matrix \((\sum (a_{ij} * X_j) </=/> b_i)\) (given as \(K * N_{organisms}\) matrix for K constraints) e.g. [0 1 1 0] for \(X_2 + X_3\) in a 4-organisms model

    • BMrhs - RHS for BMcon, K x 1 vector for K constraints

    • BMcsense - Sense of the constraint, ‘L’, ‘E’, ‘G’ for <=, =, >= (for general constraints on e.g. total carbon uptake, molecular crowding, default [])

    • MC - \(K * (N_{rxns}+N_{organisms})\) coefficient matrix, for K additional constraints

    • MCmode - \(K * (N_{rxns}+N_{organisms})\) matrix , with number 0 ~ 3

      • 0: original variable

      • 1: positive part of the variable

      • 2: negative part of the variable

      • 3: absolute value of the variable

    • MCrhs - RHS of the constraints (default all zeros if .MC is given)

    • MClhs - LHS of the constraints (default -inf if .MC is given) (parameters in the iterative algorithm, [default value])

    • GRguess [0.2] - Initial guess of the growth rate.

    • feasCrit [1] - Criteria for feasibility, 1 or 2: The algorithm tests iteratively at a given growth rate whether a feasible solution can be found.

      1. Use a threshold total biomass BMweight (see below). i.e. \(\sum X \geq BMweight\) (use it if the total biomass is known, the most common usage)

      2. Use a threshold on minimum biomass production (=specific growth rate x \(\sum biomass\), which is roughly constant over a range of growth rate if the sum of biomass is not bounded above) i.e. \(\sum X * gr \geq BMtol * BMref * GR0\) where BMref is the maximum biomass at a small growth rate GR0 and BMtol is a fraction ranging from 0 to 1

    • algorithm [1] - Algorithm to find the maximum growth rate

      1. Fzero after finding grLB and grUB with simple guessing [\(gr^T = gr * \sum X / \sum X^T\)]

      2. Simple guessing with minimum one percent step size

      3. Bisection method

    • BMweight [1] - Minimum total biomass for feasibility. Used only if feasCrit = 1. Set BMweight to a close-to-zero value to compute the wash-out dilution rate.

    • GR0 [0.001] - A small growth rate to obtain a reference value for maximum total biomass production. Used only if feasCrit = 2 or solveGR0 = true

    • BMtol [0.8] - Fractional tolerance for biomass production to check feasibility. Used only if feasCrit = 2

    • solveGR0[false] - true to solve the model at a low growth rate GR0 first to test feasibility

    • GRtol [1e-6] - Precision for the growth rate found (\(grUB - grLB < GRtol\))

    • BMtolAbs [1e-5] - Absolute tolerance for positivity of biomass

    • maxIter (1e3) - maximum nummber of iteration (parameters in the optimization model, [default value])

    • minNorm [0] - 0: No minNorm. 1: min sum of absolution flux of the final solution.

    • BMgdw [all 1s] - The gram dry weight per mmol of the biomass reaction of each organism. Maybe used to scale the biomass reactions between organisms.

    • BMobj [all 1s] - Objective coefficient for the biomass of each organism when doing the maximization at each step. (other parameters)

    • verbFlag [3] - Print level. 0, 1, 2, 3 for silence, one log per 10, 5 (default) or 1 iteration respectively

    • LPonly [false] - Return the initial LP at zero growth rate only. Calculate nothing.

    • saveModel [‘’] String, if non-empty, save the LP structure.

  • parameter – structure for solver-specific parameters. ‘param1’, value1, …: name-value pairs for solveCobraLP parameters. See solveCobraLP for details

OUTPUTS

  • sol – COBRA solution structure (Cplex structure if using ibm_cplex)

  • result – structure with the following fields:

    • GRmax: maximum specific growth rate found (/h)

    • vBM: biomass formation rate (gdw/h)

    • BM: Biomass vector at GRmax (gdw)

    • Ut: uptake fluxes (mmol/h)

    • Ex: export fluxes (mmol/h)

    • flux: flux distribution for the original model (the following ‘iter’ fields are status in each iteration:) [GR | biomass X | biomass flux (GR * X) | max. infeas. of solution])

    • iter0: stationary, no growth, gr = 0

    • iter: iterations for finding max gr

    • stat: status at the termination of the algorithm:

      • optimal: optimal growth rate found

      • maintenance: feasible at maintenance, but cannot grow

      • minimal growth: feasible at a minimal growth rate (possible only if options.solveGR0 = true)

      • infeasible: infeasible model, even with maintenance requirement only

      • LPonly: return the LP structure only. No optimization performed (only if options.LPonly = true)

      • xxx (minNorm L1-norm): in result.flux the sum of absolute fluxes is minimized. ‘xxx’ is one of the status above.

SteadyComFVA(modelCom, options, varargin)[source]

Flux variability analysis for community model at community steady-state for a range of growth rates. The function is capable of saving intermediate results and continuing from previous results if the file path is given in options.saveFVA. It also allows switch from single thread to parallel computation from intermediate results (but not the other way round).

USAGE

[minFlux, maxFlux, minFD, maxFD, GRvector, result, LP] = SteadyComFVA (modelCom, options, parameter, ‘param1’, value1, ‘param2’, value2, …)

INPUT

modelCom – A community COBRA model structure with the following fields (created using createMultipleSpeciesModel) (the first 5 fields are required, at least one of the last two is needed. Can be obtained using getMultiSpecisModelId):

  • S - Stoichiometric matrix

  • b - Right hand side

  • c - Objective coefficients

  • lb - Lower bounds

  • ub - Upper bounds

  • infoCom - structure containing community reaction info

  • indCom - the index structure corresponding to infoCom

OPTIONAL INPUTS

  • options – struct with the following possible fields:

    • optGRpercent - A vector of percentages. Perform FVA at these percents of max. growth rate (Default = [99.99])

    • optBMpercent - Only consider solutions that yield at least a certain percentage of the optimal biomass (Default = 99.99)

    • rxnNameList - List of reactions (index row vector or subset of *.rxns) for which FVA is performed. (Default = biomass reaction of each species) Or a \((N_{rxns} + N_{organism}) x K\) matrix for FVA of K linear combinations of fluxes and/or abundances e.g., [1; -2; 0] for finding the max/min of \(1 v_1 - 2 v_2 + 0 v_3\)

    • rxnFluxList - List of reactions (index vector or subset of *.rxns) whose fluxes are also returned along with the FVA result of each entry in rxnNameList (Default = biomass reaction of each species)

    • GRmax - maximum growth rate of the model (default to be found SteadyCom.m) (the two parameters below are usually determined by solving the problem during the program. Provide them only if you want to constrain the total biomass to a particular value)

    • BMmaxLB - lower bound for the total biomass (default 1)

    • BMmaxUB - upper bound for the total biomass (other parameters below)

    • saveFVA - If non-empty, become the filename to save the FVA results (default empty, not saving)

    • saveFre - save frequency. Save every (#rxns for FVA) * saveFre (default 0.1)

    • threads - for parallelization: > 1 for explicitly stating the no. of threads used, 0 or -1 for using all available threads. Default 1. (Requires Matlab parallel toolbox)

    • verbFlag - Verbose output. 1 to have waitbar, >1 to have stepwise output (default 3)

    • loadModel - (ibm_cplex only) String of filename to be loaded. If non-empty, load the cplex model (‘loadModel.mps’), basis (‘loadModel.bas’) and parameters (‘loadModel.prm’). (May add also other parameters in SteadyCom for calculating the maximum growth rate.)

  • parameter – structure for solver-specific parameters. ‘param1’, value1, … name-value pairs for solveCobraLP parameters. See solveCobraLP for details

OUTPUTS

  • minFlux – Minimum flux for each reaction

  • maxFlux – Maximum flux for each reaction

OPTIONAL OUTPUTS:
minFD: \(rxnFluxList * rxnNameList\) matrix containing the fluxes in options.rxnFluxList

corresponding to minimizing each reaction in options.rxnNameList

maxFD: \(rxnFluxList * rxnNameList\) matrix containing the fluxes in options.rxnFluxList

corresponding to maximizing each reaction in options.rxnNameList

GRvector: a vector of growth rates at which FVA has been performed result: result structure from SteadyCom LP: LP problem structure (Cplex LP object for ibm_cplex)

SteadyComPOA(modelCom, options, varargin)[source]

Analyze pairwise relationship between reaction fluxes/biomass variables for a community model at community steady-state at a given growth rate. See tutorial_SteadyCom for more details.

USAGE

[POAtable, fluxRange, Stat, GRvector] = SteadyComPOA (modelCom, options, parameters, ‘param1’, value1, ‘param2’, value2, …)

INPUT

modelCom – A community COBRA model structure with the following fields (created using createMultipleSpeciesModel) (the first 5 fields are required, at least one of the last two is needed. Can be obtained using getMultiSpecisModelId):

  • S - Stoichiometric matrix

  • b - Right hand side

  • c - Objective coefficients

  • lb - Lower bounds

  • ub - Upper bounds

  • infoCom - structure containing community reaction info

  • indCom - the index structure corresponding to infoCom

OPTIONAL INPUTS

  • options – option structure with the following fields:

    • GRmax - maximum growth rate of the model (default to be found SteadyCom.m)

    • optGRpercent - A vector of percentages. Perform FVA at these percents of max. growth rate (Default = [99.99])

    • optBMpercent - Only consider solutions that yield at least a certain percentage of the optimal biomass (Default = 99.99)

    • rxnNameList - list of reactions (IDs or .rxns) to be analyzed. Use a \((N_{rxns} + N_{organism}) x K\) matrix for POA of K linear combinations of fluxes and/or abundances (Default = biomass reaction of each organism,

    • pairList - pairs in rxnNameList to be analyzed. N_pair by 2 array of:

        • indices referring to the rxns in rxnNameList, e.g., [1 2] to analyze rxnNameList{1} vs rxnNameList{2}

        • rxn names which are members of rxnNameList, e.g., {‘EX_glc-D(e)’, ‘EX_ac(e)’}

      If not supplied, analyze all K(K-1) pairs from the K targets in rxnNameList.

    • symmetric - true to avoid running symmetric pairs (e.g. analyze pair (j,k) only if \(j > k\), total \(K(K-1)/2 pairs)\). Used only when pairList is not supplied.

    • Nstep - number of steps for fixing one flux at a value between the min. and the max. possible fluxes. Default 10. Can also be a vector indicating the fraction of intermediate value to be analyzed, e.g. [0 0.5 1] means computing at minFlux, 0.5(minFlux + maxFlux) and maxFlux

    • NstepScale - used only when Nstep is a single number.

      • -‘lin’ for a linear (uniform) scale of step size

      • -‘log’ for a log scaling of the step sizes

    • fluxRange - flux range for each entry in rxnNameList. K x 2 matrix. Defaulted to be found by SteadyComFVA.m (other parameters)

    • savePOA - filename to save the POA results (default ‘POAtmp/POA’). Must be non-empty. New folder is recommended

    • threads - for parallelization: > 1 for explicitly stating the no. of threads used, 0 or -1 for using all available threads. Default 1.

    • verbFlag - verbose output. 0 or 1.

    • loadModel - (ibm_cplex only) string of filename to be loaded. If non-empty, load the cplex model (‘loadModel.mps’), basis (‘loadModel.bas’) and parameters (‘loadModel.prm’). (May add also other parameters in SteadyCom for calculating the maximum growth rate.)

  • parameter – structure for solver-specific parameters. ‘param1’, value1, …: name-value pairs for solveCobraLP parameters. See solveCobraLP for details

OUTPUTS

  • POAtable – K x K cells. Each (i, i) -cell contains a Nstep x 1 x N_gr matrix of the fluxes at which rxnNameList{i} is fixed. Each (i, j) -cell contains a Nstep x 2 x N_gr matrix, the (p,:,q) -entry being the range of rxnNameList{j} when rxnNameList{i} is fixed at POAtable{i, i}(p, 1, q) at growth rate = GRvector(q)

  • fluxRange – K x 2 x N_gr matrix of flux range for each entry in rxnNameList

  • Stat – K x K x N_gr structure array with fields:

    • -‘cor’: the slope from linear regression between the fluxes of a pair

    • -‘r2’: the corresponding coefficient of determination (R-square)

  • GRvector – vector of growth rates being analyzed