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.0458671, -0.0125858, -0.0234423, -0.00900513…
:ACALDt => [-0.0390564, -0.00348029, -0.0158087, -0.0063688…
:ACKr => [-0.0188499, -0.0166737, -0.0197513, -0.0267942,…
:ACONTa => [6.06709, 5.90429, 6.11702, 6.16589, 5.92549, 6.…
:ACONTb => [6.06709, 5.90429, 6.11702, 6.16589, 5.92549, 6.…
:ACt2r => [-0.0188499, -0.0166737, -0.0197513, -0.0267942,…
:ADK1 => [0.0269334, 0.0387848, 0.0196545, 0.0258474, 0.0…
:AKGDH => [4.57674, 4.32262, 4.58188, 4.58158, 4.66791, 4.…
:AKGt2r => [-0.00129507, -0.00355213, -0.00414714, -0.00133…
:ALCD2x => [-0.00681074, -0.0091055, -0.00763354, -0.002636…
:ATPM => [8.42355, 8.43912, 8.41399, 8.41422, 8.41616, 8.…
:ATPS4r => [45.0973, 45.55, 44.9865, 44.9146, 44.9904, 45.1…
:BIOMASS_Ecoli_core_w_GAM => [0.865354, 0.865398, 0.865358, 0.865328, 0.86530…
:CO2t => [-22.9055, -22.9042, -22.9151, -22.931, -22.7004…
:CS => [6.06709, 5.90429, 6.11702, 6.16589, 5.92549, 6.…
:CYTBD => [43.8251, 43.8787, 43.8735, 43.9391, 43.4323, 43…
:D_LACt2 => [-0.00823089, -0.0107554, -0.0125498, -0.0075721…
:ENO => [14.774, 14.5866, 14.8127, 14.8428, 14.7534, 14.…
:ETOHt2r => [-0.00681074, -0.0091055, -0.00763354, -0.002636…
⋮ => ⋮
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 scatter plot:
[s.O2t s.CO2t]
380×2 Matrix{Float64}:
21.9125 -22.9055
21.9394 -22.9042
21.9367 -22.9151
21.9695 -22.931
21.7161 -22.7004
21.9529 -22.9207
21.9508 -22.9167
21.9859 -22.9717
21.9424 -22.9008
21.9522 -22.9015
⋮
21.9592 -22.9261
21.9986 -22.9677
21.9692 -22.9353
21.9629 -22.9382
21.9859 -22.9556
21.9677 -22.9386
22.0294 -23.0145
21.9857 -22.9522
21.966 -22.9283
This page was generated using Literate.jl.