Simulating growth with FBA#

Flux balance analysis (FBA) sets an objective, constrains the uptake rates, solves and reads the fluxes back. Every other analysis on this site is built from this loop.

Functions#

MATLAB

Python

Purpose

setParam

Reaction.bounds, Model.objective cobrapy

set bounds and the objective

solveLP

Model.optimize cobrapy

solve the LP

printFluxes

Model.summary cobrapy

show the fluxes that carry material

solveLP (minFlux)

pfba cobrapy

pick a parsimonious solution among the optima

Every Python simulation step on this page is cobrapy, marked cobrapy in the functions table. In MATLAB, solveLP needs neither the COBRA Toolbox nor anything else outside RAVEN.

Setup#

The model is yeast-GEM.yml from docs/data/, yeast-GEM v9.1.0, which has a growth objective and an aerobic glucose medium set.

Solving#

model = readYAMLmodel('yeast-GEM.yml');
sol = solveLP(model);
fprintf('objective: %s\n', model.rxns{model.c == 1});
fprintf('status:    %d\n', sol.stat);
fprintf('growth:    %.4f /h\n', sol.f);

objective: r_2111
status:    1
growth:    0.0809 /h

The objective is model.c, a vector with one entry per reaction, so model.c == 1 finds the reaction being maximised.

sol holds the whole answer. sol.f is the objective value with its natural sign: RAVEN minimises internally and negates the result, so no sign has to be undone by the caller. sol.x is the flux vector, in the order of model.rxns, and sol.stat reports how the solve ended:

stat

Meaning

1

solved to optimality

0

a feasible solution was found, but not proven optimal

-1

no feasible solution exists

-2

solved, but the second optimisation set by minFlux failed

Anything other than 1 makes the fluxes provisional at best, so check stat before reading x. A plain solve (minFlux left at its default) also returns sol.sPrice and sol.rCost, the shadow prices and reduced costs from the LP dual.

from raven_toolbox.io import read_yaml_model

model = read_yaml_model("yeast-GEM.yml")
solution = model.optimize()

print(f"objective: {model.objective.expression}")
print(f"status:    {solution.status}")
print(f"growth:    {solution.objective_value:.4f} /h")

objective: 1.0*r_2111 - 1.0*r_2111_reverse_58b69
status:    optimal
growth:    0.0809 /h

optimize maximises by default and returns a Solution whose fluxes is a pandas Series indexed by reaction id.

The objective prints as a difference of two terms because cobrapy gives every reaction a forward and a reverse variable, both non-negative, and reports the net flux as forward minus reverse. The suffix on the reverse variable is a hash, so it differs between models. Only the expression is split; r_2111 remains one reaction, and solution.fluxes['r_2111'] is the net value.

Changing the objective#

The objective is a reaction to maximise: growth, a product exchange, an ATP demand.

model = setParam(model, 'obj', 'r_2111', 1);   % growth
disp(model.rxnNames{getIndexes(model, 'r_2111', 'rxns')});

growth

setParam writes the coefficient into model.c and zeroes every other entry, so setting an objective replaces the previous one. A negative coefficient minimises that reaction.

model.objective = "r_2111"
print(model.reactions.get_by_id("r_2111").name)

growth

Assigning a reaction or its id replaces the objective. model.objective also accepts an expression, so a weighted combination of reactions is written as one expression, and model.objective_direction switches between maximising and minimising.

Constraining an uptake rate#

yeast-GEM writes its exchange reactions as metabolite <=>, so uptake is a negative flux and the lower bound limits it. A model that writes uptake as => metabolite, like smallYeast in Solvers, has the opposite sign. yeast-GEM has glucose capped at 1 mmol/gDW/h, which determines the growth rate of the first solve; ten times the glucose gives roughly ten times the growth.

idx = getIndexes(model, 'r_1714', 'rxns');     % D-glucose exchange
fprintf('shipped bounds: [%g %g]\n', model.lb(idx), model.ub(idx));

model = setParam(model, 'lb', 'r_1714', -10);
sol = solveLP(model);
fprintf('growth on 10 mmol glucose: %.4f /h\n', sol.f);

model = setParam(model, 'lb', 'r_1714', -1);   % back to the shipped medium

shipped bounds: [-1 1000]
growth on 10 mmol glucose: 0.8370 /h

glucose = model.reactions.get_by_id("r_1714")
print("shipped bounds:", glucose.bounds)

glucose.lower_bound = -10.0
faster = model.optimize()
print(f"growth on 10 mmol glucose: {faster.objective_value:.4f} /h")

glucose.lower_bound = -1.0   # back to the shipped medium

shipped bounds: (-1.0, 1000.0)
growth on 10 mmol glucose: 0.8370 /h

The upper bound of 1000 on the same reaction is the model's convention for "unbounded", and it allows secretion of glucose, which the model has no reason to do. Bounds are a pair, and only the one facing the direction of interest constrains the answer. Growth media and conditions sets all of them at once and covers the opposite sign convention that model.medium uses.

Temporary changes#

A knockout or a tighter bound is usually a question to test, not a change to keep. cobrapy's model is a context manager: changes made inside with model: are rolled back on leaving the block. MATLAB has no equivalent; copy the struct, change the copy, and let it go out of scope.

modelKO = setParam(model, 'eq', 'r_1992', 0);   % close oxygen uptake
solKO = solveLP(modelKO);
fprintf('anaerobic: %.4f /h\n', solKO.f);

% `model` itself is untouched -- MATLAB copied it on assignment
sol = solveLP(model);
fprintf('back to aerobic: %.4f /h\n', sol.f);

anaerobic: 0.0000 /h
back to aerobic: 0.0809 /h

