improve docs and docstrings

Author: stelmo

docs/src/examples/08-community-models.jl

@@ -1,9 +1,10 @@
-# # Enzyme constrained models
+# # Community FBA models
 
 using COBREXA
 
-# Here we will construct an enzyme constrained variant of the *E. coli* "core"
-# model. We will need the model, which we can download if it is not already present.
+# Here we will 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.
 
 import Downloads: download
 
@@ -15,59 +16,124 @@ import Downloads: download
 
 import JSONFBCModels
 import Tulip
+import AbstractFBCModels as A
+import ConstraintTrees as C
 
 model = load_model("e_coli_core.json")
 
-# Enzyme constrained models require parameters that are usually not used by
-# conventional constraint based models. These include reaction specific turnover
-# numbers, molar masses of enzymes, and capacity bounds.
+# Community models work by joining its members together through their exchange
+# reactions, weighted by the abundance of each microbe. These exchange reactions
+# are then linked to an environmental exchange. For more theoretical details,
+# see "Gottstein, et al, 2016, Constraint-based stoichiometric modelling from
+# single organisms to microbial communities, Journal of the Royal Society
+# Interface". 
 
-import AbstractFBCModels as A
+# ## Building a community of two *E. coli*s 
+
+# Here we will construct a simple community of two interacting microbes. To do
+# this, we need to import the models. We import the models are ConstraintTrees,
+# because it is easier to build the model explicitly than rely on an opaque
+# one-shot function.
 
-m1 = fbc_model_constraints(model)
-m2 = fbc_model_constraints(model)
+ecoli1 = fbc_model_constraints(model)
+ecoli2 = fbc_model_constraints(model)
 
+# Since the models are joined through their individual exchange reactions to an
+# environmental exchange reactionq, we need to identify all possible exchange
+# reactions in the community. Since the models are the same, this is
+# straightforward here. Additionally, we need to specify the upper and lower
+# bounds of these environmental exchange reactions.
 lbs, ubs = A.bounds(model)
-env_ex_rxns = Dict(rid => (lbs[i], ubs[i]) for (i, rid) in enumerate(A.reactions(model)) if startswith(rid, "EX_"))
 
+env_ex_rxns = Dict(
+    rid => (lbs[i], ubs[i]) for
+    (i, rid) in enumerate(A.reactions(model)) if startswith(rid, "EX_")
+)
+
+# Now we simply create an blank model that only includes environmental exchange reactions.
 
 m = build_community_environment(env_ex_rxns)
-m += :bug1^m1
-m += :bug2^m2
 
+# Next we join each member microbe to the model.
+m += :bug1^ecoli1
+m += :bug2^ecoli2
+
+# We also need to specify the abundances of each member, as this weights the
+# flux of each metabolite each member microbe can share with other members or
+# the environment.
 member_abundances = [(:bug1, 0.2), (:bug2, 0.8)]
-m *=
-    :environmental_exchange_balances^link_environmental_exchanges(
-        m,
-        [(:bug1, 0.2), (:bug2, 0.8)],
-    )
 
+m *= :environmental_exchange_balances^link_environmental_exchanges(m, member_abundances)
+
+# Finally, the most sensible community FBA simulation involves assuming the
+# growth rate of the models is the same. In this case, we simply set the growth
+# rate flux of each member to be the same.
 m *=
-    :equal_growth_rate_constraint^equal_growth_rate_constraints(
-        [(:bug1, m.bug1.fluxes.:BIOMASS_Ecoli_core_w_GAM.value), (:bug2, m.bug2.fluxes.:BIOMASS_Ecoli_core_w_GAM.value)]
-    )
+    :equal_growth_rate_constraint^equal_growth_rate_constraints([
+        (:bug1, m.bug1.fluxes.:BIOMASS_Ecoli_core_w_GAM.value),
+        (:bug2, m.bug2.fluxes.:BIOMASS_Ecoli_core_w_GAM.value),
+    ])
+
+# Since each growth rate is the same, we can pick any of the growth rates as the
+# objective for the simulation.
+m *= :objective^C.Constraint(m.bug1.fluxes.:BIOMASS_Ecoli_core_w_GAM.value)
 
+# Since the models are usually used in a mono-culture context, the glucose input
+# for each individual member is limited. We need to undo this limitation, and
+# rather rely on the constrained environmental exchange reaction (and the bounds
+# we set for it earlier).
 m.bug1.fluxes.EX_glc__D_e.bound = (-1000.0, 1000.0)
 m.bug2.fluxes.EX_glc__D_e.bound = (-1000.0, 1000.0)
-m.bug1.fluxes.CYTBD.bound = (-10.0, 10.0) # respiration limited
 
-m *= :objective^C.Constraint(m.bug1.fluxes.:BIOMASS_Ecoli_core_w_GAM.value)
-
-using Gurobi
+# We can also be interesting, and limit respiration in one of the members, to
+# see what effect this has on the community.
+m.bug1.fluxes.CYTBD.bound = (-10.0, 10.0)
 
