Table of Contents

Author(s): Ronan M.T. Fleming, Leiden University
Reviewer(s):
INTRODUCTION
TIMING
E. coli core model
MATERIALS - EQUIPMENT SETUP
PROCEDURE
Load E. coli core model
Checking the non-trivial constraints on a model
What are the default constraints on the model?
Hint: printConstraints
Example 1: Calculating growth rates
What is the growth rate of E. coli on glucose (uptake rate = 18.5 mmol/gDW/h) under aerobic and anaerobic conditions?
Hint: changeRxnBounds, changeObjective, optimizeCbModel, printFluxVector
What are the main fields to check in the FBAsolution structure?
Hint: help optimizeCbModel
What does FBAsolution.stat mean?
Example 2: Display an optimal flux vector on a metabolic map
Which reactions/pathways are in use (look at the flux vector and flux map)?
Hint: drawFlux
Example 3: Anerobic growth
What is the optimal growth rate under anaerobic conditions?
Hint: changeRxnBounds
What reactions of oxidative phosphorylation are active in anaerobic conditions?
Hint: printFluxVector drawFlux
Example 4: Growth on alternate substrates
What is the growth rate of E. coli on succinate?
Hint: changeRxnBounds
TROUBLESHOOTING
Acknowledgments
REFERENCES

In this tutorial, Flux Balance Analysis (FBA) is introduced using the E. coli core model, with functions in the COBRA Toolbox v3.0 [2].

Flux balance analysis is a solution to the optimisation problem

where c is a vector of linear objective coefficients, S is an m times n matrix of stoichiometric coefficients for m molecular species involved in n reactions. are n times 1 vectors that are the lower and upper bounds on the n times 1 variable vector of reaction rates (fluxes). The optimal objective value is is always unique, but the optimal vector is usually not unique.

In summary, the data is {c,S,l,u} and the variable being optimised is v.

< 1 hrs

A map of the E. coli core model is shown in Figure 1.

Figure 1 Map of the core E. coli metabolic network. Orange circles represent cytosolic metabolites, yellow circles represent extracellular metabolites, and the blue arrows represent reactions. Reaction name abbreviations are uppercase (blue) and metabolite name abbreviations are lowercase (rust colour). This flux map was drawn using SimPheny and edited for clarity with Adobe Illustrator.

Please ensure that all the required dependencies (e.g. , git and curl) of The COBRA Toolbox have been properly installed by following the installation guide here. Please ensure that the COBRA Toolbox has been initialised (tutorial_initialize.mlx) and verify that the pre-packaged LP and QP solvers are functional (tutorial_verify.mlx).

The most appropriate way to load a model into The COBRA Toolbox is to use the readCbModel function.

fileName = 'ecoli_core_model.mat';

if ~exist('modelOri','var')

modelOri = readCbModel(fileName);

end

%backward compatibility with primer requires relaxation of upper bound on

%ATPM

modelOri = changeRxnBounds(modelOri,'ATPM',1000,'u');

model = modelOri;

The meaning of each field in a standard model is defined in the standard COBRA model field definition.

In general, the following fields should always be present:

- S, the stoichiometric matrix
- mets, the identifiers of the metabolites
- b, Accumulation (positive) or depletion (negative) of the corresponding metabolites. 0 Indicates no concentration change.
- csense, indicator whether the b vector is a lower bound ('G'), upper bound ('L'), or hard constraint 'E' for the metabolites.
- rxns, the identifiers of the reactions
- lb, the lower bounds of the reactions
- ub, the upper bounds of the reactions
- c, the linear objective
- genes, the list of genes in your model
- rules, the Gene-protein-reaction rules in a computer readable format present in your model.
- osenseStr, the objective sense either 'max' for maximisation or 'min' for minimisation

printConstraints(model,-1000,1000)

Growth of E. coli on glucose can be simulated under aerobic conditions.

To set the maximum glucose uptake rate to 18.5 mmol gDW-1 hr-1 (millimoles per gram dry cell weight per hour, the default flux units used in the COBRA Toolbox), enter:

model = changeRxnBounds(model,'EX_glc(e)',-18.5,'l');

This changes the lower bound ('l') of the glucose exchange reaction to -18.5, a biologically realistic uptake rate. By convention, exchange reactions are written as export reactions (e.g. ‘glc[e] <==>’), so import of a metabolite is a negative flux.

To allow unlimited oxygen uptake, enter:

model = changeRxnBounds(model,'EX_o2(e)',-1000,'l');

By setting the lower bound of the oxygen uptake reaction to such a large number, it is practically unbounded. Next, to ensure that the biomass reaction is set as the objective function, enter:

