Community FBA models
using COBREXAHere we construct a community FBA model of two E. coli "core" models that can interact by exchanging selected metabolites. To do this, we will need the model, which we can download if it is not already present.
download_model(
"http://bigg.ucsd.edu/static/models/e_coli_core.json",
"e_coli_core.json",
"7bedec10576cfe935b19218dc881f3fb14f890a1871448fc19a9b4ee15b448d8",
)"e_coli_core.json"Additionally to COBREXA and the model format package, we will need a solver and a few supporting packages.
import JSONFBCModels
import HiGHS
import AbstractFBCModels.CanonicalModel as CM
import ConstraintTrees as CThe core model has an artificial bound on input glucose; here we unblock that one, and we are going to add a community-global glucose intake bound later.
ecoli = load_model("e_coli_core.json", CM.Model)
ecoli.reactions["EX_glc__D_e"].lower_bound = -1000.0
ecoli.reactions["EX_glc__D_e"].upper_bound = 1000.01000.0To create a community that is actually interesting, we need some diversity. Here we simply block a different reaction in each of the community members:
ecoli1 = deepcopy(ecoli)
ecoli1.reactions["CYTBD"].lower_bound = ecoli1.reactions["CYTBD"].upper_bound = 0.0
ecoli2 = deepcopy(ecoli)
ecoli2.reactions["FBA"].lower_bound = ecoli2.reactions["FBA"].upper_bound = 0.00.0Analysing the community
To construct the community, we have to provide identifiers for the models (these will be used in the constraint tree), and corresponding models with the abundances.
my_community = Dict("bug1" => (ecoli1, 0.2), "bug2" => (ecoli2, 0.8))Dict{String, Tuple{AbstractFBCModels.CanonicalModel.Model, Float64}} with 2 entries:
"bug2" => (Model(Dict{String, Reaction}("ACALD"=>Reaction("Acetaldehyde dehyd…
"bug1" => (Model(Dict{String, Reaction}("ACALD"=>Reaction("Acetaldehyde dehyd…The community is constructed and analysed using community_flux_balance_analysis:
solution = community_flux_balance_analysis(
my_community,
["EX_glc__D_e" => (-10.0, 0.0)],
optimizer = HiGHS.Optimizer,
)ConstraintTrees.Tree{Float64} with 6 elements:
:bug1 => ConstraintTrees.Tree{Float64}(#= 5 elements =#)
:bug2 => ConstraintTrees.Tree{Float64}(#= 5 elements =#)
:community_balance => ConstraintTrees.Tree{Float64}(#= 20 elements =#)
:community_biomass => 0.523716
:community_exchanges => ConstraintTrees.Tree{Float64}(#= 20 elements =#)
:equal_growth => ConstraintTrees.Tree{Float64}(#= 1 element =#)Investigating the solution
We can now e.g. observe the differences in individual pairs of exchanges:
C.zip(
tuple,
solution.bug1.interface.exchanges,
solution.bug2.interface.exchanges,
Tuple{Float64,Float64},
)ConstraintTrees.Tree{Tuple{Float64, Float64}} with 20 elements:
:EX_ac_e => (16.5424, 0.0)
:EX_acald_e => (0.0, 0.0)
:EX_akg_e => (0.0, 0.0)
:EX_co2_e => (-0.935723, 19.8177)
:EX_etoh_e => (15.9879, 0.0)
:EX_for_e => (35.0581, 2.52782)
:EX_fru_e => (0.0, 0.0)
:EX_fum_e => (0.0, 0.0)
:EX_glc__D_e => (-20.245, -7.43875)
:EX_gln__L_e => (0.0, 0.0)
:EX_glu__L_e => (0.0, 0.0)
:EX_h2o_e => (-13.1086, 23.6327)
:EX_h_e => (62.1062, 13.0336)
:EX_lac__D_e => (0.0, 0.0)
:EX_mal__L_e => (0.0, 0.0)
:EX_nh4_e => (-2.85572, -2.85572)
:EX_o2_e => (0.0, -20.4762)
:EX_pi_e => (-1.92659, -1.92659)
:EX_pyr_e => (0.0, 0.0)
:EX_succ_e => (0.0, 0.0)...or use screen to efficiently find out which composition is best:
screen(0.0:0.1:1.0) do ratio2
ratio1 = 1 - ratio2
res = community_flux_balance_analysis(
[("bug1" => (ecoli1, ratio1)), ("bug2" => (ecoli2, ratio2))],
["EX_glc__D_e" => (-10.0, 0.0)],
interface = :sbo, # usually more reproducible
optimizer = HiGHS.Optimizer,
)
(ratio1, ratio2) => (isnothing(res) ? nothing : res.community_biomass)
end11-element Vector{Pair{Tuple{Float64, Float64}, Float64}}:
(1.0, 0.0) => 0.2116629497353106
(0.9, 0.1) => 0.2342460462980283
(0.8, 0.2) => 0.2597197180964099
(0.7, 0.3) => 0.28867689024498494
(0.6, 0.4) => 0.3218845533955697
(0.5, 0.5) => 0.3603526987796251
(0.4, 0.6) => 0.40543883161610117
(0.30000000000000004, 0.7) => 0.459011502020198
(0.19999999999999996, 0.8) => 0.5237157737585175
(0.09999999999999998, 0.9) => 0.6034233598466012
(0.0, 1.0) => 0.7040369478590236(The result seem like the bug1 is eventually going to be completely out-grown by the other one.)
Note: interfaces of constraint systems
Internally, the community is connected via interfaces, which are small constraint trees (typically with no bounds attached) that describe parts of the constraint system that can be easily attached to other parts.
The best kind of interface to choose generally differs from model to model. COBREXA gives a few "default" choices that cover a good part of sensible metabolic modeling. For example, if the model contains SBO annotations, we can ask for the interface created using the annotated reactions:
flux_balance_constraints(ecoli, interface = :sbo).interfaceConstraintTrees.Tree{ConstraintTrees.Constraint} with 2 elements:
:biomass => ConstraintTrees.Tree{ConstraintTrees.Constraint}(#= 1 element =…
:exchanges => ConstraintTrees.Tree{ConstraintTrees.Constraint}(#= 20 elements…If there are no annotations, we can still at least detect the boundary reactions and make an interface out of them:
flux_balance_constraints(ecoli, interface = :boundary).interfaceConstraintTrees.Tree{ConstraintTrees.Constraint} with 1 element:
:boundary => ConstraintTrees.Tree{ConstraintTrees.Constraint}(#= 20 elements …The default kind of interface in community_flux_balance_analysis is :identifier_prefixes, which relies on usual prefixes of reaction names (such as EX_ for exchanges).
Even if all of these methods fail, a suitable interface yourself can be produced manually. (Additionally, we can do useful stuff, such as removing the unnecessary bounds from the exchange descriptions.)
custom_model = flux_balance_constraints(ecoli)
custom_model *= remove_bounds(
:interface^C.ConstraintTree(
:biomass => custom_model.fluxes.BIOMASS_Ecoli_core_w_GAM,
:exchanges => C.ConstraintTree(
k => v for (k, v) in custom_model.fluxes if startswith(string(k), "EX_")
),
),
)
custom_model.interface.exchangesConstraintTrees.Tree{ConstraintTrees.Constraint} with 20 elements:
:EX_ac_e => ConstraintTrees.Constraint(ConstraintTrees.LinearValue(#= ...…
:EX_acald_e => ConstraintTrees.Constraint(ConstraintTrees.LinearValue(#= ...…
:EX_akg_e => ConstraintTrees.Constraint(ConstraintTrees.LinearValue(#= ...…
:EX_co2_e => ConstraintTrees.Constraint(ConstraintTrees.LinearValue(#= ...…
:EX_etoh_e => ConstraintTrees.Constraint(ConstraintTrees.LinearValue(#= ...…
:EX_for_e => ConstraintTrees.Constraint(ConstraintTrees.LinearValue(#= ...…
:EX_fru_e => ConstraintTrees.Constraint(ConstraintTrees.LinearValue(#= ...…
:EX_fum_e => ConstraintTrees.Constraint(ConstraintTrees.LinearValue(#= ...…
:EX_glc__D_e => ConstraintTrees.Constraint(ConstraintTrees.LinearValue(#= ...…
:EX_gln__L_e => ConstraintTrees.Constraint(ConstraintTrees.LinearValue(#= ...…
:EX_glu__L_e => ConstraintTrees.Constraint(ConstraintTrees.LinearValue(#= ...…
:EX_h2o_e => ConstraintTrees.Constraint(ConstraintTrees.LinearValue(#= ...…
:EX_h_e => ConstraintTrees.Constraint(ConstraintTrees.LinearValue(#= ...…
:EX_lac__D_e => ConstraintTrees.Constraint(ConstraintTrees.LinearValue(#= ...…
:EX_mal__L_e => ConstraintTrees.Constraint(ConstraintTrees.LinearValue(#= ...…
:EX_nh4_e => ConstraintTrees.Constraint(ConstraintTrees.LinearValue(#= ...…
:EX_o2_e => ConstraintTrees.Constraint(ConstraintTrees.LinearValue(#= ...…
:EX_pi_e => ConstraintTrees.Constraint(ConstraintTrees.LinearValue(#= ...…
:EX_pyr_e => ConstraintTrees.Constraint(ConstraintTrees.LinearValue(#= ...…
:EX_succ_e => ConstraintTrees.Constraint(ConstraintTrees.LinearValue(#= ...…Connecting the community constraints manually
To connect such interfaces into a community model, simply use function interface_constraints (which is how community_flux_balance_analysis constructs the community model internally via community_flux_balance_constraints). The assembly might look roughly as follows:
custom_community = interface_constraints(
"bug1" => (
custom_model * :handicap^C.Constraint(custom_model.fluxes.CYTBD.value, 0),
0.2,
),
"bug2" =>
(custom_model * :handicap^C.Constraint(custom_model.fluxes.FBA.value, 0), 0.8),
bound = r -> r == (:exchanges, :EX_glc__D_e) ? C.Between(-10, 0) : nothing,
)
custom_community.interface.exchangesConstraintTrees.Tree{ConstraintTrees.Constraint} with 20 elements:
:EX_ac_e => ConstraintTrees.Constraint(ConstraintTrees.LinearValue(#= ...…
:EX_acald_e => ConstraintTrees.Constraint(ConstraintTrees.LinearValue(#= ...…
:EX_akg_e => ConstraintTrees.Constraint(ConstraintTrees.LinearValue(#= ...…
:EX_co2_e => ConstraintTrees.Constraint(ConstraintTrees.LinearValue(#= ...…
:EX_etoh_e => ConstraintTrees.Constraint(ConstraintTrees.LinearValue(#= ...…
:EX_for_e => ConstraintTrees.Constraint(ConstraintTrees.LinearValue(#= ...…
:EX_fru_e => ConstraintTrees.Constraint(ConstraintTrees.LinearValue(#= ...…
:EX_fum_e => ConstraintTrees.Constraint(ConstraintTrees.LinearValue(#= ...…
:EX_glc__D_e => ConstraintTrees.Constraint(ConstraintTrees.LinearValue(#= ...…
:EX_gln__L_e => ConstraintTrees.Constraint(ConstraintTrees.LinearValue(#= ...…
:EX_glu__L_e => ConstraintTrees.Constraint(ConstraintTrees.LinearValue(#= ...…
:EX_h2o_e => ConstraintTrees.Constraint(ConstraintTrees.LinearValue(#= ...…
:EX_h_e => ConstraintTrees.Constraint(ConstraintTrees.LinearValue(#= ...…
:EX_lac__D_e => ConstraintTrees.Constraint(ConstraintTrees.LinearValue(#= ...…
:EX_mal__L_e => ConstraintTrees.Constraint(ConstraintTrees.LinearValue(#= ...…
:EX_nh4_e => ConstraintTrees.Constraint(ConstraintTrees.LinearValue(#= ...…
:EX_o2_e => ConstraintTrees.Constraint(ConstraintTrees.LinearValue(#= ...…
:EX_pi_e => ConstraintTrees.Constraint(ConstraintTrees.LinearValue(#= ...…
:EX_pyr_e => ConstraintTrees.Constraint(ConstraintTrees.LinearValue(#= ...…
:EX_succ_e => ConstraintTrees.Constraint(ConstraintTrees.LinearValue(#= ...…For the model to work properly, we would need to add several other things, mainly the equal growth constraints (possibly via all_equal_constraints). community_flux_balance_constraints add these automatically, so we can equivalently just supply the constraint trees, and re-use the rest of the implementation:
custom_community = community_flux_balance_constraints(
[
"bug1" => (
custom_model * :handicap^C.Constraint(custom_model.fluxes.CYTBD.value, 0),
0.2,
),
"bug2" => (
custom_model * :handicap^C.Constraint(custom_model.fluxes.FBA.value, 0),
0.8,
),
],
["EX_glc__D_e" => (-10.0, 0.0)],
)ConstraintTrees.Tree{ConstraintTrees.Constraint} with 6 elements:
:bug1 => ConstraintTrees.Tree{ConstraintTrees.Constraint}(#= 6…
:bug2 => ConstraintTrees.Tree{ConstraintTrees.Constraint}(#= 6…
:community_balance => ConstraintTrees.Tree{ConstraintTrees.Constraint}(#= 2…
:community_biomass => ConstraintTrees.Constraint(ConstraintTrees.LinearValu…
:community_exchanges => ConstraintTrees.Tree{ConstraintTrees.Constraint}(#= 2…
:equal_growth => ConstraintTrees.Tree{ConstraintTrees.Constraint}(#= 1…This can be solved with the usual means, reaching the same result as above:
custom_solution = optimized_values(
custom_community,
objective = custom_community.community_biomass.value,
output = custom_community.community_biomass,
optimizer = HiGHS.Optimizer,
)0.5237157737585176This page was generated using Literate.jl.