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.0379445, -0.0154706, -0.0118701, -0.0184794,…
  :ACALDt                   => [-0.0293087, -0.00541103, -0.00187985, -0.004634…
  :ACKr                     => [-0.0186831, -0.0197096, -0.0195001, -0.0320747,…
  :ACONTa                   => [6.06373, 6.06492, 6.07075, 6.28881, 6.05432, 6.…
  :ACONTb                   => [6.06373, 6.06492, 6.07075, 6.28881, 6.05432, 6.…
  :ACt2r                    => [-0.0186831, -0.0197096, -0.0195001, -0.0320747,…
  :ADK1                     => [0.0308046, 0.0345466, 0.038156, 0.0340105, 0.03…
  :AKGDH                    => [4.58919, 4.54466, 4.55828, 5.15354, 4.52292, 4.…
  :AKGt2r                   => [-0.00139203, -0.00131027, -0.00125164, -0.00187…
  :ALCD2x                   => [-0.0086358, -0.0100596, -0.00999023, -0.0138451…
  :ATPM                     => [8.44619, 8.43197, 8.42848, 8.45354, 8.4358, 8.4…
  :ATPS4r                   => [45.1391, 45.1891, 45.1978, 44.8104, 45.2481, 45…
  :BIOMASS_Ecoli_core_w_GAM => [0.865344, 0.865439, 0.865445, 0.865555, 0.86545…
  :CO2t                     => [-22.9095, -22.9147, -22.9083, -22.8441, -22.916…
  :CS                       => [6.06373, 6.06492, 6.07075, 6.28881, 6.05432, 6.…
  :CYTBD                    => [43.847, 43.8982, 43.8966, 43.7736, 43.9038, 43.…
  :D_LACt2                  => [-0.00991379, -0.0101694, -0.0106216, -0.0117222…
  :ENO                      => [14.7647, 14.7511, 14.7552, 14.996, 14.7371, 14.…
  :ETOHt2r                  => [-0.0086358, -0.0100596, -0.00999023, -0.0138451…
  ⋮                         => ⋮

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.9235  -22.9095
 21.9491  -22.9147
 21.9483  -22.9083
 21.8868  -22.8441
 21.9519  -22.9167
 21.9794  -22.9532
 21.9523  -22.9178
 21.9566  -22.9232
 21.9133  -22.8774
 21.9487  -22.9139
  ⋮       
 21.9656  -22.9668
 22.0161  -23.0036
 21.9676  -22.9472
 21.9909  -22.9752
 21.9922  -22.9646
 21.9843  -22.9655
 21.9769  -22.9599
 22.0025  -23.0016
 21.9903  -22.9714

This page was generated using Literate.jl.