# 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") &&

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

Exchange directions

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.8739215069684118

Similarly 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.215263

The 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).

Screening can be easily parallelized

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