Parsimonious flux balance analysis

Here, we use parsimonious_flux_balance_analysis (pFBA) to find the optimal flux distribution in the E. coli "core" model. In essence, pFBA first uses FBA to find an optimal objective value for the model, and then minimizes the squared distance of the flux from the zero (i.e., minimizes its L2 norm). As the main benefit, this gives a unique (and possibly more realistic) solution to the model.

using COBREXA

download_model(
    "http://bigg.ucsd.edu/static/models/e_coli_core.json",
    "e_coli_core.json",
    "7bedec10576cfe935b19218dc881f3fb14f890a1871448fc19a9b4ee15b448d8",
)

"e_coli_core.json"

Notably, we need an optimizer that can solve quadratic (QP) models:

import Clarabel

import JSONFBCModels
model = load_model("e_coli_core.json")

JSONFBCModels.JSONFBCModel(#= 95 reactions, 72 metabolites =#)

The pFBA is, in its most default form, implemented in function parsimonious_flux_balance_analysis:

solution = parsimonious_flux_balance_analysis(model; optimizer = Clarabel.Optimizer)

ConstraintTrees.Tree{Float64} with 5 elements:
  :coupling               => ConstraintTrees.Tree{Float64}(#= 0 elements =#)
  :flux_stoichiometry     => ConstraintTrees.Tree{Float64}(#= 72 elements =#)
  :fluxes                 => ConstraintTrees.Tree{Float64}(#= 95 elements =#)
  :objective              => 0.873922
  :parsimonious_objective => 11414.2

solution.fluxes

ConstraintTrees.Tree{Float64} with 95 elements:
  :ACALD                    => -1.10007e-10
  :ACALDt                   => -3.18315e-11
  :ACKr                     => -4.99153e-11
  :ACONTa                   => 6.00725
  :ACONTb                   => 6.00725
  :ACt2r                    => -3.36167e-11
  :ADK1                     => -9.84745e-12
  :AKGDH                    => 5.06438
  :AKGt2r                   => -4.87392e-11
  :ALCD2x                   => -6.24263e-11
  :ATPM                     => 8.39
  :ATPS4r                   => 45.514
  :BIOMASS_Ecoli_core_w_GAM => 0.873922
  :CO2t                     => -22.8098
  :CS                       => 6.00725
  :CYTBD                    => 43.599
  :D_LACt2                  => -4.43927e-11
  :ENO                      => 14.7161
  :ETOHt2r                  => -4.18493e-11
  ⋮                         => ⋮

Using different solvers for the problem stages

Sometimes it is useful to employ a dedicated LP solver to find the solution to the original FBA problem, and then a dedicated QP solver to minimize the fluxes. We can set the optimizer and parsimonious optimizer separately using keyword arguments:

import HiGHS

solution = parsimonious_flux_balance_analysis(
    model;
    optimizer = HiGHS.Optimizer, # HiGHS is used only for LP here
    parsimonious_optimizer = Clarabel.Optimizer, # Clarabel is great for solving QPs
)

ConstraintTrees.Tree{Float64} with 5 elements:
  :coupling               => ConstraintTrees.Tree{Float64}(#= 0 elements =#)
  :flux_stoichiometry     => ConstraintTrees.Tree{Float64}(#= 72 elements =#)
  :fluxes                 => ConstraintTrees.Tree{Float64}(#= 95 elements =#)
  :objective              => 0.873922
  :parsimonious_objective => 11414.2

Using linear (L1) parsimonious constraints

For efficiency reasons, it is also possible to use a pFBA version that optimizes the L1 norm instead of the L2 one (i.e., minimizing a sum of absolute values instead of the sum of squares). In turn, the uniqueness property of the solution is lost, but we do not need a QP-capable optimizer at all:

linear_solution =
    linear_parsimonious_flux_balance_analysis(model; optimizer = HiGHS.Optimizer)

ConstraintTrees.Tree{Float64} with 7 elements:
  :coupling               => ConstraintTrees.Tree{Float64}(#= 0 elements =#)
  :flux_stoichiometry     => ConstraintTrees.Tree{Float64}(#= 72 elements =#)
  :fluxes                 => ConstraintTrees.Tree{Float64}(#= 95 elements =#)
  :fluxes_forward         => ConstraintTrees.Tree{Float64}(#= 95 elements =#)
  :fluxes_reverse         => ConstraintTrees.Tree{Float64}(#= 95 elements =#)
  :objective              => 0.873922
  :parsimonious_objective => 518.422

linear_solution.fluxes

ConstraintTrees.Tree{Float64} with 95 elements:
  :ACALD                    => 0.0
  :ACALDt                   => 0.0
  :ACKr                     => -4.65317e-14
  :ACONTa                   => 6.00725
  :ACONTb                   => 6.00725
  :ACt2r                    => -4.65317e-14
  :ADK1                     => 0.0
  :AKGDH                    => 5.06438
  :AKGt2r                   => 0.0
  :ALCD2x                   => 0.0
  :ATPM                     => 8.39
  :ATPS4r                   => 45.514
  :BIOMASS_Ecoli_core_w_GAM => 0.873922
  :CO2t                     => -22.8098
  :CS                       => 6.00725
  :CYTBD                    => 43.599
  :D_LACt2                  => 0.0
  :ENO                      => 14.7161
  :ETOHt2r                  => 0.0
  ⋮                         => ⋮

This page was generated using Literate.jl.