working community

docs/src/examples/06-thermodynamic-models.jl

@@ -112,10 +112,11 @@ m = build_max_min_driving_force_model(
     R = 8.31446261815324e-3, # kJ/K/mol
 )
 
-m *= :metabolite_ratio_constraints^log_ratio_constraints(
-    Dict("atp" => ("atp_c", "adp_c", 10.0), "nadh" => ("nadh_c", "nad_c", 0.13)),
-    m.log_metabolite_concentrations,
-)
+m *=
+    :metabolite_ratio_constraints^log_ratio_constraints(
+        Dict("atp" => ("atp_c", "adp_c", 10.0), "nadh" => ("nadh_c", "nad_c", 0.13)),
+        m.log_metabolite_concentrations,
+    )
 
 # solve the model
 mmdf_solution = optimized_constraints(

docs/src/examples/07-loopless-models.jl

@@ -33,10 +33,14 @@ sol = loopless_flux_balance_analysis(
     optimizer = GLPK.Optimizer,
     max_flux_bound = 1000.0, # needs to be an order of magnitude bigger, big M method heuristic
     strict_inequality_tolerance = 1.0, # heuristic from paper
-    modifications = [set_optimizer_attribute("tol_int", 1e-9)]
+    modifications = [set_optimizer_attribute("tol_int", 1e-9)],
 )
 
 sol.pseudo_gibbs_free_energy_reaction
 
-@test isapprox(mmdf_solution.max_min_driving_force, 0.8739215069684292, atol = TEST_TOLERANCE) #src
+@test isapprox(
+    mmdf_solution.max_min_driving_force,
+    0.8739215069684292,
+    atol = TEST_TOLERANCE,
+) #src
 

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/analysis/loopless_flux_balance.jl

@@ -0,0 +1 @@
+

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

src/builders/loopless.jl

@@ -35,65 +35,84 @@ function loopless_flux_balance_analysis(
     model;
     max_flux_bound = 1000.0, # needs to be an order of magnitude bigger, big M method heuristic
     strict_inequality_tolerance = 10.0, # heuristic from paper
-    modifications=[],
+    modifications = [],
     optimizer,
 )
-
-        m = fbc_model_constraints(model)
 
-        # find all internal reactions
-        internal_reactions = [(i, Symbol(rid)) for (i, rid) in enumerate(A.reactions(model)) if !is_boundary(model, rid)]
-        internal_reaction_ids = last.(internal_reactions)
-        internal_reaction_idxs = first.(internal_reactions)
+    m = fbc_model_constraints(model)
 
-        # add loopless variables: flux direction setters (binary) and pseudo gibbs free energy of reaction
-        m += :loopless_binary_variables^C.variables(keys=internal_reaction_ids, bounds=Ref(C.Binary))
-        m += :pseudo_gibbs_free_energy_reaction^C.variables(keys=internal_reaction_ids)
+    # find all internal reactions
+    internal_reactions = [
+        (i, Symbol(rid)) for
+        (i, rid) in enumerate(A.reactions(model)) if !is_boundary(model, rid)
+    ]
+    internal_reaction_ids = last.(internal_reactions)
+    internal_reaction_idxs = first.(internal_reactions)
 
-        # add -1000 * (1-a) ≤ v ≤ 1000 * a which need to be split into forward and backward components
-        m *= :loopless_reaction_directions^:backward^C.ConstraintTree(
+    # add loopless variables: flux direction setters (binary) and pseudo gibbs free energy of reaction
+    m +=
+        :loopless_binary_variables^C.variables(
+            keys = internal_reaction_ids,
+            bounds = Ref(C.Binary),
+        )
+    m += :pseudo_gibbs_free_energy_reaction^C.variables(keys = internal_reaction_ids)
+
+    # add -1000 * (1-a) ≤ v ≤ 1000 * a which need to be split into forward and backward components
+    m *=
+        :loopless_reaction_directions^:backward^C.ConstraintTree(
             rid => C.Constraint(
-                value = m.fluxes[rid].value + max_flux_bound * (1 - m.loopless_binary_variables[rid].value),
+                value = m.fluxes[rid].value +
+                        max_flux_bound * (1 - m.loopless_binary_variables[rid].value),
                 bound = (0, Inf),
             ) for rid in internal_reaction_ids
         )
-        m *= :loopless_reaction_directions^:forward^C.ConstraintTree(
+    m *=
+        :loopless_reaction_directions^:forward^C.ConstraintTree(
             rid => C.Constraint(
-                value = max_flux_bound * m.loopless_binary_variables[rid].value - m.fluxes[rid].value ,
+                value = max_flux_bound * m.loopless_binary_variables[rid].value -
+                        m.fluxes[rid].value,
                 bound = (0, Inf),
             ) for rid in internal_reaction_ids
         )
 
-        # add -1000*a + 1 * (1-a) ≤ Gibbs ≤ -1 * a + 1000 * (1 - a) which also need to be split
-        m *= :loopless_pseudo_gibbs_sign^:backward^C.ConstraintTree(
+    # add -1000*a + 1 * (1-a) ≤ Gibbs ≤ -1 * a + 1000 * (1 - a) which also need to be split
+    m *=
+        :loopless_pseudo_gibbs_sign^:backward^C.ConstraintTree(
             rid => C.Constraint(
-                value =  m.pseudo_gibbs_free_energy_reaction[rid].value + max_flux_bound * m.loopless_binary_variables[rid].value - strict_inequality_tolerance * (1 - m.loopless_binary_variables[rid].value),
+                value = m.pseudo_gibbs_free_energy_reaction[rid].value +
+                        max_flux_bound * m.loopless_binary_variables[rid].value -
+                        strict_inequality_tolerance *
+                        (1 - m.loopless_binary_variables[rid].value),
                 bound = (0, Inf),
             ) for rid in internal_reaction_ids
         )
-        m *= :loopless_pseudo_gibbs_sign^:forward^C.ConstraintTree(
+    m *=
+        :loopless_pseudo_gibbs_sign^:forward^C.ConstraintTree(
             rid => C.Constraint(
-                value = -strict_inequality_tolerance * m.loopless_binary_variables[rid].value + max_flux_bound * (1 - m.loopless_binary_variables[rid].value) - m.pseudo_gibbs_free_energy_reaction[rid].value,
+                value = -strict_inequality_tolerance *
+                        m.loopless_binary_variables[rid].value +
+                        max_flux_bound * (1 - m.loopless_binary_variables[rid].value) -
+                        m.pseudo_gibbs_free_energy_reaction[rid].value,
                 bound = (0, Inf),
             ) for rid in internal_reaction_ids
         )
 
-        # add N_int' * Gibbs = 0 where N_int = nullspace
-        N_int = nullspace(Array(A.stoichiometry(model)[:, internal_reaction_idxs])) # no sparse nullspace function
-        m *= :loopless_condition^C.ConstraintTree(
+    # add N_int' * Gibbs = 0 where N_int = nullspace
+    N_int = nullspace(Array(A.stoichiometry(model)[:, internal_reaction_idxs])) # no sparse nullspace function
+    m *=
+        :loopless_condition^C.ConstraintTree(
             Symbol(:nullspace_vector, i) => C.Constraint(
-            value = sum(col[j] * m.pseudo_gibbs_free_energy_reaction[internal_reaction_ids[j]].value for j in eachindex(col)),    
-            bound = 0,
+                value = sum(
+                    col[j] *
+                    m.pseudo_gibbs_free_energy_reaction[internal_reaction_ids[j]].value
+                    for j in eachindex(col)
+                ),
+                bound = 0,
             ) for (i, col) in enumerate(eachcol(N_int))
         )
 
-        # solve
-        optimized_constraints(
-            m;
-            objective = m.objective.value,
-            optimizer,
-            modifications,
-        )
+    # solve
+    optimized_constraints(m; objective = m.objective.value, optimizer, modifications)
 end
 
 export loopless_flux_balance_analysis

src/misc/utils.jl

@@ -8,6 +8,7 @@ reaction, otherwise return false. Checks if on boundary by inspecting the number
 of metabolites in the reaction stoichiometry. Boundary reactions have only one
 metabolite, e.g. an exchange reaction, or a sink/demand reaction.
 """
-is_boundary(model::A.AbstractFBCModel, rxn_id::String) = length(keys(A.reaction_stoichiometry(model, rxn_id))) == 1
+is_boundary(model::A.AbstractFBCModel, rxn_id::String) =
+    length(keys(A.reaction_stoichiometry(model, rxn_id))) == 1
 
 export is_boundary