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.0333753, -0.0128153, -0.0141539, -0.0156375,…
  :ACALDt                   => [-0.0252167, -0.00329747, -0.00303379, -0.003873…
  :ACKr                     => [-0.0143741, -0.022429, -0.0245313, -0.0221085, …
  :ACONTa                   => [6.07803, 6.07503, 6.1175, 6.13723, 6.02992, 6.1…
  :ACONTb                   => [6.07803, 6.07503, 6.1175, 6.13723, 6.02992, 6.1…
  :ACt2r                    => [-0.0143741, -0.022429, -0.0245313, -0.0221085, …
  :ADK1                     => [0.0290988, 0.028934, 0.0301241, 0.0308649, 0.01…
  :AKGDH                    => [4.5442, 4.45925, 4.54524, 4.61519, 4.45937, 4.5…
  :AKGt2r                   => [-0.00156052, -0.00119021, -0.00107673, -0.00196…
  :ALCD2x                   => [-0.00815858, -0.00951781, -0.0111201, -0.011763…
  :ATPM                     => [8.43084, 8.45818, 8.43258, 8.45021, 8.45303, 8.…
  :ATPS4r                   => [45.1287, 45.2035, 45.0854, 45.0784, 45.2153, 45…
  :BIOMASS_Ecoli_core_w_GAM => [0.865401, 0.865384, 0.865423, 0.86539, 0.865341…
  :CO2t                     => [-22.9325, -22.9449, -22.9314, -22.935, -22.9117…
  :CS                       => [6.07803, 6.07503, 6.1175, 6.13723, 6.02992, 6.1…
  :CYTBD                    => [43.8958, 43.9495, 43.9168, 43.9231, 43.8859, 43…
  :D_LACt2                  => [-0.00776911, -0.00756868, -0.00866403, -0.00847…
  :ENO                      => [14.7681, 14.753, 14.8004, 14.8189, 14.7286, 14.…
  :ETOHt2r                  => [-0.00815858, -0.00951781, -0.0111201, -0.011763…
  ⋮                         => ⋮

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.9479  -22.9325
 21.9747  -22.9449
 21.9584  -22.9314
 21.9616  -22.935
 21.943   -22.9117
 21.973   -22.9453
 21.9818  -22.9555
 21.9764  -22.9558
 22.0033  -22.9732
 21.9817  -22.9547
  ⋮       
 21.984   -22.9631
 21.9492  -22.9148
 21.9931  -22.9618
 21.9585  -22.9263
 21.9659  -22.9415
 21.9547  -22.9229
 21.9466  -22.9122
 21.958   -22.9304
 21.9521  -22.9231

This page was generated using Literate.jl.