Patch bound constructors

Project.toml

@@ -1,7 +1,7 @@
 name = "COBREXA"
 uuid = "babc4406-5200-4a30-9033-bf5ae714c842"
 authors = ["The developers of COBREXA.jl"]
-version = "2.5.0"
+version = "2.6.0"
 
 [deps]
 AbstractFBCModels = "5a4f3dfa-1789-40f8-8221-69268c29937c"

docs/src/examples/05b-enzyme-constrained-models.jl

@@ -34,6 +34,7 @@ download_model(
 # -- let's use HiGHS here:
 
 import AbstractFBCModels as A
+import ConstraintTrees as C
 import JSONFBCModels
 import HiGHS
 
@@ -377,18 +378,29 @@ ec_solution.gene_product_capacity
 # arguments as the [`enzyme_constrained_flux_balance_analysis`](@ref), but
 # automatically chooses the "fastest" reaction isozyme for each reaction
 # direction and builds the model with that.
+#
+# We additionally show how to specify more complex bounds for the capacities;
+# in this case we list all fluxes for which we have isozyme data and add a
+# realistic lower limit for the capacity.
+
+minimum_enzyme_capacity = 20.0
 
 simplified_ec_solution = simplified_enzyme_constrained_flux_balance_analysis(
     model;
     reaction_isozymes,
     gene_product_molar_masses = ecoli_core_gene_product_masses,
-    capacity = total_enzyme_capacity,
+    capacity = Dict(
+        :total => (
+            Symbol.(keys(reaction_isozymes)),
+            C.Between(minimum_enzyme_capacity, total_enzyme_capacity),
+        ),
+    ),
     optimizer = HiGHS.Optimizer,
 )
 
 # In this case, the result is the same as with the full analysis:
 
-simplified_ec_solution.capacity_limits.total_capacity
+simplified_ec_solution.capacity_limits.total
 
 # Gene product amounts are not present in the model but are reconstructed
 # nevertheless (they are uniquely determined by the flux):
@@ -425,7 +437,6 @@ ec_optimum = optimized_values(
 )
 
 # ...then creating a system constrained to near-optimal growth:
-import ConstraintTrees as C
 
 ec_system.objective.bound = C.Between(0.99 * ec_optimum, Inf)
 

src/builders/enzymes.jl

@@ -175,10 +175,11 @@ gene product IDs to numbers (such as `Dict{Symbol, Float64}`), and
 All parameter functions may return `nothing`, at which point the given
 object is considered nonexistent and is omitted from constraints.
 
-`capacity_limits` is an interable container of triples `(limit_id,
-gene_product_ids, capacity_bound)` which are converted to a constraint
+`capacity_limits` is an interable container of pairs `limit_id =>
+(gene_product_ids, capacity_bound)` which are converted to a constraint
 identified by the `limit_id` that limits the total mass of `gene_product_ids`
-(which is any iterable container) by `capacity_bound`.
+(which is any iterable container) by `capacity_bound` (which may be anything
+usable as a bound in `Constraint`s).
 """
 function enzyme_constraints(;
     fluxes_forward::C.ConstraintTree,
@@ -220,7 +221,7 @@ function enzyme_constraints(;
                        init = zero(C.LinearValue),
                    ),
                    bound,
-               ) for (id, gps, bound) in capacity_limits
+               ) for (id, (gps, bound)) in capacity_limits
            )
 end
 
@@ -271,7 +272,7 @@ Parameter functions `mass_cost_forward` and `mass_cost_reverse` take a flux ID
 the enzyme mass required to catalyze one "unit" of reaction in the forward or
 reverse direction, respectively. Returning `nothing` ignores the mass cost.
 
-`capacity_limits` is an iterable container of triples `(limit_id, flux_ids,
+`capacity_limits` is an iterable container of pairs `limit_id => (flux_ids,
 bound)`, which creates the capacity bounds over groups of fluxes (in the same
 manner as for gene products in [`enzyme_constraints`](@ref)).
 """
@@ -292,13 +293,13 @@ function simplified_enzyme_constraints(;
         id => C.Constraint(;
             value = C.sum(
                 contribution(fluxes_forward, mass_cost_forward, f) for f in fs;
-                init = zero(C.LinearValue), # TODO not type stable if LinearValueT{not Float64}
+                init = zero(C.LinearValue),
             ) + C.sum(
                 contribution(fluxes_reverse, mass_cost_reverse, f) for f in fs;
-                init = zero(C.LinearValue), # TODO not type stable if LinearValueT{not Float64}
+                init = zero(C.LinearValue),
             ),
             bound,
-        ) for (id, fs, bound) in capacity_limits
+        ) for (id, (fs, bound)) in capacity_limits
     )
 end
 

