CycleFreeFlux
CycleFreeFlux essentially defines a L1-parsimonious model which can be used to run a cycle-free FBA and FVA. In COBREXA, this is best done by reusing linear_parsimonious_flux_balance_analysis.
First, let's get a model, create a constraint tree with the model, and ask for explicitly materializing constraints for the exchanges:
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")
cs = flux_balance_constraints(model, interface = :identifier_prefixes)ConstraintTrees.Tree{ConstraintTrees.Constraint} with 5 elements:
:coupling => ConstraintTrees.Tree{ConstraintTrees.Constraint}(#= 0 …
:flux_stoichiometry => ConstraintTrees.Tree{ConstraintTrees.Constraint}(#= 72…
:fluxes => ConstraintTrees.Tree{ConstraintTrees.Constraint}(#= 95…
:interface => ConstraintTrees.Tree{ConstraintTrees.Constraint}(#= 3 …
:objective => ConstraintTrees.Constraint(ConstraintTrees.LinearValue…We will also need some existing solution of the model – CycleFreeFlux algorithm uses this one as a reference for fixing the exchange reaction flux.
some_flux =
optimized_values(cs, objective = cs.objective.value, optimizer = HiGHS.Optimizer)ConstraintTrees.Tree{Float64} with 5 elements:
:coupling => ConstraintTrees.Tree{Float64}(#= 0 elements =#)
:flux_stoichiometry => ConstraintTrees.Tree{Float64}(#= 72 elements =#)
:fluxes => ConstraintTrees.Tree{Float64}(#= 95 elements =#)
:interface => ConstraintTrees.Tree{Float64}(#= 3 elements =#)
:objective => 0.873922(Ideally, we should use a solving method that gives a more unique flux, but for this example a simple FBA optimum will do.)
With this in hand, we can start the CycleFreeFlux workflow by placing constraints on exchange reactions in a linear-parsimonious model:
import ConstraintTrees as C
cs = linear_parsimonious_flux_balance_constraints(model)
cs *=
:fixed_exchanges^C.ConstraintTree(
k => C.Constraint(cs.fluxes[k].value, relative_tolerance_bound(0.999)(v)) for
(k, v) in some_flux.interface.exchanges
)ConstraintTrees.Tree{ConstraintTrees.Constraint} with 9 elements:
:coupling => ConstraintTrees.Tree{ConstraintTrees.Constraint}…
:directional_flux_balance => ConstraintTrees.Tree{ConstraintTrees.Constraint}…
:fixed_exchanges => ConstraintTrees.Tree{ConstraintTrees.Constraint}…
:flux_stoichiometry => ConstraintTrees.Tree{ConstraintTrees.Constraint}…
:fluxes => ConstraintTrees.Tree{ConstraintTrees.Constraint}…
:fluxes_forward => ConstraintTrees.Tree{ConstraintTrees.Constraint}…
:fluxes_reverse => ConstraintTrees.Tree{ConstraintTrees.Constraint}…
:objective => ConstraintTrees.Constraint(ConstraintTrees.Linea…
:parsimonious_objective => ConstraintTrees.Constraint(ConstraintTrees.Linea…(We purposefully made the constraints a little less strict by using relative_tolerance_bound – the toy E. coli model would otherwise display no variability at all.)
Now we can get a L1-parsimonious (thus cycle-free) solution of the model with the predefined exchanges:
cycle_free_flux = parsimonious_optimized_values(
cs,
objective = cs.objective.value,
objective_value = some_flux.objective,
parsimonious_objective = cs.parsimonious_objective.value,
optimizer = HiGHS.Optimizer,
)
cycle_free_flux.fluxesConstraintTrees.Tree{Float64} with 95 elements:
:ACALD => 0.0
:ACALDt => 0.0
:ACKr => 0.0
:ACONTa => 6.00725
:ACONTb => 6.00725
:ACt2r => 0.0
:ADK1 => 0.0
:AKGDH => 5.06438
:AKGt2r => 0.0
:ALCD2x => 0.0
:ATPM => 8.39
:ATPS4r => 45.514
:BIOMASS_Ecoli_core_w_GAM => 0.873922
:CO2t => -22.8098
:CS => 6.00725
:CYTBD => 43.599
:D_LACt2 => 0.0
:ENO => 14.7161
:ETOHt2r => 0.0
⋮ => ⋮CycleFreeFVA
With this in hand, we can also run the cycle-free flux variability analysis (again with an added bit of tolerances in both the objective and parsimonious bounds):
cs.objective.bound = C.Between(cycle_free_flux.objective * 0.999, Inf)
cs.parsimonious_objective.bound =
C.Between(0, cycle_free_flux.parsimonious_objective * 1.001)
var = constraints_variability(cs, cs.fluxes, optimizer = HiGHS.Optimizer)ConstraintTrees.Tree{Tuple{Union{Nothing, Float64}, Union{Nothing, Float64}}} with 95 elements:
:ACALD => (-0.0, 0.0)
:ACALDt => (-0.0, 0.0)
:ACKr => (-0.0, 0.0)
:ACONTa => (5.92071, 6.08235)
:ACONTb => (5.92071, 6.08235)
:ACt2r => (-0.0, 0.0)
:ADK1 => (-0.0, 0.108535)
:AKGDH => (4.48889, 5.14042)
:AKGt2r => (-0.0, 9.63674e-13)
:ALCD2x => (-0.0, 0.0)
:ATPM => (8.39, 8.51247)
:ATPS4r => (45.3147, 45.6379)
:BIOMASS_Ecoli_core_w_GAM => (0.873048, 0.873922)
:CO2t => (-22.8311, -22.7898)
:CS => (5.92071, 6.08235)
:CYTBD => (43.561, 43.6426)
:D_LACt2 => (1.39223e-13, 0.0)
:ENO => (14.6209, 14.7825)
:ETOHt2r => (-0.0, -5.77316e-14)
⋮ => ⋮CycleFree sampling
Naturally, we can also run flux sampling from the above model. To implement this, we follow the implementation of flux_sample –- first we generate the warmup:
warmup = vcat(
(
transpose(v) for (_, vs) in constraints_variability(
cs,
cs.fluxes,
optimizer = HiGHS.Optimizer,
output = (_, om) -> variable_vector(om),
output_type = Vector{Float64},
) for v in vs
)...,
)190×285 Matrix{Float64}:
7.40268 -0.0 4.63538 -15.9497 5.18557 … -0.0 -0.0 -0.0 -0.0
7.3821 -0.0 4.57806 -15.927 5.24023 -0.0 -0.0 -0.0 -0.0
7.40268 -0.0 4.63538 -15.9497 5.18557 -0.0 -0.0 -0.0 -0.0
7.3821 -0.0 4.57806 -15.927 5.24023 -0.0 -0.0 -0.0 -0.0
7.40268 -0.0 4.63538 -15.9497 5.18557 -0.0 -0.0 -0.0 -0.0
7.3821 -0.0 4.57806 -15.927 5.24023 … -0.0 -0.0 -0.0 -0.0
7.3821 -0.0 4.57806 -15.927 5.24023 -0.0 -0.0 -0.0 -0.0
7.54579 -0.0 5.0732 -16.0886 4.74305 -0.0 -0.0 -0.0 -0.0
7.3821 -0.0 4.57806 -15.927 5.24023 -0.0 -0.0 -0.0 -0.0
7.54579 -0.0 5.0732 -16.0886 4.74305 -0.0 -0.0 -0.0 -0.0
⋮ ⋱
7.3821 -0.0 4.57806 -15.927 5.24023 -0.0 -0.0 -0.0 -0.0
7.40268 -0.0 4.63538 -15.9497 5.18557 -0.0 -0.0 -0.0 -0.0
7.54579 -0.0 5.0732 -16.0886 4.74305 -0.0 -0.0 -0.0 -0.0
7.54579 -0.0 5.0732 -16.0886 4.74305 -0.0 -0.0 -0.0 -0.0
7.3821 -2.27792e-13 4.57806 -15.927 5.24023 … -0.0 -0.0 -0.0 -0.0
7.54579 -0.0 5.0732 -16.0886 4.74305 -0.0 -0.0 -0.0 -0.0
7.3821 -2.27792e-13 4.57806 -15.927 5.24023 -0.0 -0.0 -0.0 -0.0
7.3821 -0.0 4.57806 -15.927 5.24023 -0.0 -0.0 -0.0 -0.0
7.54579 -0.0 5.0732 -16.0886 4.74305 -0.0 -0.0 -0.0 -0.0Next, we can run the sampling:
sample = sample_constraints(
sample_chain_achr,
cs,
start_variables = warmup,
seed = UInt(1234),
output = cs.fluxes,
n_chains = 10,
collect_iterations = collect(10:15),
)ConstraintTrees.Tree{Vector{Float64}} with 95 elements:
:ACALD => [-7.85929e-15, 3.12233e-15, -6.32648e-15, -9.255…
:ACALDt => [-3.35535e-15, 7.30498e-15, -6.93217e-16, -4.301…
:ACKr => [-2.25938e-14, -1.25307e-14, -1.87866e-14, -2.04…
:ACONTa => [5.99813, 6.00102, 5.99386, 5.9927, 5.99009, 5.9…
:ACONTb => [5.99813, 6.00102, 5.99386, 5.9927, 5.99009, 5.9…
:ACt2r => [-2.26604e-14, -1.26102e-14, -1.88598e-14, -2.05…
:ADK1 => [0.00366918, 0.00229605, 0.00255894, 0.00289447,…
:AKGDH => [4.97009, 4.92654, 4.93649, 4.92894, 4.95038, 4.…
:AKGt2r => [3.21028e-14, 4.66312e-14, 2.72193e-14, 1.34945e…
:ALCD2x => [-3.86995e-15, -3.7229e-15, -5.06453e-15, -4.348…
:ATPM => [8.391, 8.39171, 8.39097, 8.3908, 8.39019, 8.390…
:ATPS4r => [45.5143, 45.509, 45.5227, 45.5218, 45.5294, 45.…
:BIOMASS_Ecoli_core_w_GAM => [0.873127, 0.87315, 0.87316, 0.87315, 0.873184, …
:CO2t => [-22.8174, -22.8185, -22.8188, -22.8189, -22.819…
:CS => [5.99813, 6.00102, 5.99386, 5.9927, 5.99009, 5.9…
:CYTBD => [43.616, 43.6182, 43.6186, 43.619, 43.619, 43.62…
:D_LACt2 => [-2.08644e-14, -2.16803e-14, -2.07692e-14, -2.34…
:ENO => [14.6991, 14.7022, 14.6952, 14.6939, 14.6916, 14…
:ETOHt2r => [-4.90694e-15, -3.98458e-15, -5.63877e-15, -5.09…
⋮ => ⋮The results can be observed (and usually plotted) from the sample vectors, such as the one for oxygen exchange:
sample.EX_o2_e11400-element Vector{Float64}:
-21.80801279761615
-21.80909965602572
-21.80930963740098
-21.80947899277754
-21.809475488922107
-21.810421483414107
-21.80854524250612
-21.80880422807698
-21.809512140853517
-21.80952709530663
⋮
-21.810011276079326
-21.810020653415503
-21.808622072753835
-21.808747348916107
-21.809865327092407
-21.809303155329953
-21.81135072741717
-21.810153001483112
-21.81051289673462This page was generated using Literate.jl.