community_dfba

Dynamic community Flux Balance Analysis on top of cobrapy — a clean reimplementation of the dcFBA algorithms that runs on the free GLPK solver (and Gurobi/CPLEX when available), instead of cbmpy + CPLEX.

It provides the same three methods as dcFBA:

Method	Merged matrix?	Objective	Biology
DynamicJointFBA	yes (one shared matrix)	maximize total community biomass	cooperative
DynamicParallelFBA	no (each solved separately)	each organism maximizes its own biomass	competitive / selfish
EndPointFBA	yes, duplicated over N time steps	maximize community biomass at the horizon	cooperative, "clairvoyant"

Organisms interact only through a shared extracellular metabolite pool (compartment e): private (cytosolic) species/reactions are suffixed per organism (glc__D_c → glc__D_c__ecoli), extracellular species are shared, so competition and cross-feeding emerge from the merged stoichiometry.

Why this exists

cbmpy 0.8.8.1's free GLPK backend is broken (glp_set_mat_row … column index out of range, SIGABRT), and cbmpy has no interface to Gurobi/Mosek. cobrapy (via optlang) drives GLPK and Gurobi cleanly, so the dcFBA algorithms were ported onto it.

Install

Requirements are already satisfied on this machine (cobra, numpy; matplotlib

• pandas for plots). To install the package itself:

cd community_dfba
pip install --user .            # or: PYTHONPATH=. python your_script.py

Solvers

• GLPK (default, free) — works out of the box.
• Gurobi — gurobipy is installed and optlang detects it, but a WLS (Web License Service) license validates online, and outbound HTTPS to token.gurobi.com is currently blocked here (HTTP 451). Fix the network egress (whitelist *.gurobi.com / set HTTPS_PROXY) or use a node-locked gurobi.lic, then pass solver="gurobi". Only EndPointFBA QP smoothing strictly needs Gurobi/CPLEX (GLPK cannot do QP).

Note: because a broken Gurobi default makes cobra crash on model creation, importing community_dfba sets cobra's default solver to GLPK. Pass solver="gurobi" explicitly to opt in.

Quick start

import cobra; cobra.Configuration().solver = "glpk"
from cobra.io import load_model, read_sbml_model
import community_dfba as cdf

# --- load your organisms (any SBML/.xml GSMM) ---
a = read_sbml_model("bacterium_A.xml"); a.id = "A"
b = read_sbml_model("bacterium_B.xml"); b.id = "B"

# cap per-gDW uptake on the TRANSPORT reaction (not the exchange - see notes)
a.reactions.get_by_id("GLCpts").upper_bound = 10

# --- Joint FBA ---
cm = cdf.build_community_model([a, b],
        biomass_ids=["BIOMASS_A", "BIOMASS_B"], organism_ids=["A", "B"])
dj = cdf.DynamicJointFBA(cm,
        initial_biomasses={"A": 0.01, "B": 0.01},
        initial_concentrations={"glc__D_e": 15.0}, solver="glpk")
dj.simulate(dt=0.1, n=1000)

dj.get_biomasses()          # {"A":[...], "B":[...]}  biomass vs time
dj.get_metabolites()        # {"glc__D_e":[...], ...} pool vs time
dj.times                    # time points
dj.flux_series("GLCpts__A") # a reaction's flux vs time (suffixed id)
dj.relative_abundance()
dj.community_growth_rate()

DynamicParallelFBA takes the raw model list (not a merged model); EndPointFBA takes the merged model + number of steps n and is solved once with .simulate(). See examples/compare_methods.py for all three side-by-side and examples/plot_comparison.py for the figure.

Monod (Michaelis–Menten) uptake kinetics

kin = cdf.Kinetics()
# reaction id (suffixed, in the merged model), substrate species id, km, vmax
kin.add("GLCpts__A", substrate_id="glc__D_e", km=0.5, vmax=10)
dj = cdf.DynamicJointFBA(cm, {"A":0.01,"B":0.01}, {"glc__D_e":15}, kinetics=kin)

For DynamicParallelFBA, pass kinetics={"A": kin_a, "B": kin_b} (per organism, original reaction ids).

One call for any method — run_dynamic

