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.0126928, -0.0132022, -0.0956138, -0.0186338,…
  :ACALDt                   => [-0.00101728, -0.00402865, -0.0923352, -0.006125…
  :ACKr                     => [-0.0189779, -0.0136157, -0.0190046, -0.0243587,…
  :ACONTa                   => [6.06916, 5.94618, 5.99766, 6.20658, 6.06928, 6.…
  :ACONTb                   => [6.06916, 5.94618, 5.99766, 6.20658, 6.06928, 6.…
  :ACt2r                    => [-0.0189779, -0.0136157, -0.0190046, -0.0243587,…
  :ADK1                     => [0.0175932, 0.0210534, 0.0142139, 0.0433463, 0.0…
  :AKGDH                    => [4.58688, 4.19225, 4.6211, 4.93626, 4.55905, 4.9…
  :AKGt2r                   => [-0.00357699, -0.00248902, -0.000188391, -0.0026…
  :ALCD2x                   => [-0.0116756, -0.00917355, -0.00327862, -0.012507…
  :ATPM                     => [8.42164, 8.41813, 8.41175, 8.42803, 8.42801, 8.…
  :ATPS4r                   => [45.1601, 45.4086, 45.1018, 44.9614, 45.2659, 45…
  :BIOMASS_Ecoli_core_w_GAM => [0.86552, 0.865447, 0.865307, 0.865516, 0.865503…
  :CO2t                     => [-22.9006, -22.9606, -22.8425, -22.8888, -22.954…
  :CS                       => [6.06916, 5.94618, 5.99766, 6.20658, 6.06928, 6.…
  :CYTBD                    => [43.8823, 43.9801, 43.6477, 43.8429, 43.9523, 43…
  :D_LACt2                  => [-0.00986381, -0.00758744, -0.00292644, -0.01059…
  :ENO                      => [14.7553, 14.6187, 14.7462, 14.9025, 14.7434, 14…
  :ETOHt2r                  => [-0.0116756, -0.00917355, -0.00327862, -0.012507…
  ⋮                         => ⋮

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.9412  -22.9006
 21.99    -22.9606
 21.8238  -22.8425
 21.9215  -22.8888
 21.9762  -22.9543
 21.9156  -22.8821
 21.9584  -22.9414
 21.9762  -22.9533
 21.938   -22.908
 21.9466  -22.9267
  ⋮       
 22.0085  -22.9814
 21.9548  -22.9247
 21.954   -22.9261
 21.9554  -22.9257
 22.1199  -23.1073
 21.9258  -22.8969
 21.9227  -22.8716
 21.9422  -22.9135
 21.939   -22.9105

This page was generated using Literate.jl.