Parsimonious flux balance analysis (pFBA)

Parsimonious flux balance analysis attempts to find a realistic flux of a model, by trying to minimize squared sum of all fluxes while maintaining the reached optimum. COBREXA.jl implements it in function parsimonious_flux_balance_analysis (accompanied by vector- and dictionary-returning variants parsimonious_flux_balance_analysis_vec and parsimonious_flux_balance_analysis_dict).

As usual, we demonstrate the functionality on the E. coli model:

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

using COBREXA, Tulip, Clarabel

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

Because the parsimonious objective is quadratic, we need a an optimizer capable of solving quadratic programs.

As the simplest choice, we can use Clarabel.jl, but any any JuMP.jl-supported optimizer that supports quadratic programming will work.

Running of basic pFBA is perfectly analogous to running of FBA and other analyses. We add several modifications that improve the solution (using functions silence, and change_optimizer_attribute), and fix the glucose exchange (using change_constraint) in order to get a more reasonable result:

fluxes = parsimonious_flux_balance_analysis_dict(
    model,
    Clarabel.Optimizer;
    modifications = [
        silence, # optionally silence the optimizer (Clarabel is very verbose by default)
        change_constraint("R_EX_glc__D_e"; lb = -12, ub = -12), # fix glucose consumption rate
    ],
)

Dict{String, Float64} with 95 entries:
  "R_EX_fum_e"    => 5.14033e-12
  "R_ACONTb"      => 7.03277
  "R_GLNS"        => 0.270339
  "R_SUCOAS"      => -5.8921
  "R_TPI"         => 8.90908
  "R_EX_pi_e"     => -3.88931
  "R_PPC"         => 3.02966
  "R_O2t"         => 25.7859
  "R_G6PDH2r"     => 6.11782
  "R_TALA"        => 1.85013
  "R_PPCK"        => -6.28606e-13
  "R_EX_lac__D_e" => 2.08955e-11
  "R_PGL"         => 6.11782
  "R_H2Ot"        => -34.7096
  "R_GLNabc"      => -1.00922e-11
  "R_EX_co2_e"    => 27.0082
  "R_EX_gln__L_e" => 5.03705e-12
  "R_EX_nh4_e"    => -5.76498
  "R_MALt2_2"     => -1.27933e-11
  ⋮               => ⋮

Using different optimizers for linear and quadratic problems

It is quite useful to use specialized optimizers for specialized tasks in pFBA. In particular, one would usually require to get a precise solution from the linear programming part (where the precision is achievable), and trade off a little precision for vast improvements in computation time in the quadratic programming part.

In pFBA, we can use the modifications and qp_modifications parameters to switch and parametrize the solvers in the middle of the process, which allows us to implement precisely that improvement. We demonstrate the switching on a vector-returning variant of pFBA:

flux_vector = parsimonious_flux_balance_analysis_vec(
    model,
    Tulip.Optimizer; # start with Tulip
    modifications = [
        change_constraint("R_EX_glc__D_e"; lb = -12, ub = -12),
        change_optimizer_attribute("IPM_IterationsLimit", 500), # we may change Tulip-specific attributes here
    ],
    qp_modifications = [
        change_optimizer(Clarabel.Optimizer), # now switch to Clarabel (Tulip wouldn't be able to finish the computation)
        silence, # and make it quiet.
    ],
)

95-element Vector{Float64}:
 -1.6408859620599666e-11
 -5.102018392197851e-12
 -6.0394828066438044e-12
  7.032766647219417
  7.032766647214134
 -4.0464571756211966e-12
  3.588150191299765e-12
  5.892098543575038
 -8.959896920612044e-12
 -8.737191010316242e-12
  ⋮
  6.587160251335506e-12
  1.0087574266648148e-11
  5.8920996478608245
 -5.892098543573014
  1.8501299332083683
  1.2470903529652462e-11
  1.8501299332169412
  1.4684623223107778
  8.909080293451035

This page was generated using Literate.jl.