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.0132678, -0.0138165, -0.0363481, -0.0133836,…
  :ACALDt                   => [-0.00511409, -0.00282928, -0.0285296, -0.003810…
  :ACKr                     => [-0.0157661, -0.0159658, -0.0158454, -0.0171813,…
  :ACONTa                   => [6.13668, 6.10373, 6.07865, 6.02761, 6.04388, 6.…
  :ACONTb                   => [6.13668, 6.10373, 6.07865, 6.02761, 6.04388, 6.…
  :ACt2r                    => [-0.0157661, -0.0159658, -0.0158454, -0.0171813,…
  :ADK1                     => [0.0296718, 0.0238618, 0.0266923, 0.0304407, 0.0…
  :AKGDH                    => [4.66005, 4.66128, 4.58138, 4.42137, 4.31609, 4.…
  :AKGt2r                   => [-0.00195275, -0.000343564, -0.00153574, -0.0019…
  :ALCD2x                   => [-0.00815374, -0.0109872, -0.00781848, -0.009572…
  :ATPM                     => [8.43899, 8.4596, 8.43038, 8.43953, 8.42455, 8.4…
  :ATPS4r                   => [45.1244, 45.3106, 45.1615, 45.3443, 45.3256, 45…
  :BIOMASS_Ecoli_core_w_GAM => [0.865421, 0.865456, 0.865396, 0.865398, 0.86542…
  :CO2t                     => [-22.9659, -23.0116, -22.9352, -22.9659, -22.991…
  :CS                       => [6.13668, 6.10373, 6.07865, 6.02761, 6.04388, 6.…
  :CYTBD                    => [43.9819, 44.0423, 43.8906, 43.9776, 44.0171, 43…
  :D_LACt2                  => [-0.00835994, -0.00546069, -0.00725549, -0.00833…
  :ENO                      => [14.8094, 14.7721, 14.7719, 14.7009, 14.7154, 14…
  :ETOHt2r                  => [-0.00815374, -0.0109872, -0.00781848, -0.009572…
  ⋮                         => ⋮

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.991   -22.9659
 22.0212  -23.0116
 21.9453  -22.9352
 21.9888  -22.9659
 22.0085  -22.9919
 21.9851  -22.9671
 21.9877  -22.9694
 21.9405  -22.9228
 22.0211  -23.0059
 21.9696  -22.9465
  ⋮       
 21.9342  -22.9081
 22.0103  -22.9942
 21.9976  -22.9741
 21.9277  -22.9082
 21.926   -22.8818
 21.9201  -22.8937
 21.989   -22.9734
 21.9551  -22.9287
 21.9819  -22.9468

This page was generated using Literate.jl.