cameo uses 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.
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.
For example, 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), assuming that the cell is evolutionary optimized to achieve that objective, 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{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} $$Load a model.
from cameo import load_model
model = load_model('iJO1366')
In cameo, flux balance analysis can be performed with the function fba
.
from cameo import fba
%time fba_result = fba(model)
CPU times: user 256 ms, sys: 6.4 ms, total: 262 ms Wall time: 283 ms
Basically, fba
calls model.optimize()
and wraps the optimization solution in a FluxDistributionResult
object. The maximum objective values (corresponding to a maximum growth rate) can be obtained from result.objective_value
.
fba_result.data_frame
flux | |
---|---|
12DGR120tipp | 0.000000 |
12DGR140tipp | 0.000000 |
12DGR141tipp | 0.000000 |
12DGR160tipp | 0.000000 |
12DGR161tipp | 0.000000 |
12DGR180tipp | 0.000000 |
12DGR181tipp | 0.000000 |
12PPDRtex | 0.000000 |
12PPDRtpp | 0.000000 |
12PPDStex | 0.000000 |
12PPDStpp | 0.000000 |
14GLUCANabcpp | 0.000000 |
14GLUCANtexi | 0.000000 |
23CAMPtex | 0.000000 |
23CCMPtex | 0.000000 |
23CGMPtex | 0.000000 |
23CUMPtex | 0.000000 |
23DAPPAt2pp | 0.000000 |
23DAPPAtex | 0.000000 |
23PDE2pp | 0.000000 |
23PDE4pp | 0.000000 |
23PDE7pp | 0.000000 |
23PDE9pp | 0.000000 |
26DAHtex | 0.000000 |
2AGPA120tipp | 0.000000 |
2AGPA140tipp | 0.000000 |
2AGPA141tipp | 0.000000 |
2AGPA160tipp | 0.000000 |
2AGPA161tipp | 0.000000 |
2AGPA180tipp | 0.000000 |
... | ... |
VALTRS | 0.000000 |
VALabcpp | 0.000000 |
VALt2rpp | 0.000000 |
VALtex | 0.000000 |
VPAMTr | 0.000000 |
WCOS | 0.000000 |
X5PL3E | 0.000000 |
XAND | 0.000000 |
XANt2pp | 0.000000 |
XANtex | 0.000000 |
XANtpp | 0.000000 |
XMPtex | 0.000000 |
XPPT | 0.000000 |
XTSNH | 0.000000 |
XTSNt2rpp | 0.000000 |
XTSNtex | 0.000000 |
XYLI1 | 0.000000 |
XYLI2 | 0.000000 |
XYLK | 0.000000 |
XYLK2 | 0.000000 |
XYLUt2pp | 0.000000 |
XYLUtex | 0.000000 |
XYLabcpp | 0.000000 |
XYLt2pp | 0.000000 |
XYLtex | 0.000000 |
ZN2abcpp | 0.000000 |
ZN2t3pp | 0.000000 |
ZN2tpp | 0.000335 |
ZNabcpp | 0.000000 |
Zn2tex | 0.000335 |
2583 rows × 1 columns
Flux distributions can be visualized using escher :
fba_result.display_on_map("iJO1366.Central metabolism")
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 the $L_1$ norm of $\mathbf{v}$. The assumption behind pFBA is that cells try to minimize flux magnitude as well in order to keep protein costs low.
$$ \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} $$In cameo, pFBA can be performed with the function pfba
.
from cameo import pfba
%time pfba_result = pfba(model)
CPU times: user 6.46 s, sys: 77 ms, total: 6.53 s Wall time: 10.3 s
The objective_function
value is $\lvert \mathbf{v} \rvert$ ...
pfba_result.objective_value
699.0222751839516
... which is smaller than flux vector of the original FBA solution.
abs(fba_result.data_frame.flux).sum()
765.0977334751339
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 implemented a linear version of MOMA.
Simulating knockouts:
model.reactions.PGI
Reaction identifier | PGI |
Name | Glucose-6-phosphate isomerase |
Memory address | 0x011313e908 |
Stoichiometry |
g6p_c <=> f6p_c D-Glucose 6-phosphate <=> D-Fructose 6-phosphate |
GPR | b4025 |
Lower bound | -1000.0 |
Upper bound | 1000.0 |
model.reactions.PGI.knock_out()
model.reactions.PGI
Reaction identifier | PGI |
Name | Glucose-6-phosphate isomerase |
Memory address | 0x011313e908 |
Stoichiometry |
g6p_c --> f6p_c D-Glucose 6-phosphate --> D-Fructose 6-phosphate |
GPR | b4025 |
Lower bound | 0 |
Upper bound | 0 |
%time fba_knockout_result = fba(model)
fba_knockout_result[model.reactions.BIOMASS_Ec_iJO1366_core_53p95M]
CPU times: user 104 ms, sys: 4.13 ms, total: 108 ms Wall time: 124 ms
0.9761293262947268
%time pfba_knockout_result = pfba(model)
pfba_knockout_result[model.reactions.BIOMASS_Ec_iJO1366_core_53p95M]
CPU times: user 5.78 s, sys: 38.7 ms, total: 5.82 s Wall time: 7.48 s
0.9761293262947268
MOMA and ROOM rely on a reference (wild-type) flux distribution. We can use the one previously computed.
Parsimonious FBA references seem to produce better results using this methods
from cameo.flux_analysis.simulation import room, lmoma
%time lmoma_result = lmoma(model, reference=pfba_result.fluxes)
lmoma_result[model.reactions.BIOMASS_Ec_iJO1366_core_53p95M]
CPU times: user 27.4 s, sys: 201 ms, total: 27.6 s Wall time: 35.6 s
0.8724092397035721
ROOM is a difficult computational problem. If the bounds of the system are not large enough, it can take many hours to simulate. To improve the speed of the simulation and the chances of finding a solution, we increase the bounds.
for reaction in model.reactions:
if reaction.upper_bound == 1000:
reaction.upper_bound = 99999999
if reaction.lower_bound == -1000:
reaction.lower_bound = -99999999
%time room_result = room(model, reference=pfba_result.fluxes)
room_result[model.reactions.BIOMASS_Ec_iJO1366_core_53p95M]
CPU times: user 35 s, sys: 257 ms, total: 35.3 s Wall time: 44.4 s
0.9519006583451706