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 8 elements:
:coupling => 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.LinearV…
:parsimonious_objective => ConstraintTrees.Constraint(ConstraintTrees.LinearV…(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, 0.0)
: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 => (-0.0, 6.0633e-14)
:ENO => (14.6209, 14.7825)
:ETOHt2r => (-0.0, 0.0)
⋮ => ⋮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×190 Matrix{Float64}:
7.3821 -0.0 4.57806 -15.927 … -0.0 -0.0 -0.0 -0.0
7.3821 -0.0 4.57806 -15.927 -0.0 -0.0 -0.0 -0.0
7.3821 -0.0 4.57806 -15.927 -0.0 -0.0 -0.0 -0.0
7.3821 -0.0 4.57806 -15.927 -0.0 -0.0 -0.0 -0.0
7.3821 -0.0 4.57806 -15.927 -0.0 -0.0 -0.0 -0.0
7.3821 -0.0 4.57806 -15.927 … -0.0 -0.0 -0.0 -0.0
7.3821 -0.0 4.57806 -15.927 -0.0 -0.0 -0.0 -0.0
7.54579 -0.0 5.0732 -16.0886 -0.0 -0.0 6.21725e-15 -0.0
7.3821 -0.0 4.57806 -15.927 -0.0 -0.0 -0.0 -0.0
7.54579 -0.0 5.0732 -16.0886 -0.0 -0.0 6.21725e-15 -0.0
⋮ ⋱
7.3821 -0.0 4.57806 -15.927 -0.0 -0.0 -0.0 -0.0
7.40268 -0.0 4.63538 -15.9497 -0.0 -0.0 -0.0 -0.0
7.54579 -0.0 5.0732 -16.0886 -0.0 -0.0 -0.0 -0.0
7.54579 -0.0 5.0732 -16.0886 -0.0 -0.0 -0.0 -0.0
7.3821 -5.86198e-14 4.57806 -15.927 … -0.0 -0.0 -0.0 -0.0
7.54579 -0.0 5.0732 -16.0886 -0.0 -0.0 -0.0 -0.0
7.3821 -5.86198e-14 4.57806 -15.927 -0.0 -0.0 -0.0 -0.0
7.3821 -0.0 4.57806 -15.927 -0.0 -0.0 -0.0 -0.0
7.54579 -0.0 5.0732 -16.0886 -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 => [-2.86477e-15, -5.80857e-15, -6.303e-15, -4.2159…
:ACALDt => [-1.69876e-15, -2.82671e-15, -3.48804e-15, 8.023…
:ACKr => [-5.1976e-15, -4.94937e-15, -4.77809e-15, -4.416…
:ACONTa => [6.0017, 5.99204, 5.99008, 5.98695, 6.00876, 5.9…
:ACONTb => [6.0017, 5.99204, 5.99008, 5.98695, 6.00876, 5.9…
:ACt2r => [-5.26414e-15, -4.94501e-15, -4.80578e-15, -4.51…
:ADK1 => [0.0041622, 0.00376928, 0.00249937, 0.00247738, …
:AKGDH => [4.99378, 4.91115, 4.90038, 4.89263, 5.06269, 4.…
:AKGt2r => [2.03481e-15, 3.56293e-15, 4.05086e-15, 2.00581e…
:ALCD2x => [-1.39845e-15, -3.11776e-15, -2.98776e-15, -5.18…
:ATPM => [8.39128, 8.39161, 8.39152, 8.3916, 8.39043, 8.3…
:ATPS4r => [45.5167, 45.529, 45.5312, 45.5407, 45.5107, 45.…
:BIOMASS_Ecoli_core_w_GAM => [0.873185, 0.873134, 0.873139, 0.873156, 0.87326…
:CO2t => [-22.8194, -22.821, -22.8213, -22.8214, -22.8177…
:CS => [6.0017, 5.99204, 5.99008, 5.98695, 6.00876, 5.9…
:CYTBD => [43.6198, 43.6232, 43.6237, 43.624, 43.6162, 43.…
:D_LACt2 => [-2.25102e-14, -2.26823e-14, -1.8228e-14, -2.554…
:ENO => [14.7032, 14.6931, 14.6912, 14.6882, 14.7111, 14…
:ETOHt2r => [-1.55563e-15, -3.30731e-15, -3.15731e-15, -5.35…
⋮ => ⋮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.809910675462316
-21.811593820339848
-21.81186604603884
-21.81199036604807
-21.808086247677952
-21.8117797245191
-21.811417309355413
-21.811108656969214
-21.81115053554036
-21.811432268848144
⋮
-21.811507003496565
-21.812028633557663
-21.81105149470561
-21.81086295219901
-21.810914110886465
-21.811205956007473
-21.813364572555663
-21.80963333194993
-21.81300524702459This page was generated using Literate.jl.