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.
import pandas
pandas.options.display.max_rows = 100
import cobra.test
model = cobra.test.create_test_model("textbook")
model.optimize()
<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:
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
model.solution.status
'optimal'
model.solution.f
0.8739215069684305
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.
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
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
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
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.
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.
model.objective
{<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}.
# 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
{<Reaction ATPM at 0x1165107b8>: 1}
model.optimize().f
174.99999999999997
The objective function can also be changed by setting Reaction.objective_coefficient directly.
model.reactions.get_by_id("ATPM").objective_coefficient = 0.
biomass_rxn.objective_coefficient = 1.
model.objective
{<Reaction Biomass_Ecoli_core at 0x116510828>: 1.0}
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.
fva_result = cobra.flux_analysis.flux_variability_analysis(
model, model.reactions[:20])
pandas.DataFrame.from_dict(fva_result).T.round(5)
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.
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)
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 |
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
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
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.
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).
FBA_sol = model.optimize()
pFBA_sol = cobra.flux_analysis.optimize_minimal_flux(model)
These functions should give approximately the same objective value
abs(FBA_sol.f - pFBA_sol.f)
1.1102230246251565e-16