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.022002, -0.016054, -0.0131751, -0.0212693, -…
  :ACALDt                   => [-0.00848178, -0.00411332, -0.0043492, -0.006326…
  :ACKr                     => [-0.0128176, -0.0161339, -0.0128173, -0.0144995,…
  :ACONTa                   => [6.10737, 6.05675, 6.0334, 6.0977, 5.99983, 6.04…
  :ACONTb                   => [6.10737, 6.05675, 6.0334, 6.0977, 5.99983, 6.04…
  :ACt2r                    => [-0.0128176, -0.0161339, -0.0128173, -0.0144995,…
  :ADK1                     => [0.024502, 0.0366378, 0.0424833, 0.0293806, 0.03…
  :AKGDH                    => [4.62758, 4.51832, 4.47292, 4.62707, 4.42573, 4.…
  :AKGt2r                   => [-0.00325474, -0.00173299, -0.00213076, -0.00231…
  :ALCD2x                   => [-0.0135202, -0.0119407, -0.00882592, -0.0149432…
  :ATPM                     => [8.40682, 8.43451, 8.44251, 8.41195, 8.42563, 8.…
  :ATPS4r                   => [45.1231, 45.2162, 45.287, 45.129, 45.3573, 45.3…
  :BIOMASS_Ecoli_core_w_GAM => [0.865343, 0.865427, 0.865428, 0.865415, 0.86550…
  :CO2t                     => [-22.9437, -22.9371, -22.9527, -22.9297, -22.942…
  :CS                       => [6.10737, 6.05675, 6.0334, 6.0977, 5.99983, 6.04…
  :CYTBD                    => [43.9133, 43.9235, 43.9539, 43.8956, 43.935, 43.…
  :D_LACt2                  => [-0.011078, -0.00754119, -0.00779677, -0.0098178…
  :ENO                      => [14.7919, 14.737, 14.7092, 14.7829, 14.6729, 14.…
  :ETOHt2r                  => [-0.0135202, -0.0119407, -0.00882592, -0.0149432…
  ⋮                         => ⋮

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.9567  -22.9437
 21.9617  -22.9371
 21.977   -22.9527
 21.9478  -22.9297
 21.9675  -22.9422
 21.9824  -22.9521
 21.9762  -22.954
 21.9655  -22.9419
 21.9837  -22.9577
 21.9896  -22.9658
  ⋮       
 21.8773  -22.8493
 21.9206  -22.8935
 21.9653  -22.933
 21.8938  -22.8605
 21.969   -22.9377
 21.9521  -22.9307
 21.9628  -22.9351
 21.95    -22.9171
 21.9813  -22.9572

This page was generated using Literate.jl.