Flux sampling

Flux sampling gives an interesting statistical insight into the behavior of the model in the optimal feasible space, and the general "shape" of the optimal- or near-optimal set of feasible states of a given model.

For demonstration, we need the usual packages and models:

using COBREXA

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

import JSONFBCModels, HiGHS

model = load_model("e_coli_core.json")

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

Function flux_sample uses linear optimization to generate a set of warm-up points (by default, the method to generate the warm-up is basically FVA), and then runs the hit-and-run flux sampling algorithm on the near-optimal feasible space of the model:

s = flux_sample(
    model,
    optimizer = HiGHS.Optimizer,
    objective_bound = relative_tolerance_bound(0.99),
    n_chains = 2,
    collect_iterations = [10],
)

ConstraintTrees.Tree{Vector{Float64}} with 95 elements:
  :ACALD                    => [-0.027926, -0.0240032, -0.0146277, -0.013018, -…
  :ACALDt                   => [-0.0203869, -0.00567736, -0.00739976, -0.005616…
  :ACKr                     => [-0.0291238, -0.0246163, -0.0210431, -0.018874, …
  :ACONTa                   => [6.05355, 6.08472, 6.08291, 6.02655, 6.072, 5.41…
  :ACONTb                   => [6.05355, 6.08472, 6.08291, 6.02655, 6.072, 5.41…
  :ACt2r                    => [-0.0291238, -0.0246163, -0.0210431, -0.018874, …
  :ADK1                     => [0.0098783, 0.016713, 0.0223146, 0.0228432, 0.02…
  :AKGDH                    => [4.51313, 4.521, 4.54535, 4.33096, 4.50003, 2.87…
  :AKGt2r                   => [-0.00137416, -0.000961709, -0.00179263, -0.0015…
  :ALCD2x                   => [-0.00753913, -0.0183258, -0.00722792, -0.007401…
  :ATPM                     => [8.43943, 8.4519, 8.42703, 8.43573, 8.44496, 8.4…
  :ATPS4r                   => [45.3125, 45.1122, 45.1946, 45.2517, 45.2418, 46…
  :BIOMASS_Ecoli_core_w_GAM => [0.865797, 0.865481, 0.865678, 0.86549, 0.865615…
  :CO2t                     => [-22.9194, -22.9023, -22.9555, -22.9587, -22.958…
  :CS                       => [6.05355, 6.08472, 6.08291, 6.02655, 6.072, 5.41…
  :CYTBD                    => [43.8722, 43.8377, 43.9525, 43.9689, 43.9619, 44…
  :D_LACt2                  => [-0.0077662, -0.0108802, -0.0072397, -0.00678666…
  :ENO                      => [14.7547, 14.7818, 14.7634, 14.7022, 14.7411, 14…
  :ETOHt2r                  => [-0.00753913, -0.0183258, -0.00722792, -0.007401…
  ⋮                         => ⋮

The result is a tree of vectors of sampled states for each value; the order of the values in these vectors is fixed. You can thus e.g. create a good matrix for plotting the sample as 2D scatterplot:

[s.O2t s.CO2t]

380×2 Matrix{Float64}:
 21.9361  -22.9194
 21.9189  -22.9023
 21.9762  -22.9555
 21.9845  -22.9587
 21.9809  -22.9586
 22.0795  -23.0675
 21.9523  -22.9303
 21.9763  -22.9639
 21.9551  -22.933
 21.9546  -22.9323
  ⋮       
 21.9419  -22.9105
 21.939   -22.8995
 21.9646  -22.9289
 21.9623  -22.9288
 21.9524  -22.92
 21.9451  -22.9011
 21.96    -22.927
 21.9669  -22.9314
 21.9622  -22.9164

This page was generated using Literate.jl.