working community

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

@@ -23,3 +23,51 @@ model = load_model("e_coli_core.json")
 # numbers, molar masses of enzymes, and capacity bounds.
 
 import AbstractFBCModels as A
+
+m1 = fbc_model_constraints(model)
+m2 = fbc_model_constraints(model)
+
+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_"))
+
+
+m = build_community_environment(env_ex_rxns)
+m += :bug1^m1
+m += :bug2^m2
+
+member_abundances = [(:bug1, 0.2), (:bug2, 0.8)]
+m *=
+    :environmental_exchange_balances^link_environmental_exchanges(
+        m,
+        [(:bug1, 0.2), (:bug2, 0.8)],
+    )
+
+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)]
+    )
+
+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
+
+sol = optimized_constraints(m; objective = m.objective.value, optimizer=Gurobi.Optimizer)
+
+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_"))
+
+# exchange cytosolic metabolites
+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_"))

src/COBREXA.jl

@@ -37,6 +37,7 @@ include("builders/objectives.jl")
 include("builders/enzymes.jl")
 include("builders/thermodynamic.jl")
 include("builders/loopless.jl")
+include("builders/communities.jl")
 
 include("analysis/modifications.jl")
 include("analysis/flux_balance.jl")

src/builders/communities.jl

@@ -0,0 +1,51 @@
+"""
+$(TYPEDSIGNATURES)
+
+Helper function to create environmental exchange rections.
+"""
+function environment_exchange_variables(env_ex_rxns = Dict{String,Tuple{Float64,Float64}}())
+    rids = collect(keys(env_ex_rxns))
+    lbs_ubs = collect(values(env_ex_rxns))
+    C.variables(; keys = Symbol.(rids), bounds = lbs_ubs)
+end
+
+export environment_exchange_variables
+
+function build_community_environment(env_ex_rxns = Dict{String,Tuple{Float64,Float64}}())
+    C.ConstraintTree(
+        :environmental_exchange_reactions => environment_exchange_variables(env_ex_rxns),
+    )
+end
+
+export build_community_environment
+
+function link_environmental_exchanges(
+    m::C.ConstraintTree,
+    member_abundances::Vector{Tuple{Symbol,Float64}};
+    on = m.:environmental_exchange_reactions,
+    member_fluxes_id = :fluxes,
+)
+    C.ConstraintTree(
+        rid => C.Constraint(
+            value = -rxn.value + sum(
+                abundance * m[member][member_fluxes_id][rid].value for
+                (member, abundance) in member_abundances if
+                haskey(m[member][member_fluxes_id], rid);
+                init = zero(C.LinearValue),
+            ),
+            bound = 0.0,
+        ) for (rid, rxn) in on
+    )
+end
+
+export link_environmental_exchanges
+
+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
+        zip(member_biomasses[1:end-1], member_biomasses[2:end])
+    )
+end
+
+export equal_growth_rate_constraints