Edit on GitHub

Building and analysing a small community model

Here we will use COBREXA to build and analyze a small community model consisting of three E. coli mutants using the CoreModel. We will use an objective function that enforces equal growth rates.

We will first construct a community of only two mutants to illustrate the effect of the community biomass objective function. Then we will add a third member that has a lethal knockout. Due to the bounds on the exchange reactions these three models are incapable of sharing resources - hence the maximum growth rate will be zero. By changing the bounds we can allow resource sharing, "saving" the community.

Load the base model

!isfile("e_coli_core.json") &&
    download("http://bigg.ucsd.edu/static/models/e_coli_core.json", "e_coli_core.json");

using COBREXA
using Tulip

Load the models and inspect fba solutions

base_model = load_model(CoreModel, "e_coli_core.json") # base from from which the knockouts will be constructed

cytbd_knockout_model = remove_reaction(base_model, "CYTBD") # knockout the CYTBD (cytochrome oxidase) reaction
sol = flux_balance_analysis_dict(cytbd_knockout_model, Tulip.Optimizer)
sol["BIOMASS_Ecoli_core_w_GAM"] # Cytochrome oxidase knockout μ (growth rate)
0.21166294865467686
atps4r_knockout_model = remove_reaction(base_model, "ATPS4r") # knockout the ATP synthase reaction
sol = flux_balance_analysis_dict(atps4r_knockout_model, Tulip.Optimizer)
sol["BIOMASS_Ecoli_core_w_GAM"] # ATP synthase knockout μ
0.37422987390664064
eno_knockout_model = remove_reaction(base_model, "ENO") # knockout the enolase reaction
sol = flux_balance_analysis_dict(eno_knockout_model, Tulip.Optimizer)
sol["BIOMASS_Ecoli_core_w_GAM"] # Enolase knockout μ, cannot grow by itself
9.545439179585187e-16

Build a community model of the cytochrome oxidase knockout and the ATP synthase knockout models

ex_rxn_mets = Dict(
    ex_rxn => first(keys(reaction_stoichiometry(base_model, ex_rxn))) for
    ex_rxn in filter(looks_like_exchange_reaction, reactions(base_model))
) # identify exchange reactions heuristically
Dict{String, String} with 20 entries:
  "EX_for_e"    => "for_e"
  "EX_nh4_e"    => "nh4_e"
  "EX_pyr_e"    => "pyr_e"
  "EX_co2_e"    => "co2_e"
  "EX_ac_e"     => "ac_e"
  "EX_h2o_e"    => "h2o_e"
  "EX_succ_e"   => "succ_e"
  "EX_akg_e"    => "akg_e"
  "EX_fru_e"    => "fru_e"
  "EX_acald_e"  => "acald_e"
  "EX_fum_e"    => "fum_e"
  "EX_mal__L_e" => "mal__L_e"
  "EX_h_e"      => "h_e"
  "EX_pi_e"     => "pi_e"
  "EX_lac__D_e" => "lac__D_e"
  "EX_glu__L_e" => "glu__L_e"
  "EX_o2_e"     => "o2_e"
  "EX_gln__L_e" => "gln__L_e"
  "EX_glc__D_e" => "glc__D_e"
  "EX_etoh_e"   => "etoh_e"
model_names = ["cytbd_ko", "atps4r_ko"]
community_model = join_with_exchanges(
    CoreModel,
    [cytbd_knockout_model, atps4r_knockout_model],
    ex_rxn_mets;
    biomass_ids = ["BIOMASS_Ecoli_core_w_GAM", "BIOMASS_Ecoli_core_w_GAM"],
    model_names = model_names,
)
Metabolic model of type CoreModel
sparse([9, 51, 55, 64, 65, 34, 44, 59, 66, 64  …  155, 156, 157, 158, 159, 160, 161, 162, 163, 164], [1, 1, 1, 1, 1, 2, 2, 2, 2, 3  …  199, 200, 201, 202, 203, 204, 205, 206, 207, 208], [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], 166, 209)
Number of reactions: 209
Number of metabolites: 166

