ODE Growth models with package deSolve

Thomas Petzoldt

6th February 2017

Overview

Growth models
Differential equation solvers
Resource limited growth and chemostat
Optional: semicontinuous cultures

Birth, death, population growth

\[ \frac{dN}{dt} = b\cdot N - d\cdot N\\ \frac{dN}{dt} = (b-d) N = rN \]

abundance \(N\) (“number”)
birth rate \(b\), death rate \(d\)
population growth rate \(r\) can be positive or negative

Analytical solution

\[\begin{align} \frac{dN}{dt} &= rN \\ \int_0^t \frac{dN}{dt} &= r\int_0^tN \\ \int_0^t \frac{1}{N}dN &= r\int_0^t dt \\ \ln(N) &= rt + c \\ \end{align}\]

\[ N_t = N_0 e^{rt} \]

This is the exponential growth model.

Plot analytical solution as function of time

r  <- 0.5
N0 <- 10
dt <- 0.1
time <- seq(0, 10, dt)
# analytical solution
N <- N0 * exp(r * time)
plot(time, N, type="l")

Numerical simulation in a loop

N <- numeric(length(time))
N[1] <- N0
for (i in 2:length(time)) {
  N[i] <- N[i-1] + r * N[i-1] * dt
}
plot(time, N, type = "l")

This is the Euler method.

Numerical simulation with package deSolve

library(deSolve)
model <- function (time, y, parms) {
  with(as.list(c(y, parms)), {
    dN <-   r * N
    list(c(dN))
  })
}
y0     <- c(N = 0.1)
parms  <- c(r = 0.1, K = 10)
times  <- seq(0, 100, 1)
out <- ode(y0, times, model, parms) # <---
plot(out)

How is growth limited?

Unlimited
- exponential growth
Carrying capacity
- logistic growth (and others, see tomorrow)
Limiting resource
- carbon, nutrients, …
Interspecific interactions
- predation, competition, space, chemo-communication, quorum sensing

Solution of the logistic

\[\begin{align} \\ \frac{dN}{dt} &= r \cdot N \left(1-\frac{N}{K}\right) \\ \\ N_t &= \frac{K N_0 e^{rt}}{K + N_0 (e^{-r t} -1)} \end{align}\]

Numerical simulation

library(deSolve)
model <- function (time, y, parms) {
  with(as.list(c(y, parms)), {
    dN <-   r * N * (1 - N / K)
    list(c(dN))
  })
}
y0     <- c(N = 0.1)
parms  <- c(r = 0.1, K = 10)
times  <- seq(0, 100, 1)
out <- ode(y0, times, model, parms) # <---
plot(out)

Numerical solvers

… a very broad topic, please see http://desolve.r-forge.r-project.org
- links to books (e.g. Soetaert, Cash and Mazzia, 2012)
- slides, papers, conference proceedings
rk4: very popular, easy to implement, but: either inaccurate or slow
lsoda
- default, works with non-stiff and stiff systems
- automatic internal stepsize, fixed external steps
- accuracy specified by atol and rtol
important:
- atol (absolute tolerance) depends on scale of state variables, default (\(10^{-6}\))
  - should be changed if working with very small or big numbers
  - Hint: set atol = \(10^{2}\) if bacteria grow from \(10^{6}\) to \(10^{9}\)
- rtol (relative tolerance)
  - Note: R calculates with 16 significant digits, so rtol < \(10^{-16}\) does not make sense.

Resource limited phytoplankton growth

\[ f(P) = \frac{P}{k_P + P}\\ \frac{dAlg}{dt} = r \cdot f(P) \cdot Alg\\ \frac{dP}{dt} = - r \cdot 1/Y \cdot f(P) \cdot Alg \]

\(f(P)\) is a Monod function with half saturation constant \(k_P\)
\(Y\) is the yield coefficient

Numerical simulation