-sol = optimized_constraints(m; objective = m.objective.value, optimizer=Gurobi.Optimizer)
+# Finally, we can simulate the system!
+sol = optimized_constraints(
+    m;
+    objective = m.objective.value,
+    optimizer = Tulip.Optimizer,
+    modifications = [set_optimizer_attribute("IPM_IterationsLimit", 1000)],
+)
 
-Dict(k => v for (k, v) in sol.bug1.fluxes if startswith(string(k), "EX_"))
-Dict(k => v for (k, v) in sol.bug2.fluxes if startswith(string(k), "EX_"))
+@test isapprox(sol.:objective, 0.66686196344, atol = TEST_TOLERANCE) #src
 
-# exchange cytosolic metabolites
+# At the moment the members cannot really exchange any metabolites. We can
+# change this by changing their individual exchange bounds. 
 mets = [:EX_akg_e, :EX_succ_e, :EX_pyr_e, :EX_acald_e, :EX_fum_e, :EX_mal__L_e]
 for met in mets
     m.bug1.fluxes[met].bound = (-1000.0, 1000.0)
     m.bug2.fluxes[met].bound = (-1000.0, 1000.0)
 end
 
-sol = optimized_constraints(m; objective = m.objective.value, optimizer=Tulip.Optimizer, modifications=[set_optimizer_attribute("IPM_IterationsLimit", 100000)])
-
-Dict(k => v for (k, v) in sol.bug1.fluxes if startswith(string(k), "EX_"))
-Dict(k => v for (k, v) in sol.bug2.fluxes if startswith(string(k), "EX_"))
+sol = optimized_constraints(
+    m;
+    objective = m.objective.value,
+    optimizer = Tulip.Optimizer,
+    modifications = [set_optimizer_attribute("IPM_IterationsLimit", 1000)],
+)
+
+
+# We can see that by allowing the microbes to share metabolites, the growth rate
+# of the system as a whole increased! We can inspect the individual exchanges to
+# see which metabolites are being shared (pyruvate in this case).
+bug1_ex_fluxes = Dict(k => v for (k, v) in sol.bug1.fluxes if startswith(string(k), "EX_"))
+bug2_ex_fluxes = Dict(k => v for (k, v) in sol.bug2.fluxes if startswith(string(k), "EX_"))
+
+#!!! warning "Flux units"
+# The unit of the environmental exchange reactions (mmol/gDW_total_biomass/h) is
+# different to the unit of the individual species fluxes
+# (mmol/gDW_species_biomass/h). This is because the mass balance needs to take
+# into account the abundance of each species for the simulation to make sense.
+# In this specific case, look at the flux of pyruvate (EX_pyr_e). There is no
+# environmental exchange flux, so the two microbes share the metabolite.
+# However,  `bug1_ex_fluxes[:EX_pyr_e] != bug2_ex_fluxes[:EX_pyr_e]`, but rather
+# `abundance_bug1 * bug1_ex_fluxes[:EX_pyr_e] != abundance_bug2 *
+# bug2_ex_fluxes[:EX_pyr_e]`. Take care of this when comparing fluxes!
+
+@test isapprox(
+    abs(0.2 * bug1_ex_fluxes[:EX_pyr_e] + 0.8 * bug2_ex_fluxes[:EX_pyr_e]),
+    0.0,
+    atol = TEST_TOLERANCE,
+) #src

src/builders/communities.jl

@@ -11,6 +11,11 @@ end
 
 export environment_exchange_variables
 
+"""
+$(TYPEDSIGNATURES)
+
+Helper function to build a "blank" community model with only environmental exchange reactions.
+"""
 function build_community_environment(env_ex_rxns = Dict{String,Tuple{Float64,Float64}}())
     C.ConstraintTree(
         :environmental_exchange_reactions => environment_exchange_variables(env_ex_rxns),
@@ -19,6 +24,12 @@ end
 
 export build_community_environment
 
+"""
+$(TYPEDSIGNATURES)
+
+Helper function to link species specific exchange reactions to the environmental
+exchange reactions by weighting them with their abundances.
+"""
 function link_environmental_exchanges(
     m::C.ConstraintTree,
     member_abundances::Vector{Tuple{Symbol,Float64}};
@@ -40,7 +51,14 @@ end
 
 export link_environmental_exchanges
 
-function equal_growth_rate_constraints(member_biomasses::Vector{Tuple{Symbol,C.LinearValue}})
+"""
+$(TYPEDSIGNATURES)
+
+Helper function to set each species growth rate equal to each other.
+"""
+function equal_growth_rate_constraints(
+    member_biomasses::Vector{Tuple{Symbol,C.LinearValue}},
+)
     C.ConstraintTree(
         Symbol(bid1, :_, bid2) => C.Constraint(value = bval1 - bval2, bound = 0.0) for
         ((bid1, bval1), (bid2, bval2)) in