res = cdf.run_dynamic("joint",      # or "parallel" / "endpoint"
        models=[a, b], biomass_ids=["BIOMASS_A","BIOMASS_B"], organism_ids=["A","B"],
        initial_biomasses={"A":0.01,"B":0.01}, initial_concentrations={"glc__D_e":15},
        dt=0.1, n=1000, solver="glpk")   # endpoint: pass endpoint_n=<#steps>

Medium metabolite essentiality

Re-runs the simulation with each medium metabolite removed (initial conc → 0) and compares growth at time T to the baseline: fraction = growth_without / growth_with. A metabolite is essential when fraction <= growth_fraction. Growth is the biomass at time T (default: end). Scored for the whole community and/or individual members.

df = cdf.medium_essentiality(
        "joint", [a, b], ["BIOMASS_A","BIOMASS_B"], ["A","B"],
        initial_biomasses={"A":0.01,"B":0.01},
        initial_concentrations={"glc__D_e":15,"o2_e":25,"nh4_e":20,"pi_e":20},
        metabolites=None,          # None = every medium metabolite present at t0
        growth_fraction=0.5,       # essentiality threshold
        target=["community","A","B"],   # "community", a member id, or a list
        time=None,                 # time T (default = end of run)
        dt=0.5, n=200, solver="glpk")
df.to_csv("essentiality.tsv", sep="\t")

Returns a DataFrame indexed by metabolite with, per target, the columns <target>_growth, <target>_baseline, <target>_fraction, <target>_essential. Works with any method ("joint"|"parallel"|"endpoint").

Per-bacterium output (TSV)

paths = cdf.export_tsv(res, out_dir="results", prefix="joint_")
# -> results/joint_A.tsv, results/joint_B.tsv

