Flux balance analysis (FBA)

We will use flux_balance_analysis and several related functions to find the optimal flux in the E. coli "core" model.

If it is not already present, download the model and load the package:

!isfile("e_coli_core.xml") &&
    download("http://bigg.ucsd.edu/static/models/e_coli_core.xml", "e_coli_core.xml")

using COBREXA

model = load_model("e_coli_core.xml")

Metabolic model of type SBMLModel
sparse([8, 10, 21, 43, 50, 51, 8, 9, 6, 12  …  33, 66, 68, 72, 23, 26, 33, 72, 22, 33], [1, 1, 1, 1, 1, 1, 2, 2, 3, 3  …  93, 93, 93, 93, 94, 94, 94, 94, 95, 95], [-1.0, 1.0, -1.0, 1.0, -1.0, 1.0, 1.0, -1.0, -1.0, 1.0  …  1.0, -1.0, 1.0, -1.0, -1.0, 1.0, 1.0, -1.0, -1.0, 1.0], 72, 95)
Number of reactions: 95
Number of metabolites: 72

To perform any optimization-based analysis, we need to use a linear programming solver (also called an optimizer). Any of the JuMP.jl-supported optimizers will work. Here, we will demonstrate Tulip.jl and GLPK; other solvers will likely work just as well.

using Tulip

solved_model = flux_balance_analysis(model, Tulip.Optimizer)

A JuMP Model
Maximization problem with:
Variables: 95
Objective function type: JuMP.AffExpr
`JuMP.AffExpr`-in-`MathOptInterface.EqualTo{Float64}`: 72 constraints
`JuMP.AffExpr`-in-`MathOptInterface.LessThan{Float64}`: 190 constraints
Model mode: AUTOMATIC
CachingOptimizer state: ATTACHED_OPTIMIZER
Solver name: Tulip
Names registered in the model: lbs, mb, ubs, x

solved_model is now an instance of optimized JuMP model. To get the variable values out manually, we can use JuMP.value function. Flux variables are stored as vector x:

using JuMP
value.(solved_model[:x])

95-element Vector{Float64}:
  -5.62439030120323e-9
  -4.201129368838846e-9
  -3.3923987790069394e-9
   6.007249566759449
   6.00724956675944
  -3.3924030713751446e-9
   2.8591669740742613e-8
   5.064375361904729
  -3.847904111734164e-9
  -1.4232688663206453e-9
   ⋮
   2.3807910379998874e-8
   2.532925798014968e-8
 494.2414417548114
  -5.064375361904718
   1.4969838025832385
   3.373717085310436e-7
   1.4969838025832798
   1.18149814026048
   7.4773819191983355

To simplify things, there is a variant of the FBA function that does this for us automatically:

flux_balance_analysis_vec(model, Tulip.Optimizer)

95-element Vector{Float64}:
  -5.62439030120323e-9
  -4.201129368838846e-9
  -3.3923987790069394e-9
   6.007249566759449
   6.00724956675944
  -3.3924030713751446e-9
   2.8591669740742613e-8
   5.064375361904729
  -3.847904111734164e-9
  -1.4232688663206453e-9
   ⋮
   2.3807910379998874e-8
   2.532925798014968e-8
 494.2414417548114
  -5.064375361904718
   1.4969838025832385
   3.373717085310436e-7
   1.4969838025832798
   1.18149814026048
   7.4773819191983355

Likewise, there is another variant that returns the fluxes annotated by reaction names, in a dictionary:

flux_balance_analysis_dict(model, Tulip.Optimizer)

Dict{String, Float64} with 95 entries:
  "R_EX_fum_e"    => 0.0
  "R_ACONTb"      => 6.00725
  "R_GLNS"        => 0.223462
  "R_SUCOAS"      => -5.06438
  "R_TPI"         => 7.47738
  "R_EX_pi_e"     => -3.2149
  "R_PPC"         => 2.50431
  "R_O2t"         => 21.7995
  "R_G6PDH2r"     => 4.95999
  "R_TALA"        => 1.49698
  "R_PPCK"        => 5.88441e-8
  "R_EX_lac__D_e" => 2.37701e-9
  "R_PGL"         => 4.95999
  "R_H2Ot"        => -29.1758
  "R_GLNabc"      => 0.0
  "R_EX_co2_e"    => 22.8098
  "R_EX_gln__L_e" => 0.0
  "R_EX_nh4_e"    => -4.76532
  "R_MALt2_2"     => 0.0
  ⋮               => ⋮

Switching solvers is easy, and may be useful in case we need advanced functionality or performance present only in certain solvers. To switch to GLPK, we simply load the package and use a different optimizer to run the analysis:

using GLPK
flux_balance_analysis_dict(model, GLPK.Optimizer)

Dict{String, Float64} with 95 entries:
  "R_EX_fum_e"    => 0.0
  "R_ACONTb"      => 6.00725
  "R_GLNS"        => 0.223462
  "R_SUCOAS"      => -5.06438
  "R_TPI"         => 7.47738
  "R_EX_pi_e"     => -3.2149
  "R_PPC"         => 2.50431
  "R_O2t"         => 21.7995
  "R_G6PDH2r"     => 4.95998
  "R_TALA"        => 1.49698
  "R_PPCK"        => -7.10543e-15
  "R_EX_lac__D_e" => 1.83673e-14
  "R_PGL"         => 4.95998
  "R_H2Ot"        => -29.1758
  "R_GLNabc"      => 0.0
  "R_EX_co2_e"    => 22.8098
  "R_EX_gln__L_e" => 0.0
  "R_EX_nh4_e"    => -4.76532
  "R_MALt2_2"     => 0.0
  ⋮               => ⋮

To get a shortened but useful overview of what was found in the analysis, you can use flux_summary function:

flux_summary(flux_balance_analysis_dict(model, GLPK.Optimizer))

Biomass
  R_BIOMASS_Ecoli_core_w_GAM: 0.8739
Import
  R_EX_o2_e:                  -21.7995
  R_EX_glc__D_e:              -10.0
  R_EX_nh4_e:                 -4.7653
  R_EX_pi_e:                  -3.2149
Export
  R_EX_h_e:                   17.5309
  R_EX_co2_e:                 22.8098
  R_EX_h2o_e:                 29.1758

This page was generated using Literate.jl.