model <- function (time, y, parms) {
  with(as.list(c(y, parms)), {
    f <- P/(kP + P)
    dAlg <- r * f * Alg 
    dP   <- - r * 1/Y * f * Alg
    list(c(dAlg, dP))
  })
}
y      <- c(Alg = 0.1, P = 0.2)         # in mg/L
parms  <- c(r = 0.1, kP = 5e-3, Y = 41) # Y = C:P mass ratio
times  <- seq(0, 100, 1)
out <- ode(y, times, model, parms)
plot(out)

Hint: Stoichiometric calculations in R

library(marelac)

redfield(1, species="P")

##     C   H   O  N P
## 1 106 263 110 16 1

redfield(1, species="P", method="mass")

##          C        H        O        N P
## 1 41.10363 8.558477 56.82016 7.235388 1

molweight("C2H5OH")

##   C2H5OH 
## 46.06844

Reactors

Chemostat

Continuously Stirred Tank Reactor (CSTR)
open system, flow-through
continuous supply of resource
export of organisms and remaining substrate

Batch and Chemostat

Chemostat model

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

Equilibrium

\[\begin{align} \frac{dX}{dt} &= \mu X - D X &= 0\\ \frac{dS}{dt} &= D(S_0 - S) - \mu \cdot 1/Y \cdot X &= 0 \end{align}\]

Gain and loss equal \(\rightarrow\) all derivatives are zero.
\(\rightarrow\) dynamic equilibrium (German: Fließgleichgewicht), turnover can be high!

library(deSolve)
batch <- function(time, y, parms){
  with(as.list(c(y, parms)), {
    f <- r * S / (ks + S)
    dS <- - 1/Y * f * N
    dN <- f * N
    return(list(c(dS, dN)))
  })
}
y <- c(S = 10, N = 1e4)
parms <- c(r=1, ks=5, Y=1e6, S0=10)
times <- seq(0, 20, 0.1)
out <- ode(y, times, batch, parms)
plot(out)

parms <- c(r=1, ks=5, Y=1e6, S0=10, D=0.5)
times <- seq(0, 50, 0.01)
etime <- seq(5, 50, 5) # time points to triger events
eventfun <- function(t, y, parms) { # event function
  with(as.list(c(y, parms)), {
    return(c(D * S0 + (1-D) * S, (1-D) * N)) # D = dilution rate
  })
}
out <- ode(y, times, batch, parms,
           events = list(func = eventfun, time = etime))

plot(out)

library(deSolve)
batch <- function(time, y, parms){
  with(as.list(c(y, parms)), {
    f <- r * S / (ks + S)
    dS <- - 1/Y * f * N
    dN <- f * N
    return(list(c(dS, dN)))
  })
}
y <- c(S = 10, N = 1e4)
parms <- c(r=1, ks=5, Y=1e6, S0=10)
times <- seq(0, 20, 0.1)

parms <- c(r=1, ks=5, Y=1e6, S0=10, D=0.5)
times <- seq(0, 50, 0.01)
etime <- seq(5, 50, 5) # time points to triger events
eventfun <- function(t, y, parms) { # event function
  with(as.list(c(y, parms)), {
    return(c(D * S0 + (1-D) * S, (1-D) * N)) # D = dilution rate
  })
}


crit <- 0.9e7 # critical cell number that triggers dilution event
rootfun <- function (t, y, pars) return(crit - y[2])
out <- ode(y, times, batch, parms,
events = list(func = eventfun, root = TRUE), rootfun = rootfun)

plot(out)

Thanks to:

Karline Soetaert (NIOZ, NL) and
Woodrow Setzer (US EPA)

And to you.

More food …

http://desolve.r-forge.r-project.org/

http://desolve.r-forge.r-project.org/user2014/tutorial.pdf

http://desolve.r-forge.r-project.org/user2014/examples/vanderpol/vanderpol.svg

http://desolve.r-forge.r-project.org/user2014/examples/FME/fit_twocomp.svg

http://desolve.r-forge.r-project.org/user2014/examples/compiled_lorenz/compiledcode.svg

http://desolve.r-forge.r-project.org/user2014/examples/