sMOMENT
sMOMENT algorithm can be used to easily adjust the metabolic activity within the cell to respect known enzymatic parameters and enzyme mass constraints measured by proteomics and other methods.
The original description from sMOMENT is by Bekiaris, and Klamt, "Automatic construction of metabolic models with enzyme constraints.", BMC bioinformatics, 2020
Let's load some packages:
!isfile("e_coli_core.json") &&
download("http://bigg.ucsd.edu/static/models/e_coli_core.json", "e_coli_core.json")
using COBREXA, GLPK
model = load_model("e_coli_core.json")Metabolic model of type JSONModel sparse([9, 51, 55, 64, 65, 34, 44, 59, 66, 64 … 20, 22, 23, 25, 16, 17, 34, 44, 57, 59], [1, 1, 1, 1, 1, 2, 2, 2, 2, 3 … 93, 93, 94, 94, 95, 95, 95, 95, 95, 95], [1.0, 1.0, -1.0, -1.0, 1.0, -1.0, 1.0, -1.0, 1.0, 1.0 … 1.0, -1.0, 1.0, -1.0, -1.0, 1.0, -1.0, 1.0, 1.0, -1.0], 72, 95) Number of reactions: 95 Number of metabolites: 72
We will necessarily need the enzyme turnover numbers (aka "kcats") and masses of the required gene products. You do not necessarily need to know all data for the given model, but the more you have, the better the approximation will be.
For the demonstration purpose, we will generate the data randomly. In a realistic setting, you would input experimental or database-originating data here:
import Random
Random.seed!(1) # repeatability
gene_product_masses = Dict(genes(model) .=> randn(n_genes(model)) .* 10 .+ 60)Dict{String, Float64} with 137 entries:
"b4301" => 69.9065
"b1602" => 70.0268
"b4154" => 55.4119
"b3236" => 50.4756
"b1621" => 59.9869
"b1779" => 56.4588
"b3951" => 64.7369
"b1676" => 54.0343
"b3114" => 56.9763
"b1241" => 60.6193
"b2276" => 55.0068
"b1761" => 63.8676
"b3925" => 55.3065
"b3493" => 47.598
"b3733" => 52.1707
"b2926" => 68.4462
"b0979" => 45.6388
"b4015" => 56.497
"b2296" => 44.2344
⋮ => ⋮We only take the reactions that have gene products (i.e., enzymes) associated with them):
rxns = filter(
x ->
!looks_like_biomass_reaction(x) &&
!looks_like_exchange_reaction(x) &&
!isnothing(reaction_gene_association(model, x)),
reactions(model),
)69-element Vector{String}:
"PFK"
"PFL"
"PGI"
"PGK"
"PGL"
"ACALD"
"AKGt2r"
"PGM"
"PIt2r"
"ALCD2x"
⋮
"MALt2_2"
"MDH"
"ME1"
"ME2"
"NADH16"
"NADTRHD"
"NH4t"
"O2t"
"PDH"The information about each enzyme and its capabilities is stored in an Isozyme structure. For simplicity, sMOMENT ignores much of the information about the multiplicity of required gene products and other.
rxn_isozymes = Dict(
rxn => Isozyme(
Dict(vcat(reaction_gene_association(model, rxn)...) .=> 1),
randn() * 100 + 600, #forward kcat
randn() * 100 + 500, #reverse kcat
) for rxn in rxns
)Dict{String, Isozyme} with 69 entries:
"NH4t" => Isozyme(Dict("b0451"=>1, "s0001"=>1), 669.647, 627.895)
"ACALD" => Isozyme(Dict("b1241"=>1, "b0351"=>1), 592.595, 607.625)
"PTAr" => Isozyme(Dict("b2297"=>1, "b2458"=>1), 584.264, 557.713)
"PGL" => Isozyme(Dict("b0767"=>1), 540.211, 485.094)
"NADTRHD" => Isozyme(Dict("b1603"=>1, "b1602"=>1, "b3962"=>1), 612.8, 608.149)
"ALCD2x" => Isozyme(Dict("b1241"=>1, "b0356"=>1, "b1478"=>1), 462.073, 330.9…
"PGK" => Isozyme(Dict("b2926"=>1), 720.862, 424.687)
"PDH" => Isozyme(Dict("b0114"=>1, "b0115"=>1, "b0116"=>1), 669.985, 380.3…
"LDH_D" => Isozyme(Dict("b1380"=>1, "b2133"=>1), 676.15, 452.057)
"PYK" => Isozyme(Dict("b1854"=>1, "b1676"=>1), 646.569, 533.939)
"CO2t" => Isozyme(Dict("s0001"=>1), 599.569, 371.856)
"PIt2r" => Isozyme(Dict("b2987"=>1, "b3493"=>1), 460.427, 476.199)
"ME1" => Isozyme(Dict("b1479"=>1), 871.006, 523.277)
"MALt2_2" => Isozyme(Dict("b3528"=>1), 524.693, 464.676)
"CS" => Isozyme(Dict("b0720"=>1), 634.028, 438.943)
"PGM" => Isozyme(Dict("b4395"=>1, "b3612"=>1, "b0755"=>1), 433.345, 483.6…
"TKT1" => Isozyme(Dict("b2935"=>1, "b2465"=>1), 615.688, 426.719)
"ATPS4r" => Isozyme(Dict("b3734"=>1, "b3737"=>1, "b3732"=>1, "b3736"=>1, "b3…
"ACONTa" => Isozyme(Dict("b0118"=>1, "b1276"=>1), 379.032, 588.815)
⋮ => ⋮In case some of the reactions are missing in rxn_isozymes, sMOMENT simply ignores them.
Once the data is gathered, we create a model that wraps the original model with additional sMOMENT structure:
smoment_model =
model |> with_smoment(
reaction_isozyme = rxn_isozymes,
gene_product_molar_mass = gene_product_masses,
total_enzyme_capacity = 50.0,
)Metabolic model of type SMomentModel sparse([9, 51, 55, 64, 65, 34, 44, 59, 66, 64 … 23, 25, 23, 25, 16, 17, 34, 44, 57, 59], [1, 1, 1, 1, 1, 2, 2, 2, 2, 3 … 129, 129, 130, 130, 131, 131, 131, 131, 131, 131], [1.0, 1.0, -1.0, -1.0, 1.0, -1.0, 1.0, -1.0, 1.0, -1.0 … -1.0, 1.0, 1.0, -1.0, -1.0, 1.0, -1.0, 1.0, 1.0, -1.0], 72, 131) Number of reactions: 131 Number of metabolites: 72
(You could alternatively use the make_smoment_model to create the structure more manually; but with_smoment is easier to use e.g. with screen.)
In turn, you should have a complete model with unidirectional reactions and additional coupling, as specified by the sMOMENT method:
[stoichiometry(smoment_model); coupling(smoment_model)]73×131 SparseArrays.SparseMatrixCSC{Float64, Int64} with 599 stored entries:
⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠰⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠣⠀⠀⠀⠀⠀⠀⠈⠶⠶⠄⠀⡀⡀⠀⠀⠀
⠂⠀⠓⠂⠒⠂⠘⠩⠓⠃⠛⠚⠂⠘⠘⡊⠓⠀⠂⠀⠀⠀⠀⠐⢁⠀⠀⠑⠈⠓⠒⠚⠓⠊⠡⣖⠓⠒⠀⠀
⠀⠀⠶⠀⠰⠆⠀⠀⠀⠰⠀⢰⠀⠀⠀⠀⠀⠀⡆⠀⠀⠀⠀⠀⠈⠀⠀⠀⠈⢡⡶⣠⣤⢠⢠⠶⠹⢿⣆⠰
⠀⣉⠀⢠⠀⠀⠀⠀⠀⠀⠀⢀⠀⠀⠐⠠⠄⠀⠀⠀⠀⠀⠀⠀⠀⠆⠀⠀⠀⠀⠉⠀⠀⠀⠀⠀⠀⠀⠐⠆
⠀⠀⠆⠀⣒⠀⡘⠀⢃⠃⠒⠛⠀⠈⠀⠈⠁⠒⠀⠀⠀⠀⠀⠀⠀⢀⠀⠂⠈⠠⠞⠒⠀⠐⠀⠀⠈⠁⠀⠀
⠈⠀⠀⠀⠀⠀⠀⠀⢀⠁⠀⠩⠁⣠⡔⠀⠐⠀⠠⠄⢀⠀⠀⠀⠀⠀⠁⠀⠊⠀⠈⠀⠀⢀⠀⠉⠀⠙⠀⠈
⠠⠀⠭⠀⠈⠛⠀⠀⠀⢀⠠⠤⠀⠠⢠⠀⡤⣬⠉⠁⠐⠀⠀⠀⠀⠀⡀⠀⠄⠀⠀⠀⠀⠀⠀⠄⠄⠀⠀⠠
⠄⠤⢠⡄⠀⠀⠎⠉⠄⠤⠴⠦⠀⠂⠀⠀⠀⠤⠐⠒⠀⡀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠤⠤⠀⠠⠀⠀⠀⠀⠀
⢁⠉⣀⠀⠀⠀⠑⠂⠃⢙⢉⣉⠶⢐⠀⠀⠀⣉⠀⢀⡀⠠⠀⠀⠀⠀⢀⠀⠀⠀⠀⠉⠀⠐⠐⠀⡀⠒⠀⢐
⢤⠄⠀⠀⠈⠁⠀⠀⠀⠀⠀⠨⠀⠀⠀⠀⠘⠨⠅⠭⠀⠐⡀⠀⠀⠀⠠⢄⠠⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀
⣠⣄⣀⣀⣀⣀⣀⣀⡀⣀⣀⡠⣀⣀⣀⣀⣐⣠⣤⣤⡄⠀⠑⠀⠀⠀⢠⣀⣚⣢⣤⣀⣀⣀⣀⣀⣀⣀⣀⣀the type (SMomentModel) is a model wrapper – it is a thin additional layer that just adds the necessary sMOMENT-relevant information atop the original model, which is unmodified. That makes the process very efficient and suitable for large-scale data processing. You can still access the original "base" model hidden in the SMomentModel using unwrap_model.
Other than that, the SMomentModel is a model type like any other, and you can run any analysis you want on it, such as FBA:
flux_balance_analysis_dict(smoment_model, GLPK.Optimizer)Dict{String, Float64} with 95 entries:
"ACALD" => 7.21514e-15
"PTAr" => 4.6195
"ALCD2x" => 7.21514e-15
"PDH" => 5.64771
"PYK" => 2.26508
"CO2t" => -5.61213
"EX_nh4_e" => -1.16158
"MALt2_2" => 0.0
"CS" => 0.229832
"PGM" => -6.97218
"TKT1" => 0.0769017
"EX_mal__L_e" => 0.0
"ACONTa" => 0.229832
"EX_pi_e" => -0.783652
"GLNS" => 0.0544703
"ICL" => -7.26379e-16
"EX_o2_e" => -5.36585
"FBA" => 3.65917
"EX_gln__L_e" => 7.27645e-17
⋮ => ⋮(Notice that the total reaction fluxes are reported despite the fact that reactions are indeed split in the model! The underlying mechanism is provided by reaction_flux accessor.)
Variability of the sMOMENT model can be explored as such:
flux_variability_analysis(smoment_model, GLPK.Optimizer, bounds = gamma_bounds(0.95))95×2 Matrix{Float64}:
2.45267 5.27497
0.0 8.65379
1.02326 5.40892
-10.4873 -5.50558
5.32907e-15 2.38384
-1.94964 0.0
-0.364104 0.0
-10.1846 -5.20283
0.74447 0.783652
-1.94964 0.0
⋮
0.0 0.0
-1.3876 0.949426
0.0 0.82907
0.0 1.3876
7.28026 11.5838
0.0 1.43239
1.1035 1.41386
3.64013 6.0218
-3.55271e-15 6.98001...and a sMOMENT model sample can be obtained as usual with sampling:
(
affine_hit_and_run(
smoment_model,
warmup_from_variability(smoment_model, GLPK.Optimizer),
)' * reaction_flux(smoment_model)
)1310×95 Matrix{Float64}:
3.10685 1.10726 2.59042 -5.80754 … 0.219219 2.2906 2.48831
2.69583 0.741132 2.27535 -4.99543 0.194943 2.32345 2.12898
3.04288 1.09169 2.50473 -5.66116 0.232269 2.44912 2.28091
3.3285 1.15695 2.8487 -6.0796 0.2824 2.19665 2.46634
2.97507 1.03514 2.41419 -5.2681 0.177603 2.5011 1.82948
2.98859 1.04911 2.55938 -5.56872 … 0.226673 2.38212 2.25291
3.10703 1.03917 2.60557 -5.70924 0.218012 2.34946 2.3762
3.02855 1.08782 2.4878 -5.61031 0.231883 2.45454 2.23157
2.96803 0.912577 2.43025 -5.58785 0.247873 2.35008 2.34231
2.90577 0.980347 2.37952 -5.39473 0.21198 2.41175 2.10716
⋮ ⋱
2.99797 0.773903 2.58985 -5.53668 0.174147 2.24359 2.29044
2.67869 0.842503 2.15491 -5.03359 0.262696 2.50289 1.9988
2.90926 0.666214 2.41727 -5.28405 0.264136 2.38063 2.24734
2.96929 0.818245 2.5056 -5.42826 0.241253 2.33236 2.20057
2.93624 0.912764 2.44179 -5.41595 … 0.227794 2.37043 2.06027
2.98602 0.893283 2.51572 -5.51114 0.218293 2.3529 2.17738
2.26609 0.749011 1.9116 -4.52937 0.275108 2.36304 1.85706
2.9427 0.840576 2.48082 -5.43956 0.201142 2.32729 2.17416
2.88146 0.823753 2.26378 -5.2138 0.28943 2.47918 2.15513This page was generated using Literate.jl.