{ "cells": [ { "cell_type": "markdown", "metadata": {}, "source": [ "# Simulate models\n", "\n", "**cameo** uses the model data structures defined by [cobrapy](https://opencobra.github.io/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." ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "## Primer: Constraint-Based Modeling\n", "\n", "Constraint-based modeling is a powerful modeling framework for analyzing metabolism on the genome scale ([McCloskey et al., 2013](http://www.ncbi.nlm.nih.gov/pubmed/23632383)). 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:\n", "\n", "$$\n", "\\begin{align}\n", "\\mathbf{S} \\mathbf{v} = 0\\,,\n", "\\end{align}\n", "$$\n", "\n", "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](https://en.wikipedia.org/wiki/Linear_programming).\n", "\n", "\n", "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):\n", "\n", "$$\n", "\\begin{align}\n", " Max ~ & ~ Z_{obj} = \\mathbf{c}^{T} \\mathbf{v}\\\\\n", " \\text{s.t.}~ & ~ \\mathbf{S} \\mathbf{v} = 0 \\\\\n", " ~ & ~ \\mathbf{v}_{lb} \\leq \\mathbf{v} \\leq \\mathbf{v}_{ub} \\,.\n", "\\end{align}\n", "$$" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "## Flux Balance Analysis\n", "\n", "Load a model." ] }, { "cell_type": "code", "execution_count": 1, "metadata": {}, "outputs": [], "source": [ "from cameo import load_model\n", "model = load_model('iJO1366')" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "In cameo, flux balance analysis can be performed with the function `fba`." ] }, { "cell_type": "code", "execution_count": 2, "metadata": {}, "outputs": [ { "name": "stdout", "output_type": "stream", "text": [ "CPU times: user 256 ms, sys: 6.4 ms, total: 262 ms\n", "Wall time: 283 ms\n" ] } ], "source": [ "from cameo import fba\n", "%time fba_result = fba(model)" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "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`." ] }, { "cell_type": "code", "execution_count": 3, "metadata": {}, "outputs": [ { "data": { "text/html": [ "
\n", "\n", "\n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", "
flux
12DGR120tipp0.000000
12DGR140tipp0.000000
12DGR141tipp0.000000
12DGR160tipp0.000000
12DGR161tipp0.000000
12DGR180tipp0.000000
12DGR181tipp0.000000
12PPDRtex0.000000
12PPDRtpp0.000000
12PPDStex0.000000
12PPDStpp0.000000
14GLUCANabcpp0.000000
14GLUCANtexi0.000000
23CAMPtex0.000000
23CCMPtex0.000000
23CGMPtex0.000000
23CUMPtex0.000000
23DAPPAt2pp0.000000
23DAPPAtex0.000000
23PDE2pp0.000000
23PDE4pp0.000000
23PDE7pp0.000000
23PDE9pp0.000000
26DAHtex0.000000
2AGPA120tipp0.000000
2AGPA140tipp0.000000
2AGPA141tipp0.000000
2AGPA160tipp0.000000
2AGPA161tipp0.000000
2AGPA180tipp0.000000
......
VALTRS0.000000
VALabcpp0.000000
VALt2rpp0.000000
VALtex0.000000
VPAMTr0.000000
WCOS0.000000
X5PL3E0.000000
XAND0.000000
XANt2pp0.000000
XANtex0.000000
XANtpp0.000000
XMPtex0.000000
XPPT0.000000
XTSNH0.000000
XTSNt2rpp0.000000
XTSNtex0.000000
XYLI10.000000
XYLI20.000000
XYLK0.000000
XYLK20.000000
XYLUt2pp0.000000
XYLUtex0.000000
XYLabcpp0.000000
XYLt2pp0.000000
XYLtex0.000000
ZN2abcpp0.000000
ZN2t3pp0.000000
ZN2tpp0.000335
ZNabcpp0.000000
Zn2tex0.000335
\n", "

2583 rows × 1 columns

\n", "
" ], "text/plain": [ " flux\n", "12DGR120tipp 0.000000\n", "12DGR140tipp 0.000000\n", "12DGR141tipp 0.000000\n", "12DGR160tipp 0.000000\n", "12DGR161tipp 0.000000\n", "12DGR180tipp 0.000000\n", "12DGR181tipp 0.000000\n", "12PPDRtex 0.000000\n", "12PPDRtpp 0.000000\n", "12PPDStex 0.000000\n", "12PPDStpp 0.000000\n", "14GLUCANabcpp 0.000000\n", "14GLUCANtexi 0.000000\n", "23CAMPtex 0.000000\n", "23CCMPtex 0.000000\n", "23CGMPtex 0.000000\n", "23CUMPtex 0.000000\n", "23DAPPAt2pp 0.000000\n", "23DAPPAtex 0.000000\n", "23PDE2pp 0.000000\n", "23PDE4pp 0.000000\n", "23PDE7pp 0.000000\n", "23PDE9pp 0.000000\n", "26DAHtex 0.000000\n", "2AGPA120tipp 0.000000\n", "2AGPA140tipp 0.000000\n", "2AGPA141tipp 0.000000\n", "2AGPA160tipp 0.000000\n", "2AGPA161tipp 0.000000\n", "2AGPA180tipp 0.000000\n", "... ...\n", "VALTRS 0.000000\n", "VALabcpp 0.000000\n", "VALt2rpp 0.000000\n", "VALtex 0.000000\n", "VPAMTr 0.000000\n", "WCOS 0.000000\n", "X5PL3E 0.000000\n", "XAND 0.000000\n", "XANt2pp 0.000000\n", "XANtex 0.000000\n", "XANtpp 0.000000\n", "XMPtex 0.000000\n", "XPPT 0.000000\n", "XTSNH 0.000000\n", "XTSNt2rpp 0.000000\n", "XTSNtex 0.000000\n", "XYLI1 0.000000\n", "XYLI2 0.000000\n", "XYLK 0.000000\n", "XYLK2 0.000000\n", "XYLUt2pp 0.000000\n", "XYLUtex 0.000000\n", "XYLabcpp 0.000000\n", "XYLt2pp 0.000000\n", "XYLtex 0.000000\n", "ZN2abcpp 0.000000\n", "ZN2t3pp 0.000000\n", "ZN2tpp 0.000335\n", "ZNabcpp 0.000000\n", "Zn2tex 0.000335\n", "\n", "[2583 rows x 1 columns]" ] }, "execution_count": 3, "metadata": {}, "output_type": "execute_result" } ], "source": [ "fba_result.data_frame" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "Flux distributions can be visualized using [*escher*](https://escher.github.io) :" ] }, { "cell_type": "code", "execution_count": 4, "metadata": {}, "outputs": [ { "data": { "text/html": [ "\n", "\n", "\n", " \n", " \n", "\n", " \n", " \n", " \n", " \n", "\n", "\n", "\n", "
\n", "\n", " \n" ], "text/plain": [ "" ] }, "metadata": {}, "output_type": "display_data" } ], "source": [ "fba_result.display_on_map(\"iJO1366.Central metabolism\")" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "## Parsimonious Flux Balance Analysis\n", "\n", "Parsimonious flux balance analysis ([Lewis et al., 2010](http://www.ncbi.nlm.nih.gov/pubmed/20664636)), 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.\n", "\n", "$$\n", "\\begin{align}\n", " Max ~ & ~ \\lvert \\mathbf{v} \\rvert\\\\\n", " \\text{s.t.}~ & ~ \\mathbf{S} \\mathbf{v} = 0 \\\\\n", " & ~ \\mathbf{c}^{T} \\mathbf{v} \\ge Z_{obj} \\\\\n", " ~ & ~ \\mathbf{v}_{lb} \\leq \\mathbf{v} \\leq \\mathbf{v}_{ub} \\,.\n", "\\end{align}\n", "$$\n", "\n", "In cameo, pFBA can be performed with the function `pfba`." ] }, { "cell_type": "code", "execution_count": 5, "metadata": { "scrolled": true }, "outputs": [ { "name": "stdout", "output_type": "stream", "text": [ "CPU times: user 6.46 s, sys: 77 ms, total: 6.53 s\n", "Wall time: 10.3 s\n" ] } ], "source": [ "from cameo import pfba\n", "%time pfba_result = pfba(model)" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "The `objective_function` value is $\\lvert \\mathbf{v} \\rvert$ ..." ] }, { "cell_type": "code", "execution_count": 6, "metadata": {}, "outputs": [ { "data": { "text/plain": [ "699.0222751839516" ] }, "execution_count": 6, "metadata": {}, "output_type": "execute_result" } ], "source": [ "pfba_result.objective_value" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "... which is smaller than flux vector of the original FBA solution." ] }, { "cell_type": "code", "execution_count": 7, "metadata": {}, "outputs": [ { "data": { "text/plain": [ "765.0977334751339" ] }, "execution_count": 7, "metadata": {}, "output_type": "execute_result" } ], "source": [ "abs(fba_result.data_frame.flux).sum()" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "Setp 2: Simulate knockouts phenotypes\n", "-----------------------------------\n", "\n", "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.\n", "\n", "*******************************************\n", "Simulating knockouts:\n", "\n", "* Manipulate the bounds of the reaction (or use the shorthand method knock_out)" ] }, { "cell_type": "code", "execution_count": 8, "metadata": {}, "outputs": [ { "data": { "text/html": [ "\n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", "
Reaction identifierPGI
NameGlucose-6-phosphate isomerase
Memory address0x011313e908
Stoichiometry\n", "

g6p_c <=> f6p_c

\n", "

D-Glucose 6-phosphate <=> D-Fructose 6-phosphate

\n", "
GPRb4025
Lower bound-1000.0
Upper bound1000.0
\n", " " ], "text/plain": [ "" ] }, "execution_count": 8, "metadata": {}, "output_type": "execute_result" } ], "source": [ "model.reactions.PGI" ] }, { "cell_type": "code", "execution_count": 9, "metadata": {}, "outputs": [ { "data": { "text/html": [ "\n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", "
Reaction identifierPGI
NameGlucose-6-phosphate isomerase
Memory address0x011313e908
Stoichiometry\n", "

g6p_c --> f6p_c

\n", "

D-Glucose 6-phosphate --> D-Fructose 6-phosphate

\n", "
GPRb4025
Lower bound0
Upper bound0
\n", " " ], "text/plain": [ "" ] }, "execution_count": 9, "metadata": {}, "output_type": "execute_result" } ], "source": [ "model.reactions.PGI.knock_out()\n", "model.reactions.PGI" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "* Simulate using different methods:" ] }, { "cell_type": "code", "execution_count": 10, "metadata": {}, "outputs": [ { "name": "stdout", "output_type": "stream", "text": [ "CPU times: user 104 ms, sys: 4.13 ms, total: 108 ms\n", "Wall time: 124 ms\n" ] }, { "data": { "text/plain": [ "0.9761293262947268" ] }, "execution_count": 10, "metadata": {}, "output_type": "execute_result" } ], "source": [ "%time fba_knockout_result = fba(model)\n", "fba_knockout_result[model.reactions.BIOMASS_Ec_iJO1366_core_53p95M]" ] }, { "cell_type": "code", "execution_count": 11, "metadata": {}, "outputs": [ { "name": "stdout", "output_type": "stream", "text": [ "CPU times: user 5.78 s, sys: 38.7 ms, total: 5.82 s\n", "Wall time: 7.48 s\n" ] }, { "data": { "text/plain": [ "0.9761293262947268" ] }, "execution_count": 11, "metadata": {}, "output_type": "execute_result" } ], "source": [ "%time pfba_knockout_result = pfba(model)\n", "pfba_knockout_result[model.reactions.BIOMASS_Ec_iJO1366_core_53p95M]" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "MOMA and ROOM rely on a reference (wild-type) flux distribution. We can use the one previously computed.\n", "\n", "**Parsimonious FBA references seem to produce better results using this methods**" ] }, { "cell_type": "code", "execution_count": 12, "metadata": {}, "outputs": [], "source": [ "from cameo.flux_analysis.simulation import room, lmoma" ] }, { "cell_type": "code", "execution_count": 13, "metadata": {}, "outputs": [ { "name": "stdout", "output_type": "stream", "text": [ "CPU times: user 27.4 s, sys: 201 ms, total: 27.6 s\n", "Wall time: 35.6 s\n" ] }, { "data": { "text/plain": [ "0.8724092397035721" ] }, "execution_count": 13, "metadata": {}, "output_type": "execute_result" } ], "source": [ "%time lmoma_result = lmoma(model, reference=pfba_result.fluxes)\n", "lmoma_result[model.reactions.BIOMASS_Ec_iJO1366_core_53p95M]" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "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." ] }, { "cell_type": "code", "execution_count": 14, "metadata": {}, "outputs": [], "source": [ "for reaction in model.reactions:\n", " if reaction.upper_bound == 1000:\n", " reaction.upper_bound = 99999999\n", " if reaction.lower_bound == -1000:\n", " reaction.lower_bound = -99999999" ] }, { "cell_type": "code", "execution_count": 15, "metadata": {}, "outputs": [ { "name": "stdout", "output_type": "stream", "text": [ "CPU times: user 35 s, sys: 257 ms, total: 35.3 s\n", "Wall time: 44.4 s\n" ] }, { "data": { "text/plain": [ "0.9519006583451706" ] }, "execution_count": 15, "metadata": {}, "output_type": "execute_result" } ], "source": [ "%time room_result = room(model, reference=pfba_result.fluxes)\n", "room_result[model.reactions.BIOMASS_Ec_iJO1366_core_53p95M]" ] } ], "metadata": { "anaconda-cloud": {}, "kernelspec": { "display_name": "Python 3", "language": "python", "name": "python3" }, "language_info": { "codemirror_mode": { "name": "ipython", "version": 3 }, "file_extension": ".py", "mimetype": "text/x-python", "name": "python", "nbconvert_exporter": "python", "pygments_lexer": "ipython3", "version": "3.9.2" } }, "nbformat": 4, "nbformat_minor": 4 }