From e2e70586d1203d9329808ca3b85847f894956121 Mon Sep 17 00:00:00 2001 From: "St. Elmo Wilken" Date: Fri, 2 Feb 2024 20:22:10 +0100 Subject: [PATCH 1/8] implement algo --- src/builders/enzymes.jl | 52 +++++++++++++++++++++++++ src/frontend/enzymes.jl | 86 +++++++++++++++++++++++++++++++++++++++++ 2 files changed, 138 insertions(+) diff --git a/src/builders/enzymes.jl b/src/builders/enzymes.jl index a3cb20d9..81070c43 100644 --- a/src/builders/enzymes.jl +++ b/src/builders/enzymes.jl @@ -239,3 +239,55 @@ enzyme_constraints(; ) export enzyme_constraints + +""" +$(TYPEDSIGNATURES) + +Returns a constraint tree, built from `constraints`, with simplified enzyme +constraints added based on reactions in `fluxes`. + +Inputs: +- `fastest_isozyme_forward` and `fastest_isozyme_reverse` are functions that + take a reaction ID and return the smallest MW/kcat for all isozymes that can + catalyze the reaction. +- `capacity_limits` + +Returns simplified enzyme constraints using the variables. Only reactions in +`reaction_isozymes`, which is a mapping of reaction identifiers to +[`Isozyme`](@ref) descriptions are included in the constraints. +""" +function simplified_enzyme_constraints( + constraints::C.ConstraintTree; + fluxes = constraints.fluxes, + fastest_isozyme_forward = _ -> 0.0, + fastest_isozyme_reverse = _ -> 0.0, + capacity_limits = [], +) + + # allocate variables for everything (nb. += wouldn't associate right here) + constraints = + constraints + + sign_split_variables(fluxes; positive = :fluxes_forward, negative = :fluxes_reverse) + + # connect all parts with constraints + constraints * + :directional_flux_balance^sign_split_constraints( + positive = constraints.fluxes_forward, + negative = constraints.fluxes_reverse, + signed = constraints.fluxes, + ) * + :capacity_bounds^C.ConstraintTree( + Symbol(id) => C.Constraint( + value = sum( + C.value(constraints.fluxes_forward[rid]) * fastest_isozyme_forward(rid) + for rid in flxs + ) + sum( + C.value(constraints.fluxes_reverse[rid]) * fastest_isozyme_reverse(rid) + for rid in flxs + ), + bound = C.Between(0, limit), + ) for (id, flxs, limit) in capacity_limits + ) +end + +export simplified_enzyme_constraints diff --git a/src/frontend/enzymes.jl b/src/frontend/enzymes.jl index e44963f7..e623208f 100644 --- a/src/frontend/enzymes.jl +++ b/src/frontend/enzymes.jl @@ -168,3 +168,89 @@ enzyme_constrained_flux_balance_analysis(model::A.AbstractFBCModel; kwargs...) = ) export enzyme_constrained_flux_balance_analysis + +""" +$(TYPEDSIGNATURES) + +Construct a simplified enzyme-constrained flux-balance constraint system. The +model is parameterized by `reaction_isozymes`, which is a mapping of reaction +identifiers to [`Isozyme`](@ref) descriptions. + +For each gene product in `reaction_isozymes`, a corresponding entry in +`gene_product_molar_masses` must be present, which is a mapping of gene products +to their molar masses. Internally, the cheapest/fastest isozyme (minimum of +MW/kcat for each isozyme) is used in the capacity bound. + +`capacity` may be a single number, which sets the limit of protein required for +all the fluxes in `reaction_isozymes` (mass/mass DW units). Alternatively, +`capacity` may be a vector of identifier-fluxes-limit triples that make a +constraint (identified by the given identifier) that limits the enzyme +requirements of the listed fluxes to the given limit (mass/mass DW units). + +## Note +[`simplified_enzyme_constrained_flux_balance_constraints`](@ref) and +[`enzyme_constrained_flux_balance_constraints`](@ref) differ in how the capacity +bound(s) are formulated. For the former, fluxes are used, but for the latter, +gene products are used directly. +""" +function simplified_enzyme_constrained_flux_balance_constraints( + model; + reaction_isozymes::Dict{String,Dict{String,Isozyme}}, + gene_product_molar_masses::Dict{String,Float64}, + capacity::Union{Vector{Tuple{String,Vector{String},Float64}},Float64}, +) + + # setup: get fastest isozyme for each reaction (flux * MW/kcat = isozyme + # mass fraction required to catalyze reaction) + enzyme_speeds = Dict( + Symbol(rid) => (; + forward = minimum( + sum( + stoich * gene_product_molar_masses[gid] for + (gid, stoich) in iso.gene_product_stoichiometry + ) / iso.kcat_forward for iso in values(isos) + ), + reverse = minimum( + sum( + stoich * gene_product_molar_masses[gid] for + (gid, stoich) in iso.gene_product_stoichiometry + ) / iso.kcat_reverse for iso in values(isos) + ), + ) for (rid, isos) in reaction_isozymes + ) + + # fails ungracefully if rid is not in enzyme_speeds - prevents footguns + fastest_isozyme_forward(rid) = enzyme_speeds[rid].forward + fastest_isozyme_reverse(rid) = enzyme_speeds[rid].reverse + + capacity_limits = + capacity isa Real ? [("totalcapacity", keys(enzyme_speeds), capacity)] : capacity + + c = simplified_enzyme_constraints( + flux_balance_constraints(model); + fastest_isozyme_forward, + fastest_isozyme_reverse, + capacity_limits, + ) +end + +export simplified_enzyme_constrained_flux_balance_constraints + +""" +$(TYPEDSIGNATURES) + +Perform the enzyme-constrained flux balance analysis on the `model` and return the solved constraint system. + +Arguments are forwarded to +[`simplified_enzyme_constrained_flux_balance_constraints`](@ref); solver configuration +arguments are forwarded to [`optimized_values`](@ref). +""" +simplified_enzyme_constrained_flux_balance_analysis(model::A.AbstractFBCModel; kwargs...) = + frontend_optimized_values( + simplified_enzyme_constrained_flux_balance_constraints, + model; + objective = x -> x.objective.value, + kwargs..., + ) + +export simplified_enzyme_constrained_flux_balance_analysis From 88a3b37c138a4dda3ddc298fe9f12b49c66a45c8 Mon Sep 17 00:00:00 2001 From: "St. Elmo Wilken" Date: Fri, 2 Feb 2024 20:22:27 +0100 Subject: [PATCH 2/8] add smoment example --- .../examples/05-enzyme-constrained-models.jl | 19 +++++++++++++ src/builders/enzymes.jl | 7 ++--- src/frontend/enzymes.jl | 27 +++++++++++-------- 3 files changed, 37 insertions(+), 16 deletions(-) diff --git a/docs/src/examples/05-enzyme-constrained-models.jl b/docs/src/examples/05-enzyme-constrained-models.jl index 58c8cce9..5e8aa873 100644 --- a/docs/src/examples/05-enzyme-constrained-models.jl +++ b/docs/src/examples/05-enzyme-constrained-models.jl @@ -342,3 +342,22 @@ ec_solution = enzyme_constrained_flux_balance_analysis( 0.011875920383431717, #src atol = TEST_TOLERANCE, #src ) #src + +# ## Running a simplified enzyme constrained model +# We can also simplify the normal enzyme constrained model, by only considering +# the fastest isozyme, hence obviating the need to create variables for each +# isozyme, and the actual gene products. + +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, + optimizer = GLPK.Optimizer, +) + +@test isapprox( + ec_solution.objective, + simplified_ec_solution.objective, + atol = TEST_TOLERANCE, +) #src diff --git a/src/builders/enzymes.jl b/src/builders/enzymes.jl index 81070c43..e077842e 100644 --- a/src/builders/enzymes.jl +++ b/src/builders/enzymes.jl @@ -250,11 +250,8 @@ Inputs: - `fastest_isozyme_forward` and `fastest_isozyme_reverse` are functions that take a reaction ID and return the smallest MW/kcat for all isozymes that can catalyze the reaction. -- `capacity_limits` - -Returns simplified enzyme constraints using the variables. Only reactions in -`reaction_isozymes`, which is a mapping of reaction identifiers to -[`Isozyme`](@ref) descriptions are included in the constraints. +- `capacity_limits` is an iterable container of triples `(capacity_id, + reaction_ids, capacity_bound)`, which creates capacity bound(s). """ function simplified_enzyme_constraints( constraints::C.ConstraintTree; diff --git a/src/frontend/enzymes.jl b/src/frontend/enzymes.jl index e623208f..b5c9510c 100644 --- a/src/frontend/enzymes.jl +++ b/src/frontend/enzymes.jl @@ -153,7 +153,8 @@ export enzyme_constrained_flux_balance_constraints """ $(TYPEDSIGNATURES) -Perform the enzyme-constrained flux balance analysis on the `model` and return the solved constraint system. +Perform the enzyme-constrained flux balance analysis on the `model` and return +the solved constraint system. Arguments are forwarded to [`enzyme_constrained_flux_balance_constraints`](@ref); solver configuration @@ -172,14 +173,17 @@ export enzyme_constrained_flux_balance_analysis """ $(TYPEDSIGNATURES) -Construct a simplified enzyme-constrained flux-balance constraint system. The -model is parameterized by `reaction_isozymes`, which is a mapping of reaction -identifiers to [`Isozyme`](@ref) descriptions. +Construct a simplified enzyme-constrained flux-balance constraint system. Based +on the SMOMENT algorithm described in *Bekiaris, P.S., Klamt, S. Automatic +construction of metabolic models with enzyme constraints. BMC Bioinformatics 21, +19 (2020). https://doi.org/10.1186/s12859-019-3329-9*. -For each gene product in `reaction_isozymes`, a corresponding entry in -`gene_product_molar_masses` must be present, which is a mapping of gene products -to their molar masses. Internally, the cheapest/fastest isozyme (minimum of -MW/kcat for each isozyme) is used in the capacity bound. +The model is parameterized by `reaction_isozymes`, which is a mapping of +reaction identifiers to [`Isozyme`](@ref) descriptions. For each gene product in +`reaction_isozymes`, a corresponding entry in `gene_product_molar_masses` must +be present, which is a mapping of gene products to their molar masses. +Internally, the cheapest/fastest isozyme (minimum of MW/kcat for each isozyme) +is used in the capacity bound. `capacity` may be a single number, which sets the limit of protein required for all the fluxes in `reaction_isozymes` (mass/mass DW units). Alternatively, @@ -239,11 +243,12 @@ export simplified_enzyme_constrained_flux_balance_constraints """ $(TYPEDSIGNATURES) -Perform the enzyme-constrained flux balance analysis on the `model` and return the solved constraint system. +Perform the enzyme-constrained flux balance analysis on the `model` and return +the solved constraint system. Arguments are forwarded to -[`simplified_enzyme_constrained_flux_balance_constraints`](@ref); solver configuration -arguments are forwarded to [`optimized_values`](@ref). +[`simplified_enzyme_constrained_flux_balance_constraints`](@ref); solver +configuration arguments are forwarded to [`optimized_values`](@ref). """ simplified_enzyme_constrained_flux_balance_analysis(model::A.AbstractFBCModel; kwargs...) = frontend_optimized_values( From b0aeab0852c924f01aaa63bcbdbe1494afc92ef2 Mon Sep 17 00:00:00 2001 From: stelmo Date: Fri, 2 Feb 2024 19:33:18 +0000 Subject: [PATCH 3/8] automatic formatting triggered by @stelmo on PR #12 --- src/builders/enzymes.jl | 2 +- src/frontend/enzymes.jl | 4 ++-- 2 files changed, 3 insertions(+), 3 deletions(-) diff --git a/src/builders/enzymes.jl b/src/builders/enzymes.jl index e077842e..76730efd 100644 --- a/src/builders/enzymes.jl +++ b/src/builders/enzymes.jl @@ -272,7 +272,7 @@ function simplified_enzyme_constraints( positive = constraints.fluxes_forward, negative = constraints.fluxes_reverse, signed = constraints.fluxes, - ) * + ) * :capacity_bounds^C.ConstraintTree( Symbol(id) => C.Constraint( value = sum( diff --git a/src/frontend/enzymes.jl b/src/frontend/enzymes.jl index b5c9510c..a334f5d3 100644 --- a/src/frontend/enzymes.jl +++ b/src/frontend/enzymes.jl @@ -183,13 +183,13 @@ reaction identifiers to [`Isozyme`](@ref) descriptions. For each gene product in `reaction_isozymes`, a corresponding entry in `gene_product_molar_masses` must be present, which is a mapping of gene products to their molar masses. Internally, the cheapest/fastest isozyme (minimum of MW/kcat for each isozyme) -is used in the capacity bound. +is used in the capacity bound. `capacity` may be a single number, which sets the limit of protein required for all the fluxes in `reaction_isozymes` (mass/mass DW units). Alternatively, `capacity` may be a vector of identifier-fluxes-limit triples that make a constraint (identified by the given identifier) that limits the enzyme -requirements of the listed fluxes to the given limit (mass/mass DW units). +requirements of the listed fluxes to the given limit (mass/mass DW units). ## Note [`simplified_enzyme_constrained_flux_balance_constraints`](@ref) and From 03ddeaeb5863d1eefaf9a47b6bb0dcc8e418e265 Mon Sep 17 00:00:00 2001 From: "St. Elmo Wilken" Date: Fri, 2 Feb 2024 20:37:31 +0100 Subject: [PATCH 4/8] make docs render --- docs/src/examples/05-enzyme-constrained-models.jl | 8 ++++---- 1 file changed, 4 insertions(+), 4 deletions(-) diff --git a/docs/src/examples/05-enzyme-constrained-models.jl b/docs/src/examples/05-enzyme-constrained-models.jl index 5e8aa873..866eab1e 100644 --- a/docs/src/examples/05-enzyme-constrained-models.jl +++ b/docs/src/examples/05-enzyme-constrained-models.jl @@ -356,8 +356,8 @@ simplified_ec_solution = simplified_enzyme_constrained_flux_balance_analysis( optimizer = GLPK.Optimizer, ) -@test isapprox( - ec_solution.objective, - simplified_ec_solution.objective, - atol = TEST_TOLERANCE, +@test isapprox( #src + ec_solution.objective, #src + simplified_ec_solution.objective, #src + atol = TEST_TOLERANCE, #src ) #src From 06021f1d22326fdd3e1c0c57627b5548afcb2336 Mon Sep 17 00:00:00 2001 From: "St. Elmo Wilken" Date: Fri, 2 Feb 2024 21:38:28 +0100 Subject: [PATCH 5/8] add interface kwargs to frontend --- src/frontend/enzymes.jl | 8 ++++++-- 1 file changed, 6 insertions(+), 2 deletions(-) diff --git a/src/frontend/enzymes.jl b/src/frontend/enzymes.jl index a334f5d3..4b0f5cc4 100644 --- a/src/frontend/enzymes.jl +++ b/src/frontend/enzymes.jl @@ -191,6 +191,8 @@ all the fluxes in `reaction_isozymes` (mass/mass DW units). Alternatively, constraint (identified by the given identifier) that limits the enzyme requirements of the listed fluxes to the given limit (mass/mass DW units). +`interface` and `interface_name` are forwarded to [`flux_balance_constraints`](@ref). + ## Note [`simplified_enzyme_constrained_flux_balance_constraints`](@ref) and [`enzyme_constrained_flux_balance_constraints`](@ref) differ in how the capacity @@ -202,6 +204,8 @@ function simplified_enzyme_constrained_flux_balance_constraints( reaction_isozymes::Dict{String,Dict{String,Isozyme}}, gene_product_molar_masses::Dict{String,Float64}, capacity::Union{Vector{Tuple{String,Vector{String},Float64}},Float64}, + interface::Maybe{Symbol} = nothing, + interface_name = :interface, ) # setup: get fastest isozyme for each reaction (flux * MW/kcat = isozyme @@ -230,8 +234,8 @@ function simplified_enzyme_constrained_flux_balance_constraints( capacity_limits = capacity isa Real ? [("totalcapacity", keys(enzyme_speeds), capacity)] : capacity - c = simplified_enzyme_constraints( - flux_balance_constraints(model); + simplified_enzyme_constraints( + flux_balance_constraints(model; interface, interface_name); fastest_isozyme_forward, fastest_isozyme_reverse, capacity_limits, From a4b6ccee54e482cd2a887a156a603b36da51587c Mon Sep 17 00:00:00 2001 From: Mirek Kratochvil Date: Wed, 6 Mar 2024 21:19:23 +0100 Subject: [PATCH 6/8] working version, debug needed --- .../examples/05-enzyme-constrained-models.jl | 38 +++- src/builders/enzymes.jl | 184 ++++++++++-------- src/frontend/enzymes.jl | 133 +++++++------ 3 files changed, 208 insertions(+), 147 deletions(-) diff --git a/docs/src/examples/05-enzyme-constrained-models.jl b/docs/src/examples/05-enzyme-constrained-models.jl index 866eab1e..5be13186 100644 --- a/docs/src/examples/05-enzyme-constrained-models.jl +++ b/docs/src/examples/05-enzyme-constrained-models.jl @@ -329,6 +329,21 @@ ec_solution = enzyme_constrained_flux_balance_analysis( optimizer = GLPK.Optimizer, ) +# We can notice that the objective function is a little lower than with +# unconstrained E. coli core: + +ec_solution.objective + +# One can also observe many interesting thing, e.g. the amount of gene product +# material required for the system to run: + +ec_solution.gene_product_amount + +# The total amount of required gene product mass is, by default, present as +# `total_capacity`: + +ec_solution.gene_product_capacity + #src these values should be unique (glucose transporter is the only way to get carbon into the system) @test isapprox(ec_solution.objective, 0.706993382849705, atol = TEST_TOLERANCE) #src @test isapprox( #src @@ -343,10 +358,16 @@ ec_solution = enzyme_constrained_flux_balance_analysis( atol = TEST_TOLERANCE, #src ) #src -# ## Running a simplified enzyme constrained model -# We can also simplify the normal enzyme constrained model, by only considering -# the fastest isozyme, hence obviating the need to create variables for each -# isozyme, and the actual gene products. +# ## Simplified models +# +# Because most reactions typically only use a single isozyme, we may also use a +# simplified representation of the problem where this fact is reflected, saving +# the variable allocation for the isozymes. +# +# [`simplified_enzyme_constrained_flux_balance_analysis`](@ref) takes similar +# 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. simplified_ec_solution = simplified_enzyme_constrained_flux_balance_analysis( model; @@ -356,6 +377,15 @@ simplified_ec_solution = simplified_enzyme_constrained_flux_balance_analysis( optimizer = GLPK.Optimizer, ) +# In this case, the result is the same as with the full analysis: + +simplified_ec_solution.capacity_limits.total_capacity + +# Gene product amounts are not present in the model but are reconstructed +# nevertheless (they are uniquely determined by the flux: + +simplified_ec_solution.gene_product_amounts + @test isapprox( #src ec_solution.objective, #src simplified_ec_solution.objective, #src diff --git a/src/builders/enzymes.jl b/src/builders/enzymes.jl index 76730efd..858261af 100644 --- a/src/builders/enzymes.jl +++ b/src/builders/enzymes.jl @@ -69,9 +69,9 @@ export isozyme_flux_constraints """ $(TYPEDSIGNATURES) -A constraint tree that binds the isozyme amounts to gene product amounts -accordingly to their multiplicities (aka. stoichiometries, protein units, ...) -given by `isozyme_stoichiometry`. +A constraint tree that computes the gene product amounts from given isozyme +amounts their multiplicities (aka. stoichiometries, protein units, ...) given +by `isozyme_stoichiometry`. Values in ConstraintTree `gene_product_amounts` should describe the gene product allocations. Allocation for the isozyme is ignored if the gene product @@ -85,17 +85,13 @@ and reverse-direction amounts at once. To only use a single tree, use an uni-tuple: `isozyme_amounts = tuple(my_tree)`. Parameter function `isozyme_stoichiometry` gets called with a reaction and -isozyme ID as given by the isozyme amount trees. It should return `nothing` in +isozyme IDs as given by the isozyme amount trees. It should return `nothing` in case there's no information -- in such case, the isozyme is *not going to be included* in the calculation of gene product mass. """ -function gene_product_isozyme_constraints( - gene_product_amounts::C.ConstraintTree, - isozyme_amounts, - isozyme_stoichiometry, -) +function isozyme_gene_product_amount_constraints(isozyme_amounts, isozyme_stoichiometry) res = C.ConstraintTree() - # This needs to invert the stoichiometry mapping, + # This needs to "invert" the stoichiometry mapping, # so we patch up a fresh new CT in place. for iss in isozyme_amounts for (rid, is) in iss @@ -103,12 +99,10 @@ function gene_product_isozyme_constraints( gpstoi = isozyme_stoichiometry(rid, iid) isnothing(gpstoi) && continue for (gp, stoi) in gpstoi - haskey(gene_product_amounts, gp) || continue if haskey(res, gp) res[gp].value += i.value * stoi else - res[gp] = - C.Constraint(i.value * stoi - gene_product_amounts[gp].value, 0) + res[gp] = C.Constraint(i.value * stoi, 0) end end end @@ -117,7 +111,7 @@ function gene_product_isozyme_constraints( res end -export gene_product_isozyme_constraints +export isozyme_gene_product_amount_constraints """ $(TYPEDSIGNATURES) @@ -126,8 +120,6 @@ Returns a constraint tree with enzyme capacity constraints, added for reactions in `fluxes_forward` and `fluxes_reverse`. This is used to construct the constraint system in [`enzyme_constrained_flux_balance_constraints`](@ref). -Parameter `gene_ids` is an iterable collection of gene identifiers (as `Symbol`s). - Parameter function `isozyme_ids` takes a reaction ID and returns `nothing` if the reaction does not have isozymes associated with it, or an iterable container of all the isozyme IDs for that reaction (as `Symbol`s). @@ -136,10 +128,6 @@ Parameters `isozyme_forward_ids` and `isozyme_reverse_ids` can be used to fine-tune the generated isozymes in either direction; both default to `isozyme_ids`. -Parameter function `gene_product_bound` can be used to generate a bound on the -gene product amount variables (the `Symbol` key is passed as a parameter). By -default, all amounts are bounded to be non-negative. - The keys in the output constraint tree can be customized by setting `isozyme_forward_amounts_name`, `isozyme_reverse_amounts_name` and `gene_product_amounts_name`. @@ -147,14 +135,11 @@ The keys in the output constraint tree can be customized by setting enzyme_variables(; fluxes_forward::C.ConstraintTree, fluxes_reverse::C.ConstraintTree, - gene_ids = Symbol[], isozyme_ids = _ -> nothing, isozyme_forward_ids = isozyme_ids, isozyme_reverse_ids = isozyme_ids, - gene_product_bound = _ -> C.Between(0, Inf), isozyme_forward_amounts_name = :isozyme_forward_amounts, isozyme_reverse_amounts_name = :isozyme_reverse_amounts, - gene_product_amounts_name = :gene_product_amounts, ) = isozyme_forward_amounts_name^isozyme_amount_variables( Symbol.(keys(fluxes_forward)), @@ -163,10 +148,6 @@ enzyme_variables(; isozyme_reverse_amounts_name^isozyme_amount_variables( Symbol.(keys(fluxes_reverse)), isozyme_reverse_ids, - ) + - gene_product_amounts_name^C.variables( - keys = gene_ids, - bounds = gene_product_bound.(gene_ids), ) export enzyme_variables @@ -196,12 +177,11 @@ 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`. """ -enzyme_constraints(; +function enzyme_constraints(; fluxes_forward::C.ConstraintTree, fluxes_reverse::C.ConstraintTree, isozyme_forward_amounts::C.ConstraintTree, isozyme_reverse_amounts::C.ConstraintTree, - gene_product_amounts::C.ConstraintTree, kcat_forward = (_, _) -> nothing, kcat_reverse = (_, _) -> nothing, isozyme_gene_product_stoichiometry = (_, _) -> nothing, @@ -209,81 +189,113 @@ enzyme_constraints(; capacity_limits = [], isozyme_flux_forward_balance_name = :isozyme_flux_forward_balance, isozyme_flux_reverse_balance_name = :isozyme_flux_reverse_balance, - gene_product_isozyme_balance_name = :gene_product_isozyme_balance, + gene_product_amounts_name = :gene_product_amounts, gene_product_capacity_name = :gene_product_capacity, -) = - isozyme_flux_forward_balance_name^isozyme_flux_constraints( - isozyme_forward_amounts, - fluxes_forward, - kcat_forward, - ) * - isozyme_flux_reverse_balance_name^isozyme_flux_constraints( - isozyme_reverse_amounts, - fluxes_reverse, - kcat_reverse, - ) * - gene_product_isozyme_balance_name^gene_product_isozyme_constraints( - gene_product_amounts, +) + gene_product_amounts = isozyme_gene_product_amount_constraints( (isozyme_forward_amounts, isozyme_reverse_amounts), isozyme_gene_product_stoichiometry, - ) * - gene_product_capacity_name^C.ConstraintTree( - Symbol(id) => C.Constraint(; - value = sum( - gene_product_amounts[gp].value * gpmm for - (gp, gpmm) in ((gp, gene_product_molar_mass(gp)) for gp in gps) if - !isnothing(gpmm) - ), - bound, - ) for (id, gps, bound) in capacity_limits ) + return isozyme_flux_forward_balance_name^isozyme_flux_constraints( + isozyme_forward_amounts, + fluxes_forward, + kcat_forward, + ) * + isozyme_flux_reverse_balance_name^isozyme_flux_constraints( + isozyme_reverse_amounts, + fluxes_reverse, + kcat_reverse, + ) * + gene_product_amounts_name^gene_product_amounts * + gene_product_capacity_name^C.ConstraintTree( + id => C.Constraint(; + value = sum( + gene_product_amounts[gp].value * gpmm for + (gp, gpmm) in ((gp, gene_product_molar_mass(gp)) for gp in gps) if + !isnothing(gpmm) && haskey(gene_product_amounts, gp) + ), + bound, + ) for (id, gps, bound) in capacity_limits + ) +end + export enzyme_constraints """ $(TYPEDSIGNATURES) -Returns a constraint tree, built from `constraints`, with simplified enzyme -constraints added based on reactions in `fluxes`. +Like [`isozyme_gene_product_amount_constraints`](@ref), but works with the +"simplified" view where each reaction has an uniquely determined catalytic +isozyme, as with [`simplified_enzyme_constraints`](@ref). -Inputs: -- `fastest_isozyme_forward` and `fastest_isozyme_reverse` are functions that - take a reaction ID and return the smallest MW/kcat for all isozymes that can - catalyze the reaction. -- `capacity_limits` is an iterable container of triples `(capacity_id, - reaction_ids, capacity_bound)`, which creates capacity bound(s). +As the main difference, `isozyme_amounts` is not supposed to contain a subtree +of isozymes for each reaction, and `isozyme_stoichiometry` accordingly expects +only a single identifier parameter for the reaction ID. """ -function simplified_enzyme_constraints( - constraints::C.ConstraintTree; - fluxes = constraints.fluxes, - fastest_isozyme_forward = _ -> 0.0, - fastest_isozyme_reverse = _ -> 0.0, +function simplified_isozyme_gene_product_amount_constraints( + isozyme_amounts_stoichiometry_kcat, +) + res = C.ConstraintTree() + for (iss, stoikcatfn) in isozyme_amounts + for (rid, i) in iss + x = stoikcatfn(rid) + isnothing(x) && continue + gpstoi, kcat = x + isnothing(kcat) && continue + for (gp, stoi) in gpstoi + if haskey(res, gp) + res[gp].value += i.value * stoi / kcat + else + res[gp] = C.Constraint(i.value * stoi / kcat, 0) + end + end + end + end + res +end + +export simplified_isozyme_gene_product_amount_constraints + +""" +$(TYPEDSIGNATURES) + +Build a constraint system that bounds fluxes according to their enzyme mass +requirements, with respect to per-reaction enzyme mass costs. + +Parameter functions `mass_cost_forward` and `mass_cost_reverse` take a flux ID +(corresponding to a flux in `fluxes_forward` and `fluxes_reverse`) and return +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, +bound)`, which creates the capacity bounds over groups of fluxes (in the same +manner as for gene products in [`enzyme_constraints`](@ref)). +""" +function simplified_enzyme_constraints(; + fluxes_forward::C.ConstraintTree, + fluxes_reverse::C.ConstraintTree, + mass_cost_forward = _ -> nothing, + mass_cost_reverse = _ -> nothing, capacity_limits = [], ) + function contribution(fl, cost, rid) + c = cost(rid) + isnothing(c) && !haskey(fl, rid) && return zero(C.LinearValue) + return c * fl[rid] + end - # allocate variables for everything (nb. += wouldn't associate right here) - constraints = - constraints + - sign_split_variables(fluxes; positive = :fluxes_forward, negative = :fluxes_reverse) - - # connect all parts with constraints - constraints * - :directional_flux_balance^sign_split_constraints( - positive = constraints.fluxes_forward, - negative = constraints.fluxes_reverse, - signed = constraints.fluxes, - ) * - :capacity_bounds^C.ConstraintTree( - Symbol(id) => C.Constraint( + return C.ConstraintTree( + id => C.Constraint(; value = sum( - C.value(constraints.fluxes_forward[rid]) * fastest_isozyme_forward(rid) - for rid in flxs + contribution(fluxes_forward, mass_cost_forward, i) for i in fluxes; + init = zero(C.LinearValue), ) + sum( - C.value(constraints.fluxes_reverse[rid]) * fastest_isozyme_reverse(rid) - for rid in flxs + contribution(fluxes_reverse, mass_cost_reverse, i) for i in fluxes; + init = zero(C.LinearValue), ), - bound = C.Between(0, limit), - ) for (id, flxs, limit) in capacity_limits + bound = bound, + ) for (id, fluxes, bound) in capacity_limits ) end diff --git a/src/frontend/enzymes.jl b/src/frontend/enzymes.jl index 4b0f5cc4..4646c893 100644 --- a/src/frontend/enzymes.jl +++ b/src/frontend/enzymes.jl @@ -105,7 +105,6 @@ function enzyme_constrained_flux_balance_constraints( gene_product_molar_masses[k] else throw(DomainError(k, "missing a required gene product molar mass")) - # TODO: maybe better default to 0 and ignore the issue ? end end @@ -121,7 +120,6 @@ function enzyme_constrained_flux_balance_constraints( constraints += enzyme_variables(; fluxes_forward = constraints.fluxes_forward, fluxes_reverse = constraints.fluxes_reverse, - gene_ids, isozyme_forward_ids, isozyme_reverse_ids, ) @@ -137,14 +135,15 @@ function enzyme_constrained_flux_balance_constraints( fluxes_reverse = constraints.fluxes_reverse, isozyme_forward_amounts = constraints.isozyme_forward_amounts, isozyme_reverse_amounts = constraints.isozyme_reverse_amounts, - gene_product_amounts = constraints.gene_product_amounts, kcat_forward, kcat_reverse, isozyme_gene_product_stoichiometry, gene_product_molar_mass, capacity_limits = capacity isa Real ? - [(:total_capacity, gene_ids, capacity)] : - [(k, Symbol.(gs), cap) for (k, gs, cap) in capacity], + [(:total_capacity, gene_ids, C.Between(0, capacity))] : + [ + (Symbol(k), Symbol.(gs), C.Between(0, cap)) for (k, gs, cap) in capacity + ], ) end @@ -173,31 +172,16 @@ export enzyme_constrained_flux_balance_analysis """ $(TYPEDSIGNATURES) -Construct a simplified enzyme-constrained flux-balance constraint system. Based -on the SMOMENT algorithm described in *Bekiaris, P.S., Klamt, S. Automatic -construction of metabolic models with enzyme constraints. BMC Bioinformatics 21, -19 (2020). https://doi.org/10.1186/s12859-019-3329-9*. - -The model is parameterized by `reaction_isozymes`, which is a mapping of -reaction identifiers to [`Isozyme`](@ref) descriptions. For each gene product in -`reaction_isozymes`, a corresponding entry in `gene_product_molar_masses` must -be present, which is a mapping of gene products to their molar masses. -Internally, the cheapest/fastest isozyme (minimum of MW/kcat for each isozyme) -is used in the capacity bound. - -`capacity` may be a single number, which sets the limit of protein required for -all the fluxes in `reaction_isozymes` (mass/mass DW units). Alternatively, -`capacity` may be a vector of identifier-fluxes-limit triples that make a -constraint (identified by the given identifier) that limits the enzyme -requirements of the listed fluxes to the given limit (mass/mass DW units). - -`interface` and `interface_name` are forwarded to [`flux_balance_constraints`](@ref). - -## Note -[`simplified_enzyme_constrained_flux_balance_constraints`](@ref) and -[`enzyme_constrained_flux_balance_constraints`](@ref) differ in how the capacity -bound(s) are formulated. For the former, fluxes are used, but for the latter, -gene products are used directly. +Like [`enzyme_constrained_flux_balance_constraints`](@ref), but automatically +selects a single "fastest" isozyme for each reaction direction. In turn, the +system requires much less variables in the constraint system description, and +usually solves more efficiently, for the price of possibly finding suboptimal +solutions. The method follows the SMOMENT algorithm described in *Bekiaris, +P.S., Klamt, S. Automatic construction of metabolic models with enzyme +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). """ function simplified_enzyme_constrained_flux_balance_constraints( model; @@ -208,38 +192,73 @@ function simplified_enzyme_constrained_flux_balance_constraints( interface_name = :interface, ) - # setup: get fastest isozyme for each reaction (flux * MW/kcat = isozyme - # mass fraction required to catalyze reaction) - enzyme_speeds = Dict( - Symbol(rid) => (; - forward = minimum( - sum( - stoich * gene_product_molar_masses[gid] for - (gid, stoich) in iso.gene_product_stoichiometry - ) / iso.kcat_forward for iso in values(isos) - ), - reverse = minimum( - sum( - stoich * gene_product_molar_masses[gid] for - (gid, stoich) in iso.gene_product_stoichiometry - ) / iso.kcat_reverse for iso in values(isos) - ), - ) for (rid, isos) in reaction_isozymes + isozyme_mass(i) = sum( + ( + ( + haskey(gene_product_molar_masses, gid) ? + coeff * gene_product_molar_masses[gid] : + throw(DomainError(gid, "missing a required gene product molar mass")) + ) for (gid, stoich) in i.gene_product_stoichiometry + ), + init = 0.0, ) - # fails ungracefully if rid is not in enzyme_speeds - prevents footguns - fastest_isozyme_forward(rid) = enzyme_speeds[rid].forward - fastest_isozyme_reverse(rid) = enzyme_speeds[rid].reverse + min_isozyme_cost_forward = Dict( + Symbol(rid) => argmin(last, (i, isozyme_mass(i) / i.kcat_forward) for i in is) + for (rid, is) in reaction_isozymes if !isnothing(is.kcat_forward) + ) + min_isozyme_cost_reverse = Dict( + Symbol(rid) => argmin(last, (i, isozyme_mass(i) / i.kcat_reverse) for i in is) + for (rid, is) in reaction_isozymes if !isnothing(is.kcat_reverse) + ) - capacity_limits = - capacity isa Real ? [("totalcapacity", keys(enzyme_speeds), capacity)] : capacity + constraints = flux_balance_constraints(model; interface, interface_name) - simplified_enzyme_constraints( - flux_balance_constraints(model; interface, interface_name); - fastest_isozyme_forward, - fastest_isozyme_reverse, - capacity_limits, + constraints += sign_split_variables( + constraints.fluxes, + positive = :fluxes_forward, + negative = :fluxes_reverse, ) + + return constraints * + sign_split_constraints(; + positive = constraints.fluxes_forward, + negative = constraints.fluxes_reverse, + signed = constraints.fluxes, + ) * + :capacity_limits^simplified_enzyme_constraints( + constraints.fluxes_forward, + constraints.fluxes_reverse, + mass_cost_forward = rid -> + maybemap(last, get(min_cost_forward, rid, nothing)), + mass_cost_reverse = rid -> + maybemap(last, get(min_cost_reverse, rid, nothing)), + capacity_limits = capacity isa Real ? + [( + :total_capacity, + Symbol.(a.genes(model)), + C.Between(0, capacity), + )] : + [ + (Symbol(k), Symbol.(gs), C.Between(0, cap)) for (k, gs, cap) in capacity + ], + ) * + :gene_product_amounts^simplified_isozyme_gene_product_amount_constraints( + ( + constraints.fluxes_forward, + rid -> maybemap( + i -> (i.gene_product_stoichiometry, i.kcat_forward), + get(min_isozyme_cost_forward, rid, nothing), + ), + ), + ( + constraints.fluxes_reverse, + rid -> maybemap( + i -> (i.gene_product_stoichiometry, i.kcat_reverse), + get(min_isozyme_cost_reverse, rid, nothing), + ), + ), + ) end export simplified_enzyme_constrained_flux_balance_constraints From 67617e4d69c1a7a142db1daa179e2325cfbf4f7c Mon Sep 17 00:00:00 2001 From: Mirek Kratochvil Date: Wed, 6 Mar 2024 23:09:01 +0100 Subject: [PATCH 7/8] fix quite a few errors, simplify the gene product view in the model --- .../examples/05-enzyme-constrained-models.jl | 2 +- src/builders/enzymes.jl | 27 +++++----- src/frontend/enzymes.jl | 54 +++++++++++++------ 3 files changed, 52 insertions(+), 31 deletions(-) diff --git a/docs/src/examples/05-enzyme-constrained-models.jl b/docs/src/examples/05-enzyme-constrained-models.jl index 5be13186..4d6e398e 100644 --- a/docs/src/examples/05-enzyme-constrained-models.jl +++ b/docs/src/examples/05-enzyme-constrained-models.jl @@ -337,7 +337,7 @@ ec_solution.objective # One can also observe many interesting thing, e.g. the amount of gene product # material required for the system to run: -ec_solution.gene_product_amount +ec_solution.gene_product_amounts # The total amount of required gene product mass is, by default, present as # `total_capacity`: diff --git a/src/builders/enzymes.jl b/src/builders/enzymes.jl index 858261af..ca521479 100644 --- a/src/builders/enzymes.jl +++ b/src/builders/enzymes.jl @@ -102,7 +102,7 @@ function isozyme_gene_product_amount_constraints(isozyme_amounts, isozyme_stoich if haskey(res, gp) res[gp].value += i.value * stoi else - res[gp] = C.Constraint(i.value * stoi, 0) + res[gp] = C.Constraint(i.value * stoi) end end end @@ -229,25 +229,24 @@ Like [`isozyme_gene_product_amount_constraints`](@ref), but works with the "simplified" view where each reaction has an uniquely determined catalytic isozyme, as with [`simplified_enzyme_constraints`](@ref). -As the main difference, `isozyme_amounts` is not supposed to contain a subtree -of isozymes for each reaction, and `isozyme_stoichiometry` accordingly expects -only a single identifier parameter for the reaction ID. +As the main difference, the arguments are tuples that contain first the +constraint tree without the "isozyme" layer (i.e., fluxes), and second a +function that returns the gene product stoichiometry and the turnover number +(again in a tuple) for the given flux identifier. """ -function simplified_isozyme_gene_product_amount_constraints( - isozyme_amounts_stoichiometry_kcat, -) +function simplified_isozyme_gene_product_amount_constraints(x...) res = C.ConstraintTree() - for (iss, stoikcatfn) in isozyme_amounts + for (iss, stoikcatfn) in x for (rid, i) in iss x = stoikcatfn(rid) isnothing(x) && continue - gpstoi, kcat = x + (gpstoi, kcat) = x isnothing(kcat) && continue for (gp, stoi) in gpstoi if haskey(res, gp) res[gp].value += i.value * stoi / kcat else - res[gp] = C.Constraint(i.value * stoi / kcat, 0) + res[gp] = C.Constraint(i.value * stoi / kcat) end end end @@ -279,10 +278,10 @@ function simplified_enzyme_constraints(; mass_cost_reverse = _ -> nothing, capacity_limits = [], ) - function contribution(fl, cost, rid) - c = cost(rid) - isnothing(c) && !haskey(fl, rid) && return zero(C.LinearValue) - return c * fl[rid] + function contribution(fl, cost, id) + c = cost(id) + isnothing(c) && !haskey(fl, id) && return zero(C.LinearValue) + return c * fl[id] end return C.ConstraintTree( diff --git a/src/frontend/enzymes.jl b/src/frontend/enzymes.jl index 4646c893..7cd7dc9a 100644 --- a/src/frontend/enzymes.jl +++ b/src/frontend/enzymes.jl @@ -125,7 +125,7 @@ function enzyme_constrained_flux_balance_constraints( ) return constraints * - sign_split_constraints(; + :directional_flux_balance^sign_split_constraints(; positive = constraints.fluxes_forward, negative = constraints.fluxes_reverse, signed = constraints.fluxes, @@ -191,25 +191,35 @@ function simplified_enzyme_constrained_flux_balance_constraints( interface::Maybe{Symbol} = nothing, interface_name = :interface, ) - + # TODO this deserves a rewrite once more + # prepare some data isozyme_mass(i) = sum( ( ( haskey(gene_product_molar_masses, gid) ? coeff * gene_product_molar_masses[gid] : throw(DomainError(gid, "missing a required gene product molar mass")) - ) for (gid, stoich) in i.gene_product_stoichiometry + ) for (gid, coeff) in i.gene_product_stoichiometry ), init = 0.0, ) + # here I wish we had proper failure handling and patternmatcing in list comprehensions min_isozyme_cost_forward = Dict( - Symbol(rid) => argmin(last, (i, isozyme_mass(i) / i.kcat_forward) for i in is) - for (rid, is) in reaction_isozymes if !isnothing(is.kcat_forward) + Symbol(rid) => argmin( + last, + (i, isozyme_mass(i) / i.kcat_forward) for + (_, i) in is if !isnothing(i.kcat_forward) + ) for (rid, is) in reaction_isozymes if + any(!isnothing(i.kcat_forward) for (_, i) in is) ) min_isozyme_cost_reverse = Dict( - Symbol(rid) => argmin(last, (i, isozyme_mass(i) / i.kcat_reverse) for i in is) - for (rid, is) in reaction_isozymes if !isnothing(is.kcat_reverse) + Symbol(rid) => argmin( + last, + (i, isozyme_mass(i) / i.kcat_reverse) for + (_, i) in is if !isnothing(i.kcat_reverse) + ) for (rid, is) in reaction_isozymes if + any(!isnothing(i.kcat_forward) for (_, i) in is) ) constraints = flux_balance_constraints(model; interface, interface_name) @@ -221,22 +231,22 @@ function simplified_enzyme_constrained_flux_balance_constraints( ) return constraints * - sign_split_constraints(; + :directional_flux_balance^sign_split_constraints(; positive = constraints.fluxes_forward, negative = constraints.fluxes_reverse, signed = constraints.fluxes, ) * - :capacity_limits^simplified_enzyme_constraints( - constraints.fluxes_forward, - constraints.fluxes_reverse, + :capacity_limits^simplified_enzyme_constraints(; + fluxes_forward = constraints.fluxes_forward, + fluxes_reverse = constraints.fluxes_reverse, mass_cost_forward = rid -> - maybemap(last, get(min_cost_forward, rid, nothing)), + maybemap(last, get(min_isozyme_cost_forward, rid, nothing)), mass_cost_reverse = rid -> - maybemap(last, get(min_cost_reverse, rid, nothing)), + maybemap(last, get(min_isozyme_cost_reverse, rid, nothing)), capacity_limits = capacity isa Real ? [( :total_capacity, - Symbol.(a.genes(model)), + Symbol.(A.genes(model)), C.Between(0, capacity), )] : [ @@ -247,14 +257,26 @@ function simplified_enzyme_constrained_flux_balance_constraints( ( constraints.fluxes_forward, rid -> maybemap( - i -> (i.gene_product_stoichiometry, i.kcat_forward), + ic -> ( + Dict( + Symbol(k) => v for + (k, v) in first(ic).gene_product_stoichiometry + ), + first(ic).kcat_forward, + ), get(min_isozyme_cost_forward, rid, nothing), ), ), ( constraints.fluxes_reverse, rid -> maybemap( - i -> (i.gene_product_stoichiometry, i.kcat_reverse), + ic -> ( + Dict( + Symbol(k) => v for + (k, v) in first(ic).gene_product_stoichiometry + ), + first(ic).kcat_reverse, + ), get(min_isozyme_cost_reverse, rid, nothing), ), ), From a631fb86b3a2fc2846ce6a8fd688cf057acc36ee Mon Sep 17 00:00:00 2001 From: Mirek Kratochvil Date: Thu, 7 Mar 2024 09:12:15 +0100 Subject: [PATCH 8/8] fix more stuff, add comments --- src/builders/enzymes.jl | 12 ++++++------ src/frontend/enzymes.jl | 28 ++++++++++++++++++++++------ 2 files changed, 28 insertions(+), 12 deletions(-) diff --git a/src/builders/enzymes.jl b/src/builders/enzymes.jl index ca521479..12d69518 100644 --- a/src/builders/enzymes.jl +++ b/src/builders/enzymes.jl @@ -280,21 +280,21 @@ function simplified_enzyme_constraints(; ) function contribution(fl, cost, id) c = cost(id) - isnothing(c) && !haskey(fl, id) && return zero(C.LinearValue) - return c * fl[id] + (isnothing(c) || !haskey(fl, id)) && return zero(C.LinearValue) + return c * fl[id].value end return C.ConstraintTree( id => C.Constraint(; value = sum( - contribution(fluxes_forward, mass_cost_forward, i) for i in fluxes; + contribution(fluxes_forward, mass_cost_forward, f) for f in fs; init = zero(C.LinearValue), ) + sum( - contribution(fluxes_reverse, mass_cost_reverse, i) for i in fluxes; + contribution(fluxes_reverse, mass_cost_reverse, f) for f in fs; init = zero(C.LinearValue), ), - bound = bound, - ) for (id, fluxes, bound) in capacity_limits + bound, + ) for (id, fs, bound) in capacity_limits ) end diff --git a/src/frontend/enzymes.jl b/src/frontend/enzymes.jl index 7cd7dc9a..dfc0fb1e 100644 --- a/src/frontend/enzymes.jl +++ b/src/frontend/enzymes.jl @@ -181,7 +181,10 @@ P.S., Klamt, S. Automatic construction of metabolic models with enzyme 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). +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. """ function simplified_enzyme_constrained_flux_balance_constraints( model; @@ -191,8 +194,14 @@ function simplified_enzyme_constrained_flux_balance_constraints( interface::Maybe{Symbol} = nothing, interface_name = :interface, ) - # TODO this deserves a rewrite once more - # prepare some data + # TODO this deserves a rewrite once more -- Isozyme struct is hiding a bit + # too much of uncertainty for the code of this thing to be short and + # concise...maybe we should have an isozyme with only one kcat which is + # never missing so that the nothing-handling code and duplicate getters can + # disappear? + + # prepare getters and data + isozyme_mass(i) = sum( ( ( @@ -204,7 +213,6 @@ function simplified_enzyme_constrained_flux_balance_constraints( init = 0.0, ) - # here I wish we had proper failure handling and patternmatcing in list comprehensions min_isozyme_cost_forward = Dict( Symbol(rid) => argmin( last, @@ -222,14 +230,19 @@ function simplified_enzyme_constrained_flux_balance_constraints( any(!isnothing(i.kcat_forward) for (_, i) in is) ) + # allocate the model and variables + 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, ) + # connect everything with constraints + return constraints * :directional_flux_balance^sign_split_constraints(; positive = constraints.fluxes_forward, @@ -246,11 +259,11 @@ function simplified_enzyme_constrained_flux_balance_constraints( capacity_limits = capacity isa Real ? [( :total_capacity, - Symbol.(A.genes(model)), + keys(constraints.fluxes), C.Between(0, capacity), )] : [ - (Symbol(k), Symbol.(gs), C.Between(0, cap)) for (k, gs, cap) in capacity + (Symbol(k), Symbol.(fs), C.Between(0, cap)) for (k, fs, cap) in capacity ], ) * :gene_product_amounts^simplified_isozyme_gene_product_amount_constraints( @@ -281,6 +294,9 @@ function simplified_enzyme_constrained_flux_balance_constraints( ), ), ) + # TODO add some generic functionality that allows people to add additional + # capacity group bounds over gene product amounts. That code can be shared + # with enzyme_constraints. end export simplified_enzyme_constrained_flux_balance_constraints