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

4.1. Running FBA

In [2]:
model.optimize()
Out[2]:
<Solution 0.87 at 0x10ddd0080>

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

4.1.1. 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    AKGDH    akg_c + coa_c + nad_c --> co2_c + nadh_c + succoa_c
 13.1%    5.1      MDH              mal__L_c + nad_c <=> h_c + nadh_c + oaa_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

4.2. 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 0x116510828>: 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 0x1165107b8>: 1}
In [11]:
model.optimize().f
Out[11]:
174.99999999999997

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 0x116510828>: 1.0}

4.3. 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.round(5)
Out[13]:
maximum minimum
ACALD 0.00000 0.00000
ACALDt -0.00000 0.00000
ACKr -0.00000 0.00000
ACONTa 6.00725 6.00725
ACONTb 6.00725 6.00725
ACt2r 0.00000 0.00000
ADK1 -0.00000 0.00000
AKGDH 5.06438 5.06438
AKGt2r 0.00000 0.00000
ALCD2x 0.00000 0.00000
ATPM 8.39000 8.39000
ATPS4r 45.51401 45.51401
Biomass_Ecoli_core 0.87392 0.87392
CO2t -22.80983 -22.80983
CS 6.00725 6.00725
CYTBD 43.59899 43.59899
D_LACt2 0.00000 0.00000
ENO 14.71614 14.71614
ETOHt2r 0.00000 0.00000
EX_ac_e -0.00000 0.00000

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.round(5)
Out[14]:
maximum minimum
ACALD 0.00000 -2.54237
ACALDt -0.00000 -2.54237
ACKr -0.00000 -3.81356
ACONTa 8.89452 0.84859
ACONTb 8.89452 0.84859
ACt2r 0.00000 -3.81356
ADK1 17.16100 0.00000
AKGDH 8.04593 0.00000
AKGt2r 0.00000 -1.43008
ALCD2x 0.00000 -2.21432
ATPM 25.55100 8.39000
ATPS4r 59.38106 34.82562
Biomass_Ecoli_core 0.87392 0.78653
CO2t -15.20653 -26.52885
CS 8.89452 0.84859
CYTBD 51.23909 35.98486
D_LACt2 0.00000 -2.14512
ENO 16.73252 8.68659
ETOHt2r 0.00000 -2.21432
EX_ac_e 3.81356 0.00000

4.3.1. 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.optimize()
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.874
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
 85.0%   9.76 ± 0.24   GLCpts                     glc__D_e + pep_c --> g6p_c + pyr_c
 15.0%   6.13 ± 6.13      PYK                  adp_c + h_c + pep_c --> atp_c + pyr_c

CONSUMING REACTIONS -- Pyruvate
-------------------------------
  %             FLUX   RXN ID                        REACTION
 78.9%  11.34 ± 7.43      PDH     coa_c + nad_c + pyr_c --> accoa_c + co2_c + nadh_c
 21.1%   0.85 ± 0.02 Bioma...   1.496 3pg_c + 3.7478 accoa_c + 59.81 atp_c + 0.36...

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.

4.4. 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]:
1.1102230246251565e-16