Split unidirectional reactions in models
By default, the constraint system produced by e.g. flux_balance_constraints assumes a single variable for each reaction flux that may be both positive and negative (depending on the reaction). This example explains several ways to "split" such bidirectional fluxes into unidirectional "forward" and "reverse" parts. This is useful in modeling of capacity constraints (such system can be found e.g. in enzyme_constrained_flux_balance_analysis and linear_parsimonious_flux_balance_analysis) and many other purposes.
Here we show how to create such system for the toy E. coli model:
using COBREXA
download_model(
"http://bigg.ucsd.edu/static/models/e_coli_core.json",
"e_coli_core.json",
"7bedec10576cfe935b19218dc881f3fb14f890a1871448fc19a9b4ee15b448d8",
)
import JSONFBCModels
import HiGHS
import ConstraintTrees as C
model = load_model("e_coli_core.json")JSONFBCModels.JSONFBCModel(#= 95 reactions, 72 metabolites =#)We will only work with the constraint representation of the model. The fluxes there are bidirectional:
c = flux_balance_constraints(model);
c.fluxesConstraintTrees.Tree{ConstraintTrees.Constraint} with 95 elements:
:ACALD => ConstraintTrees.Constraint(ConstraintTrees.Linea…
:ACALDt => ConstraintTrees.Constraint(ConstraintTrees.Linea…
:ACKr => ConstraintTrees.Constraint(ConstraintTrees.Linea…
:ACONTa => ConstraintTrees.Constraint(ConstraintTrees.Linea…
:ACONTb => ConstraintTrees.Constraint(ConstraintTrees.Linea…
:ACt2r => ConstraintTrees.Constraint(ConstraintTrees.Linea…
:ADK1 => ConstraintTrees.Constraint(ConstraintTrees.Linea…
:AKGDH => ConstraintTrees.Constraint(ConstraintTrees.Linea…
:AKGt2r => ConstraintTrees.Constraint(ConstraintTrees.Linea…
:ALCD2x => ConstraintTrees.Constraint(ConstraintTrees.Linea…
:ATPM => ConstraintTrees.Constraint(ConstraintTrees.Linea…
:ATPS4r => ConstraintTrees.Constraint(ConstraintTrees.Linea…
:BIOMASS_Ecoli_core_w_GAM => ConstraintTrees.Constraint(ConstraintTrees.Linea…
:CO2t => ConstraintTrees.Constraint(ConstraintTrees.Linea…
:CS => ConstraintTrees.Constraint(ConstraintTrees.Linea…
:CYTBD => ConstraintTrees.Constraint(ConstraintTrees.Linea…
:D_LACt2 => ConstraintTrees.Constraint(ConstraintTrees.Linea…
:ENO => ConstraintTrees.Constraint(ConstraintTrees.Linea…
:ETOHt2r => ConstraintTrees.Constraint(ConstraintTrees.Linea…
⋮ => ⋮As the simplest approach, we can allocate 2 sets of variables for the forward and reverse fluxes via sign_split_variables and connect them to the fluxes with sign_split_constraints. These functions ensure that the bounds of the unidirectional fluxes are within the expectable limit, and, respectively, that the original fluxes are constrained to match the sum of the split directions:
c += sign_split_variables(c.fluxes, positive = :fluxes_forward, negative = :fluxes_reverse);
c *=
:directional_flux_balance^sign_split_constraints(
positive = c.fluxes_forward,
negative = c.fluxes_reverse,
signed = c.fluxes,
)ConstraintTrees.Tree{ConstraintTrees.Constraint} with 7 elements:
:coupling => ConstraintTrees.Tree{ConstraintTrees.Constraint}…
:directional_flux_balance => 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…We can solve this system as usual and obvserve the results as usual
x = optimized_values(c, objective = c.objective.value, optimizer = HiGHS.Optimizer)
C.zip(tuple, x.fluxes, x.fluxes_forward, x.fluxes_reverse, Tuple)ConstraintTrees.Tree{Tuple} with 95 elements:
:ACALD => (0.0, 0.0, 0.0)
:ACALDt => (0.0, 1000.0, 1000.0)
:ACKr => (0.0, 1000.0, 1000.0)
:ACONTa => (6.00725, 6.00725, 0.0)
:ACONTb => (6.00725, 6.00725, 0.0)
:ACt2r => (0.0, 1000.0, 1000.0)
:ADK1 => (0.0, 0.0, 0.0)
:AKGDH => (5.06438, 5.06438, 0.0)
:AKGt2r => (0.0, 1000.0, 1000.0)
:ALCD2x => (0.0, 1000.0, 1000.0)
:ATPM => (8.39, 8.39, 0.0)
:ATPS4r => (45.514, 45.514, 0.0)
:BIOMASS_Ecoli_core_w_GAM => (0.873922, 0.873922, 0.0)
:CO2t => (-22.8098, 0.0, 22.8098)
:CS => (6.00725, 6.00725, 0.0)
:CYTBD => (43.599, 43.599, 0.0)
:D_LACt2 => (0.0, 1000.0, 1000.0)
:ENO => (14.7161, 14.7161, 0.0)
:ETOHt2r => (0.0, 1000.0, 1000.0)
⋮ => ⋮Simplifying the system using asymmetric construction
If used naively, the above construction uses 3 variables and many constraints for each flux, which is not quite efficient. To ameliorate the usage of solver resources, one may construct a slightly simpler (but asymmetric) system that only uses 2 variables:
c2 = flux_balance_constraints(model);
c2 += :fluxes_forward^unsigned_positive_contribution_variables(c2.fluxes);
c2 *=
:fluxes_reverse^unsigned_negative_contribution_constraints(c2.fluxes, c2.fluxes_forward);This way, only forward fluxes are allocated as variables, and reverse fluxes are "computed" as linearly dependent values. Additionally, since the bounds on the forward and reverse fluxes completely subsume the original bounds on fluxes, one can further simplify the system by removing the original bounds:
c2.fluxes = remove_bounds(c2.fluxes)ConstraintTrees.Tree{ConstraintTrees.Constraint} with 95 elements:
:ACALD => ConstraintTrees.Constraint(ConstraintTrees.Linea…
:ACALDt => ConstraintTrees.Constraint(ConstraintTrees.Linea…
:ACKr => ConstraintTrees.Constraint(ConstraintTrees.Linea…
:ACONTa => ConstraintTrees.Constraint(ConstraintTrees.Linea…
:ACONTb => ConstraintTrees.Constraint(ConstraintTrees.Linea…
:ACt2r => ConstraintTrees.Constraint(ConstraintTrees.Linea…
:ADK1 => ConstraintTrees.Constraint(ConstraintTrees.Linea…
:AKGDH => ConstraintTrees.Constraint(ConstraintTrees.Linea…
:AKGt2r => ConstraintTrees.Constraint(ConstraintTrees.Linea…
:ALCD2x => ConstraintTrees.Constraint(ConstraintTrees.Linea…
:ATPM => ConstraintTrees.Constraint(ConstraintTrees.Linea…
:ATPS4r => ConstraintTrees.Constraint(ConstraintTrees.Linea…
:BIOMASS_Ecoli_core_w_GAM => ConstraintTrees.Constraint(ConstraintTrees.Linea…
:CO2t => ConstraintTrees.Constraint(ConstraintTrees.Linea…
:CS => ConstraintTrees.Constraint(ConstraintTrees.Linea…
:CYTBD => ConstraintTrees.Constraint(ConstraintTrees.Linea…
:D_LACt2 => ConstraintTrees.Constraint(ConstraintTrees.Linea…
:ENO => ConstraintTrees.Constraint(ConstraintTrees.Linea…
:ETOHt2r => ConstraintTrees.Constraint(ConstraintTrees.Linea…
⋮ => ⋮The system solves just like the "symmetric" one:
x2 = optimized_values(c2, objective = c2.objective.value, optimizer = HiGHS.Optimizer)ConstraintTrees.Tree{Float64} with 6 elements:
:coupling => ConstraintTrees.Tree{Float64}(#= 0 elements =#)
:flux_stoichiometry => ConstraintTrees.Tree{Float64}(#= 72 elements =#)
:fluxes => ConstraintTrees.Tree{Float64}(#= 95 elements =#)
:fluxes_forward => ConstraintTrees.Tree{Float64}(#= 95 elements =#)
:fluxes_reverse => ConstraintTrees.Tree{Float64}(#= 95 elements =#)
:objective => 0.873922We can also compare the raw variable counts:
(C.var_count(c), C.var_count(c2))(285, 190)Simplifying the system by removing original variables
If one can assume that the original system is just allocated variables with no other semantics attached, one can reduce the variable and constraint count even in the "nicer" symmetric case from above.
In particular, it is possible to substitute a combination of forward and reverse flux for the bidirectional variables, which allows them to be pruned out of the system together with their original associated bounds:
subst_vals = [C.variable(; idx).value for idx = 1:C.var_count(c)]
c.fluxes = C.zip(c.fluxes, c.fluxes_forward, c.fluxes_reverse) do f, p, n
(var_idx,) = f.value.idxs
subst_value = p.value - n.value
subst_vals[var_idx] = subst_value
C.Constraint(subst_value) # the bidirectional bound is dropped here
end
c = C.prune_variables(C.substitute(c, subst_vals))ConstraintTrees.Tree{ConstraintTrees.Constraint} with 7 elements:
:coupling => ConstraintTrees.Tree{ConstraintTrees.Constraint}…
:directional_flux_balance => 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…The variable count is now reduced, and the system again solves just as the original:
C.var_count(c)190x = optimized_values(c, objective = c.objective.value, optimizer = HiGHS.Optimizer);
x.objective0.8739215069684304The bidirectional flux values are computed transparently in the result:
x.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
⋮ => ⋮This page was generated using Literate.jl.