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.027926, -0.0240032, -0.0146277, -0.013018, -…
:ACALDt => [-0.0203869, -0.00567736, -0.00739976, -0.005616…
:ACKr => [-0.0291238, -0.0246163, -0.0210431, -0.018874, …
:ACONTa => [6.05355, 6.08472, 6.08291, 6.02655, 6.072, 5.41…
:ACONTb => [6.05355, 6.08472, 6.08291, 6.02655, 6.072, 5.41…
:ACt2r => [-0.0291238, -0.0246163, -0.0210431, -0.018874, …
:ADK1 => [0.0098783, 0.016713, 0.0223146, 0.0228432, 0.02…
:AKGDH => [4.51313, 4.521, 4.54535, 4.33096, 4.50003, 2.87…
:AKGt2r => [-0.00137416, -0.000961709, -0.00179263, -0.0015…
:ALCD2x => [-0.00753913, -0.0183258, -0.00722792, -0.007401…
:ATPM => [8.43943, 8.4519, 8.42703, 8.43573, 8.44496, 8.4…
:ATPS4r => [45.3125, 45.1122, 45.1946, 45.2517, 45.2418, 46…
:BIOMASS_Ecoli_core_w_GAM => [0.865797, 0.865481, 0.865678, 0.86549, 0.865615…
:CO2t => [-22.9194, -22.9023, -22.9555, -22.9587, -22.958…
:CS => [6.05355, 6.08472, 6.08291, 6.02655, 6.072, 5.41…
:CYTBD => [43.8722, 43.8377, 43.9525, 43.9689, 43.9619, 44…
:D_LACt2 => [-0.0077662, -0.0108802, -0.0072397, -0.00678666…
:ENO => [14.7547, 14.7818, 14.7634, 14.7022, 14.7411, 14…
:ETOHt2r => [-0.00753913, -0.0183258, -0.00722792, -0.007401…
⋮ => ⋮
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.9361 -22.9194
21.9189 -22.9023
21.9762 -22.9555
21.9845 -22.9587
21.9809 -22.9586
22.0795 -23.0675
21.9523 -22.9303
21.9763 -22.9639
21.9551 -22.933
21.9546 -22.9323
⋮
21.9419 -22.9105
21.939 -22.8995
21.9646 -22.9289
21.9623 -22.9288
21.9524 -22.92
21.9451 -22.9011
21.96 -22.927
21.9669 -22.9314
21.9622 -22.9164
This page was generated using Literate.jl.