Simulating with FBA

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.


In [1]:
import pandas
pandas.options.display.max_rows = 100

import cobra.test
model = cobra.test.create_test_model("textbook")

Running FBA


In [2]:
model.optimize()


Out[2]:
<Solution 0.87 at 0x10f309b10>

The Model.optimize() function will return a Solution object, which will also be stored at model.solution. A solution object has several attributes:

  • f: the objective value
  • status: the status from the linear programming solver
  • x_dict: a dictionary of {reaction_id: flux_value} (also called "primal")
  • x: a list for x_dict
  • y_dict: a dictionary of {metabolite_id: dual_value}.
  • y: a list for y_dict

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


In [3]:
model.solution.status


Out[3]:
'optimal'

In [4]:
model.solution.f


Out[4]:
0.8739215069684304

Analyzing FBA solutions

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.


In [5]:
model.summary()


IN FLUXES                     OUT FLUXES                    OBJECTIVES          
o2_e       -21.80             h2o_e    29.18                Biomass_Ecoli_core    0.874
glc__D_e   -10.00             co2_e    22.81                                    
nh4_e       -4.77             h_e      17.53                                    
pi_e        -3.21                                                               

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


In [6]:
model.metabolites.nadh_c.summary()


PRODUCING REACTIONS -- Nicotinamide adenine dinucleotide - reduced
------------------------------------------------------------------
  %      FLUX   RXN ID                        REACTION                       
 41.6%     16     GAPD        g3p_c + nad_c + pi_c <=> 13dpg_c + h_c + nadh_c
 24.1%    9.3      PDH     coa_c + nad_c + pyr_c --> accoa_c + co2_c + nadh_c
 13.1%    5.1      MDH              mal__L_c + nad_c <=> h_c + nadh_c + oaa_c
 13.1%    5.1    AKGDH    akg_c + coa_c + nad_c --> co2_c + nadh_c + succoa_c
  8.0%    3.1 Bioma...   1.496 3pg_c + 3.7478 accoa_c + 59.81 atp_c + 0.36...

CONSUMING REACTIONS -- Nicotinamide adenine dinucleotide - reduced
------------------------------------------------------------------
  %      FLUX   RXN ID                        REACTION                       
100.0%    -39   NADH16   4.0 h_c + nadh_c + q8_c --> 3.0 h_e + nad_c + q8h2_c

Or to get a sense of the main energy production and consumption reactions


In [7]:
model.metabolites.atp_c.summary()


PRODUCING REACTIONS -- ATP
--------------------------
  %      FLUX   RXN ID                        REACTION                       
 66.6%     46   ATPS4r     adp_c + 4.0 h_e + pi_c <=> atp_c + h2o_c + 3.0 h_c
 23.4%     16      PGK                      3pg_c + atp_c <=> 13dpg_c + adp_c
  7.4%    5.1   SUCOAS     atp_c + coa_c + succ_c <=> adp_c + pi_c + succoa_c
  2.6%    1.8      PYK                  adp_c + h_c + pep_c --> atp_c + pyr_c

CONSUMING REACTIONS -- ATP
--------------------------
  %      FLUX   RXN ID                        REACTION                       
 76.5%    -52 Bioma...   1.496 3pg_c + 3.7478 accoa_c + 59.81 atp_c + 0.36...
 12.3%   -8.4     ATPM                   atp_c + h2o_c --> adp_c + h_c + pi_c
 10.9%   -7.5      PFK                  atp_c + f6p_c --> adp_c + fdp_c + h_c

Changing the Objectives

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.


In [8]:
biomass_rxn = model.reactions.get_by_id("Biomass_Ecoli_core")

Currently in the model, there is only one objective reaction (the biomass reaction), with an objective coefficient of 1.


In [9]:
model.objective


Out[9]:
{<Reaction Biomass_Ecoli_core at 0x10f085bd0>: 1.0}

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}.


In [10]:
# change the objective to ATPM
model.objective = "ATPM"

# The upper bound should be 1000, so that we get
# the actual optimal value
model.reactions.get_by_id("ATPM").upper_bound = 1000.
model.objective


Out[10]:
{<Reaction ATPM at 0x10f0b1dd0>: 1}

In [11]:
model.optimize().f


Out[11]:
175.0

The objective function can also be changed by setting Reaction.objective_coefficient directly.


In [12]:
model.reactions.get_by_id("ATPM").objective_coefficient = 0.
biomass_rxn.objective_coefficient = 1.