src/frontend/enzymes.jl

@@ -51,6 +51,78 @@ const Isozyme = IsozymeT{Float64}
 
 export Isozyme
 
+
+"""
+$(TYPEDSIGNATURES)
+
+Expand the `capacity` argument as given to
+[`enzyme_constrained_flux_balance_constraints`](@ref) and
+[`simplified_enzyme_constrained_flux_balance_constraints`](@ref) into a form
+accepted by [`enzyme_constraints`](@ref) and
+[`simplified_enzyme_constraints`](@ref) (respectively).
+
+By default, `Bound`s are kept intact, `Real` values are converted to a fixed
+interval between a zero and the value. All other values are assumed to be lists
+of capacities. (See [`expand_enzyme_capacity_bound`](@ref) for translation of
+actual bounds).
+
+Overloading this function (or [`expand_enzyme_capacity_bound`](@ref)) gives a
+way to simplify the interface of the functions by accomodating custom capacity
+types.
+
+The second argument is provided to this function as a list of full scope of the
+capacities it can work with, by default "all capacities".
+"""
+expand_enzyme_capacity(x, all) =
+    return [:total_capacity => (all, expand_enzyme_capacity_bound(x))]
+
+"""
+$(TYPEDSIGNATURES)
+
+Overload of [`expand_enzyme_capacity`](@ref) for all `Dict`-like iterables.
+"""
+expand_enzyme_capacity(x::Union{Vector{Pair},Dict}, _) =
+    return expand_enzyme_capacity_iterable(x)
+
+"""
+$(TYPEDSIGNATURES)
+
+Overload of [`expand_enzyme_capacity`](@ref) that provides compatibility with
+the earlier capacity specifications (using triples instead of pairs).
+"""
+expand_enzyme_capacity(x::Vector{<:Tuple}, _) =
+    return expand_enzyme_capacity_iterable(id => (grp, cap) for (id, grp, cap) in x)
+
+export expand_enzyme_capacity
+
+"""
+$(TYPEDSIGNATURES)
+
+Internal helper for implementation of [`expand_enzyme_capacity`](@ref) over
+iterable `Dict`-like objects.
+"""
+expand_enzyme_capacity_iterable(x) =
+    return [id => (grp, expand_enzyme_capacity_bound(cap)) for (id, (grp, cap)) in x]
+
+"""
+$(TYPEDSIGNATURES)
+
+Expand a single capacity bound for use in enzyme-constrained models.
+Overloading this function provides additional ways to interpret the capacity
+specifications. Typically used via [`expand_enzyme_capacity`](@ref).
+"""
+expand_enzyme_capacity_bound(x::Real) = return (zero(x), x)
+
+"""
+$(TYPEDSIGNATURES)
+
+By default, [`expand_enzyme_capacity_bound`](@ref) leaves all `Bound`s intact.
+This overload implements this property.
+"""
+expand_enzyme_capacity_bound(x::C.Bound) = return x
+
+export expand_enzyme_capacity_bound
+
 """
 $(TYPEDSIGNATURES)
 
@@ -70,7 +142,9 @@ material is limited by `capacity`.
 `capacity` may be a single number, which sets the mass limit for "all described
 enzymes". Alternatively, `capacity` may be a vector of identifier-genes-limit
 triples that together form a constraint (identified by the given identifier)
-that limits the total sum of the listed genes to the given limit.
+that limits the total sum of the listed genes to the given limit. The
+interpretation of `capacity` is implemented (and can be extended) via
+[`expand_enzyme_capacity`](@ref).
 
 `interface` and `interface_name` are forwarded to
 [`flux_balance_constraints`](@ref).
@@ -79,7 +153,7 @@ function enzyme_constrained_flux_balance_constraints(
     model::A.AbstractFBCModel;
     reaction_isozymes::Dict{String,Dict{String,IsozymeT{R}}},
     gene_product_molar_masses::Dict{String,Float64},
-    capacity::Union{Vector{Tuple{String,Vector{String},R}},R},
+    capacity,
     interface::Maybe{Symbol} = nothing,
     interface_name = :interface,
 ) where {R<:Real}
@@ -137,11 +211,7 @@ function enzyme_constrained_flux_balance_constraints(
         kcat_reverse,
         isozyme_gene_product_stoichiometry,
         gene_product_molar_mass,
-        capacity_limits = capacity isa Real ?
-                          [(:total_capacity, gene_ids, (zero(capacity), capacity))] :
-                          [
-            (Symbol(k), Symbol.(gs), (zero(cap), cap)) for (k, gs, cap) in capacity
-        ],
+        capacity_limits = expand_enzyme_capacity(capacity, gene_ids),
     )
 end
 
@@ -180,15 +250,15 @@ constraints. BMC Bioinformatics 21, 19 (2020).
 https://doi.org/10.1186/s12859-019-3329-9*.
 
 Arguments are as with [`enzyme_constrained_flux_balance_constraints`](@ref),
-with a major difference in `capacity` handling: the identifier lists (2nd
-elements of the triples given in the list) are not identifiers of gene
-products, but identifiers of reactions.
+with a major difference in `capacity` handling: the identifier lists contain
+reactions identifiers (i.e., keys of `fluxes` in the constraint tree), instead
+of the gene product identifiers.
 """
 function simplified_enzyme_constrained_flux_balance_constraints(
     model;
     reaction_isozymes::Dict{String,Dict{String,IsozymeT{R}}},
     gene_product_molar_masses::Dict{String,Float64},
-    capacity::Union{Vector{Tuple{String,Vector{String},R}},R},
+    capacity,
     interface::Maybe{Symbol} = nothing,
     interface_name = :interface,
 ) where {R<:Real}
