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.124328, -0.0116282, -0.0127713, -0.0154784, …
  :ACALDt                   => [-0.12031, -0.00484284, -0.00355829, -0.00489226…
  :ACKr                     => [-0.00809812, -0.0208583, -0.0162745, -0.014721,…
  :ACONTa                   => [6.01066, 5.70407, 5.99526, 6.15381, 6.1755, 5.9…
  :ACONTb                   => [6.01066, 5.70407, 5.99526, 6.15381, 6.1755, 5.9…
  :ACt2r                    => [-0.00809812, -0.0208583, -0.0162745, -0.014721,…
  :ADK1                     => [0.0249071, 0.0137628, 0.0592348, 0.0310834, 0.0…
  :AKGDH                    => [4.69007, 3.76438, 4.22458, 4.64432, 4.76718, 4.…
  :AKGt2r                   => [-0.000518366, -0.00249488, -9.42458e-5, -0.0014…
  :ALCD2x                   => [-0.00401802, -0.0067854, -0.00921299, -0.010586…
  :ATPM                     => [8.41584, 8.41199, 8.42288, 8.42691, 8.43086, 8.…
  :ATPS4r                   => [45.0615, 45.8243, 45.3144, 44.9985, 44.9965, 45…
  :BIOMASS_Ecoli_core_w_GAM => [0.865324, 0.865418, 0.865355, 0.865561, 0.86525…
  :CO2t                     => [-22.8404, -22.9889, -22.9858, -22.9602, -22.927…
  :CS                       => [6.01066, 5.70407, 5.99526, 6.15381, 6.1755, 5.9…
  :CYTBD                    => [43.5927, 44.0029, 44.018, 43.9628, 43.9085, 43.…
  :D_LACt2                  => [-0.00341631, -0.00456997, -0.0066076, -0.010702…
  :ENO                      => [14.7744, 14.379, 14.6634, 14.831, 14.8608, 14.6…
  :ETOHt2r                  => [-0.00401802, -0.0067854, -0.00921299, -0.010586…
  ⋮                         => ⋮

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.7963  -22.8404
 22.0015  -22.9889
 22.009   -22.9858
 21.9814  -22.9602
 21.9543  -22.9275
 21.9599  -22.9297
 21.9317  -22.9
 22.0079  -22.9852
 21.9601  -22.9364
 22.0288  -23.0135
  ⋮       
 21.9795  -22.9481
 21.9787  -22.9475
 21.9605  -22.9301
 21.9937  -22.9659
 21.9597  -22.9293
 21.9708  -22.9399
 21.9813  -22.9459
 21.9069  -22.8809
 21.9495  -22.9241

This page was generated using Literate.jl.