diff --git a/Project.toml b/Project.toml index 4f18eb48..e1d25a3e 100644 --- a/Project.toml +++ b/Project.toml @@ -1,7 +1,7 @@ name = "COBREXA" uuid = "babc4406-5200-4a30-9033-bf5ae714c842" authors = ["The developers of COBREXA.jl"] -version = "2.3.1" +version = "2.4.0" [deps] AbstractFBCModels = "5a4f3dfa-1789-40f8-8221-69268c29937c" diff --git a/docs/src/examples/03d-unidirectional.jl b/docs/src/examples/03d-unidirectional.jl new file mode 100644 index 00000000..5ce31b77 --- /dev/null +++ b/docs/src/examples/03d-unidirectional.jl @@ -0,0 +1,151 @@ + +# Copyright (c) 2021-2024, University of Luxembourg #src +# Copyright (c) 2021-2024, Heinrich-Heine University Duesseldorf #src +# #src +# Licensed under the Apache License, Version 2.0 (the "License"); #src +# you may not use this file except in compliance with the License. #src +# You may obtain a copy of the License at #src +# #src +# http://www.apache.org/licenses/LICENSE-2.0 #src +# #src +# Unless required by applicable law or agreed to in writing, software #src +# distributed under the License is distributed on an "AS IS" BASIS, #src +# WITHOUT WARRANTIES OR CONDITIONS OF ANY KIND, either express or implied. #src +# See the License for the specific language governing permissions and #src +# limitations under the License. #src + +# # Split unidirectional reactions in models +# +# By default, the constraint system produced by e.g. +# [`flux_balance_constraints`](@ref) 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`](@ref) and +# [`linear_parsimonious_flux_balance_analysis`](@ref)) 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") + +# We will only work with the constraint representation of the model. The fluxes +# there are bidirectional: +c = flux_balance_constraints(model); +c.fluxes + +# As the simplest approach, we can allocate 2 sets of variables for the forward +# and reverse fluxes via [`sign_split_variables`](@ref) and connect them to the +# fluxes with [`sign_split_constraints`](@ref). 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, + ) + +# 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) + +@test all( #src + isapprox(0, atol = TEST_TOLERANCE), #src + values( #src + C.zip( #src + (f, p, n) -> (p - n - f), #src + x.fluxes, #src + x.fluxes_forward, #src + x.fluxes_reverse, #src + Float64, #src + ), #src + ), #src +) #src + +# ## 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) + +# The system solves just like the "symmetric" one: +x2 = optimized_values(c2, objective = c2.objective.value, optimizer = HiGHS.Optimizer) + +@test isapprox(x.objective, x2.objective, atol = TEST_TOLERANCE) #src + +# We can also compare the raw variable counts: + +(C.var_count(c), C.var_count(c2)) + +@test C.var_count(c) == 285 #src +@test C.var_count(c2) == 190 #src + +# ## 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)) + +# The variable count is now reduced, and the system again solves just as the +# original: + +C.var_count(c) +@test C.var_count(c) == 190 #src + +# + +x = optimized_values(c, objective = c.objective.value, optimizer = HiGHS.Optimizer); +x.objective + +@test isapprox(x.objective, x2.objective, atol = TEST_TOLERANCE) #src + +# The bidirectional flux values are computed transparently in the result: + +x.fluxes diff --git a/src/builders/unsigned.jl b/src/builders/unsigned.jl index 6ff7e666..28833dd1 100644 --- a/src/builders/unsigned.jl +++ b/src/builders/unsigned.jl @@ -63,6 +63,44 @@ export unsigned_negative_contribution_variables """ $(TYPEDSIGNATURES) +A constraint tree that connects negative unsigned variable contributions to +signed ones, while acting as positive contributions. +""" +unsigned_positive_contribution_constraints( + cs::C.ConstraintTree, + negative::C.ConstraintTree, +) = + C.zip(cs, negative) do c, n + C.Constraint( + c.value + n.value, + positive_bound_contribution(something(c.bound, C.Between(-Inf, Inf))), + ) + end + +export unsigned_positive_contribution_constraints + +""" +$(TYPEDSIGNATURES) + +A constraint tree that connects positive unsigned variable contributions to +signed ones, while acting as negative contributions. +""" +unsigned_negative_contribution_constraints( + cs::C.ConstraintTree, + positive::C.ConstraintTree, +) = + C.zip(cs, positive) do c, p + C.Constraint( + p.value - c.value, + positive_bound_contribution(-something(c.bound, C.Between(-Inf, Inf))), + ) + end + +export unsigned_negative_contribution_constraints + +""" +$(TYPEDSIGNATURES) + Shortcut for making a pair of named variable groups created by [`unsigned_positive_contribution_variables`](@ref) and [`unsigned_negative_contribution_variables`](@ref), in subtrees named by diff --git a/src/frontend/enzymes.jl b/src/frontend/enzymes.jl index 62ab5c5f..6fdcbcd8 100644 --- a/src/frontend/enzymes.jl +++ b/src/frontend/enzymes.jl @@ -104,11 +104,13 @@ function enzyme_constrained_flux_balance_constraints( # allocate all variables and build the system constraints = flux_balance_constraints(model; interface, interface_name) - constraints += sign_split_variables( - constraints.fluxes, - positive = :fluxes_forward, - negative = :fluxes_reverse, - ) + constraints += + :fluxes_forward^unsigned_positive_contribution_variables(constraints.fluxes) + constraints *= + :fluxes_reverse^unsigned_negative_contribution_constraints( + constraints.fluxes, + constraints.fluxes_forward, + ) constraints += enzyme_variables(; fluxes_forward = constraints.fluxes_forward, @@ -117,27 +119,21 @@ function enzyme_constrained_flux_balance_constraints( isozyme_reverse_ids, ) - return constraints * - :directional_flux_balance^sign_split_constraints(; - positive = constraints.fluxes_forward, - negative = constraints.fluxes_reverse, - signed = constraints.fluxes, - ) * - enzyme_constraints(; - fluxes_forward = constraints.fluxes_forward, - fluxes_reverse = constraints.fluxes_reverse, - isozyme_forward_amounts = constraints.isozyme_forward_amounts, - isozyme_reverse_amounts = constraints.isozyme_reverse_amounts, - kcat_forward, - kcat_reverse, - isozyme_gene_product_stoichiometry, - gene_product_molar_mass, - capacity_limits = capacity isa Real ? - [(:total_capacity, gene_ids, C.Between(0, capacity))] : - [ - (Symbol(k), Symbol.(gs), C.Between(0, cap)) for (k, gs, cap) in capacity - ], - ) + return constraints * enzyme_constraints(; + fluxes_forward = constraints.fluxes_forward, + fluxes_reverse = constraints.fluxes_reverse, + isozyme_forward_amounts = constraints.isozyme_forward_amounts, + isozyme_reverse_amounts = constraints.isozyme_reverse_amounts, + kcat_forward, + kcat_reverse, + isozyme_gene_product_stoichiometry, + gene_product_molar_mass, + capacity_limits = capacity isa Real ? + [(:total_capacity, gene_ids, C.Between(0, capacity))] : + [ + (Symbol(k), Symbol.(gs), C.Between(0, cap)) for (k, gs, cap) in capacity + ], + ) end export enzyme_constrained_flux_balance_constraints @@ -227,21 +223,17 @@ function simplified_enzyme_constrained_flux_balance_constraints( constraints = flux_balance_constraints(model; interface, interface_name) - # TODO consider only splitting the reactions that we have to split - constraints += sign_split_variables( - constraints.fluxes, - positive = :fluxes_forward, - negative = :fluxes_reverse, - ) + constraints += + :fluxes_forward^unsigned_positive_contribution_variables(constraints.fluxes) + constraints *= + :fluxes_reverse^unsigned_negative_contribution_constraints( + constraints.fluxes, + constraints.fluxes_forward, + ) # connect everything with constraints return constraints * - :directional_flux_balance^sign_split_constraints(; - positive = constraints.fluxes_forward, - negative = constraints.fluxes_reverse, - signed = constraints.fluxes, - ) * :capacity_limits^simplified_enzyme_constraints(; fluxes_forward = constraints.fluxes_forward, fluxes_reverse = constraints.fluxes_reverse, diff --git a/src/frontend/moma.jl b/src/frontend/moma.jl index d5aa89f2..051f5d00 100644 --- a/src/frontend/moma.jl +++ b/src/frontend/moma.jl @@ -130,24 +130,24 @@ function linear_metabolic_adjustment_minimization_constraints( ) constraints = flux_balance_constraints(model) - difference = C.zip(constraints.fluxes, reference_fluxes) do orig, ref + # `difference` actually doesn't need to go to the CT, but we include it + # anyway to calm the curiosity of good neighbors. + constraints *= :reference_diff^C.zip(constraints.fluxes, reference_fluxes) do orig, ref C.Constraint(orig.value - ref) end - difference_split_variables = - C.variables(keys = collect(keys(difference)), bounds = C.Between(0, Inf)) - constraints += :reference_positive_diff^difference_split_variables - constraints += :reference_negative_diff^difference_split_variables + # split the reference_diff into positive and negative contributions + constraints += + :reference_positive_diff^unsigned_positive_contribution_variables( + constraints.reference_diff, + ) + constraints *= + :reference_negative_diff^unsigned_negative_contribution_constraints( + constraints.reference_diff, + constraints.reference_positive_diff, + ) - # `difference` actually doesn't need to go to the CT, but we include it - # anyway to calm the curiosity of good neighbors. constraints * - :reference_diff^difference * - :reference_directional_diff_balance^sign_split_constraints( - positive = constraints.reference_positive_diff, - negative = constraints.reference_negative_diff, - signed = difference, - ) * :linear_minimal_adjustment_objective^C.Constraint( sum_value(constraints.reference_positive_diff, constraints.reference_negative_diff), ) diff --git a/src/frontend/parsimonious.jl b/src/frontend/parsimonious.jl index 6f041673..04016d68 100644 --- a/src/frontend/parsimonious.jl +++ b/src/frontend/parsimonious.jl @@ -73,19 +73,17 @@ Keyword arguments are forwarded to [`flux_balance_constraints`](@ref). """ function linear_parsimonious_flux_balance_constraints(model::A.AbstractFBCModel; kwargs...) constraints = flux_balance_constraints(model; kwargs...) - constraints = - constraints + - :fluxes_forward^unsigned_positive_contribution_variables(constraints.fluxes) + + constraints += :fluxes_reverse^unsigned_negative_contribution_variables(constraints.fluxes) + constraints *= + :fluxes_forward^unsigned_positive_contribution_constraints( + constraints.fluxes, + constraints.fluxes_reverse, + ) return constraints * - :directional_flux_balance^sign_split_constraints( - positive = constraints.fluxes_forward, - negative = constraints.fluxes_reverse, - signed = constraints.fluxes, - ) * :parsimonious_objective^C.Constraint( - sum_value(constraints.fluxes_forward, constraints.fluxes_reverse), - ) + sum_value(constraints.fluxes_forward, constraints.fluxes_reverse), + ) end export linear_parsimonious_flux_balance_constraints