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.0122753, -0.0118483, -0.00925479, -0.010066,…
  :ACALDt                   => [-0.00590139, -0.00340045, -0.00450486, -0.00305…
  :ACKr                     => [-0.0198697, -0.0209617, -0.017076, -0.0188654, …
  :ACONTa                   => [5.97863, 6.07626, 6.16074, 6.08592, 5.97739, 6.…
  :ACONTb                   => [5.97863, 6.07626, 6.16074, 6.08592, 5.97739, 6.…
  :ACt2r                    => [-0.0198697, -0.0209617, -0.017076, -0.0188654, …
  :ADK1                     => [0.0479306, 0.0324859, 0.0293246, 0.0290303, 0.0…
  :AKGDH                    => [4.27207, 4.51675, 4.63767, 4.46307, 4.47642, 4.…
  :AKGt2r                   => [-0.00104636, -0.00224753, -0.00271545, -0.00240…
  :ALCD2x                   => [-0.00637387, -0.00844781, -0.00474994, -0.00701…
  :ATPM                     => [8.4322, 8.42634, 8.43923, 8.4159, 8.40586, 8.43…
  :ATPS4r                   => [45.403, 45.1762, 45.0207, 45.1239, 45.1952, 45.…
  :BIOMASS_Ecoli_core_w_GAM => [0.865418, 0.865501, 0.865513, 0.8655, 0.865536,…
  :CO2t                     => [-22.9648, -22.9489, -22.9782, -22.9643, -22.850…
  :CS                       => [5.97863, 6.07626, 6.16074, 6.08592, 5.97739, 6.…
  :CYTBD                    => [43.974, 43.9459, 44.0005, 43.9795, 43.7381, 43.…
  :D_LACt2                  => [-0.00808285, -0.00783509, -0.00683705, -0.00768…
  :ENO                      => [14.6561, 14.7553, 14.8277, 14.759, 14.7128, 14.…
  :ETOHt2r                  => [-0.00637387, -0.00844781, -0.00474994, -0.00701…
  ⋮                         => ⋮

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.987   -22.9648
 21.973   -22.9489
 22.0002  -22.9782
 21.9897  -22.9643
 21.8691  -22.8504
 21.9496  -22.9241
 21.9926  -22.9695
 21.9834  -22.9603
 21.9735  -22.9495
 21.9354  -22.9112
  ⋮       
 21.9584  -22.9429
 22.0014  -22.9774
 21.9629  -22.9446
 22.0641  -23.0442
 21.9657  -22.9358
 21.9893  -22.9722
 21.9645  -22.9471
 21.8948  -22.8798
 21.9562  -22.9397

This page was generated using Literate.jl.