model = changeObjective(model,'Biomass_Ecoli_core_N(w/GAM)-Nmet2');

The default behaviour of the function optimizeCbModel is to conduct flux balance analysis on a model structure. This function is the most powerful function in the COBRA toolbox and has many optional inputs and outputs but for a straightforward flux balance analysis, the following fields of a model are the required input:

% model: (the following fields are required - others can be supplied)

%

% * S - `m x 1` Stoichiometric matrix

% * c - `n x 1` Linear objective coefficients

% * lb - `n x 1` Lower bounds

% * ub - `n x 1` Upper bounds

To perform FBA with maximization of the biomass reaction as the objective, enter:

FBAsolution = optimizeCbModel(model,'max')

There are many outputs fields in the FBAsolution structure, but the key ones for any FBA solution are

% solution: solution object:

% * f - Objective value

% * v - Reaction rates (Optimal primal variable, legacy FBAsolution.x)

% * stat - Solver status in standardized form

FBAsolution.f gives the value of the objective function, which should be 1.6531 mmol gDW-1 hr-1

FBAsolution.f

ans = 1.6531

This is the same as

model.c'*FBAsolution.v

ans = 1.6531

This means that the model predicts a growth rate of 1.6531 hr-1. Inspection of the flux distribution vector FBAsolution.v shows that there is high flux in the glycolysis, pentose phosphate, TCA cycle, and oxidative phosphorylation pathways, and that no organic by-products are secreted (Figure 2a).

FBAsolution.v gives the optimal flux vector

FBAsolution.v

Inspection of the flux distribution is more convenient with the printFluxVector function

fluxData = FBAsolution.v;

nonZeroFlag = 1;

printFluxVector(model, fluxData, nonZeroFlag)

outputFormatOK = changeCbMapOutput('matlab');

map=readCbMap('ecoli_core_map');

options.zeroFluxWidth = 0.1;

options.rxnDirMultiplier = 10;

drawFlux(map, model, FBAsolution.v, options);

Growth of E. coli on glucose can be simulated under anaerobic conditions.

Next, the same simulation is performed under anaerobic conditions. With the same model:

model = changeRxnBounds(model,'EX_o2(e)',0,'l');

The lower bound of the oxygen exchange reaction is now 0, so oxygen may not enter the system. When optimizeCbModel is used as before, the resulting growth rate is now much lower, 0.4706 hr-1.

FBAsolution2 = optimizeCbModel(model,'max');

FBAsolution2.f

ans = 0.4706

The flux distribution shows that oxidative phosphorylation is not used in these conditions, and that acetate, formate, and ethanol are produced by fermentation pathways (Figure 2b). Inspection of both flux vectors for comparison:

fluxData = [FBAsolution.v,FBAsolution2.v];

nonZeroFlag = 1;

excFlag = 1;

printFluxVector(model, fluxData, nonZeroFlag, excFlag)

drawFlux(map, model, FBAsolution2.v, options);

Figure 2 Flux vectors computed by FBA can be visualized on network maps. In these two examples, the thick blue arrows represent reactions carrying flux, and the thin black arrows represent unused reactions. These maps show the state of the E. coli core model with maximum growth rate as the objective (Z) under aerobic (a) and anaerobic (b) conditions. Reactions that are in use have thick blue arrows, while reactions that carry 0 flux have thin black arrows. The metabolic pathways shown in these maps are glycolysis (Glyc), pentose phosphate pathway (PPP), TCA cycle (TCA), oxidative phosphorylation (OxP), anaplerotic reactions (Ana), and fermentation pathways (Ferm).

Just as FBA was used to calculate growth rates of E. coli on glucose, it can also be used to simulate growth on other substrates. The core E. coli model contains exchange reactions for 13 different organic compounds, each of which can be used as the sole carbon source under aerobic conditions.

Before trying out new boundary conditions on a model, make sure one's starting point is appropriate.

model = modelOri;

For example, to simulate growth on succinate instead of glucose, first use the changeRxnBounds function to set the lower bound of the glucose exchange reaction (EX_glc(e)) to 0.

model = changeRxnBounds(model,'EX_glc(e)',0,'l');

Then use changeRxnBounds to set the lower bound of the succinate exchange reaction (EX_succ(e)) to -20 mmol gDW-1 hr-1 (an arbitrary uptake rate).

model = changeRxnBounds(model,'EX_succ(e)',-20,'l');

As in the glucose examples, make sure that the biomass reaction is set as the objective (the function checkObjective can be used to identify the objective reaction(s)),

checkObjective(model);