diff --git a/Project.toml b/Project.toml index 8c9a4c95..865849f0 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.5.0" +version = "2.6.0" [deps] AbstractFBCModels = "5a4f3dfa-1789-40f8-8221-69268c29937c" diff --git a/docs/src/examples/05b-enzyme-constrained-models.jl b/docs/src/examples/05b-enzyme-constrained-models.jl index 54baacef..67c83a8b 100644 --- a/docs/src/examples/05b-enzyme-constrained-models.jl +++ b/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) diff --git a/src/builders/enzymes.jl b/src/builders/enzymes.jl index 86f2a8aa..e4b1ce8c 100644 --- a/src/builders/enzymes.jl +++ b/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 diff --git a/src/frontend/enzymes.jl b/src/frontend/enzymes.jl index 0313e211..1025f13b 100644 --- a/src/frontend/enzymes.jl +++ b/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( ( diff --git a/test/misc.jl b/test/misc.jl index 5994f46c..62217c4f 100644 --- a/test/misc.jl +++ b/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 diff --git a/test/runtests.jl b/test/runtests.jl index e7e8cd16..21da5a14 100644 --- a/test/runtests.jl +++ b/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