@@ -250,15 +320,7 @@ function simplified_enzyme_constrained_flux_balance_constraints(
                    maybemap(last, get(min_isozyme_cost_forward, rid, nothing)),
                mass_cost_reverse = rid ->
                    maybemap(last, get(min_isozyme_cost_reverse, rid, nothing)),
-               capacity_limits = capacity isa Real ?
-                                 [(
-                   :total_capacity,
-                   keys(constraints.fluxes),
-                   (zero(capacity), capacity),
-               )] :
-                                 [
-                   (Symbol(k), Symbol.(fs), (zero(cap), cap)) for (k, fs, cap) in capacity
-               ],
+               capacity_limits = expand_enzyme_capacity(capacity, keys(constraints.fluxes)),
            ) *
            :gene_product_amounts^simplified_isozyme_gene_product_amount_constraints(
                (

test/misc.jl

@@ -18,6 +18,8 @@
 # documentation. If you want to add tests here, first consider actually
 # documenting the functionality in `docs/src/examples/` instead.
 
+import ConstraintTrees as C
+
 @testset "Switch bound" begin
     x = Switch(5, 10)
     y = -(((1 + 0.5 * (((1 - x) + 1) * 4)) / 2) - 2.5)
@@ -124,3 +126,21 @@ end
         universal_stoichiometry = stoi,
     )
 end
+
+@testset "Enzyme capacity expansion compat & corner cases" begin
+    x = C.EqualTo(123.0)
+    all = [:ident]
+    y = expand_enzyme_capacity(x, all)
+    @test length(y) == 1
+    (_, (ks, v)) = y[1]
+    # these things should not be touched, thus triple =
+    @test ks === all
+    @test v === x
+
+    x = expand_enzyme_capacity([(:test, [:ident], 123)], [:defa, :ults])
+    @test length(x) == 1
+    (i, (ks, v)) = x[1]
+    @test i == :test
+    @test ks == [:ident]
+    @test v == (0, 123)
+end

test/runtests.jl

@@ -56,10 +56,10 @@ end
 # documentation, which doesn't get erased to improve the test caching.
 
 @testset "COBREXA test suite" begin
-    run_doc_examples()
     @testset "Miscellaneous tests" begin
         run_test_file("misc.jl")
     end
+    run_doc_examples()
     run_test_file("aqua.jl")
 end