Skip to content

Commit

Permalink
fxip
Browse files Browse the repository at this point in the history
  • Loading branch information
@stelmo @exaexa
stelmo authored and exaexa committed Feb 2, 2024
1 parent 173c0aa commit 8e1db92
Show file tree
Hide file tree
Showing 2 changed files with 106 additions and 75 deletions.
94 changes: 89 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,88 @@ function gene_product_isozyme_constraints(
end

export gene_product_isozyme_constraints

"""
$(TYPEDSIGNATURES)
Returns a constraint tree with enzyme constrained added to it. Splits reactions into forward and backward
Function inputs:
- `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,
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(Symbol(rid)))],
rid -> Symbol.(isozyme_ids(rid)),
)

# 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,
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[Symbol(gp)].value *
gene_product_molar_mass[gp] for gp in gps
),
bound = C.Between(0, limit),
) for (id, gps, limit) in capacity_limits
)
constraints
end
87 changes: 17 additions & 70 deletions src/frontend/enzymes.jl
View file
Original file line number Diff line number Diff line change
Expand Up @@ -54,76 +54,23 @@ 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)])),
)

# 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
)
# these functions should not fail gracefully to prevent user footguns
isozyme_ids(rid) = haskey(reaction_isozymes, rid) ? keys(reaction_isozymes[rid]) : nothing
kcat_forward(rid, iso_id) = reaction_isozymes[rid][iso_id].kcat_forward
kcat_reverse(rid, iso_id) = reaction_isozymes[rid][iso_id].kcat_reverse
isozyme_subunit_stoichiometry(rid, iso_id) = reaction_isozymes[rid][iso_id].gene_product_stoichiometry
gene_product_molar_mass(gid) = gene_product_molar_masses[gid]

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

enzyme_constraints(
flux_balance_constraints(model);
isozyme_ids,
kcat_forward,
kcat_reverse,
isozyme_subunit_stoichiometry,
gene_product_molar_mass,
capacity_limits,
)
end

Expand Down

0 comments on commit 8e1db92

Please sign in to comment.