Each file: one row per time point; columns time, biomass, conc_<metabolite> … (shared pool), flux_<reaction> … (that organism's reactions). A flux describes the interval starting at that time, so the final row is blank for flux columns. cdf.per_organism_frames(res) returns the same data as in-memory pandas DataFrames instead of files.

AGORA models + a diet (gut community)

Load a folder of AGORA/SBML models (biomass reaction auto-detected) and an AGORA/VMH diet file, and run with the diet as fixed per-gDW uptake bounds:

models, biomass_ids, organism_ids = cdf.load_models("models")          # 12 GSMMs
diet = cdf.load_diet("medium/western_agora2.txt")   # {"ac[e]": 10.0, ...} (metabolite -> mmol/gDW/h)

res = cdf.run_dynamic("parallel", models, biomass_ids, organism_ids,
                      initial_biomasses={o: 0.01 for o in organism_ids},
                      initial_concentrations={},     # pool starts empty
                      diet=diet, dt=0.2, n=60)
res.relative_abundance();  cdf.export_tsv(res, "results")

# essentiality of diet metabolites (removes each from the diet supply)
df = cdf.medium_essentiality("parallel", models, biomass_ids, organism_ids,
        initial_biomasses={o:0.01 for o in organism_ids},
        diet=diet, growth_fraction=0.5, target="community")

Two ways to apply a medium:

• Finite pool (batch) — pass the medium as initial_concentrations (mmol amounts) with a per-gDW cap max_uptake. The pool depletes; uptake = min(max_uptake·biomass, pool/dt); growth plateaus when the limiting nutrient runs out. This is what the pairwise simulator uses. (The max_uptake cap is required for genome-scale models — AGORA transporters are wide open, so without it growth runs away.)
• Sustained supply (diet=) — each diet metabolite is available at up to diet[m]·biomass (mmol/gDW/h), supplied continuously and not depleted; non-diet metabolites come only from the shared cross-feeding pool.

Either way, exchange reaction bounds are recomputed every step from the current concentration (lower_bound = −min(cap, pool/dt) etc.) — setting a bound on the exchange reaction by hand has no effect; feed the medium through initial_concentrations/diet instead.

Scale/perf notes: Parallel scales to all 12 models; Joint is one ~13k-reaction LP per step (slower but fine); EndPoint is impractical at this scale on GLPK (needs Gurobi). AGORA models on the full rich diet grow fast (high initial µ) — this is an AGORA calibration property; scale the diet or add ATP maintenance to calibrate absolute rates. Relative abundance, cross-feeding, and essentiality are meaningful regardless. See examples/run_gut_community.py.

Preprocessing a model directory

Validate a folder of models and cache the info community generation needs:

models_df, shared_df, meta = cdf.preprocess_models("models")   # writes metaData/*.tsv

Checks each model is usable ("current": loads, has a biomass objective + an extracellular compartment + exchanges, and grows), extracts the biomass reaction (from the objective), and the shared extracellular metabolites (present in ≥2 models) with their average exchange flux. Writes metaData/Preprocess_result.tsv (one row per model: validity, counts, biomass reaction, single-model growth) and metaData/Preprocess_shared_metabolites.tsv. cdf.load_preprocess("metaData") reads them back.

Pairwise simulator + interaction inference

Run Joint + Parallel on every unique pair, plot each bacterium's absolute biomass and its relative abundance over time (absolute is drawn next to the monoculture "alone" curve), and classify the interaction:

summary = cdf.run_pairwise("models", n_steps=40,
             medium_file="medium/western_agora2.txt",  # finite pool (initial conc)
             max_uptake=1.0,                  # per-gDW uptake cap (mmol/gDW/h)
             proportions=(0.5, 0.5),          # initial biomass split of the pair
             out_dir="pairwise_out")

For every pair it runs the two organisms alone (monoculture) and together, then compares final biomass — sign(biomass_in_pair − biomass_alone) for each — to classify:

effect (A, B)	interaction
(+, +)	mutualism
(−, −)	competition
(+, −)	parasitism / exploitation
(+, 0)	commensalism
(−, 0)	amensalism
(0, 0)	neutralism

Initial biomass split defaults to 50/50 (proportions=(0.5, 0.5), total_biomass=0.02 → 0.01 gDW/L each).

Outputs, per pair, two figures (same filename in each folder):

• pairwise_out/plots_absolute/ — absolute biomass (solid = in pair, dashed = alone)
• pairwise_out/plots_relative/ — relative abundance (community composition; the fractions sum to 1 at every time point)

plus pairwise_interactions.tsv (solo/pair biomass + interaction per method) and per-method interaction-type matrices (examples/run_pairwise.py). Absolute mass — not relative abundance — is what reveals the interaction: two organisms can sit at 50/50 relative abundance while both grow far less than they would alone (competition). The relative-abundance plot is complementary — it shows how the composition shifts, not whether either partner was helped or harmed.

cdf.plot_biomass(result, solo=...), cdf.plot_relative_abundance(result) and cdf.classify_interaction(...) are also usable standalone.

Using your own models — practical notes

• initial_concentrations keys are extracellular metabolite ids (e.g. glc__D_e), in mmol/L; this defines the medium and which nutrient is limiting.
• Uptake caps: the code overwrites each exchange's lower bound every step to -concentration/dt (can't take up more than is present). So to cap a per-gDW uptake rate, set the bound on the internal transport reaction (e.g. GLCpts) or use Kinetics — not the EX_ reaction, whose bound is recomputed each step.
• Biomass reaction ids are the original ids of each model; after merging they become "<bm>__<organism>" (see cm.biomass_reactions).

Assumptions (same as dcFBA)

Quasi–steady-state metabolism each step; growth = objective; reaction capacity scales linearly with biomass; forward-Euler integration (Parallel shrinks dt adaptively if a metabolite would go negative); one well-mixed extracellular compartment; batch culture (no in/outflow unless you add it via the deviate hook); interaction only through shared metabolites. Joint/EndPoint maximize the community sum (can give winner-take-all splits between organisms); Parallel is the realistic model for independent competitors.

Package layout

community_dfba/
  community.py   build_community_model + CommunityModel (the merge)
  joint.py       DynamicJointFBA
  parallel.py    DynamicParallelFBA
  endpoint.py    EndPointFBA (+ optional set_qp smoothing)
  kinetics.py    Monod kinetics container
  runner.py      run_dynamic() - one call for any method
  essentiality.py medium_essentiality() - metabolite essentiality screen
  export.py      export_tsv() / per_organism_frames() - per-bacterium output
examples/
  test_merge.py              merge correctness check
  test_joint_parallel.py     Joint + Parallel validation
  test_endpoint.py           EndPointFBA validation
  test_essentiality_export.py essentiality + TSV export validation
  compare_methods.py         all three on one community
  plot_comparison.py         saves comparison.png