Growth media analysis
Nutrient availability is a major driving factor for growth of microorganisms and energy production in cells. Here, we demonstrate two main ways to examine the nutrient consumption with COBREXA.jl: Simulating deficiency of nutrients, and finding the minimal flux of nutrients required to support certain model output.
As always, we work on the toy model of E. coli:
using COBREXA, GLPK
!isfile("e_coli_core.xml") &&
download("http://bigg.ucsd.edu/static/models/e_coli_core.xml", "e_coli_core.xml")
model = load_model(StandardModel, "e_coli_core.xml")Metabolic model of type StandardModel sparse([8, 10, 21, 43, 50, 51, 8, 9, 6, 12 … 33, 66, 68, 72, 23, 26, 33, 72, 22, 33], [1, 1, 1, 1, 1, 1, 2, 2, 3, 3 … 93, 93, 93, 93, 94, 94, 94, 94, 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
What nutrients does my model need to grow?
The models usually ingest nutrients through exchange reactions. By changing the bounds on the exchange reactions, you can limit the intake of the nutrients and thus simulate the nutrient deficiency. If applied programatically to multiple exchanges, this can give you a good overview of what nutrients impact the model most.
To check the viability of a single nutrient, you can simply change a bound on a selected exchange reaction and simulate the model with a limited amount.
biomass = "R_BIOMASS_Ecoli_core_w_GAM"
model_limited = change_bound(model, "R_EX_glc__D_e", lower = -1.0)Metabolic model of type StandardModel sparse([8, 10, 21, 43, 50, 51, 8, 9, 6, 12 … 33, 66, 68, 72, 23, 26, 33, 72, 22, 33], [1, 1, 1, 1, 1, 1, 2, 2, 3, 3 … 93, 93, 93, 93, 94, 94, 94, 94, 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
By a convention, the direction of exchange reaction usually goes from the model into the environment, representing the "production". Limiting the intake thus happens by disabling the "negative production", i.e., placing a lower bound.
original_production = flux_balance_analysis_dict(model, GLPK.Optimizer)[biomass]
limited_production = flux_balance_analysis_dict(model_limited, GLPK.Optimizer)[biomass]
original_production, limited_production(0.8739215069684118, 0.0483849722968142)
Function flux_summary can help with quickly spotting what has changed:
flux_summary(flux_balance_analysis_dict(model_limited, GLPK.Optimizer))Biomass R_BIOMASS_Ecoli_core_w_GAM: 0.0484 Import R_EX_o2_e: -3.885 R_EX_glc__D_e: -1.0 R_EX_nh4_e: -0.2638 R_EX_pi_e: -0.178 Export R_EX_h_e: 0.9706 R_EX_co2_e: 3.941 R_EX_h2o_e: 4.2934
Similarly, you can check that the model can survive without oxygen, at the cost of switching the metabolism to ethanol production:
flux_summary(
flux_balance_analysis_dict(
change_bound(model, "R_EX_o2_e", lower = 0.0),
GLPK.Optimizer,
),
)Biomass R_BIOMASS_Ecoli_core_w_GAM: 0.2117 Import R_EX_glc__D_e: -10.0 R_EX_h2o_e: -7.1158 R_EX_nh4_e: -1.1542 R_EX_pi_e: -0.7786 R_EX_co2_e: -0.3782 Export R_EX_etoh_e: 8.2795 R_EX_ac_e: 8.5036 R_EX_for_e: 17.8047 R_EX_h_e: 30.5542
The effect of all nutrients on the metabolism can be scanned using screen. The change_bound function is, for this purpose, packed in a variant specified with_changed_bound:
exchanges = filter(looks_like_exchange_reaction, reactions(model))
exchanges .=> screen(
model,
variants = [[with_changed_bound(exchange, lower = 0.0)] for exchange in exchanges],
analysis = m -> begin
res = flux_balance_analysis_dict(m, GLPK.Optimizer)
isnothing(res) ? nothing : res[biomass]
end,
)20-element Vector{Pair{String}}:
"R_EX_ac_e" => 0.8739215069684118
"R_EX_acald_e" => 0.8739215069684118
"R_EX_akg_e" => 0.8739215069684118
"R_EX_co2_e" => 0.8739215069684118
"R_EX_etoh_e" => 0.8739215069684118
"R_EX_for_e" => 0.8739215069684118
"R_EX_fru_e" => 0.8739215069684118
"R_EX_fum_e" => 0.8739215069684118
"R_EX_glc__D_e" => nothing
"R_EX_gln__L_e" => 0.8739215069684118
"R_EX_glu__L_e" => 0.8739215069684118
"R_EX_h2o_e" => 0.8739215069684118
"R_EX_h_e" => 0.8739215069684293
"R_EX_lac__D_e" => 0.8739215069684118
"R_EX_mal__L_e" => 0.8739215069684118
"R_EX_nh4_e" => 2.9250893503872995e-16
"R_EX_o2_e" => 0.21166294973531524
"R_EX_pi_e" => -1.340006737455969e-15
"R_EX_pyr_e" => 0.8739215069684118
"R_EX_succ_e" => 0.8739215069684118Similarly to gene knockouts, you can also examine the effect of combined nutrient deficiencies. To obtain a more interesting result, we may examine the effect of slight deficiencies of pairs of intake metabolites. For simplicity, we show the result only on a small subset of the exchanges:
selected_exchanges = [
"R_EX_pi_e",
"R_EX_gln__L_e",
"R_EX_nh4_e",
"R_EX_pyr_e",
"R_EX_fru_e",
"R_EX_glu__L_e",
"R_EX_glc__D_e",
"R_EX_o2_e",
]
screen(
model,
variants = [
[with_changed_bounds([e1, e2], lower = [-1.0, -0.1])] for e1 in selected_exchanges,
e2 in selected_exchanges
],
analysis = m -> begin
res = flux_balance_analysis_dict(m, GLPK.Optimizer)
isnothing(res) ? nothing : res[biomass]
end,
)8×8 Matrix{Union{Nothing, Float64}}:
0.0271835 0.271835 0.0183392 … 0.271835 nothing 0.215263
0.0271835 0.880924 0.385123 0.946959 0.0248973 0.222998
0.0271835 0.22007 0.0183392 0.201731 nothing 0.183392
0.0271835 0.915298 0.0183392 0.915171 nothing 0.22357
0.0271835 0.972588 0.0183392 0.972461 0.0576796 0.245722
0.0271835 0.949672 0.201731 … 0.880796 0.028908 0.227731
0.0271835 0.0549658 0.0183392 0.0552895 nothing nothing
0.0271835 0.250863 0.0183392 0.250983 nothing 0.215263The result describes combinations of nutrient deficiencies – the nutrient that corresponds to the row is mildly deficient (limited to uptake 1.0), and the one that corresponds to the column is severely limited (to uptake 0.1).
To speed up larger analyses, remember that execution of screen can be parallelized to gain speedup. Parallelization in screen is optimized to avoid unnecessary data transfers that may occur when using trivial pmap.
What is the minimal flux of nutrients for my model to grow?
You can compute the minimal flux (i.e., mass per time) of required nutrients by constraining the model growth to a desired lower bound, and then optimize the model with an objective that minimizes intake of all exchanges (i.e., given the directionality convention of the exchanges, actually maximizes the flux through all exchange reactions along their direction).
model_with_bounded_production = change_bound(model, biomass, lower = 0.1) #minimum required growth
minimal_intake_production = flux_balance_analysis_dict(
model_with_bounded_production,
GLPK.Optimizer,
modifications = [change_objective(exchanges)],
);Metabolite "cost" data may be supplemented using the weights argument of change_objective, to reflect e.g. the different molar masses or energetic values of different nutrients.
In our simple case, we obtain the following minimal required intake:
flux_summary(minimal_intake_production)Biomass R_BIOMASS_Ecoli_core_w_GAM: 0.1 Import R_EX_o2_e: -36.9074 R_EX_glc__D_e: -10.0 R_EX_nh4_e: -0.5453 R_EX_pi_e: -0.3679 Export R_EX_co2_e: 18.3015 R_EX_h2o_e: 19.0299 R_EX_for_e: 37.443 R_EX_h_e: 39.449
This page was generated using Literate.jl.