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 0x7f1bb8789550>
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.8739215069684909
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 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
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 0x7f1b92f03470>: 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 0x7f1b92f036d8>: 1}
In [11]:
model.optimize().f
Out[11]:
175.00000000002336
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 0x7f1b92f03470>: 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.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% 6.13 ± 6.13 PYK adp_c + h_c + pep_c --> atp_c + pyr_c
50.0% 9.76 ± 0.24 GLCpts glc__D_e + pep_c --> g6p_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.
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]:
6.072919944699606e-14