model.objective


Out[12]:
{<Reaction Biomass_Ecoli_core at 0x10f085bd0>: 1.0}

Running FVA

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.


In [13]:
fva_result = cobra.flux_analysis.flux_variability_analysis(
    model, model.reactions[:20])
pandas.DataFrame.from_dict(fva_result).T


Out[13]:
maximum minimum
ACALD 0.000000 0.000000
ACALDt 0.000000 0.000000
ACKr 0.000000 0.000000
ACONTa 6.007250 6.007250
ACONTb 6.007250 6.007250
ACt2r 0.000000 0.000000
ADK1 0.000000 0.000000
AKGDH 5.064376 5.064376
AKGt2r 0.000000 0.000000
ALCD2x 0.000000 0.000000
ATPM 8.390000 8.390000
ATPS4r 45.514010 45.514010
Biomass_Ecoli_core 0.873922 0.873922
CO2t -22.809833 -22.809833
CS 6.007250 6.007250
CYTBD 43.598985 43.598985
D_LACt2 0.000000 0.000000
ENO 14.716140 14.716140
ETOHt2r 0.000000 0.000000
EX_ac_e 0.000000 0.000000

Setting parameter fraction_of_optimium=0.90 would give the flux ranges for reactions at 90% optimality.


In [14]:
fva_result = cobra.flux_analysis.flux_variability_analysis(
    model, model.reactions[:20], fraction_of_optimum=0.9)
pandas.DataFrame.from_dict(fva_result).T


Out[14]:
maximum minimum
ACALD 0.000000 -2.542370
ACALDt 0.000000 -2.542370
ACKr 0.000000 -3.813556
ACONTa 8.894520 0.848587
ACONTb 8.894520 0.848587
ACt2r 0.000000 -3.813556
ADK1 17.161000 0.000000
AKGDH 8.045934 0.000000
AKGt2r 0.000000 -1.430083
ALCD2x 0.000000 -2.214323
ATPM 25.551000 8.390000
ATPS4r 59.381062 34.825618
Biomass_Ecoli_core 0.873922 0.786529
CO2t -15.206526 -26.528850
CS 8.894520 0.848587
CYTBD 51.239087 35.984865
D_LACt2 0.000000 -2.145125
ENO 16.732521 8.686588
ETOHt2r 0.000000 -2.214323
EX_ac_e 3.813556 0.000000

Running FVA in summary methods

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


In [15]:
model.summary(fva=0.95)


IN FLUXES                     OUT FLUXES                    OBJECTIVES          
o2_e        -21.80 ± 1.91     h2o_e       27.86 ± 2.86      Biomass_Ecoli_core    0.000
glc__D_e     -9.76 ± 0.24     co2_e       21.81 ± 2.86                          
nh4_e        -4.84 ± 0.32     h_e         19.51 ± 2.86                          
pi_e         -3.13 ± 0.08     for_e        2.86 ± 2.86                          
                              ac_e         0.95 ± 0.95                          
                              acald_e      0.64 ± 0.64                          
                              pyr_e        0.64 ± 0.64                          
                              etoh_e       0.55 ± 0.55                          
                              lac__D_e     0.54 ± 0.54                          
                              succ_e       0.42 ± 0.42                          
                              akg_e        0.36 ± 0.36                          
                              glu__L_e     0.32 ± 0.32                          

Similarly, variability in metabolite mass balances can also be checked with flux variability analysis


In [16]:
model.metabolites.pyr_c.summary(fva=0.95)


PRODUCING REACTIONS -- Pyruvate
-------------------------------
  %             FLUX   RXN ID                        REACTION                       
 50.0%   9.76 ± 0.24   GLCpts                     glc__D_e + pep_c --> g6p_c + pyr_c
 50.0%   6.13 ± 6.13      PYK                  adp_c + h_c + pep_c --> atp_c + pyr_c

CONSUMING REACTIONS -- Pyruvate
-------------------------------
  %             FLUX   RXN ID                        REACTION                       
100.0%  11.34 ± 7.43      PDH     coa_c + nad_c + pyr_c --> accoa_c + co2_c + nadh_c

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.

Running pFBA

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).


In [17]:
FBA_sol = model.optimize()
pFBA_sol = cobra.flux_analysis.optimize_minimal_flux(model)

These functions should give approximately the same objective value


In [18]:
abs(FBA_sol.f - pFBA_sol.f)


Out[18]:
0.0