Commit 1dd5d7b7 authored by Miroslav Kratochvil's avatar Miroslav Kratochvil
Browse files

many small smoment fixes that should make it just work

parent f4067a14
Pipeline #54970 failed with stages
in 9 minutes and 38 seconds
......@@ -33,7 +33,8 @@ function make_smoment_model(
push!(coupling_row_reaction, i)
mw = sum(
gene_product_molar_mass(gid) * ps for (gid, ps) in isozyme.stoichiometry
gene_product_molar_mass(gid) * ps for
(gid, ps) in isozyme.gene_product_count
)
if min(lbs[i], ubs[i]) < 0 && isozyme.kcat_reverse > _constants.tolerance
......@@ -51,7 +52,7 @@ function make_smoment_model(
)
end
if min(lbs[i], ubs[i]) > 0 && isozyme.kcat_forward > _constants.tolerance
if max(lbs[i], ubs[i]) > 0 && isozyme.kcat_forward > _constants.tolerance
# reaction can run forward
push!(
columns,
......
......@@ -15,20 +15,48 @@ end
# TODO fix the docstring
"""
mutable struct SMomentModel <: ModelWrapper
Construct an enzyme capacity constrained model see `Bekiaris, Pavlos Stephanos,
and Steffen Klamt. "Automatic construction of metabolic models with enzyme
constraints." BMC bioinformatics, 2020.` for implementation details.
Note, `"§"` is reserved for internal use as a delimiter, no reaction id should
contain that character. Also note, SMOMENT assumes that each reaction only has a
single enzyme (one GRR) associated with it. It is required that a model be
modified to ensure that this condition is met. For ease-of-use,
[`remove_slow_isozymes!`](@ref) is supplied to effect this. Currently only
`modifications` that change attributes of the `optimizer` are supported.
"""
mutable struct SMomentModel <: ModelWrapper
struct SMomentModel <: ModelWrapper
Construct an enzyme-capacity constrained model using sMOMENT algorithm, as
described by *Bekiaris, Pavlos Stephanos, and Steffen Klamt, "Automatic
construction of metabolic models with enzyme constraints" BMC bioinformatics,
2020*.
Use [`make_smoment_model`](@ref) or [`with_smoment`](@ref) to construct the
models.
The model is constructed as follows:
- stoichiometry of the original model is retained as much as possible, but
enzymatic reations are split into forward and reverse parts (marked by a
suffix like `...#forward` and `...#reverse`),
- sums of forward and reverse reaction pair fluxes are constrained accordingly
to the original model,
- stoichiometry is expanded by a virtual metabolite "enzyme capacity" which is
consumed by all enzymatic reactions at a rate given by enzyme mass divided by
the corresponding kcat,
- the total consumption of the enzyme capacity is constrained by a fixed
maximum.
The `SMomentModel` structure contains a worked-out representation of the
optimization problem atop a wrapped [`MetabolicModel`](@ref), in particular the
separation of certain reactions into unidirectional forward and reverse parts,
the grouping of these reactions together into virtual "arm" reactions constrained
by bounds from the inner model, an "enzyme capacity" required for each
reaction, and the value of the maximum capacity constraint.
In the structure, field `columns` describes the correspondence of stoichiometry
columns to the stoichiometry and data of the internal wrapped model; field
`coupling_row_reaction` maps the generated coupling constraints to reaction
indexes in the wrapped model, and `total_enzyme_capacity` is the total bound on
the enzyme capacity consumption as specified in sMOMENT algorithm.
This implementation allows easy access to fluxes from the split reactions
(available in `reactions(model)`), while the original "simple" reactions from
the wrapped model are retained as [`fluxes`](@ref). All additional constraints
are implemented using [`coupling`](@ref) and [`coupling_bounds`](@ref).
Original coupling is retained.
"""
struct SMomentModel <: ModelWrapper
columns::Vector{_smoment_column}
coupling_row_reaction::Vector{Int}
total_enzyme_capacity::Float64
......@@ -91,7 +119,7 @@ Get the mapping of the reaction rates in [`SMomentModel`](@ref) to the original
fluxes in the wrapped model.
"""
reaction_flux(model::SMomentModel) =
reaction_flux(model.inner)' * _smoment_column_reactions(model)
reaction_flux(model.inner) * _smoment_column_reactions(model)
"""
coupling(model::SMomentModel)
......@@ -101,7 +129,7 @@ the wrapped model, coupling for split reactions, and the coupling for the total
enzyme capacity.
"""
coupling(model::SMomentModel) = vcat(
coupling(model.inner),
coupling(model.inner) * _smoment_column_reactions(model),
_smoment_reaction_coupling(model),
[col.capacity_required for col in model.columns]',
)
......@@ -122,7 +150,7 @@ The coupling bounds for [`SMomentModel`](@ref) (refer to [`coupling`](@ref) for
details).
"""
function coupling_bounds(model::SMomentModel)
(ilb, iub) = n_coupling_constraints(model.inner)
(ilb, iub) = coupling_bounds(model.inner)
(rlb, rub) = _smoment_reaction_coupling_bounds(model)
return (vcat(ilb, rlb, 0), vcat(iub, rub, model.total_enzyme_capacity))
return (vcat(ilb, rlb, [0.0]), vcat(iub, rub, [model.total_enzyme_capacity]))
end
......@@ -34,39 +34,6 @@ function get_genes_with_kcats(rid_isozymes::Dict{String,Vector{Isozyme}})
return unique(gids)
end
"""
remove_slow_isozymes!(
model::StandardModel,
rid_isozymes = Dict{String, Vector{Isozyme}}(),
)
Remove all but the fastest isozymes from `rid_isozymes`. Use the largest kcat
(for, rev) for these calculations. Modifies `rid_isozymes` in place.
"""
function remove_slow_isozymes!(
model::StandardModel,
rid_isozymes = Dict{String,Vector{Isozyme}}(),
)
for (rid, isozymes) in rid_isozymes
kcat_effs = Float64[]
for isozyme in isozymes
gid_stoich = isozyme.stoichiometry
kcats = isozyme.kcats
push!(
kcat_effs,
dot(
[stoich for stoich in values(gid_stoich)],
[model.genes[gid].molar_mass for gid in keys(gid_stoich)],
) / maximum(kcats),
)
end
idx = argmin(kcat_effs)
rid_isozymes[rid] = [rid_isozymes[rid][idx]]
end
return nothing
end
"""
remove_low_expressed_isozymes!(
model::StandardModel,
......
......@@ -42,8 +42,7 @@ _smoment_reaction_coupling(model::SMomentModel) = sparse(
Internal helper for determining the number of required couplings to account for
"arm" reactions.
"""
_smoment_n_reaction_couplings(model::SMomentModel) =
length(model.coupling_row_reaction)
_smoment_n_reaction_couplings(model::SMomentModel) = length(model.coupling_row_reaction)
"""
_smoment_reaction_coupling_bounds(model::SMomentModel)
......@@ -53,5 +52,37 @@ values are taken from the "original" inner model.
"""
_smoment_reaction_coupling_bounds(model::SMomentModel) =
let (lbs, ubs) = bounds(model.inner)
(lbs[coupling_row_reaction], ubs[coupling_row_reaction])
(lbs[model.coupling_row_reaction], ubs[model.coupling_row_reaction])
end
"""
smoment_isozyme_speed(isozyme::Isozyme, gene_product_molar_mass)
Compute a "score" for picking the most viable isozyme for
[`make_smoment_model`](@ref), based on maximum kcat divided by relative mass of
the isozyme. This is used because sMOMENT algorithm can not handle multiple
isozymes for one reaction.
"""
smoment_isozyme_speed(isozyme::Isozyme, gene_product_molar_mass) =
max(isozyme.kcat_forward, isozyme.kcat_reverse) / sum(
count * gene_product_molar_mass(gene) for
(gene, count) in isozyme.gene_product_count
)
"""
smoment_isozyme_speed(gene_product_molar_mass::Function)
A piping- and argmax-friendly overload of [`smoment_isozyme_speed`](@ref).
# Example
```
gene_mass_function = gid -> 1.234
best_isozyme_for_smoment = argmax(
smoment_isozyme_speed(gene_mass_function),
my_isozyme_vector,
)
```
"""
smoment_isozyme_speed(gene_product_molar_mass::Function) =
isozyme -> smoment_isozyme_speed(isozyme, gene_product_molar_mass)
@testset "SMOMENT" begin
model = load_model(StandardModel, model_paths["e_coli_core.json"])
total_protein_mass = 100 # mg/gdW
#: construct isozymes from model
rid_isozymes = Dict{String,Vector{Isozyme}}()
for (rid, kcats) in ecoli_core_reaction_kcats
grrs = reaction_gene_association(model, rid)
rid_isozymes[rid] = [
Isozyme(
Dict(grrs[i] .=> ecoli_core_protein_stoichiometry[rid][i]),
(kcats[i][1], kcats[i][2]),
) for i = 1:length(grrs)
]
end
get_gene_product_mass = gid -> get(ecoli_core_protein_masses, gid, 0.0)
#: add molar mass to genes in model
for (gid, g) in model.genes
model.genes[gid].molar_mass = get(ecoli_core_protein_masses, gid, nothing)
end
get_reaction_isozyme =
rid ->
haskey(ecoli_core_reaction_kcats, rid) ?
argmax(
smoment_isozyme_speed(get_gene_product_mass),
Isozyme(
Dict(grr .=> ecoli_core_protein_stoichiometry[rid][i]),
ecoli_core_reaction_kcats[rid][i]...,
) for (i, grr) in enumerate(reaction_gene_association(model, rid))
) : nothing
remove_slow_isozymes!(model, rid_isozymes)
smm = make_smomentmodel(model; rid_isozymes, enzyme_capacity = total_protein_mass)
change_bounds!(
smm,
["EX_glc__D_e", "GLCpts"];
lower = [-1000.0, -1.0],
upper = [nothing, 12.0],
)
smoment_model =
model |>
with_changed_bounds(
["EX_glc__D_e", "GLCpts"],
lower = [-1000.0, -1.0],
upper = [nothing, 12.0],
) |>
with_smoment(
reaction_isozymes = get_reaction_isozyme,
gene_product_molar_mass = get_gene_product_mass,
total_enzyme_capacity = 100.0,
)
rxn_fluxes = flux_balance_analysis_dict(
smm,
smoment_model,
Tulip.Optimizer;
modifications = [change_optimizer_attribute("IPM_IterationsLimit", 1000)],
)
......
......@@ -73,7 +73,7 @@ test_toyModel() = CoreModel(
["m1[c]", "m3[c]", "m2[c]", "m1[e]", "m3[e]", "biomass[c]"],
)
const reaction_standard_gibbs_free_energies = Dict(
const reaction_standard_gibbs_free_energies = Dict{String,Float64}(
#=
ΔᵣG⁰ data from Equilibrator using the E. coli core model's reactions
To generate this data manually, go to https://equilibrator.weizmann.ac.il/ and
......@@ -155,7 +155,7 @@ const reaction_standard_gibbs_free_energies = Dict(
"FUM" => -3.424133018702122,
)
const ecoli_core_protein_masses = Dict(
const ecoli_core_protein_masses = Dict{String,Float64}(
#=
Data downloaded from Uniprot for E. coli K12,
gene mass in kDa. To obtain these data yourself, go to
......@@ -300,7 +300,7 @@ const ecoli_core_protein_masses = Dict(
"b2279" => 10.845,
)
const ecoli_core_protein_stoichiometry = Dict(
const ecoli_core_protein_stoichiometry = Dict{String,Vector{Vector{Float64}}}(
#=
Data made up, each isozyme is assumed to be composed of
only one subunit each.
......@@ -381,7 +381,7 @@ const ecoli_core_protein_stoichiometry = Dict(
"FUM" => [[1.0], [1.0], [1.0]],
)
const ecoli_core_reaction_kcats = Dict(
const ecoli_core_reaction_kcats = Dict{String,Vector{Tuple{Float64,Float64}}}(
#=
Data taken from Heckmann, David, et al. "Machine learning applied to enzyme
turnover numbers reveals protein structural correlates and improves metabolic
......@@ -389,164 +389,164 @@ const ecoli_core_reaction_kcats = Dict(
kcats are the same, and each isozyme has the same kcat.
=#
"ACALD" =>
[[568.1130792316333, 568.1130792316333], [568.856126503717, 568.856126503717]],
[(568.1130792316333, 568.1130792316333), (568.856126503717, 568.856126503717)],
"PTAr" => [
[1171.9703624351055, 1171.9703624351055],
[1173.7231032615289, 1173.7231032615289],
(1171.9703624351055, 1171.9703624351055),
(1173.7231032615289, 1173.7231032615289),
],
"ALCD2x" => [
[75.9547881894345, 75.9547881894345],
[75.96334310351442, 75.96334310351442],
[76.1472359297987, 76.1472359297987],
(75.9547881894345, 75.9547881894345),
(75.96334310351442, 75.96334310351442),
(76.1472359297987, 76.1472359297987),
],
"PDH" => [[529.7610874857239, 529.7610874857239]],
"PDH" => [(529.7610874857239, 529.7610874857239)],
"PYK" => [
[422.0226052080562, 422.0226052080562],
[422.1332899347833, 422.1332899347833],
(422.0226052080562, 422.0226052080562),
(422.1332899347833, 422.1332899347833),
],
"MALt2_2" => [[234.03664660088714, 234.03664660088714]],
"CS" => [[113.29607453875758, 113.29607453875758]],
"MALt2_2" => [(234.03664660088714, 234.03664660088714)],
"CS" => [(113.29607453875758, 113.29607453875758)],
"PGM" => [
[681.4234715886669, 681.4234715886669],
[681.6540601244343, 681.6540601244343],
[680.5234799168278, 680.5234799168278],
(681.4234715886669, 681.4234715886669),
(681.6540601244343, 681.6540601244343),
(680.5234799168278, 680.5234799168278),
],
"TKT1" => [
[311.16139580671637, 311.16139580671637],
[311.20967965149947, 311.20967965149947],
(311.16139580671637, 311.16139580671637),
(311.20967965149947, 311.20967965149947),
],
"ACONTa" => [
[191.02308213992006, 191.02308213992006],
[191.03458045697235, 191.03458045697235],
(191.02308213992006, 191.02308213992006),
(191.03458045697235, 191.03458045697235),
],
"GLNS" => [
[89.83860937287024, 89.83860937287024],
[89.82177852142014, 89.82177852142014],
(89.83860937287024, 89.83860937287024),
(89.82177852142014, 89.82177852142014),
],
"ICL" => [[17.45922330097792, 17.45922330097792]],
"ICL" => [(17.45922330097792, 17.45922330097792)],
"FBA" => [
[373.425646787578, 373.425646787578],
[372.74936053215833, 372.74936053215833],
[372.88627228768166, 372.88627228768166],
(373.425646787578, 373.425646787578),
(372.74936053215833, 372.74936053215833),
(372.88627228768166, 372.88627228768166),
],
"FORt2" => [
[233.93045260179326, 233.93045260179326],
[233.84804009142908, 233.84804009142908],
(233.93045260179326, 233.93045260179326),
(233.84804009142908, 233.84804009142908),
],
"G6PDH2r" => [[589.3761070080022, 589.3761070080022]],
"AKGDH" => [[264.48071159327156, 264.48071159327156]],
"G6PDH2r" => [(589.3761070080022, 589.3761070080022)],
"AKGDH" => [(264.48071159327156, 264.48071159327156)],
"TKT2" => [
[467.4226876901618, 467.4226876901618],
[468.1440593542596, 468.1440593542596],
(467.4226876901618, 467.4226876901618),
(468.1440593542596, 468.1440593542596),
],
"FRD7" => [[90.20637824912605, 90.20637824912605]],
"SUCOAS" => [[18.494387648707622, 18.494387648707622]],
"FRD7" => [(90.20637824912605, 90.20637824912605)],
"SUCOAS" => [(18.494387648707622, 18.494387648707622)],
"FBP" => [
[568.5346256470805, 568.5346256470805],
[567.6367759041788, 567.6367759041788],
(568.5346256470805, 568.5346256470805),
(567.6367759041788, 567.6367759041788),
],
"ICDHyr" => [[39.62446791678959, 39.62446791678959]],
"AKGt2r" => [[234.99097804446805, 234.99097804446805]],
"GLUSy" => [[33.262997317319055, 33.262997317319055]],
"TPI" => [[698.301904211076, 698.301904211076]],
"ICDHyr" => [(39.62446791678959, 39.62446791678959)],
"AKGt2r" => [(234.99097804446805, 234.99097804446805)],
"GLUSy" => [(33.262997317319055, 33.262997317319055)],
"TPI" => [(698.301904211076, 698.301904211076)],
"FORt" => [
[234.38391855848187, 234.38391855848187],
[234.34725576182922, 234.34725576182922],
(234.38391855848187, 234.38391855848187),
(234.34725576182922, 234.34725576182922),
],
"ACONTb" => [
[159.74612206327865, 159.74612206327865],
[159.81975755249232, 159.81975755249232],
(159.74612206327865, 159.74612206327865),
(159.81975755249232, 159.81975755249232),
],
"GLNabc" => [[233.80358131677775, 233.80358131677775]],
"GLNabc" => [(233.80358131677775, 233.80358131677775)],
"RPE" => [
[1772.4850826683305, 1772.4850826683305],
[1768.8536177485582, 1768.8536177485582],
(1772.4850826683305, 1772.4850826683305),
(1768.8536177485582, 1768.8536177485582),
],
"ACKr" => [
[554.611547307207, 554.611547307207],
[555.112707891257, 555.112707891257],
[555.2464368932744, 555.2464368932744],
(554.611547307207, 554.611547307207),
(555.112707891257, 555.112707891257),
(555.2464368932744, 555.2464368932744),
],
"THD2" => [[24.739139801185537, 24.739139801185537]],
"THD2" => [(24.739139801185537, 24.739139801185537)],
"PFL" => [
[96.56316095411077, 96.56316095411077],
[96.65024313036014, 96.65024313036014],
[96.60761818004025, 96.60761818004025],
[96.49541118899961, 96.49541118899961],
(96.56316095411077, 96.56316095411077),
(96.65024313036014, 96.65024313036014),
(96.60761818004025, 96.60761818004025),
(96.49541118899961, 96.49541118899961),
],
"RPI" => [
[51.771578021074234, 51.771578021074234],
[51.81603467243345, 51.81603467243345],
(51.771578021074234, 51.771578021074234),
(51.81603467243345, 51.81603467243345),
],
"D_LACt2" => [
[233.51709131524734, 233.51709131524734],
[233.83187606098016, 233.83187606098016],
(233.51709131524734, 233.51709131524734),
(233.83187606098016, 233.83187606098016),
],
"TALA" => [
[109.05210545422884, 109.05210545422884],
[109.04246437049026, 109.04246437049026],
(109.05210545422884, 109.05210545422884),
(109.04246437049026, 109.04246437049026),
],
"PPCK" => [[218.4287805666016, 218.4287805666016]],
"PGL" => [[2120.4297518987964, 2120.4297518987964]],
"PPCK" => [(218.4287805666016, 218.4287805666016)],
"PGL" => [(2120.4297518987964, 2120.4297518987964)],
"NADTRHD" => [
[186.99387360624777, 186.99387360624777],
[187.16629305266423, 187.16629305266423],
(186.99387360624777, 186.99387360624777),
(187.16629305266423, 187.16629305266423),
],
"PGK" => [[57.641966636896335, 57.641966636896335]],
"PGK" => [(57.641966636896335, 57.641966636896335)],
"LDH_D" => [
[31.11118891764946, 31.11118891764946],
[31.12493425054357, 31.12493425054357],
(31.11118891764946, 31.11118891764946),
(31.12493425054357, 31.12493425054357),
],
"ME1" => [[487.0161203971232, 487.0161203971232]],
"ME1" => [(487.0161203971232, 487.0161203971232)],
"PIt2r" => [
[233.8651331835765, 233.8651331835765],
[234.27374798581067, 234.27374798581067],
(233.8651331835765, 233.8651331835765),
(234.27374798581067, 234.27374798581067),
],
"ATPS4r" => [
[7120.878030435999, 7120.878030435999],
[7116.751386037507, 7116.751386037507],
(7120.878030435999, 7120.878030435999),
(7116.751386037507, 7116.751386037507),
],
"GLCpts" => [
[233.9009878400008, 233.9009878400008],
[233.66656882114864, 233.66656882114864],
[233.66893882934883, 233.66893882934883],
(233.9009878400008, 233.9009878400008),
(233.66656882114864, 233.66656882114864),
(233.66893882934883, 233.66893882934883),
],
"GLUDy" => [[105.32811069172409, 105.32811069172409]],
"GLUDy" => [(105.32811069172409, 105.32811069172409)],
"CYTBD" => [
[153.18512795009505, 153.18512795009505],
[153.2429537682265, 153.2429537682265],
],
"FUMt2_2" => [[234.37495609395967, 234.37495609395967]],
"FRUpts2" => [[234.1933863380989, 234.1933863380989]],
"GAPD" => [[128.76795529111456, 128.76795529111456]],
"PPC" => [[165.52424516841342, 165.52424516841342]],
"NADH16" => [[971.7487306963936, 971.7487306963936]],
(153.18512795009505, 153.18512795009505),
(153.2429537682265, 153.2429537682265),
],
"FUMt2_2" => [(234.37495609395967, 234.37495609395967)],
"FRUpts2" => [(234.1933863380989, 234.1933863380989)],
"GAPD" => [(128.76795529111456, 128.76795529111456)],
"PPC" => [(165.52424516841342, 165.52424516841342)],
"NADH16" => [(971.7487306963936, 971.7487306963936)],
"PFK" => [
[1000.4626204522712, 1000.4626204522712],
[1000.5875517343595, 1000.5875517343595],
],
"MDH" => [[25.931655783969283, 25.931655783969283]],
"PGI" => [[468.11833198138834, 468.11833198138834]],
"ME2" => [[443.0973626307168, 443.0973626307168]],
"GND" => [[240.1252264230952, 240.1252264230952]],
"SUCCt2_2" => [[234.18109388303225, 234.18109388303225]],
(1000.4626204522712, 1000.4626204522712),
(1000.5875517343595, 1000.5875517343595),
],
"MDH" => [(25.931655783969283, 25.931655783969283)],
"PGI" => [(468.11833198138834, 468.11833198138834)],
"ME2" => [(443.0973626307168, 443.0973626307168)],
"GND" => [(240.1252264230952, 240.1252264230952)],
"SUCCt2_2" => [(234.18109388303225, 234.18109388303225)],
"GLUN" => [
[44.76358496525738, 44.76358496525738],
[44.84850207360875, 44.84850207360875],
[44.76185250415503, 44.76185250415503],
(44.76358496525738, 44.76358496525738),
(44.84850207360875, 44.84850207360875),
(44.76185250415503, 44.76185250415503),
],
"ADK1" => [[111.64869652600649, 111.64869652600649]],
"SUCDi" => [[680.3193833053011, 680.3193833053011]],
"ENO" => [[209.35855069219886, 209.35855069219886]],
"ADK1" => [(111.64869652600649, 111.64869652600649)],
"SUCDi" => [(680.3193833053011, 680.3193833053011)],
"ENO" => [(209.35855069219886, 209.35855069219886)],
"MALS" => [
[252.7540503869977, 252.7540503869977],
[252.2359738678874, 252.2359738678874],
(252.7540503869977, 252.7540503869977),
(252.2359738678874, 252.2359738678874),
],
"GLUt2r" => [[234.22890837451837, 234.22890837451837]],
"PPS" => [[706.1455885214322, 706.1455885214322]],
"GLUt2r" => [(234.22890837451837, 234.22890837451837)],
"PPS" => [(706.1455885214322, 706.1455885214322)],
"FUM" => [
[1576.8372583425075, 1576.8372583425075],
[1576.233088455828, 1576.233088455828],
[1575.9638204848736, 1575.9638204848736],
(1576.8372583425075, 1576.8372583425075),
(1576.233088455828, 1576.233088455828),
(1575.9638204848736, 1575.9638204848736),
],
)
Supports Markdown
0% or .
You are about to add 0 people to the discussion. Proceed with caution.
Finish editing this message first!
Please register or to comment