Steady community composition models

using COBREXA

In the community construction example, we have constructed a 2-member community of interacting E. coli with fixed abundances. Here we show how to explore the community abundances in a more dynamic way using community_composition_balance_constraints, which follows the methodology established by the SteadyCom method.

In short, we start with a similar community of knockouts, and determine what are the feasible ranges for the community abundances.

As usual, we start with the toy E. coli model and a few packages:

download_model(
    "http://bigg.ucsd.edu/static/models/e_coli_core.json",
    "e_coli_core.json",
    "7bedec10576cfe935b19218dc881f3fb14f890a1871448fc19a9b4ee15b448d8",
)

import JSONFBCModels
import HiGHS
import AbstractFBCModels.CanonicalModel as CM
import ConstraintTrees as C

[ Info: using cached `e_coli_core.json'

We remove the artificial limit of glucose intake (so that we can move it to the community level 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.0

1000.0

We create a few reaction knockouts:

knockouts = ["CYTBD", "FBA"];
ko_ecolis = Dict(begin
    m = deepcopy(ecoli)
    m.reactions[r].lower_bound = m.reactions[r].upper_bound = 0
    Symbol(r) => m
end for r in knockouts)

Dict{Symbol, AbstractFBCModels.CanonicalModel.Model} with 2 entries:
  :CYTBD => Model(Dict{String, Reaction}("ACALD"=>Reaction("Acetaldehyde dehydr…
  :FBA   => Model(Dict{String, Reaction}("ACALD"=>Reaction("Acetaldehyde dehydr…

Finding feasible community compositions

Now, let's use community_composition_variability_analysis to determine the feasible range of abundances if the community if these 2 reaction knockouts is forced to grow 0.5 g/gDW/h:

res = community_composition_variability_analysis(
    ko_ecolis,
    0.5,
    ["EX_glc__D_e" => (-10.0, 0.0)],
    optimizer = HiGHS.Optimizer,
)

ConstraintTrees.Tree{Tuple{Union{Nothing, Float64}, Union{Nothing, Float64}}} with 2 elements:
  :CYTBD => (-0.0, 0.234389)
  :FBA   => (0.765611, 1.0)

It seems like the FBA-knockout has to be present for the community to be feasible at this rate.

Finding the community composition at optimum growth

The variable community composition allows us to also scan for optimal growth (via simple interval splitting), determining the community composition near this optimum:

res = community_composition_balance_analysis(
    ko_ecolis,
    1.0,
    ["EX_glc__D_e" => (-10.0, 0.0)],
    optimizer = HiGHS.Optimizer,
    tolerance = 1e-4,
)

ConstraintTrees.Tree{Float64} with 8 elements:
  :abundances                 => ConstraintTrees.Tree{Float64}(#= 2 elements =#)
  :biomass_constraints        => ConstraintTrees.Tree{Float64}(#= 2 elements =#)
  :community                  => ConstraintTrees.Tree{Float64}(#= 2 elements =#)
  :community_balance          => ConstraintTrees.Tree{Float64}(#= 20 elements =…
  :community_exchanges        => ConstraintTrees.Tree{Float64}(#= 20 elements =…
  :diluted_constraints        => ConstraintTrees.Tree{Float64}(#= 2 elements =#)
  :growth                     => 0.703979
  :total_abundance_constraint => 1.0

The achieved growth is reported in the result tree:

res.growth

0.7039794921875

We can also have a look at the abundances at optimum. In this case, it looks like CYTBD knockout will be practically missing (within the tolerance requested above, it is technically zero):

res.abundances

ConstraintTrees.Tree{Float64} with 2 elements:
  :CYTBD => 4.90992e-5
  :FBA   => 0.999951

This page was generated using Literate.jl.