Differential Equation Model of a Chemostat

Introduction

The example shows a chemostat model (cf. Novick and Szilard 1950; Herbert, Elsworth, and Telling 1956) by example of heterotrophic or autotrophic organisms (bacteria, algae, \(X\)), where growth is limited by a single nutrient (e.g. phosphorus, \(P\)) via a Monod equation.

Implementation

The system is solved numerically with solvers from package deSolve (Soetaert, Petzoldt, and Setzer 2010b), while equilibria at infinite time are estimated with package rootSolve (Soetaert 2009). The implementation follows the standard formulation of these packages, see package documentation of deSolve and rootSolve, and Soetaert, Petzoldt, and Setzer (2010b) or Soetaert, Petzoldt, and Setzer (2010a) for details. The code is written in the R programming language for statistical computing (R Core Team 2018).

library("deSolve")
library("rootSolve")

chemostat <- function(time, init, parms) {
  with(as.list(c(init, parms)), {
    mu   <- mumax * P/(kp + P)  # Monod equation
    dX <- mu * X - D * X
    dP   <-  D *(P0 - P) - 1/Y * mu * X
    list(c(dX, dP), mu=mu)
   })
}
parms <- c(
  mumax = 0.5,    # 1/d
  kp    = 0.01,   # half saturation constant (mg/L)
  Y     = 41,     # yield coefficient (stoichiometric C:P ratio)
  D     = 0.1,    # 1/d
  P0    = 0.05    # P in inflow (mg/L)
)
times <- seq(0, 40, 0.1)  # (d)
init  <- c(X=0.01, P=0.05) # Phytoplankton C and Phosphorus P (mg/L)

Dynamic simulation

A dynamic simulation can then be performed by function ode, using the default algorithm lsoda with automatic integration step size. The result (out) is then a matrix-like object of class deSolve that is supported by a generic (i.e. object oriented) plot function.

out <- ode(init, times, chemostat, parms)
plot(out)

The plot function allows also to show several scenarios simultanaeously:

init  <- c(X=1.0, P=0.05)
p1 <- p2 <- p3 <- parms
p1["D"] <- 0; p2["D"] <- 0.3; p3["D"] <- 0.5
out <- ode(init, times, chemostat, parms)
out1 <- ode(init, times, chemostat, p1)
out2 <- ode(init, times, chemostat, p2)
out3 <- ode(init, times, chemostat, p3)

plot(out, out1, out2, out3, which=c("X", "P"))

Steady State

The equilibrium can be aproximated numerically with package rootSolve. We use a loop here for simplicity, even if list-based approaches (e.g. lapply) may appear as more elegant for the intermediate and advanced R user.

state <- data.frame(
  D = seq(0, 0.6, length.out = 100),
  X = 0,
  S = 0
)

init  <- c(X=0.01, P=0.05) 

for (i in 1:nrow(state)) {
  parms["D"] <- state$D[i]
  times <- c(0, Inf)
  out <- runsteady(init, times, chemostat, parms)
  state[i, 2:3] <- out$y
}

par(mfrow = c(1, 3))
plot(S ~ D, data = state, type = "l")
plot(X ~ D, data = state, type = "l")
plot(S * X ~ D, data = state, type = "l")

Analytical Solution

The same equilibrium can be reproduced by the analytical solution of chemostat equations.

D <- seq(0, 0.6, length.out = 100)
mumax = 0.5;  kp    = 0.01; Y     = 41; P0    = 0.05

P <- D * kp / (mumax - D)
P <- ifelse(D > (mumax * P0)/(kp + P0), P0, P)

X <- Y * (P0 - P)

par(mfrow=c(1,3))
plot(D, P, type="l")
plot(D, X, type="l")
plot(D, X*P, type="l")

Outlook

The chemostat equations can serve as a starting point for further extensions, e.g. the production of specific substances in bioreactors, or additional trophical levels. The equations are the fundamental for the description of other technical or non-technical flow-through systems, up to waste water treatment plants, or lakes.

References

Herbert, Denis, R Elsworth, and RC Telling. 1956. “The Continuous Culture of Bacteria; a Theoretical and Experimental Study.” Microbiology 14 (3). Microbiology Society: 601–22.

Novick, Aaron, and Leo Szilard. 1950. “Description of the Chemostat.” Science 112 (2920). American Association for the Advancement of Science: 715–16.

R Core Team. 2018. R: A Language and Environment for Statistical Computing. Vienna, Austria: R Foundation for Statistical Computing. https://www.R-project.org/.

Soetaert, Karline. 2009. rootSolve: Nonlinear Root Finding, Equilibrium and Steady-State Analysis of Ordinary Differential Equations. R Package 1.6. https://cran.r-project.org/package=rootSolve.

Soetaert, Karline, Thomas Petzoldt, and R. Woodrow Setzer. 2010a. “Solving Differential Equations in R.” The R Journal 2 (2): 5–15. https://journal.r-project.org/archive/2010-2/RJournal_2010-2_Soetaert~et~al.pdf.

———. 2010b. “Solving Differential Equations in R: Package deSolve.” Journal of Statistical Software 33 (9): 1–25. doi:10.18637/jss.v033.i09.