From f888b9d96a61ce760e0433d08c48fd2f3fd8d9ce Mon Sep 17 00:00:00 2001 From: "St. Elmo Wilken" Date: Sat, 27 Jan 2024 12:09:48 +0100 Subject: [PATCH] implement smoment algo --- .../examples/05-enzyme-constrained-models.jl | 13 ++ src/builders/enzymes.jl | 209 ++++++++++++++++++ src/frontend/enzymes.jl | 141 +++++------- 3 files changed, 276 insertions(+), 87 deletions(-) diff --git a/docs/src/examples/05-enzyme-constrained-models.jl b/docs/src/examples/05-enzyme-constrained-models.jl index a978b7358..ce6d9a2f3 100644 --- a/docs/src/examples/05-enzyme-constrained-models.jl +++ b/docs/src/examples/05-enzyme-constrained-models.jl @@ -133,3 +133,16 @@ ec_solution = enzyme_constrained_flux_balance_analysis( 0.011875920383431717, #src atol = TEST_TOLERANCE, #src ) #src + +### Running a simplified enzyme constrained model + +ec_solution2 = simplified_enzyme_constrained_flux_balance_analysis( + model; + reaction_isozymes, + gene_product_molar_masses, + capacity = total_enzyme_capacity, + optimizer = Tulip.Optimizer, + settings = [set_optimizer_attribute("IPM_IterationsLimit", 10_000)], +) + +@test isapprox(ec_solution.objective, ec_solution2.objective, atol = TEST_TOLERANCE) #src diff --git a/src/builders/enzymes.jl b/src/builders/enzymes.jl index e0c8ed828..d4a664943 100644 --- a/src/builders/enzymes.jl +++ b/src/builders/enzymes.jl @@ -14,6 +14,24 @@ # See the License for the specific language governing permissions and # limitations under the License. +""" +$(TYPEDEF) + +A simple struct storing information about the isozyme composition, including +subunit stoichiometry and turnover numbers. Use with +[`enzyme_constrained_flux_balance_analysis`](@ref). + +# Fields +$(TYPEDFIELDS) +""" +Base.@kwdef mutable struct Isozyme + gene_product_stoichiometry::Dict{String,Float64} + kcat_forward::Maybe{Float64} = nothing + kcat_reverse::Maybe{Float64} = nothing +end + +export Isozyme + """ $(TYPEDSIGNATURES) @@ -113,3 +131,194 @@ function gene_product_isozyme_constraints( end export gene_product_isozyme_constraints + +""" +$(TYPEDSIGNATURES) + +Adds enzyme-constraints to `constraints`. Requires `reaction_isozymes`, which is +a mapping of reaction identifiers to [`Isozyme`](@ref) descriptions, and +`gene_product_molar_masses` which is a mapping of gene products to their molar +masses. + +`capacity` may be a single number, which sets the limit for all the gene +products contained in `reaction_isozymes`. Alternatively, `capacity` may be a +vector of identifier-genes-limit triples that make a constraint (identified by +the given identifier) that limits the listed genes to the given limit. +""" +function enzyme_constraints( + constraints::C.ConstraintTree; + fluxes = constraints.fluxes, + reaction_isozymes::Dict{String,Dict{String,Isozyme}}, + gene_product_molar_masses::Dict{String,Float64}, + capacity::Union{Vector{Tuple{String,Vector{String},Float64}},Float64}, +) + + gene_ids = unique([ + Symbol(gid) for isos_dict in values(reaction_isozymes) for + iso in values(isos_dict) for gid in keys(iso.gene_product_stoichiometry) + ]) + + # might be nice to omit some conditionally (e.g. slash the direction if one + # kcat is nothing) + isozyme_amounts = isozyme_amount_variables( + Symbol.(keys(reaction_isozymes)), + rid -> Symbol.(keys(reaction_isozymes[string(rid)])), + ) + + # allocate variables for everything (nb. += wouldn't associate right here) + constraints = + constraints + + :fluxes_forward^unsigned_positive_contribution_variables(fluxes) + + :fluxes_reverse^unsigned_negative_contribution_variables(fluxes) + + :isozyme_forward_amounts^isozyme_amounts + + :isozyme_reverse_amounts^isozyme_amounts + + :gene_product_amounts^C.variables(keys = gene_ids, bounds = C.Between(0, Inf)) + + # connect all parts with constraints + constraints = + constraints * + :directional_flux_balance^sign_split_constraints( + positive = constraints.fluxes_forward, + negative = constraints.fluxes_reverse, + signed = fluxes, + ) * + :isozyme_flux_forward_balance^isozyme_flux_constraints( + constraints.isozyme_forward_amounts, + constraints.fluxes_forward, + (rid, isozyme) -> maybemap( + x -> x.kcat_forward, + maybeget(reaction_isozymes, string(rid), string(isozyme)), + ), + ) * + :isozyme_flux_reverse_balance^isozyme_flux_constraints( + constraints.isozyme_reverse_amounts, + constraints.fluxes_reverse, + (rid, isozyme) -> maybemap( + x -> x.kcat_reverse, + maybeget(reaction_isozymes, string(rid), string(isozyme)), + ), + ) * + :gene_product_isozyme_balance^gene_product_isozyme_constraints( + constraints.gene_product_amounts, + (constraints.isozyme_forward_amounts, constraints.isozyme_reverse_amounts), + (rid, isozyme) -> maybemap( + x -> [(Symbol(k), v) for (k, v) in x.gene_product_stoichiometry], + maybeget(reaction_isozymes, string(rid), string(isozyme)), + ), + ) * + :gene_product_capacity^( + capacity isa Float64 ? + C.Constraint( + value = sum( + gpa.value * gene_product_molar_masses[String(gp)] for + (gp, gpa) in constraints.gene_product_amounts + ), + bound = C.Between(0, capacity), + ) : + C.ConstraintTree( + Symbol(id) => C.Constraint( + value = sum( + constraints.gene_product_amounts[Symbol(gp)].value * + gene_product_molar_masses[gp] for gp in gps + ), + bound = C.Between(0, limit), + ) for (id, gps, limit) in capacity_limits + ) + ) + + constraints +end + +export enzyme_constraints + +""" +$(TYPEDSIGNATURES) + +Adds simplified enzyme constraints to `constraints`. Requires +`reaction_isozymes`, which is a mapping of reaction identifiers to +[`Isozyme`](@ref) descriptions, and `gene_product_molar_masses` which is a +mapping of gene products to their molar masses. Internally, the cheapest and +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`. 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. + +Note: [`simplified_enzyme_constraints`](@ref) and +[`enzyme_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_constraints( + constraints::C.ConstraintTree; + fluxes = constraints.fluxes, + reaction_isozymes::Dict{String,Dict{String,Isozyme}}, + gene_product_molar_masses::Dict{String,Float64}, + capacity::Union{Vector{Tuple{String,Vector{String},Float64}},Float64}, +) + # get fastest isozyme for each reaction (flux * MW/kcat = isozyme mass required) + 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 + ) + + # allocate variables for everything (nb. += wouldn't associate right here) + constraints = + constraints + + :fluxes_forward^unsigned_positive_contribution_variables(fluxes) + + :fluxes_reverse^unsigned_negative_contribution_variables(fluxes) + + # connect all parts with constraints + constraints = + constraints * + :directional_flux_balance^sign_split_constraints( + positive = constraints.fluxes_forward, + negative = constraints.fluxes_reverse, + signed = fluxes, + ) * + :capacity^( + capacity isa Float64 ? + C.Constraint( + value = sum( + C.value(c) * enzyme_speeds[rid].forward for + (rid, c) in constraints.fluxes_forward if haskey(enzyme_speeds, rid) + ) + sum( + C.value(c) * enzyme_speeds[rid].reverse for + (rid, c) in constraints.fluxes_reverse if haskey(enzyme_speeds, rid) + ), + bound = C.Between(0, capacity), + ) : + C.ConstraintTree( + Symbol(id) => C.Constraint( + value = sum( + C.value(get(constraints.fluxes_forward, rid, 0)) * + enzyme_speeds[rid].forward for + rid in flxs if haskey(enzyme_speeds, rid) + ) + sum( + C.value(get(constraints.fluxes_reverse, rid´, 0)) * + enzyme_speeds[rid].reverse for + (rid, c) in flxs if haskey(enzyme_speeds, rid) + ), + bound = C.Between(0, limit), + ) for (id, flxs, limit) in capacity_limits + ) + ) + + constraints + +end diff --git a/src/frontend/enzymes.jl b/src/frontend/enzymes.jl index 14bd8896e..ecf17b2ea 100644 --- a/src/frontend/enzymes.jl +++ b/src/frontend/enzymes.jl @@ -14,24 +14,6 @@ # See the License for the specific language governing permissions and # limitations under the License. -""" -$(TYPEDEF) - -A simple struct storing information about the isozyme composition, including -subunit stoichiometry and turnover numbers. Use with -[`enzyme_constrained_flux_balance_analysis`](@ref). - -# Fields -$(TYPEDFIELDS) -""" -Base.@kwdef mutable struct Isozyme - gene_product_stoichiometry::Dict{String,Float64} - kcat_forward::Maybe{Float64} = nothing - kcat_reverse::Maybe{Float64} = nothing -end - -export Isozyme - """ $(TYPEDSIGNATURES) @@ -47,6 +29,8 @@ the mass of enzymes in the whole model. enzymes". Alternatively, `capacity` may be a vector of identifier-genes-limit triples that make a constraint (identified by the given identifier) that limits the listed genes to the given limit. + +Uses [`enzyme_constraints`](@ref) internally. """ function enzyme_constrained_flux_balance_analysis( model::A.AbstractFBCModel; @@ -58,77 +42,13 @@ function enzyme_constrained_flux_balance_analysis( ) constraints = flux_balance_constraints(model) - # might be nice to omit some conditionally (e.g. slash the direction if one - # kcat is nothing) - isozyme_amounts = isozyme_amount_variables( - Symbol.(keys(reaction_isozymes)), - rid -> Symbol.(keys(reaction_isozymes[string(rid)])), + constraints = enzyme_constraints( + constraints; + reaction_isozymes, + gene_product_molar_masses, + capacity, ) - # allocate variables for everything (nb. += wouldn't associate right here) - constraints = - constraints + - :fluxes_forward^unsigned_positive_contribution_variables(constraints.fluxes) + - :fluxes_reverse^unsigned_negative_contribution_variables(constraints.fluxes) + - :isozyme_forward_amounts^isozyme_amounts + - :isozyme_reverse_amounts^isozyme_amounts + - :gene_product_amounts^C.variables( - keys = Symbol.(A.genes(model)), - bounds = C.Between(0, Inf), - ) - - # connect all parts with constraints - constraints = - constraints * - :directional_flux_balance^sign_split_constraints( - positive = constraints.fluxes_forward, - negative = constraints.fluxes_reverse, - signed = constraints.fluxes, - ) * - :isozyme_flux_forward_balance^isozyme_flux_constraints( - constraints.isozyme_forward_amounts, - constraints.fluxes_forward, - (rid, isozyme) -> maybemap( - x -> x.kcat_forward, - maybeget(reaction_isozymes, string(rid), string(isozyme)), - ), - ) * - :isozyme_flux_reverse_balance^isozyme_flux_constraints( - constraints.isozyme_reverse_amounts, - constraints.fluxes_reverse, - (rid, isozyme) -> maybemap( - x -> x.kcat_reverse, - maybeget(reaction_isozymes, string(rid), string(isozyme)), - ), - ) * - :gene_product_isozyme_balance^gene_product_isozyme_constraints( - constraints.gene_product_amounts, - (constraints.isozyme_forward_amounts, constraints.isozyme_reverse_amounts), - (rid, isozyme) -> maybemap( - x -> [(Symbol(k), v) for (k, v) in x.gene_product_stoichiometry], - maybeget(reaction_isozymes, string(rid), string(isozyme)), - ), - ) * - :gene_product_capacity^( - capacity isa Float64 ? - C.Constraint( - value = sum( - gpa.value * gene_product_molar_masses[String(gp)] for - (gp, gpa) in constraints.gene_product_amounts - ), - bound = C.Between(0, capacity), - ) : - C.ConstraintTree( - Symbol(id) => C.Constraint( - value = sum( - constraints.gene_product_amounts[Symbol(gp)].value * - gene_product_molar_masses[gp] for gp in gps - ), - bound = C.Between(0, limit), - ) for (id, gps, limit) in capacity_limits - ) - ) - optimized_constraints( constraints; objective = constraints.objective.value, @@ -138,3 +58,50 @@ function enzyme_constrained_flux_balance_analysis( end export enzyme_constrained_flux_balance_analysis + + + +""" +$(TYPEDSIGNATURES) + +Run a basic simplified enzyme-constrained flux balance analysis on `model`. +Requires `reaction_isozymes`, which is a mapping of reaction identifiers to +[`Isozyme`](@ref) descriptions, and `gene_product_molar_masses` which is a +mapping of gene products to their molar masses. Internally, the cheapest and +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`. 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. + +Uses [`simplified_enzyme_constraints`](@ref) internally. +""" +function simplified_enzyme_constrained_flux_balance_analysis( + model::A.AbstractFBCModel; + reaction_isozymes::Dict{String,Dict{String,Isozyme}}, + gene_product_molar_masses::Dict{String,Float64}, + capacity::Union{Vector{Tuple{String,Vector{String},Float64}},Float64}, + optimizer, + settings = [], +) + constraints = flux_balance_constraints(model) + + constraints = simplified_enzyme_constraints( + constraints; + reaction_isozymes, + gene_product_molar_masses, + capacity, + ) + + optimized_constraints( + constraints; + objective = constraints.objective.value, + optimizer, + settings, + ) +end + +export simplified_enzyme_constrained_flux_balance_analysis