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.
[1]:
import cobra.test
model = cobra.test.create_test_model("textbook")
4.1. Running FBA¶
[2]:
solution = model.optimize()
print(solution)
<Solution 0.874 at 0x112eb3d30>
The Model.optimize() function will return a Solution object. A solution object has several attributes:
objective_value
: the objective valuestatus
: the status from the linear programming solverfluxes
: 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.
[3]:
solution.objective_value
[3]:
0.8739215069684307
The solvers that can be used with cobrapy are so fast that for many small to mid-size models computing the solution can be even faster than it takes to collect the values from the solver and convert to them python objects. With model.optimize
, we gather values for all reactions and metabolites and that can take a significant amount of time if done repeatedly. If we are only interested in the flux value of a single reaction or the objective, it is faster to instead use model.slim_optimize
which only does the optimization and returns the objective value leaving it up to you to fetch other values that you may need.
[4]:
%%time
model.optimize().objective_value
CPU times: user 3.84 ms, sys: 672 µs, total: 4.51 ms
Wall time: 6.16 ms
[4]:
0.8739215069684307
[5]:
%%time
model.slim_optimize()
CPU times: user 229 µs, sys: 19 µs, total: 248 µs
Wall time: 257 µs
[5]:
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.
[6]:
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
[7]:
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
[8]:
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.
[9]:
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.
[10]:
from cobra.util.solver import linear_reaction_coefficients
linear_reaction_coefficients(model)
[10]:
{<Reaction Biomass_Ecoli_core at 0x112eab4a8>: 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}
.
[11]:
# 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)
[11]:
{<Reaction ATPM at 0x112eab470>: 1.0}
[12]:
model.optimize().objective_value
[12]:
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.
[13]:
from cobra.flux_analysis import flux_variability_analysis
[14]:
flux_variability_analysis(model, model.reactions[:10])
[14]:
maximum | minimum | |
---|---|---|
ACALD | -2.208811e-30 | -5.247085e-14 |
ACALDt | 0.000000e+00 | -5.247085e-14 |
ACKr | 0.000000e+00 | -8.024953e-14 |
ACONTa | 2.000000e+01 | 2.000000e+01 |
ACONTb | 2.000000e+01 | 2.000000e+01 |
ACt2r | 0.000000e+00 | -8.024953e-14 |
ADK1 | 3.410605e-13 | 0.000000e+00 |
AKGDH | 2.000000e+01 | 2.000000e+01 |
AKGt2r | 0.000000e+00 | -2.902643e-14 |
ALCD2x | 0.000000e+00 | -4.547474e-14 |
Setting parameter fraction_of_optimium=0.90
would give the flux ranges for reactions at 90% optimality.
[15]:
cobra.flux_analysis.flux_variability_analysis(
model, model.reactions[:10], fraction_of_optimum=0.9)
[15]:
maximum | minimum | |
---|---|---|
ACALD | 0.000000e+00 | -2.692308 |
ACALDt | 0.000000e+00 | -2.692308 |
ACKr | 6.635712e-30 | -4.117647 |
ACONTa | 2.000000e+01 | 8.461538 |
ACONTb | 2.000000e+01 | 8.461538 |
ACt2r | 0.000000e+00 | -4.117647 |
ADK1 | 1.750000e+01 | 0.000000 |
AKGDH | 2.000000e+01 | 2.500000 |
AKGt2r | 2.651196e-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.
[16]:
loop_reactions = [model.reactions.FRD7, model.reactions.SUCDi]
flux_variability_analysis(model, reaction_list=loop_reactions, loopless=False)
[16]:
maximum | minimum | |
---|---|---|
FRD7 | 980.0 | 0.0 |
SUCDi | 1000.0 | 20.0 |
[17]:
flux_variability_analysis(model, reaction_list=loop_reactions, loopless=True)
[17]:
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
[18]:
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]
pi_e 0 [0, 0.171] 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]
Similarly, variability in metabolite mass balances can also be checked with flux variability analysis.
[19]:
model.metabolites.pyr_c.summary(fva=0.95)
PRODUCING REACTIONS -- Pyruvate (pyr_c)
---------------------------------------
% FLUX RANGE RXN ID REACTION
---- ------ ------------ ---------- ----------------------------------------
50% 10 [1.25, 18.8] PYK adp_c + h_c + pep_c --> atp_c + pyr_c
50% 10 [9.5, 10] GLCpts glc__D_e + pep_c --> g6p_c + pyr_c
0% 0 [0, 8.75] ME1 mal__L_c + nad_c --> co2_c + nadh_c +...
0% 0 [0, 8.75] ME2 mal__L_c + nadp_c --> co2_c + nadph_c...
CONSUMING REACTIONS -- Pyruvate (pyr_c)
---------------------------------------
% FLUX RANGE RXN ID REACTION
---- ------ ------------ ---------- ----------------------------------------
100% 20 [13, 28.8] PDH coa_c + nad_c + pyr_c --> accoa_c + c...
0% 0 [0, 8.75] PPS atp_c + h2o_c + pyr_c --> amp_c + 2.0...
0% 0 [0, 5.83] PFL coa_c + pyr_c --> accoa_c + for_c
0% 0 [0, 1.35] PYRt2 h_e + pyr_e <=> h_c + pyr_c
0% 0 [0, 1.13] LDH_D lac__D_c + nad_c <=> h_c + nadh_c + p...
0% 0 [0, 0.132] Biomass... 1.496 3pg_c + 3.7478 accoa_c + 59.81 ...
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).
[20]:
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.
[21]:
abs(fba_solution.fluxes["Biomass_Ecoli_core"] - pfba_solution.fluxes[
"Biomass_Ecoli_core"])
[21]:
7.7715611723760958e-16
4.5. Running geometric FBA¶
Geometric FBA finds a unique optimal flux distribution which is central to the range of possible fluxes. For more details on geometric FBA, please see K Smallbone, E Simeonidis (2009).
[22]:
geometric_fba_sol = cobra.flux_analysis.geometric_fba(model)
geometric_fba_sol
[22]:
fluxes | reduced_costs | |
---|---|---|
ACALD | 1.785214e-14 | 0.0 |
ACALDt | 1.785214e-14 | 0.0 |
ACKr | 0.000000e+00 | 0.0 |
ACONTa | 6.007250e+00 | 0.0 |
ACONTb | 6.007250e+00 | 0.0 |
... | ... | ... |
TALA | 1.496984e+00 | 0.0 |
THD2 | 1.522652e-14 | 0.0 |
TKT1 | 1.496984e+00 | 0.0 |
TKT2 | 1.181498e+00 | 0.0 |
TPI | 7.477382e+00 | 0.0 |
95 rows × 2 columns