Skip to content

Commit

Permalink
functionalize enzyme constraint builder
Browse files Browse the repository at this point in the history
  • Loading branch information
@stelmo
stelmo committed Feb 2, 2024
1 parent 173c0aa commit 3b89b47
Show file tree
Hide file tree
Showing 3 changed files with 118 additions and 75 deletions.
2 changes: 1 addition & 1 deletion docs/src/examples/05-enzyme-constrained-models.jl
View file
Original file line number Diff line number Diff line change
Expand Up @@ -331,7 +331,7 @@ ec_solution = enzyme_constrained_flux_balance_analysis(

#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(ec_solution.gene_product_capacity, 50.0, atol = TEST_TOLERANCE) #src
@test isapprox(ec_solution.gene_product_capacity.totalcapacity, 50.0, atol = TEST_TOLERANCE) #src
@test isapprox(ec_solution.fluxes.EX_glc__D_e, -10, atol = TEST_TOLERANCE) #src
@test isapprox( #src
ec_solution.gene_product_amounts.b2417, #src
Expand Down
97 changes: 92 additions & 5 deletions src/builders/enzymes.jl
View file
Original file line number Diff line number Diff line change
Expand Up @@ -47,7 +47,7 @@ For practical purposes, both fluxes and isozymes are here considered to be
unidirectional, i.e., one would typically apply this twice to constraint both
"forward" and "reverse" fluxes.
Function `kcat` should retugn the kcat value for a given reaction and isozyme
Function `kcat` should return the kcat value for a given reaction and isozyme
(IDs of which respectively form the 2 parameters for each call).
"""
function isozyme_flux_constraints(
Expand Down Expand Up @@ -81,8 +81,7 @@ reverse-direction amounts at once, but it is possible to use a single tree,
e.g., in a uni-tuple: `tuple(my_tree)`.
`isozyme_stoichiometry` gets called with a reaction and isozyme ID as given by
the isozyme amount trees. It may return `nothing` in case there's no
information.
the isozyme amount trees.
"""
function gene_product_isozyme_constraints(
gene_product_amounts::C.ConstraintTree,
Expand All @@ -96,9 +95,9 @@ function gene_product_isozyme_constraints(
for (rid, is) in iss
for (iid, i) in is
gpstoi = isozyme_stoichiometry(rid, iid)
isnothing(gpstoi) && continue
isnothing(gpstoi) && continue # TODO: possible footgun
for (gp, stoi) in gpstoi
haskey(gene_product_amounts, gp) || continue
haskey(gene_product_amounts, gp) || continue # TODO: possible footgun
if haskey(res, gp)
res[gp].value += i.value * stoi
else
Expand All @@ -113,3 +112,91 @@ function gene_product_isozyme_constraints(
end

export gene_product_isozyme_constraints

"""
$(TYPEDSIGNATURES)
Returns a constraint tree, built from `constraints`, with enzyme constraints
added based on reactions in `fluxes`.
Inputs:
- `gene_ids` is a list of gene product IDs.
- `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.
- `kcat_forward` and `kcat_reverse` takes `(rid, iso_id)` and returns the
relevant enzyme turnover number, where `rid` is a reaction ID, and `iso_id` is
the isozyme ID of that reaction.
- `isozyme_subunit_stoichiometry` takes `(rid, iso_id)` and returns the gene
product subunit stoichiometry of the associated isozyme, or `nothing` if it
does not exist.
- `gene_product_molar_mass` takes `gid` and returns the molar mass of the gene
product with ID `gid`.
- `capacity_limits` is an iterable container of triples `(capacity_id,
gene_product_ids, capacity_bound)`, which creates capacity bound(s).
"""
function enzyme_constraints(
constraints::C.ConstraintTree;
fluxes = constraints.fluxes,
gene_ids = Symbol[],
isozyme_ids = _ -> nothing,
kcat_forward = (_, _) -> 0.0,
kcat_reverse = (_, _) -> 0.0,
isozyme_subunit_stoichiometry = (_, _) -> nothing,
gene_product_molar_mass = _ -> 0.0,
capacity_limits = [],
)
# TODO: would be nice if directionally constrained

# might be nice to omit some conditionally (e.g. slash the direction if one
# kcat is nothing)
isozyme_amounts = isozyme_amount_variables(
[rid for rid in keys(fluxes) if !isnothing(isozyme_ids(rid))],
rid -> isozyme_ids(rid),
)

# allocate variables for everything (nb. += wouldn't associate right here)
constraints =
constraints +
sign_split_variables(
constraints.fluxes;
positive = :fluxes_forward,
negative = :fluxes_reverse,
) +
: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 *
: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,
kcat_forward,
) *
:isozyme_flux_reverse_balance^isozyme_flux_constraints(
constraints.isozyme_reverse_amounts,
constraints.fluxes_reverse,
kcat_reverse,
) *
:gene_product_isozyme_balance^gene_product_isozyme_constraints(
constraints.gene_product_amounts,
(constraints.isozyme_forward_amounts, constraints.isozyme_reverse_amounts),
isozyme_subunit_stoichiometry,
) *
:gene_product_capacity^C.ConstraintTree(
Symbol(id) => C.Constraint(
value = sum(
constraints.gene_product_amounts[gp].value * gene_product_molar_mass(gp) for gp in gps
),
bound = C.Between(0, limit),
) for (id, gps, limit) in capacity_limits
)
end

export enzyme_constraints
94 changes: 25 additions & 69 deletions src/frontend/enzymes.jl
View file
Original file line number Diff line number Diff line change
Expand Up @@ -54,76 +54,32 @@ function enzyme_constrained_flux_balance_constraints(
gene_product_molar_masses::Dict{String,Float64},
capacity::Union{Vector{Tuple{String,Vector{String},Float64}},Float64},
)
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)])),
# these functions should not fail gracefully to prevent user footguns
# they take strings and but cts use symbols internally, hence lots of conversions are required - kind of annoying
isozyme_ids(rid) =
haskey(reaction_isozymes, String(rid)) ?
Symbol.(keys(reaction_isozymes[String(rid)])) : nothing
kcat_forward(rid, iso_id) = reaction_isozymes[String(rid)][String(iso_id)].kcat_forward
kcat_reverse(rid, iso_id) = reaction_isozymes[String(rid)][String(iso_id)].kcat_reverse
isozyme_subunit_stoichiometry(rid, iso_id) = Dict(
Symbol(k) => v for (k, v) in
reaction_isozymes[String(rid)][String(iso_id)].gene_product_stoichiometry
)

# 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 *
: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
)
gene_product_molar_mass(gid) = gene_product_molar_masses[String(gid)]

capacity_limits =
capacity isa Real ? [("totalcapacity", Symbol.(A.genes(model)), capacity)] :
capacity

enzyme_constraints(
flux_balance_constraints(model);
gene_ids = Symbol.(A.genes(model)),
isozyme_ids,
kcat_forward,
kcat_reverse,
isozyme_subunit_stoichiometry,
gene_product_molar_mass,
capacity_limits,
)
end

Expand Down

0 comments on commit 3b89b47

Please sign in to comment.