'eq' sets both bounds at once, so a single call closes the reaction in both directions. MATLAB structs are value types, so modelKO is a copy from the moment it is assigned and model cannot be reached through it.

with model:
    model.reactions.get_by_id("r_1992").bounds = (0, 0)   # oxygen exchange
    print(f"anaerobic: {model.optimize().objective_value:.4f} /h")

print(f"back to aerobic: {model.optimize().objective_value:.4f} /h")

anaerobic: -0.0000 /h
back to aerobic: 0.0809 /h

Without with, the bound change would persist and every later result of these examples would be an anaerobic one. The block undoes the edits cobrapy records: bounds, the objective, reactions and metabolites added or removed, stoichiometry changed with add_metabolites, gene rules and gene knockouts. Plain attribute edits are not recorded and survive the block: names, formulas, charges, subsystem, the model id and anything written into annotation or notes. Make those on model.copy().

The anaerobic growth rate is zero. Closing the oxygen exchange is not enough to make yeast-GEM grow fermentatively; growth also needs sterol and fatty-acid uptake and a different biomass composition, which a condition provides (see Growth media and conditions). The minus sign on that zero is solver noise and is not a negative growth rate. Compare the value against a tolerance, not against 0.

Inspecting the fluxes#

An FBA solution has thousands of numbers, most of them zero. Both toolboxes provide a filtered view.

exchanges = {'r_1714', 'r_1992', 'r_1672', 'r_1761'};
labels = {'glucose', 'oxygen', 'CO2', 'ethanol'};
for i = 1:numel(exchanges)
    idx = getIndexes(model, exchanges{i}, 'rxns');
    fprintf('%8s: %9.4f\n', labels{i}, sol.x(idx));
end

 glucose:   -1.0000
  oxygen:   -2.3592
     CO2:    2.6597
 ethanol:    0.0000

printFluxes(model, sol.x) does the same filtering without a list of ids: by default it prints every exchange reaction whose absolute flux is at least 1e-8. Its cutOffFlux raises that floor, onlyExchange set to false includes internal reactions, metaboliteList restricts the print-out to reactions touching named metabolites, and outputString controls the columns, with %eqn, %lower and %upper available alongside %flux.

solution = model.optimize()
exchanges = {
    "r_1714": "glucose",
    "r_1992": "oxygen",
    "r_1672": "CO2",
    "r_1761": "ethanol",
}
for rxn_id, label in exchanges.items():
    print(f"{label:>8}: {solution.fluxes[rxn_id]:9.4f}")

 glucose:   -1.0000
  oxygen:   -2.3592
     CO2:    2.6597
 ethanol:    0.0000

model.summary() prints the same picture (uptake, secretion and the objective) as a table, and solution.fluxes is a pandas Series, so ordinary filtering works: solution.fluxes[solution.fluxes.abs() > 1e-6].

In yeast-GEM's orientation the signs say which way material moves: negative through an exchange is uptake, positive is secretion. Glucose and oxygen enter, CO2 leaves, and ethanol is zero because the model has enough oxygen to respire everything it takes up.

Parsimonious solution#

An FBA optimum is rarely unique: many flux distributions reach the same growth rate, and a plain solve returns an arbitrary one, often with internal loops that carry flux without contributing to the objective. Parsimonious FBA (pFBA) keeps the objective at its optimum and then minimises the total flux, on the reasoning that a cell does not run reactions it gains nothing from.

solPars = solveLP(model, 'minFlux', 1);   % minimise sum(abs(fluxes))
fprintf('growth:     %.4f /h\n', solPars.x(getIndexes(model, 'r_2111', 'rxns')));
fprintf('total flux: %.1f\n', sum(abs(solPars.x)));

growth:     0.0809 /h
total flux: 100.6

minFlux selects the second optimisation. 1 minimises the sum of absolute fluxes, which requires one further LP. 3 minimises the number of active reactions, which is a mixed-integer problem: the result is easier to read as a pathway, and the solve is far slower. Leaving minFlux at 0 skips the second solve and is the only setting that reports shadow prices.

from cobra.flux_analysis import pfba

parsimonious = pfba(model)
print(f"growth:     {parsimonious.fluxes['r_2111']:.4f} /h")
print(f"total flux: {parsimonious.fluxes.abs().sum():.1f}")

growth:     0.0809 /h
total flux: 100.6

pfba is the equivalent of minFlux set to 1. It takes fraction_of_optimum, which relaxes the objective before minimising flux, so 0.95 asks for the leanest distribution that still reaches 95 % of the optimum.

The growth rate is unchanged and the total flux is the smallest that achieves it. The number of active reactions varies between solves: pFBA pins the total flux and leaves open which reactions carry it, so a test on that count fails intermittently.

Troubleshooting

  • The status is infeasible: something is over-constrained, most often an uptake that is closed, or a bound set to the wrong sign. See Growth media and conditions.

  • Growth is zero but the solve succeeded: the model is feasible and the optimum is zero, because a nutrient is missing or a gap blocks the biomass pseudoreaction.

  • The same model gives different flux distributions: this is expected, as alternative optima exist. Use pFBA, or compare ranges with flux variability analysis (FVA), not one solution.

  • Reading sol.f as a negated objective: RAVEN minimises -c'x internally and negates the result before returning, so sol.f is the growth rate, not its negative.

  • A model whose coefficients span many orders of magnitude: ill-scaled stoichiometry makes the solver report an optimum that shifts between solvers. solveLP takes params.maxRatio, which splits any reaction whose coefficients span more than that ratio through auxiliary metabolites before solving; the feasible region is unchanged.

See also#