from pandas import options
options.display.max_rows = 8
from cameo import load_model
model = load_model("iJO1366")

Simulating models with

computer aided metabolic engineering and optimization

cameo uses and extends the model data structures defined by cobrapy, our favorite COnstraints-Based Reconstruction and Analysis tool for Python. cameo is thus 100% compatible with cobrapy. For efficiency reasons, however, cameo implements its own simulation methods that take advantage of a more advanced solver interface.

Primer: Constraint-Based Modeling

Constraint-based modeling is a powerful modeling framework for analyzing metabolism on the genome scale (McCloskey et al., 2013). For a model that encompasses \(n\) reactions that involve \(m\) metabolites, \(\mathbf{S}\) is a matrix of dimension \(m \times n\) that encodes the stoichiometry of the metabolic reaction system; it is usually referred to as stoichiometric matrix. Assuming that the system is in a steady state—the concentration of metabolites are constant—the system of flux-balances can be formulated as

\[\begin{align} \mathbf{S} \mathbf{v} = 0\,, \end{align}\]

where \(\mathbf{v}\) is the vector of flux rates. With the addition of a biologically meaningful objective, flux capacity constraints, information about the reversibility of reactions under physiological conditions, an optimization problem can be formulated that can easily be solved using linear programming.

, e.g., maximimization of biomass production,Given the maximization of growth rate as one potential biological objective \(v_{biomass}\), i.e., the flux of an artificial reaction that consumes biomass components in empirically determined proportions, and assuming that the cell is evolutionary optimized to achieve that objective, and incorporating knowledge about reaction reversibility, uptake and secretion rates, and maximum flux capacities in the form of lower and uppers bounds (\(\mathbf{v}_{lb}\) and \(\mathbf{v}_{ub}\)) on the flux variables \(\mathbf{v}\), one can formulate and solve an optimization problem to identify an optimal set of flux rates using flux balance analysis (FBA):

\[\begin{split}\begin{align} Max ~ & ~ Z_{obj} = \mathbf{c}^{T} \mathbf{v}\\ \text{s.t.}~ & ~ \mathbf{S} \mathbf{v} = 0 \\ ~ & ~ \mathbf{v}_{lb} \leq \mathbf{v} \leq \mathbf{v}_{ub} \,. \end{align}\end{split}\]

Flux Balance Analysis

In cameo, flux balance analysis can be performed with the function fba.

from cameo import fba
fba_result = fba(model)

Basically, fba calls model.solve() and wraps the optimization solution in a FluxDistributionResult object. The maximum objective values (corresponding to a maximum growth rate) can obtained throug result.objective_value.

fba_result.objective_value

0.9823718127269799

Parsimonious Flux Balance Analysis

Parsimonious flux balance analysis (Lewis et al., 2010), a variant of FBA, performs FBA in in a first step to determine the maximum objective value \(Z_{obj}\), fixes it in form of an additional model constraint (\(\mathbf{c}^{T} \mathbf{v} \ge Z_{obj}\)), and then minimizes in a second optimization the \(L_1\) norm of \(\mathbf{v}\). The assumption behind the pFBA is that cells try to minimize flux magnitude as well in order to keep the costs of protein low.

\[\begin{split}\begin{align} Max ~ & ~ \lvert \mathbf{v} \rvert\\ \text{s.t.}~ & ~ \mathbf{S} \mathbf{v} = 0 \\ & ~ \mathbf{c}^{T} \mathbf{v} \ge Z_{obj} \\ ~ & ~ \mathbf{v}_{lb} \leq \mathbf{v} \leq \mathbf{v}_{ub} \,. \end{align}\end{split}\]

In cameo, pFBA can be performed with the function pfba.

from cameo import pfba
pfba_result = pfba(model)

The objective_function value is \(\lvert \mathbf{v} \rvert\) ...

pfba_result.objective_value

699.0222751839377

... whis is significantly smaller than flux vector of the original FBA solution.

abs(fba_result.data_frame.flux).sum()

764.91487969777245

Setp 2: Simulate knockouts phenotypes

Although PFBA and FBA can be used to simulate the effect of knockouts, other methods have been proven more valuable for that task: MOMA and ROOM. In cameo we implement a linear version of MOMA.

Simulating knockouts:

  • Manipulate the bounds of the reaction (or use the shorthand method knock_out)
model.reactions.PGI
IdPGI
Stoichiometryg6p_c <=> f6p_c
Lower bound-999999.000000
Upper bound999999.000000 
model.reactions.PGI.knock_out()
model.reactions.PGI
IdPGI
Stoichiometryg6p_c --> f6p_c
Lower bound0.000000
Upper bound0.000000
  • Simulate using different methods:
%time
fba_knockout_result = simulation.fba(model)
fba_knockout_result[model.objective]

CPU times: user 2 µs, sys: 0 ns, total: 2 µs
Wall time: 5.01 µs

0.905983

pfba_knockout_result = simulation.pfba(model)
pfba_knockout_result[model.objective]

0.905983

MOMA and ROOM relly on a reference (wild-type) flux distribution and we can use the one previously computed.

Parsimonious FBA references seem to produce better results using this methods

lmoma_result["2 * EX_glc_lp_e_rp_"]

-18.7358

%time
lmoma_result = simulation.lmoma(model, reference=pfba_result.fluxes)
lmoma_result[model.objective]

CPU times: user 2 µs, sys: 1 µs, total: 3 µs
Wall time: 5.01 µs

0.791393

%time
room_result = simulation.room(model, reference=pfba_result.fluxes)
room_result[model.objective]

CPU times: user 2 µs, sys: 1 µs, total: 3 µs
Wall time: 5.01 µs

0.887440

room_result

<cameo.core.result.FluxDistributionResult at 0x10aa75b50>