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.0458671, -0.0125858, -0.0234423, -0.00900513…
  :ACALDt                   => [-0.0390564, -0.00348029, -0.0158087, -0.0063688…
  :ACKr                     => [-0.0188499, -0.0166737, -0.0197513, -0.0267942,…
  :ACONTa                   => [6.06709, 5.90429, 6.11702, 6.16589, 5.92549, 6.…
  :ACONTb                   => [6.06709, 5.90429, 6.11702, 6.16589, 5.92549, 6.…
  :ACt2r                    => [-0.0188499, -0.0166737, -0.0197513, -0.0267942,…
  :ADK1                     => [0.0269334, 0.0387848, 0.0196545, 0.0258474, 0.0…
  :AKGDH                    => [4.57674, 4.32262, 4.58188, 4.58158, 4.66791, 4.…
  :AKGt2r                   => [-0.00129507, -0.00355213, -0.00414714, -0.00133…
  :ALCD2x                   => [-0.00681074, -0.0091055, -0.00763354, -0.002636…
  :ATPM                     => [8.42355, 8.43912, 8.41399, 8.41422, 8.41616, 8.…
  :ATPS4r                   => [45.0973, 45.55, 44.9865, 44.9146, 44.9904, 45.1…
  :BIOMASS_Ecoli_core_w_GAM => [0.865354, 0.865398, 0.865358, 0.865328, 0.86530…
  :CO2t                     => [-22.9055, -22.9042, -22.9151, -22.931, -22.7004…
  :CS                       => [6.06709, 5.90429, 6.11702, 6.16589, 5.92549, 6.…
  :CYTBD                    => [43.8251, 43.8787, 43.8735, 43.9391, 43.4323, 43…
  :D_LACt2                  => [-0.00823089, -0.0107554, -0.0125498, -0.0075721…
  :ENO                      => [14.774, 14.5866, 14.8127, 14.8428, 14.7534, 14.…
  :ETOHt2r                  => [-0.00681074, -0.0091055, -0.00763354, -0.002636…
  ⋮                         => ⋮

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.9125  -22.9055
 21.9394  -22.9042
 21.9367  -22.9151
 21.9695  -22.931
 21.7161  -22.7004
 21.9529  -22.9207
 21.9508  -22.9167
 21.9859  -22.9717
 21.9424  -22.9008
 21.9522  -22.9015
  ⋮       
 21.9592  -22.9261
 21.9986  -22.9677
 21.9692  -22.9353
 21.9629  -22.9382
 21.9859  -22.9556
 21.9677  -22.9386
 22.0294  -23.0145
 21.9857  -22.9522
 21.966   -22.9283

This page was generated using Literate.jl.