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.00812914, -0.0098026, -0.0186774, -0.0089507…
  :ACALDt                   => [-0.0007826, -0.00236039, -0.0114581, -0.0061269…
  :ACKr                     => [-0.0116891, -0.0186744, -0.015783, -0.0132327, …
  :ACONTa                   => [6.04314, 5.9706, 6.03092, 5.66849, 6.12118, 6.0…
  :ACONTb                   => [6.04314, 5.9706, 6.03092, 5.66849, 6.12118, 6.0…
  :ACt2r                    => [-0.0116891, -0.0186744, -0.015783, -0.0132327, …
  :ADK1                     => [0.0358674, 0.0395008, 0.0320222, 0.0418153, 0.0…
  :AKGDH                    => [4.5028, 4.2879, 4.43608, 3.58027, 4.56863, 4.45…
  :AKGt2r                   => [-0.00250927, -0.00174955, -0.0016111, -0.001675…
  :ALCD2x                   => [-0.00734654, -0.00744221, -0.00721934, -0.00282…
  :ATPM                     => [8.45021, 8.44118, 8.44268, 8.44661, 8.46325, 8.…
  :ATPS4r                   => [45.3739, 45.4317, 45.2539, 46.0575, 45.2494, 45…
  :BIOMASS_Ecoli_core_w_GAM => [0.865436, 0.865391, 0.865387, 0.865314, 0.86546…
  :CO2t                     => [-22.9413, -22.9631, -22.9482, -23.0339, -22.926…
  :CS                       => [6.04314, 5.9706, 6.03092, 5.66849, 6.12118, 6.0…
  :CYTBD                    => [43.9515, 43.9855, 43.9387, 44.0886, 43.9276, 43…
  :D_LACt2                  => [-0.00860253, -0.0075182, -0.0107408, -0.0098557…
  :ENO                      => [14.7136, 14.6442, 14.7121, 14.3318, 14.7946, 14…
  :ETOHt2r                  => [-0.00734654, -0.00744221, -0.00721934, -0.00282…
  ⋮                         => ⋮

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.9758  -22.9413
 21.9927  -22.9631
 21.9694  -22.9482
 22.0443  -23.0339
 21.9638  -22.926
 21.9843  -22.9597
 21.9741  -22.944
 21.9563  -22.9252
 22.0078  -22.9937
 21.9995  -22.972
  ⋮       
 21.9696  -22.9408
 21.9576  -22.922
 21.9458  -22.9086
 21.9652  -22.9302
 21.9249  -22.8862
 21.9535  -22.9191
 21.9553  -22.9183
 21.9458  -22.9105
 21.9541  -22.9196

This page was generated using Literate.jl.