Simulate Bacterial and Phage Population Dynamics in Continuous Culture
by Stephen T. Abedon Ph.D. (abedon.1@osu.edu)
phage.org | phage-therapy.org | biologyaspoetry.org | abedon.phage.org | google scholar
Jump to: 🔬 Simulator | 📖 Background | ⚙️ Methodology | 📘 Abedon 2009 | 🧮 More Calculators
chemostat.phage.org · Abedon’s Books
How can I improve this page? contact: chemostat@phage.org
The chemostat is a device for maintaining continuous culture of microorganisms under steady-state conditions. Fresh medium containing a limiting nutrient flows into the vessel at a constant dilution rate D, while an equal volume of culture flows out, washing cells and phages out of the system at the same rate. At steady state, bacterial growth rate equals the dilution rate: μ = D.
The chemostat is invaluable as an experimental and theoretical tool because it allows precise control of microbial growth rate and population density through adjustable dilution rate and nutrient supply. It has been used extensively to study microbial evolution, competition, and phage-bacteria coevolution under well-defined conditions.
Bacterial growth rate in the chemostat depends on the concentration of the limiting nutrient S according to the Monod equation:
where μmax is the maximum specific growth rate and Ks is the half-saturation constant (the nutrient concentration at which growth rate is half its maximum). At steady state in the phage-free chemostat, the nutrient concentration settles at:
and bacterial density at:
where Y is the yield coefficient (cells produced per unit nutrient consumed) and S₀ is the incoming nutrient concentration. When the dilution rate approaches μmax, S* → S₀ and N* → 0 (washout).
The foundational theoretical work on phages in chemostats was developed by Levin, Stewart, and Chao [1] and expanded by Campbell [2] and others. These models demonstrate that phages can stably coexist with bacteria in the chemostat when the phage can sustain an adequate adsorption rate relative to washout.
Phage persistence in the chemostat requires that the phage net growth rate exceed its loss by dilution. This produces a threshold condition on the adsorption rate constant k, burst size β, and the bacterial density supported in the phage-free chemostat. When phages do persist, they typically depress the bacterial population density below its phage-free steady state — sometimes dramatically — while the nutrient concentration rises (because fewer bacteria are consuming nutrient). This is sometimes informally described as "killing the winner": the dominant bacterial competitor is suppressed, potentially allowing other strains to persist.
The simplest phage-in-chemostat models collapse all intracellular phage replication into an instantaneous lysis rate δ (ODE model). This is mathematically convenient but biologically approximate: real phage lytic cycles have an explicit latent period τ during which infected cells neither grow nor release virions.
Incorporating the latent period explicitly via a delay-differential equation (DDE) model can produce qualitatively different dynamics compared to the ODE model [3,4]. The delay means that phages released at time t were produced by infections that occurred at time t − τ. Since bacteria and phages in the chemostat at time t − τ may have been at different densities, this creates a lag in the phage response to bacterial density changes. The DDE model can exhibit sustained oscillations — limit cycles — that the ODE model damps to a steady state. These oscillations become more prominent when the latent period is long relative to characteristic ecological timescales.
The side-by-side comparison in this calculator allows direct visualization of the differences between these two modeling approaches for any given set of biological parameters.
When phages are absent, both models reduce to the standard two-dimensional chemostat system of ordinary differential equations:
where D is the dilution rate, S₀ is the incoming nutrient concentration, Y is the yield coefficient, μmax is the maximum specific growth rate, and Ks is the half-saturation constant.
The ODE model adds an infected-cell compartment I and a free phage compartment P. Adsorption creates infected cells at rate k·N·P, infected cells lyse at rate δ (the instantaneous lysis rate), and each lysis event releases β phages:
Note that infected cells are washed out at rate D just like uninfected cells. The total lysis rate is δ·I, and each lysis event produces β free phages net (the parental phage particle is consumed on adsorption, accounted for in the −k·N·P term of the phage equation).
The DDE model replaces the infected-cell compartment with an explicit latent period delay τ. Phages released at time t come from cells infected at time t − τ, discounted by dilution-rate washout of the infected cell over the latent period:
The factor e−D·τ accounts for the probability that an infected cell is not washed out during the latent period. When D·τ is small this factor approaches 1 and the two models converge; when D·τ is large, washout during the latent period substantially reduces effective phage reproduction.
When phage-resistant bacteria R are introduced at time tR, a fifth equation is added to both models. Resistant bacteria share the nutrient dynamics with sensitive bacteria but carry no phage interaction terms:
The fitness cost factor fc provides a convenient shorthand: setting fc < 1 scales μmax,R = fc × μmax downward and Ks,R = Ks / fc upward, reflecting the typical observation that phage resistance carries a metabolic cost that reduces both maximum growth rate and nutrient affinity. The separate μmax,R and Ks,R inputs allow independent fine-tuning.
The six modeling modifications described in Abedon (2009) [1] are implemented as optional adjustments that overlay the base models. All default to neutral (no effect). When enabled, they modify the equations as follows:
The B&L fit values (accessible via the "Load B&L fit values" button) are those used in Abedon (2009) [1] Fig. 9E/F to approximate the first ~200 hours of the Bohannan & Lenski (1997) [2] chemostat. All six refinements are active: tlag = 5 hr, dP = 0.0031 min⁻¹, fY = 0.75, fτ = 1.333 (τ → 48 min), fβ = 0.75 (β → 60), fk = 0.75, fμ = 0.75. Wall population input is not active. See the 📘 Abedon 2009 tab for the full parameter table and discussion.
The ODE model (and the bacteria-only model) are integrated using an adaptive-step RK4 method. Each step is evaluated twice — once as a single step of size h, and once as two half-steps of size h/2. The difference provides a local error estimate. If the error exceeds the tolerance (relative + absolute, tol = 10⁻⁵), the step is rejected and h is reduced; if well below tolerance, h is increased. Step size is bounded between 10⁻⁷ hr and 0.02 hr. When phage-resistant bacteria are introduced at tR, the ODE simulation is split into two segments — 0 to tR and tR to tmax — with R₀ added to the state at the boundary, ensuring the discontinuity is handled cleanly.
The DDE model uses a fixed-step RK4 with a circular history buffer, because adaptive step sizing would require variable-length delay-buffer interpolation. The fixed step is set to min(τeff/40, 0.002 hr). The delay terms N(t−τ) and P(t−τ) are retrieved from the buffer at a fixed lag of round(τeff/h) steps; initial conditions are used for t < τeff. All populations are floored at zero after each step.
For the phage-free chemostat, the steady states are:
For the phage model, the phage-positive interior equilibrium satisfies μ(S**) = D + k·P** with additional conditions; the calculator reports final simulated values rather than solving the algebraic system analytically, since the DDE model and the Abedon (2009) refinements have no simple closed-form equilibria.
All internal calculations use hours (hr) as the time unit and µg/mL as the nutrient concentration unit. The rate/time parameter display toggle converts values between per-minute and per-hour for convenience; simulation duration is always in hours. The Bohannan & Lenski (1997) preset loads parameters converted from their original per-minute values.
This chapter takes as its central case study the 0.5 mg glucose/L chemostat experiment reported by Bohannan and Lenski [1] — specifically their Figure 3B, which paired a long-term experimental record of bacterial and phage population density with a mechanistic model. The central argument is that obtaining genuine quantitative precision in modeling phage-bacterial chemostat dynamics requires more thorough experimental characterization of the system's components — its "deconstruction" — before and alongside the main experiment, rather than relying on post hoc parameter adjustment afterward.
The chapter is explicitly positioned toward the "precision" end of the realism–generality–precision tradeoff in model building. Rather than asking what general principles apply to predator–prey systems, it asks: given this specific phage, this specific host, and these specific conditions, can we predict what actually happens quantitatively?
Comparing the Bohannan and Lenski [1] simulation against their actual chemostat data reveals five systematic discrepancies that the standard model fails to capture:
Following Schrag and Mittler, bacterial populations adhering to chemostat vessel walls can serve as a phage-protected reservoir, continuously seeding the planktonic phase with susceptible cells. Even a modest input of 0.1 cells/mL/min from wall populations substantially stabilizes the simulated chemostat and brings bacterial dynamics qualitatively closer to the observed pattern. Whether wall populations actually contributed to the Bohannan–Lenski [1] chemostat remains unknown — but the chapter argues that this possibility should be empirically assessed before publication of any long-term chemostat experiment.
One of the chapter's most distinctive findings concerns the rate at which free phage densities decline during troughs. By measuring the slopes of log-transformed phage decline from the actual chemostat data and comparing them against the slope expected from dilution alone, the author finds that free phages are disappearing roughly twice as fast as outflow can account for. This pattern is not unique to the Bohannan–Lenski [1] system: a systematic survey of phage decline rates across more than a dozen published phage-bacterial chemostat studies finds a mean excess decline ratio of approximately 1.7-fold, suggesting the phenomenon may be widespread. Candidate mechanisms include adsorption to bacterial debris, phage instability in the specific growth medium, or undetected variation in actual flow rates — but the mechanism remains unresolved. When this excess decay is incorporated into the model, the degree of phage–bacterial "nestling" increases substantially and the chemostat approaches damped rather than sustained oscillations.
Bohannan and Lenski [1] initiated their chemostat with stationary-phase bacteria, creating an initial physiological lag before exponential growth begins. Modeling a lag of approximately five hours brings the early phage and bacterial trajectories into reasonable agreement with the actual data. However, even with this lag, simulated bacteria do not reach the peak densities observed experimentally, suggesting additional factors are at work.
The parameter values used by Bohannan and Lenski [1] were derived from experiments conducted at nutrient densities far higher (300 mg glucose/L) than those in their chemostat (0.5 mg glucose/L). Under nutrient limitation, phage T4 infections of smaller, slower-growing bacteria may yield reduced burst sizes, longer latent periods, and lower adsorption rate constants. Cumulatively reducing these parameters — to approximately 75% of their standard values — in combination with the lag phase and excess phage decay, produces a simulation (the chapter's Fig. 9F) that tracks the actual chemostat data more closely than any of the individual modifications alone.
Phage T4, the organism used by Bohannan and Lenski [1], displays lysis inhibition [2] — a phenotype in which secondary adsorption of an already-infected bacterium triggers an extended latent period and enlarged burst size. At the phage densities observed in the 0.5 mg/L chemostat, rough calculations suggest that a substantial fraction of infected bacteria (~20–40%) may be undergoing lysis inhibition at peak phage densities. This could simultaneously slow early phage population growth (when few secondary adsorptions occur) and amplify peak phage production, producing greater instability overall. The chapter notes that the 13-fold increase in phage equilibrium densities observed when glucose is raised from 0.1 to 0.5 mg/L — far exceeding the 3-fold increase the standard model predicts — is consistent with lysis inhibition playing a density-dependent role. Because lysis inhibition is difficult to model mechanistically, the chapter's recommendation is direct: avoid phages that display lysis inhibition (or lysogeny) in long-term chemostat experiments unless lysis inhibition itself is the object of study.
The chapter's constructive contribution goes beyond diagnosing discrepancies. It outlines a minimum set of short-term control experiments — collectively termed "chemostat deconstruction" — that should be performed and reported alongside any long-term phage-bacterial chemostat study:
The chapter closes by noting that these controls are not merely desirable but necessary if the goal is quantitative prediction rather than qualitative description. The Bohannan and Lenski [1] experiment itself ran for 600 hours — a substantial investment — and the chapter argues that investing that kind of time in a long-term chemostat experiment without first characterizing its components is a missed opportunity for the kind of precision that could connect laboratory phage ecology to real-world applications such as phage therapy.
The table below shows the complete parameter set used in the Abedon (2009) [3] spreadsheet model (their Fig. 9F equivalent), which was the closest approximation to the first ~200 hours of the Bohannan & Lenski (1997) [1] chemostat. The 📘 Load B&L fit values button in the Simulator tab loads all of these values simultaneously, overriding both the main B&L (1997) preset parameters and the Abedon (2009) refinements.
| Parameter | Symbol | B&L (1997) value | Abedon (2009) Fig. 9F | Units | Notes |
|---|---|---|---|---|---|
| Chemostat & Bacterial Parameters | |||||
| Flow rate | F | 0.00333 min⁻¹ | 0.00333 min⁻¹ | = 0.2 hr⁻¹ | Unchanged |
| Nutrient supply | Sr | 0.5 mg/L | 0.5 mg/L | = 0.5 µg/mL | Unchanged |
| Max growth rate | μmax | 0.01296 min⁻¹ | 0.01296 min⁻¹ | = 0.7726 hr⁻¹ | Unchanged |
| Half-saturation | Ks | 0.0727 mg/L | 0.0727 mg/L | = 0.0727 µg/mL | Unchanged |
| Reciprocal yield | Sd | 2×10⁻¹² g/cell | 1.5×10⁻¹² g/cell | = fY = 0.75 | ⚙️ Modified (Abedon 2009) |
| Initial bacteria | N₀ | 10⁴ /mL | 2×10⁴ /mL | /mL | ⚙️ Boosted (24-hr stationary phase inoculum) |
| Phage Parameters | |||||
| Adsorption rate | k | 5×10⁻⁹ mL/min | 3.75×10⁻⁹ mL/min | = 2.25×10⁻⁷ mL/hr | ⚙️ fk = 0.75 (panels D–F) |
| Burst size | β | 80 | 60 | phages/cell | ⚙️ fβ = 0.75 (panels C–F) |
| Latent period | τ | 36 min | 48 min | = 0.8 hr | ⚙️ fτ = 1.333 (panels B–F) |
| Initial phage | P₀ | 10⁶ /mL | 10⁶ /mL | /mL | Unchanged |
| Abedon (2009) Refinements — all panels B–F share the first three | |||||
| Bacterial lag phase | tlag | 0 | 5 hr (300 min) | hr | ⚙️ All panels B–F. Inoculum was 24-hr stationary-phase bacteria |
| Excess phage decay | dP | 0 | 0.0031 min⁻¹ | min⁻¹ | ⚙️ All panels B–F. Mechanism unknown; ≈ outflow rate |
| Yield factor | fY | 1.0 | 0.75 | — | ⚙️ All panels B–F. Sd = 1.5×10⁻¹² g/cell |
| Latent period factor | fτ | 1.0 | 1.333 | — | ⚙️ All panels B–F. τ → 48 min (zero-min rise) |
| Burst size factor | fβ | 1.0 | 0.75 | — | ⚙️ Panels C–F. β → 60 |
| Adsorption factor | fk | 1.0 | 0.75 | — | ⚙️ Panels D–F. k → 3.75×10⁻⁹ mL/min |
| Growth rate factor | fμ | 1.0 | 0.75 | — | ⚙️ Panels E–F. μmax → 0.00972 min⁻¹ |
| Wall population input | Nw | 0 | 0 | /mL/min | Not active in Fig. 9F |
The B&L (1997) [1] preset in this calculator reproduces the standard model without refinements. The 📘 Load B&L fit values button applies all Fig. 9E/F parameters simultaneously. The chapter's analysis suggests that no single modification resolves all five discrepancies simultaneously — the cumulative effect of all six active adjustments (panels E/F) is required for reasonable agreement with the first ~200 hours of the actual chemostat.