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 cobra.test
model = cobra.test.create_test_model("textbook")

4.1. Running FBA

In [2]:
solution = model.optimize()
print(solution)
<Solution 0.874 at 0x10efd5cc0>

The Model.optimize() function will return a Solution object. A solution object has several attributes:

  • objective_value: the objective value
  • status: the status from the linear programming solver
  • fluxes: a pandas series with flux indexed by reaction identifier. The flux for a reaction variable is the difference of the primal values for the forward and reverse reaction variables.
  • shadow_prices: a pandas series with shadow price indexed by the metabolite identifier.

For example, after the last call to model.optimize(), if the optimization succeeds it’s status will be optimal. In case the model is infeasible an error is raised.

In [3]:
solution.objective_value
Out[3]:
0.8739215069684307

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 [4]:
model.summary()
IN FLUXES        OUT FLUXES    OBJECTIVES
---------------  ------------  ----------------------
o2_e      21.8   h2o_e  29.2   Biomass_Ecol...  0.874
glc__D_e  10     co2_e  22.8
nh4_e      4.77  h_e    17.5
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 [5]:
model.metabolites.nadh_c.summary()
PRODUCING REACTIONS -- Nicotinamide adenine dinucleotide - reduced (nadh_c)
---------------------------------------------------------------------------
%       FLUX  RXN ID      REACTION
----  ------  ----------  --------------------------------------------------
42%    16     GAPD        g3p_c + nad_c + pi_c <=> 13dpg_c + h_c + nadh_c
24%     9.28  PDH         coa_c + nad_c + pyr_c --> accoa_c + co2_c + nadh_c
13%     5.06  AKGDH       akg_c + coa_c + nad_c --> co2_c + nadh_c + succ...
13%     5.06  MDH         mal__L_c + nad_c <=> h_c + nadh_c + oaa_c
8%      3.1   Biomass...  1.496 3pg_c + 3.7478 accoa_c + 59.81 atp_c + 0....

CONSUMING REACTIONS -- Nicotinamide adenine dinucleotide - reduced (nadh_c)
---------------------------------------------------------------------------
%       FLUX  RXN ID      REACTION
----  ------  ----------  --------------------------------------------------
100%   38.5   NADH16      4.0 h_c + nadh_c + q8_c --> 3.0 h_e + nad_c + q...

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

In [6]:
model.metabolites.atp_c.summary()
PRODUCING REACTIONS -- ATP (atp_c)
----------------------------------
%      FLUX  RXN ID      REACTION
---  ------  ----------  --------------------------------------------------
67%  45.5    ATPS4r      adp_c + 4.0 h_e + pi_c <=> atp_c + h2o_c + 3.0 h_c
23%  16      PGK         3pg_c + atp_c <=> 13dpg_c + adp_c
7%    5.06   SUCOAS      atp_c + coa_c + succ_c <=> adp_c + pi_c + succoa_c
3%    1.76   PYK         adp_c + h_c + pep_c --> atp_c + pyr_c

CONSUMING REACTIONS -- ATP (atp_c)
----------------------------------
%      FLUX  RXN ID      REACTION
---  ------  ----------  --------------------------------------------------
76%  52.3    Biomass...  1.496 3pg_c + 3.7478 accoa_c + 59.81 atp_c + 0....
12%   8.39   ATPM        atp_c + h2o_c --> adp_c + h_c + pi_c
11%   7.48   PFK         atp_c + f6p_c --> adp_c + fdp_c + h_c
0%    0.223  GLNS        atp_c + glu__L_c + nh4_c --> adp_c + gln__L_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 [7]:
biomass_rxn = model.reactions.get_by_id("Biomass_Ecoli_core")

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

