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.0142433, -0.0233536, -0.0122654, -0.0162743,…
  :ACALDt                   => [-0.00347733, -0.00915342, -0.00150051, -0.00498…
  :ACKr                     => [-0.0145782, -0.0113297, -0.015535, -0.0127172, …
  :ACONTa                   => [6.20856, 6.33897, 6.13551, 6.16939, 6.22453, 6.…
  :ACONTb                   => [6.20856, 6.33897, 6.13551, 6.16939, 6.22453, 6.…
  :ACt2r                    => [-0.0145782, -0.0113297, -0.015535, -0.0127172, …
  :ADK1                     => [0.0229999, 0.04163, 0.0255975, 0.0244223, 0.017…
  :AKGDH                    => [4.79196, 5.23481, 4.59487, 4.67697, 4.82674, 4.…
  :AKGt2r                   => [-0.00189919, -0.00192017, -0.00158513, -0.00163…
  :ALCD2x                   => [-0.010766, -0.0142002, -0.0107649, -0.011294, -…
  :ATPM                     => [8.43396, 8.42513, 8.43447, 8.42909, 8.43828, 8.…
  :ATPS4r                   => [44.9355, 44.6973, 45.0723, 44.9954, 44.8902, 44…
  :BIOMASS_Ecoli_core_w_GAM => [0.865491, 0.865513, 0.865484, 0.865455, 0.86544…
  :CO2t                     => [-22.9553, -22.9067, -22.9611, -22.9603, -22.941…
  :CS                       => [6.20856, 6.33897, 6.13551, 6.16939, 6.22453, 6.…
  :CYTBD                    => [43.9585, 43.8659, 43.9719, 43.9638, 43.9285, 43…
  :D_LACt2                  => [-0.00864784, -0.0127269, -0.0089709, -0.0091644…
  :ENO                      => [14.8869, 15.0263, 14.8107, 14.8464, 14.91, 14.8…
  :ETOHt2r                  => [-0.010766, -0.0142002, -0.0107649, -0.011294, -…
  ⋮                         => ⋮

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.9792  -22.9553
 21.9329  -22.9067
 21.9859  -22.9611
 21.9819  -22.9603
 21.9642  -22.9415
 21.9729  -22.9492
 21.9882  -22.9661
 21.9863  -22.9622
 22.0293  -23.0088
 21.9793  -22.9571
  ⋮       
 21.9465  -22.9154
 21.9441  -22.9152
 21.9414  -22.9095
 21.8586  -22.8109
 21.9227  -22.8899
 21.9238  -22.8887
 21.9192  -22.8809
 21.9475  -22.9199
 21.9536  -22.9334

This page was generated using Literate.jl.