Enzyme-constrained communities
This example demonstrates a simple way to create enzyme-constrained community models. We use the genome-scale iML1515 E. coli model to simulate a community of 2 auxotrophic organisms (both of which exhibit no growth in isolation due to the lack of the amino acid), and observe how they can save each other by supplying themselves amino-acids. Such analysis easily extends to more auxotrophes and even communities of several different species.
The simulations are, very roughly, replicating the logic of the experimental work by Mee, Michael T., et al. "Syntrophic exchange in synthetic microbial communities." Proceedings of the National Academy of Sciences 111.20 (2014).
As usual, we start by loading packages and downloading models:
using COBREXA
import AbstractFBCModels as A
import JSONFBCModels
import ConstraintTrees as C
import HiGHS
download_model(
"http://bigg.ucsd.edu/static/models/iML1515.json",
"iML1515.json",
"b0f9199f048779bb08a14dfa6c09ec56d35b8750d2f99681980d0f098355fbf5",
)"iML1515.json"Collecting data and parameters
Enzyme-constrained models require parameters for protein molar masses and reaction turnover numbers (kcats). COBREXA supplies prepared example data for the iML1515 model; in this section we summarize the loading of the data into Julia structures from the used format. Other formats will work just as well.
The loading is hidden by default for brevity:
Loading the model parameters
import CSV
data_dir = joinpath(dirname(pathof(COBREXA)), "..", "docs", "src", "examples", "data");
e_coli_gp_mass = Dict{String,Float64}(
x.gene_product => x.mass for
x in CSV.File(joinpath(data_dir, "e_coli_gp_mass.tsv"), delim = '\t')
);
kcat_scale = 3600 / 1e3;
e_coli_rxn_kcat_isozyme = Dict{String,Isozyme}(
x.reaction => Isozyme(
kcat_forward = x.kcat * kcat_scale,
kcat_reverse = x.kcat * kcat_scale,
gene_product_stoichiometry = Dict(),
) for x in CSV.File(joinpath(data_dir, "e_coli_reaction_kcat.tsv"), delim = '\t')
);
e_coli_rxn_isozymes = Dict{String,Dict{String,Isozyme}}();
for x in CSV.File(joinpath(data_dir, "e_coli_isozyme_gp.tsv"), delim = '\t')
haskey(e_coli_rxn_kcat_isozyme, x.reaction) || continue
rxn = get!(e_coli_rxn_isozymes, x.reaction, Dict{String,Isozyme}())
iso = get!(rxn, x.isozyme, deepcopy(e_coli_rxn_kcat_isozyme[x.reaction]))
iso.gene_product_stoichiometry[x.gene_product] = x.stoichiometry
end;In the end, we have gene product weight data (just like in the enzyme-constrained model example):
e_coli_gp_massDict{String, Float64} with 1517 entries:
"b1329" => 59.9
"b3236" => 32.337
"b0688" => 58.361
"b2052" => 36.141
"b0832" => 33.238
"b0586" => 141.991
"b2245" => 28.916
"b1759" => 15.046
"b3772" => 56.195
"b1692" => 54.58
"b3850" => 21.226
"b1006" => 45.557
"b3428" => 93.173
"b0635" => 70.857
"b1387" => 73.003
"b1588" => 89.987
"b2541" => 28.5
"b4478" => 42.523
"b2143" => 31.54
⋮ => ⋮... as well as isozyme data with kcats:
e_coli_rxn_isozymesDict{String, Dict{String, Isozyme}} with 2266 entries:
"PACOAT" => Dict("iso1"=>Isozyme(Dict("b1396"=>4.0), 33.9551, 33.9551))
"Zn2tex" => Dict("iso3"=>Isozyme(Dict("b1377"=>3.0), 648.0, 648.0), "i…
"GUI1" => Dict("iso1"=>Isozyme(Dict("b3092"=>1.0), 21.4406, 21.4406))
"DXYLK" => Dict("iso1"=>Isozyme(Dict("b3564"=>2.0), 28.7528, 28.7528))
"CBL1tonex" => Dict("iso1"=>Isozyme(Dict("b3005"=>1.0, "b3006"=>4.0, "b39…
"FE3DCITtonex" => Dict("iso1"=>Isozyme(Dict("b3005"=>3.0, "b3006"=>6.0, "b42…
"FACOAL180t2pp" => Dict("iso1"=>Isozyme(Dict("b1701"=>2.0), 28.4896, 28.4896)…
"METSOXR1" => Dict("iso3"=>Isozyme(Dict("b3781"=>1.0, "b3551"=>1.0), 115…
"LIPOtex" => Dict("iso3"=>Isozyme(Dict("b1377"=>3.0), 648.0, 648.0), "i…
"NTD11" => Dict("iso1"=>Isozyme(Dict("b2744"=>1.0), 24.6906, 24.6906)…
"GLUNpp" => Dict("iso1"=>Isozyme(Dict("b2957"=>4.0), 18.2572, 18.2572))
"ORNDC" => Dict("iso1"=>Isozyme(Dict("b0693"=>1.0), 130.25, 130.25), …
"ALAGLUE" => Dict("iso1"=>Isozyme(Dict("b1325"=>1.0), 24.1219, 24.1219))
"UAGCVT" => Dict("iso1"=>Isozyme(Dict("b3189"=>1.0), 11.88, 11.88))
"I2FE2ST" => Dict("iso1"=>Isozyme(Dict("b2529"=>1.0, "b2528"=>1.0), 11.…
"6D6SPA" => Dict("iso1"=>Isozyme(Dict("b3881"=>1.0), 149.289, 149.289))
"BMOCOS" => Dict("iso1"=>Isozyme(Dict("b0827"=>2.0), 0.072, 0.072))
"GLCURt2rpp" => Dict("iso1"=>Isozyme(Dict("b3093"=>1.0), 648.0, 648.0), "i…
"PGSA141" => Dict("iso1"=>Isozyme(Dict("b1912"=>1.0), 40.7639, 40.7639))
⋮ => ⋮Model assembly
For simplicity, we will work with the "canonical" Julia-structured view of the iML1515:
wt_model = load_model("iML1515.json", A.CanonicalModel.Model)AbstractFBCModels.CanonicalModel.Model(
reactions = Dict{String, AbstractFBCModels.CanonicalModel.Reaction}("PACOAT" …
metabolites = Dict{String, AbstractFBCModels.CanonicalModel.Metabolite}("hphh…
genes = Dict{String, AbstractFBCModels.CanonicalModel.Gene}("b1329" => Abstra…
couplings = Dict{String, AbstractFBCModels.CanonicalModel.Coupling}(),
)
As the usual quirk, we loosen the lower bound on glucose intake that is required for plain FBA:
wt_model.reactions["EX_glc__D_e"].lower_bound = -1000.0;Additionally we allow the models isoleucine and methionine uptake:
wt_model.reactions["EX_ile__L_e"].lower_bound = -1000.0;
wt_model.reactions["EX_met__L_e"].lower_bound = -1000.0;...and for good manners, we also remove the biomass annotation from the biomass reaction that we are not interested in:
wt_model.reactions["BIOMASS_Ec_iML1515_WT_75p37M"].annotations["sbo"] = [];Let's create these two knockouts– one incapable of producing isoleucine:
ile_model = deepcopy(wt_model)
delete!(ile_model.reactions, "THRD_L");...and another one without the reaction that is required for producing methionine:
met_model = deepcopy(wt_model)
delete!(met_model.reactions, "HSST");For brevity, let's make a shortcut that creates enzyme-constrained FBA system from the model together with a proper interface for community building:
ecfba_constraints(m, capacity) = enzyme_constrained_flux_balance_constraints(
m,
reaction_isozymes = e_coli_rxn_isozymes,
gene_product_molar_masses = e_coli_gp_mass,
interface = :sbo;
capacity,
)ecfba_constraints (generic function with 1 method)We can now create the community by creating each model's constraint tree with an interface with enzyme_constrained_flux_balance_constraints, and connecting them via community_flux_balance_constraints. We have to pick the model abundances for cFBA, so we pick 1:1 abundance ratio. We also have to pick the capacities for the enzyme-constrained models (these will be properly diluted by the community FBA formulation), and specify that the community is not allowed to exchange either of our two selected amino acids externally (the individual models might cheat the auxotrophe community setting by consuming these).
community_constraints = community_flux_balance_constraints(
[
"ile_ko" => (ecfba_constraints(ile_model, 100.0), 0.5),
"met_ko" => (ecfba_constraints(met_model, 100.0), 0.5),
],
["EX_ile__L_e" => 0.0, "EX_met__L_e" => 0.0],
)ConstraintTrees.ConstraintTree with 6 elements:
:community_balance => ConstraintTrees.ConstraintTree(#= 331 elements =#)
:community_biomass => ConstraintTrees.Constraint(ConstraintTrees.LinearValu…
:community_exchanges => ConstraintTrees.ConstraintTree(#= 331 elements =#)
:equal_growth => ConstraintTrees.ConstraintTree(#= 1 element =#)
:ile_ko => ConstraintTrees.ConstraintTree(#= 13 elements =#)
:met_ko => ConstraintTrees.ConstraintTree(#= 13 elements =#)Simulating the community
Since the community constraints created above form a completely normal optimization problem, we can optimize them as usual via optimized_values; picking the community_biomass value as an objective:
res = optimized_values(
community_constraints,
objective = community_constraints.community_biomass.value,
optimizer = HiGHS.Optimizer,
)ConstraintTrees.Tree{Float64} with 6 elements:
:community_balance => ConstraintTrees.Tree{Float64}(#= 331 elements =#)
:community_biomass => 0.230488
:community_exchanges => ConstraintTrees.Tree{Float64}(#= 331 elements =#)
:equal_growth => ConstraintTrees.Tree{Float64}(#= 1 element =#)
:ile_ko => ConstraintTrees.Tree{Float64}(#= 13 elements =#)
:met_ko => ConstraintTrees.Tree{Float64}(#= 13 elements =#)We can observe that the community indeed grows, although not as quickly as the WT model normally would:
res.community_biomass0.2304879809596633One may also observe the "global" community exchanges:
sort(collect(res.community_exchanges), by = last)331-element Vector{Pair{Symbol, Union{Float64, ConstraintTrees.Tree{Float64}}}}:
:EX_glc__D_e => -11.45938376478452
:EX_o2_e => -10.09733248430452
:EX_nh4_e => -2.4892506028859476
:EX_pi_e => -0.22233008936155146
:EX_so4_e => -0.05804148336526241
:EX_k_e => -0.04498964046745957
:EX_mg2_e => -0.0019994832348250793
:EX_fe2_e => -0.0019022173068601014
:EX_fe3_e => -0.001799650155333051
:EX_ca2_e => -0.0011996899408950475
⋮
:EX_xylu__L_e => 0.0
:EX_meoh_e => 4.6097596191932663e-7
:EX_dxylnt_e => 0.00010279763950800903
:EX_co2_e => 0.45749610377308664
:EX_for_e => 0.9776511883414105
:EX_pyr_e => 1.4998199654895021
:EX_5dglcn_e => 8.893365393706606
:EX_h_e => 13.490554529094084
:EX_h2o_e => 16.93633191576877Appropriately, we can check that the individual community members exchange the expected amino acids (the individual values are scaled to the individual members' biomasses; in the community view these values would be halved by the 0.5 abundances):
[
res.met_ko.fluxes.EX_ile__L_e res.met_ko.fluxes.EX_met__L_e
res.ile_ko.fluxes.EX_ile__L_e res.ile_ko.fluxes.EX_met__L_e
]2×2 Matrix{Float64}:
0.0669634 -0.0354746
-0.0669634 0.0354746We can see that isoleucine is indeed moving into the isoleucine knockout, and methionine into the methionine knockout. (The signs follow the usual exchange convention where negative values mean uptake and positive values mean secretion.)
Finally, one might be interested in finding the optimal community composition for the auxotrophes. The example on community building describes a common way to find the optimal abundance ratios via screening.
This page was generated using Literate.jl.