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.0624972, -0.014811, -0.0293074, -0.0121536, …
  :ACALDt                   => [-0.0563532, -0.00686748, -0.0210267, -0.0047627…
  :ACKr                     => [-0.0133212, -0.0222815, -0.0129283, -0.0122634,…
  :ACONTa                   => [6.03023, 6.14183, 6.06512, 5.93585, 6.06407, 6.…
  :ACONTb                   => [6.03023, 6.14183, 6.06512, 5.93585, 6.06407, 6.…
  :ACt2r                    => [-0.0133212, -0.0222815, -0.0129283, -0.0122634,…
  :ADK1                     => [0.03136, 0.0325063, 0.0345115, 0.0340099, 0.037…
  :AKGDH                    => [4.57749, 4.69788, 4.5449, 4.20853, 4.46754, 4.5…
  :AKGt2r                   => [-0.000802057, -0.00143212, -0.000958101, -0.001…
  :ALCD2x                   => [-0.00614399, -0.00794351, -0.00828073, -0.00739…
  :ATPM                     => [8.42835, 8.44628, 8.42567, 8.42549, 8.4312, 8.4…
  :ATPS4r                   => [45.1761, 45.0705, 45.1892, 45.468, 45.2682, 45.…
  :BIOMASS_Ecoli_core_w_GAM => [0.865442, 0.865608, 0.865565, 0.865617, 0.86556…
  :CO2t                     => [-22.9028, -22.9398, -22.9287, -22.9506, -22.950…
  :CS                       => [6.03023, 6.14183, 6.06512, 5.93585, 6.06407, 6.…
  :CYTBD                    => [43.794, 43.9326, 43.9118, 43.9827, 43.9789, 43.…
  :D_LACt2                  => [-0.00595511, -0.0091377, -0.00613819, -0.007121…
  :ENO                      => [14.7436, 14.8221, 14.7492, 14.6029, 14.7286, 14…
  :ETOHt2r                  => [-0.00614399, -0.00794351, -0.00828073, -0.00739…
  ⋮                         => ⋮

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 scatter plot:

[s.O2t s.CO2t]

380×2 Matrix{Float64}:
 21.897   -22.9028
 21.9663  -22.9398
 21.9559  -22.9287
 21.9913  -22.9506
 21.9894  -22.9504
 21.9665  -22.929
 21.9823  -22.9413
 21.9204  -22.8768
 22.0033  -22.9755
 21.9636  -22.9255
  ⋮       
 21.9326  -22.8867
 21.9677  -22.9269
 21.9566  -22.9141
 21.9657  -22.9283
 21.9748  -22.941
 21.9563  -22.9245
 21.9852  -22.9535
 21.9648  -22.9288
 21.9477  -22.9032

This page was generated using Literate.jl.