{ "cells": [ { "cell_type": "markdown", "metadata": {}, "source": [ "# Simulating with FBA\n", "\n", "Simulations using flux balance analysis can be solved using Model.optimize(). This will maximize or minimize (maximizing is the default) flux through the objective reactions." ] }, { "cell_type": "code", "execution_count": 1, "metadata": { "collapsed": false }, "outputs": [], "source": [ "import pandas\n", "pandas.options.display.max_rows = 100\n", "\n", "import cobra.test\n", "model = cobra.test.create_test_model(\"textbook\")" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "## Running FBA" ] }, { "cell_type": "code", "execution_count": 2, "metadata": { "collapsed": false }, "outputs": [ { "data": { "text/plain": [ "" ] }, "execution_count": 2, "metadata": {}, "output_type": "execute_result" } ], "source": [ "model.optimize()" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "The Model.optimize() function will return a Solution object, which will also be stored at model.solution. A solution object has several attributes:\n", "\n", " - f: the objective value\n", " - status: the status from the linear programming solver\n", " - x_dict: a dictionary of {reaction_id: flux_value} (also called \"primal\")\n", " - x: a list for x_dict\n", " - y_dict: a dictionary of {metabolite_id: dual_value}.\n", " - y: a list for y_dict" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "For example, after the last call to model.optimize(), the status should be 'optimal' if the solver returned no errors, and f should be the objective value" ] }, { "cell_type": "code", "execution_count": 3, "metadata": { "collapsed": false }, "outputs": [ { "data": { "text/plain": [ "'optimal'" ] }, "execution_count": 3, "metadata": {}, "output_type": "execute_result" } ], "source": [ "model.solution.status" ] }, { "cell_type": "code", "execution_count": 4, "metadata": { "collapsed": false }, "outputs": [ { "data": { "text/plain": [ "0.8739215069684305" ] }, "execution_count": 4, "metadata": {}, "output_type": "execute_result" } ], "source": [ "model.solution.f" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "### Analyzing FBA solutions\n", "Models solved using FBA can be further analyzed by using summary methods, which output printed text to give a quick representation of model behavior. Calling the summary method on the entire model displays information on the input and output behavior of the model, along with the optimized objective." ] }, { "cell_type": "code", "execution_count": 5, "metadata": { "collapsed": false }, "outputs": [ { "name": "stdout", "output_type": "stream", "text": [ "IN FLUXES OUT FLUXES OBJECTIVES \n", "o2_e -21.80 h2o_e 29.18 Biomass_Ecoli_core 0.874\n", "glc__D_e -10.00 co2_e 22.81 \n", "nh4_e -4.77 h_e 17.53 \n", "pi_e -3.21 \n" ] } ], "source": [ "model.summary()" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "In addition, the input-output behavior of individual metabolites can also be inspected using summary methods. For instance, the following commands can be used to examine the overall redox balance of the model" ] }, { "cell_type": "code", "execution_count": 6, "metadata": { "collapsed": false }, "outputs": [ { "name": "stdout", "output_type": "stream", "text": [ "PRODUCING REACTIONS -- Nicotinamide adenine dinucleotide - reduced\n", "------------------------------------------------------------------\n", " % FLUX RXN ID REACTION \n", " 41.6% 16 GAPD g3p_c + nad_c + pi_c <=> 13dpg_c + h_c + nadh_c\n", " 24.1% 9.3 PDH coa_c + nad_c + pyr_c --> accoa_c + co2_c + nadh_c\n", " 13.1% 5.1 AKGDH akg_c + coa_c + nad_c --> co2_c + nadh_c + succoa_c\n", " 13.1% 5.1 MDH mal__L_c + nad_c <=> h_c + nadh_c + oaa_c\n", " 8.0% 3.1 Bioma... 1.496 3pg_c + 3.7478 accoa_c + 59.81 atp_c + 0.36...\n", "\n", "CONSUMING REACTIONS -- Nicotinamide adenine dinucleotide - reduced\n", "------------------------------------------------------------------\n", " % FLUX RXN ID REACTION \n", "100.0% -39 NADH16 4.0 h_c + nadh_c + q8_c --> 3.0 h_e + nad_c + q8h2_c\n" ] } ], "source": [ "model.metabolites.nadh_c.summary()" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "Or to get a sense of the main energy production and consumption reactions" ] }, { "cell_type": "code", "execution_count": 7, "metadata": { "collapsed": false }, "outputs": [ { "name": "stdout", "output_type": "stream", "text": [ "PRODUCING REACTIONS -- ATP\n", "--------------------------\n", " % FLUX RXN ID REACTION \n", " 66.6% 46 ATPS4r adp_c + 4.0 h_e + pi_c <=> atp_c + h2o_c + 3.0 h_c\n", " 23.4% 16 PGK 3pg_c + atp_c <=> 13dpg_c + adp_c\n", " 7.4% 5.1 SUCOAS atp_c + coa_c + succ_c <=> adp_c + pi_c + succoa_c\n", " 2.6% 1.8 PYK adp_c + h_c + pep_c --> atp_c + pyr_c\n", "\n", "CONSUMING REACTIONS -- ATP\n", "--------------------------\n", " % FLUX RXN ID REACTION \n", " 76.5% -52 Bioma... 1.496 3pg_c + 3.7478 accoa_c + 59.81 atp_c + 0.36...\n", " 12.3% -8.4 ATPM atp_c + h2o_c --> adp_c + h_c + pi_c\n", " 10.9% -7.5 PFK atp_c + f6p_c --> adp_c + fdp_c + h_c\n" ] } ], "source": [ "model.metabolites.atp_c.summary()" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "## Changing the Objectives\n", "\n", "The objective function is determined from the objective_coefficient attribute of the objective reaction(s). Generally, a \"biomass\" function which describes the composition of metabolites which make up a cell is used." ] }, { "cell_type": "code", "execution_count": 8, "metadata": { "collapsed": false }, "outputs": [], "source": [ "biomass_rxn = model.reactions.get_by_id(\"Biomass_Ecoli_core\")" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "Currently in the model, there is only one objective reaction (the biomass reaction), with an objective coefficient of 1." ] }, { "cell_type": "code", "execution_count": 9, "metadata": { "collapsed": false }, "outputs": [ { "data": { "text/plain": [ "{: 1.0}" ] }, "execution_count": 9, "metadata": {}, "output_type": "execute_result" } ], "source": [ "model.objective" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "The objective function can be changed by assigning Model.objective, which can be a reaction object (or just it's name), or a dict of {Reaction: objective_coefficient}." ] }, { "cell_type": "code", "execution_count": 10, "metadata": { "collapsed": false }, "outputs": [ { "data": { "text/plain": [ "{: 1}" ] }, "execution_count": 10, "metadata": {}, "output_type": "execute_result" } ], "source": [ "# change the objective to ATPM\n", "model.objective = \"ATPM\"\n", "\n", "# The upper bound should be 1000, so that we get\n", "# the actual optimal value\n", "model.reactions.get_by_id(\"ATPM\").upper_bound = 1000.\n", "model.objective" ] }, { "cell_type": "code", "execution_count": 11, "metadata": { "collapsed": false }, "outputs": [ { "data": { "text/plain": [ "174.99999999999997" ] }, "execution_count": 11, "metadata": {}, "output_type": "execute_result" } ], "source": [ "model.optimize().f" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "The objective function can also be changed by setting Reaction.objective_coefficient directly." ] }, { "cell_type": "code", "execution_count": 12, "metadata": { "collapsed": false }, "outputs": [ { "data": { "text/plain": [ "{: 1.0}" ] }, "execution_count": 12, "metadata": {}, "output_type": "execute_result" } ], "source": [ "model.reactions.get_by_id(\"ATPM\").objective_coefficient = 0.\n", "biomass_rxn.objective_coefficient = 1.\n", "\n", "model.objective" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "## Running FVA\n", "\n", "FBA will not give always give unique solution, because multiple flux states can achieve the same optimum. FVA (or flux variability analysis) finds the ranges of each metabolic flux at the optimum." ] }, { "cell_type": "code", "execution_count": 13, "metadata": { "collapsed": false }, "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", "
maximumminimum
ACALD0.000000.00000
ACALDt-0.000000.00000
ACKr-0.000000.00000
ACONTa6.007256.00725
ACONTb6.007256.00725
ACt2r0.000000.00000
ADK1-0.000000.00000
AKGDH5.064385.06438
AKGt2r0.000000.00000
ALCD2x0.000000.00000
ATPM8.390008.39000
ATPS4r45.5140145.51401
Biomass_Ecoli_core0.873920.87392
CO2t-22.80983-22.80983
CS6.007256.00725
CYTBD43.5989943.59899
D_LACt20.000000.00000
ENO14.7161414.71614
ETOHt2r0.000000.00000
EX_ac_e-0.000000.00000
\n", "
" ], "text/plain": [ " maximum minimum\n", "ACALD 0.00000 0.00000\n", "ACALDt -0.00000 0.00000\n", "ACKr -0.00000 0.00000\n", "ACONTa 6.00725 6.00725\n", "ACONTb 6.00725 6.00725\n", "ACt2r 0.00000 0.00000\n", "ADK1 -0.00000 0.00000\n", "AKGDH 5.06438 5.06438\n", "AKGt2r 0.00000 0.00000\n", "ALCD2x 0.00000 0.00000\n", "ATPM 8.39000 8.39000\n", "ATPS4r 45.51401 45.51401\n", "Biomass_Ecoli_core 0.87392 0.87392\n", "CO2t -22.80983 -22.80983\n", "CS 6.00725 6.00725\n", "CYTBD 43.59899 43.59899\n", "D_LACt2 0.00000 0.00000\n", "ENO 14.71614 14.71614\n", "ETOHt2r 0.00000 0.00000\n", "EX_ac_e -0.00000 0.00000" ] }, "execution_count": 13, "metadata": {}, "output_type": "execute_result" } ], "source": [ "fva_result = cobra.flux_analysis.flux_variability_analysis(\n", " model, model.reactions[:20])\n", "pandas.DataFrame.from_dict(fva_result).T.round(5)" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "Setting parameter fraction_of_optimium=0.90 would give the flux ranges for reactions at 90% optimality." ] }, { "cell_type": "code", "execution_count": 14, "metadata": { "collapsed": false }, "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", "
maximumminimum
ACALD0.00000-2.54237
ACALDt-0.00000-2.54237
ACKr-0.00000-3.81356
ACONTa8.894520.84859
ACONTb8.894520.84859
ACt2r0.00000-3.81356
ADK117.161000.00000
AKGDH8.045930.00000
AKGt2r0.00000-1.43008
ALCD2x0.00000-2.21432
ATPM25.551008.39000
ATPS4r59.3810634.82562
Biomass_Ecoli_core0.873920.78653
CO2t-15.20653-26.52885
CS8.894520.84859
CYTBD51.2390935.98486
D_LACt20.00000-2.14512
ENO16.732528.68659
ETOHt2r0.00000-2.21432
EX_ac_e3.813560.00000
\n", "
" ], "text/plain": [ " maximum minimum\n", "ACALD 0.00000 -2.54237\n", "ACALDt -0.00000 -2.54237\n", "ACKr -0.00000 -3.81356\n", "ACONTa 8.89452 0.84859\n", "ACONTb 8.89452 0.84859\n", "ACt2r 0.00000 -3.81356\n", "ADK1 17.16100 0.00000\n", "AKGDH 8.04593 0.00000\n", "AKGt2r 0.00000 -1.43008\n", "ALCD2x 0.00000 -2.21432\n", "ATPM 25.55100 8.39000\n", "ATPS4r 59.38106 34.82562\n", "Biomass_Ecoli_core 0.87392 0.78653\n", "CO2t -15.20653 -26.52885\n", "CS 8.89452 0.84859\n", "CYTBD 51.23909 35.98486\n", "D_LACt2 0.00000 -2.14512\n", "ENO 16.73252 8.68659\n", "ETOHt2r 0.00000 -2.21432\n", "EX_ac_e 3.81356 0.00000" ] }, "execution_count": 14, "metadata": {}, "output_type": "execute_result" } ], "source": [ "fva_result = cobra.flux_analysis.flux_variability_analysis(\n", " model, model.reactions[:20], fraction_of_optimum=0.9)\n", "pandas.DataFrame.from_dict(fva_result).T.round(5)" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "### Running FVA in summary methods\n", "Flux variability analysis can also be embedded in calls to summary methods. For instance, the expected variability in substrate consumption and product formation can be quickly found by" ] }, { "cell_type": "code", "execution_count": 15, "metadata": { "collapsed": false }, "outputs": [ { "name": "stdout", "output_type": "stream", "text": [ "IN FLUXES OUT FLUXES OBJECTIVES \n", "o2_e -21.80 ± 1.91 h2o_e 27.86 ± 2.86 Biomass_Ecoli_core 0.874\n", "glc__D_e -9.76 ± 0.24 co2_e 21.81 ± 2.86 \n", "nh4_e -4.84 ± 0.32 h_e 19.51 ± 2.86 \n", "pi_e -3.13 ± 0.08 for_e 2.86 ± 2.86 \n", " ac_e 0.95 ± 0.95 \n", " acald_e 0.64 ± 0.64 \n", " pyr_e 0.64 ± 0.64 \n", " etoh_e 0.55 ± 0.55 \n", " lac__D_e 0.54 ± 0.54 \n", " succ_e 0.42 ± 0.42 \n", " akg_e 0.36 ± 0.36 \n", " glu__L_e 0.32 ± 0.32 \n" ] } ], "source": [ "model.optimize()\n", "model.summary(fva=0.95)" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "Similarly, variability in metabolite mass balances can also be checked with flux variability analysis" ] }, { "cell_type": "code", "execution_count": 16, "metadata": { "collapsed": false }, "outputs": [ { "name": "stdout", "output_type": "stream", "text": [ "PRODUCING REACTIONS -- Pyruvate\n", "-------------------------------\n", " % FLUX RXN ID REACTION \n", " 85.0% 9.76 ± 0.24 GLCpts glc__D_e + pep_c --> g6p_c + pyr_c\n", " 15.0% 6.13 ± 6.13 PYK adp_c + h_c + pep_c --> atp_c + pyr_c\n", "\n", "CONSUMING REACTIONS -- Pyruvate\n", "-------------------------------\n", " % FLUX RXN ID REACTION \n", " 78.9% 11.34 ± 7.43 PDH coa_c + nad_c + pyr_c --> accoa_c + co2_c + nadh_c\n", " 21.1% 0.85 ± 0.02 Bioma... 1.496 3pg_c + 3.7478 accoa_c + 59.81 atp_c + 0.36...\n" ] } ], "source": [ "model.metabolites.pyr_c.summary(fva=0.95)" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "In these summary methods, the values are reported as a the center point +/- the range of the FVA solution, calculated from the maximum and minimum values." ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "## Running pFBA\n", "\n", "Parsimonious FBA (often written pFBA) finds a flux distribution which gives the optimal growth rate, but minimizes the total sum of flux. This involves solving two sequential linear programs, but is handled transparently by cobrapy. For more details on pFBA, please see [Lewis et al. (2010)](http://dx.doi.org/10.1038/msb.2010.47)." ] }, { "cell_type": "code", "execution_count": 17, "metadata": { "collapsed": false }, "outputs": [], "source": [ "FBA_sol = model.optimize()\n", "pFBA_sol = cobra.flux_analysis.optimize_minimal_flux(model)" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "These functions should give approximately the same objective value" ] }, { "cell_type": "code", "execution_count": 18, "metadata": { "collapsed": false }, "outputs": [ { "data": { "text/plain": [ "1.1102230246251565e-16" ] }, "execution_count": 18, "metadata": {}, "output_type": "execute_result" } ], "source": [ "abs(FBA_sol.f - pFBA_sol.f)" ] } ], "metadata": { "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.4.4" } }, "nbformat": 4, "nbformat_minor": 0 }