In [8]:
from cobra.util.solver import linear_reaction_coefficients
linear_reaction_coefficients(model)
Out[8]:
{<Reaction Biomass_Ecoli_core at 0x10efce518>: 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 [9]:
# 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.
linear_reaction_coefficients(model)
Out[9]:
{<Reaction ATPM at 0x10efce4e0>: 1.0}
In [10]:
model.optimize().objective_value
Out[10]:
174.99999999999966

We can also have more complicated objectives including quadratic terms.

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 [11]:
from cobra.flux_analysis import flux_variability_analysis
In [12]:
flux_variability_analysis(model, model.reactions[:10])
Out[12]:
maximum minimum
ACALD 5.243759e-30 -4.809828e-14
ACALDt 0.000000e+00 -4.809828e-14
ACKr -9.253692e-30 -7.356207e-14
ACONTa 2.000000e+01 2.000000e+01
ACONTb 2.000000e+01 2.000000e+01
ACt2r 0.000000e+00 -7.356207e-14
ADK1 3.126388e-13 0.000000e+00
AKGDH 2.000000e+01 2.000000e+01
AKGt2r 0.000000e+00 -2.660756e-14
ALCD2x 0.000000e+00 -4.168517e-14

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

In [13]:
cobra.flux_analysis.flux_variability_analysis(
    model, model.reactions[:10], fraction_of_optimum=0.9)
Out[13]:
maximum minimum
ACALD 6.957398e-16 -2.692308
ACALDt 0.000000e+00 -2.692308
ACKr 1.919072e-15 -4.117647
ACONTa 2.000000e+01 8.461538
ACONTb 2.000000e+01 8.461538
ACt2r 1.714115e-15 -4.117647
ADK1 1.750000e+01 0.000000
AKGDH 2.000000e+01 2.500000
AKGt2r 9.070758e-16 -1.489362
ALCD2x 0.000000e+00 -2.333333

The standard FVA may contain loops, i.e. high absolute flux values that only can be high if they are allowed to participate in loops (a mathematical artifact that cannot happen in vivo). Use the loopless argument to avoid such loops. Below, we can see that FRD7 and SUCDi reactions can participate in loops but that this is avoided when using the looplesss FVA.

In [14]:
loop_reactions = [model.reactions.FRD7, model.reactions.SUCDi]
flux_variability_analysis(model, reaction_list=loop_reactions, loopless=False)
Out[14]:
maximum minimum
FRD7 980.0 0.0
SUCDi 1000.0 20.0
In [15]:
flux_variability_analysis(model, reaction_list=loop_reactions, loopless=True)
Out[15]:
maximum minimum
FRD7 0.0 0.0
SUCDi 20.0 20.0

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 [16]:
model.optimize()
model.summary(fva=0.95)
IN FLUXES                     OUT FLUXES                     OBJECTIVES
----------------------------  -----------------------------  ------------
id          Flux  Range       id          Flux  Range        ATPM  175
--------  ------  ----------  --------  ------  -----------
o2_e          60  [55.9, 60]  co2_e         60  [54.2, 60]
glc__D_e      10  [9.5, 10]   h2o_e         60  [54.2, 60]
nh4_e          0  [0, 0.673]  for_e          0  [0, 5.83]
                              h_e            0  [0, 5.83]
                              ac_e           0  [0, 2.06]
                              acald_e        0  [0, 1.35]
                              pyr_e          0  [0, 1.35]
                              etoh_e         0  [0, 1.17]
                              lac__D_e       0  [0, 1.13]
                              succ_e         0  [0, 0.875]
                              akg_e          0  [0, 0.745]
                              glu__L_e       0  [0, 0.673]
                              pi_e           0  [0, -0.171]

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

In [17]:
model.metabolites.pyr_c.summary(fva=0.95)
PRODUCING REACTIONS -- Pyruvate (pyr_c)
---------------------------------------
%       FLUX  RANGE       RXN ID      REACTION
----  ------  ----------  ----------  ----------------------------------------
52%    10     [0, -5.83]  GLCpts      glc__D_e + pep_c --> g6p_c + pyr_c
48%     9.12  [0, -5.83]  PYK         adp_c + h_c + pep_c --> atp_c + pyr_c

CONSUMING REACTIONS -- Pyruvate (pyr_c)
---------------------------------------
%       FLUX  RANGE       RXN ID      REACTION
----  ------  ----------  ----------  ----------------------------------------
100%   19.1   [0, 5.83]   PDH         coa_c + nad_c + pyr_c --> accoa_c + c...
0%      0     [0, 5.83]   Biomass...  1.496 3pg_c + 3.7478 accoa_c + 59.81 ...
0%      0     [0, 5.83]   FRUpts2     fru_e + pep_c --> f6p_c + pyr_c
0%      0     [0, 5.83]   LDH_D       lac__D_c + nad_c <=> h_c + nadh_c + p...
0%      0     [0, 5.83]   ME1         mal__L_c + nad_c --> co2_c + nadh_c +...
0%      0     [0, 5.83]   ME2         mal__L_c + nadp_c --> co2_c + nadph_c...
0%      0     [0, 5.83]   PFL         coa_c + pyr_c --> accoa_c + for_c
0%      0     [0, 5.83]   PPS         atp_c + h2o_c + pyr_c --> amp_c + 2.0...
0%      0     [0, 5.83]   PYRt2       h_e + pyr_e <=> h_c + pyr_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.

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 [18]:
model.objective = 'Biomass_Ecoli_core'
fba_solution = model.optimize()
pfba_solution = cobra.flux_analysis.pfba(model)

These functions should give approximately the same objective value.

In [19]:
abs(fba_solution.fluxes["Biomass_Ecoli_core"] - pfba_solution.fluxes[
    "Biomass_Ecoli_core"])
Out[19]:
2.2204460492503131e-16