Set exchange reaction bounds of community model based on the bounds of the individual models

env_ex_rxn_idxs = indexin(keys(ex_rxn_mets), reactions(community_model)) # identify the global (environmental exchange reactions)
cytbd_ex_rxn_idxs = indexin(keys(ex_rxn_mets), reactions(cytbd_knockout_model)) # identify the indices of the corresponding exchange reactions in the original models
atps4r_ex_rxn_idxs = indexin(keys(ex_rxn_mets), reactions(atps4r_knockout_model))
20-element Vector{Union{Nothing, Int64}}:
 48
 58
 61
 46
 43
 55
 62
 45
 49
 44
 50
 57
 54
 60
 56
 53
 59
 52
 51
 47

In case some exchange reactions are not present in both models, set environmental exchange bound to the sum of the individual exchange bounds

for (env_ex, m2_ex, m1_ex) in zip(env_ex_rxn_idxs, cytbd_ex_rxn_idxs, atps4r_ex_rxn_idxs)
    m2lb = isnothing(m2_ex) ? 0.0 : atps4r_knockout_model.xl[m2_ex]
    m2ub = isnothing(m2_ex) ? 0.0 : atps4r_knockout_model.xu[m2_ex]
    m1lb = isnothing(m1_ex) ? 0.0 : cytbd_knockout_model.xl[m1_ex]
    m1ub = isnothing(m1_ex) ? 0.0 : cytbd_knockout_model.xu[m1_ex]
    change_bounds!(community_model, [env_ex]; lower = [m1lb + m2lb], upper = [m1ub + m2ub])
end

Add objective function to community model`

biomass_ids_weights = Dict(model_names .* "_BIOMASS_Ecoli_core_w_GAM" .=> 1.0)
update_community_objective!(community_model, "community_biomass", biomass_ids_weights)

Perform community FBA

d = flux_balance_analysis_dict(
    community_model,
    Tulip.Optimizer;
    modifications = [change_optimizer_attribute("IPM_IterationsLimit", 1000)],
)
d["community_biomass"] # community μ
0.21166294961691232

Notice, the growth rate is limited to the slowest organism as per the objective function

Add the enolase knockout to the community model

community_model = add_model_with_exchanges(
    community_model,
    eno_knockout_model,
    ex_rxn_mets;
    model_name = "eno_ko",
    biomass_id = "BIOMASS_Ecoli_core_w_GAM",
)

push!(model_names, "eno_ko")
biomass_ids_weights = Dict(model_names .* "_BIOMASS_Ecoli_core_w_GAM" .=> 1.0)
update_community_objective!(community_model, "community_biomass", biomass_ids_weights)

d = flux_balance_analysis_dict(
    community_model,
    Tulip.Optimizer;
    modifications = [change_optimizer_attribute("IPM_IterationsLimit", 1000)],
)
println("Community μ = ", d["community_biomass"])
Community μ = 1.6243153515245851e-15

Notice that the high communal growth rate is 0, due to the enolase knockout. The reason for this behaviour: enolase is a central reaction in glycolysis - without it the organism cannot access the lower glycolysis pathways or the TCA cycle, hence the model predicts no growth for the knockout, and hence no growth for the system since they all have to have the same growth rate.

Allow the mutants to rescue each other by sharing pyruvate

pyr_exs = model_names .* "_EX_pyr_e"
change_bounds!(community_model, pyr_exs; lower = fill(-1000.0, 3), upper = fill(1000.0, 3))

d = flux_balance_analysis_dict(
    community_model,
    Tulip.Optimizer;
    modifications = [change_optimizer_attribute("IPM_IterationsLimit", 1000)],
)
d["community_biomass"] # community μ
0.2360711079620927

Notice that the growth rate is now above 0! Nutrient sharing saved the day!

This page was generated using Literate.jl.