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.0158589, -0.016158, -0.0102801, -0.0140256, …
  :ACALDt                   => [-0.00881421, -0.00728834, -0.00237504, -0.00807…
  :ACKr                     => [-0.0139521, -0.0220105, -0.0150958, -0.0111963,…
  :ACONTa                   => [6.10535, 6.47902, 5.94617, 5.99034, 6.08856, 6.…
  :ACONTb                   => [6.10535, 6.47902, 5.94617, 5.99034, 6.08856, 6.…
  :ACt2r                    => [-0.0139521, -0.0220105, -0.0150958, -0.0111963,…
  :ADK1                     => [0.0637421, 0.0512532, 0.0350986, 0.051915, 0.04…
  :AKGDH                    => [4.6362, 5.4739, 4.3828, 4.37782, 4.56038, 4.817…
  :AKGt2r                   => [-0.00254333, -0.0026107, -0.000545288, -0.00268…
  :ALCD2x                   => [-0.00704465, -0.00886963, -0.00790507, -0.00594…
  :ATPM                     => [8.43641, 8.46866, 8.42731, 8.43804, 8.43835, 8.…
  :ATPS4r                   => [45.2151, 44.4444, 45.3787, 45.4199, 45.1578, 45…
  :BIOMASS_Ecoli_core_w_GAM => [0.865442, 0.865416, 0.865733, 0.865457, 0.86541…
  :CO2t                     => [-22.9466, -22.8917, -22.9103, -22.9613, -22.939…
  :CS                       => [6.10535, 6.47902, 5.94617, 5.99034, 6.08856, 6.…
  :CYTBD                    => [43.924, 43.8714, 43.9092, 43.9617, 43.9391, 43.…
  :D_LACt2                  => [-0.0076368, -0.00729615, -0.00967269, -0.007175…
  :ENO                      => [14.7849, 15.1647, 14.6243, 14.6631, 14.7653, 14…
  :ETOHt2r                  => [-0.00704465, -0.00886963, -0.00790507, -0.00594…
  ⋮                         => ⋮

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.962   -22.9466
 21.9357  -22.8917
 21.9546  -22.9103
 21.9809  -22.9613
 21.9696  -22.939
 21.95    -22.9
 21.9912  -22.9742
 21.9736  -22.9474
 21.9679  -22.9399
 21.9573  -22.9427
  ⋮       
 21.9431  -22.9054
 21.9692  -22.9355
 21.9487  -22.9134
 21.9788  -22.9465
 21.9837  -22.9545
 21.9674  -22.9314
 21.953   -22.919
 22.0026  -22.9746
 21.9853  -22.9545

This page was generated using Literate.jl.