Dynamical Systems in Biology

Exercise 1: Lotka-Volterra

The ODEs of the Lotka-Volterra model, $\dot{\vec x} = f(\vec x)$, are given by: $$\begin{align} \dot{x}(t) & = a \cdot x(t) - b \cdot x(t) \cdot y(t), \\ \dot{y}(t) & = c \cdot x(t) \cdot y(t) - d \cdot y(t). \end{align}$$ where $\mathbb{R} \ni t \mapsto x(t) \in \mathbb{R}$ and $\mathbb{R} \ni t \mapsto y(t) \in \mathbb{R}$, and all parameters $(a, b, c, d) \in \mathbb{R}_{+}^4$ are positive.

• Write a function euler_LV <- function(ini, pars, dt, tmax) integrating the LV system via the Euler method for a set of initial values ini, parameters pars, stepsize dt for $t\in[0,{\tt tmax}]$. The function should return a data.frame with columns time, x, y. Reminder: $\vec x_{n+1} = \vec x_n + f(\vec x_n)\cdot dt$.

# Lotka-Volterra euler method
euler_LV <- function(ini, pars, dt, tmax) {
  with(as.list(pars),{
  N <- tmax/dt
  x <- y <- numeric(N+1)
  x[1] <- ini[1]
  y[1] <- ini[2]
  
  for(i in 1:N){
    dx <- a*x[i] - b*x[i]*y[i]
    dy <- c*x[i]*y[i] - d*y[i]
    x[i+1] <- x[i] + dt*dx
    y[i+1] <- y[i] + dt*dy
  }
  times <- seq(from = 0, by = dt, length.out = N+1)
  return(data.frame(time = times, x = x, y = y))
  })
}

• Plot solutions of the LV system using ggplot2() in configuration and phase space for $a = 2/3$, $b = 4/30$, $c = 1/10$, $d = 1$ and a set of initial values $(x(t=0), y(t=0)) = \{ (10,15), (15,10)\}$. To do so, first arrange the output of euler_LV in a config or phase space data.frame with columns:
  • time, value, species, ini, where species is a character variable labelling the states x and y.
  • x, y, ini, where ini is a character or factor variable to label the initial value.

library(ggplot2)
theme_set(theme_bw())

pars <- c(a = 2/3, b = 4/30, c = 1/10 , d = 1)

# initial values: x<y
ini <- c(x = 10, y = 15)
out1 <- euler_LV(ini, pars, 1e-3, 50)

# initial values: x>y
ini <- c(x = 15, y = 10)
out2 <- euler_LV(ini, pars, 1e-3, 50)

# convert the outputs to data.frames
out2config <- function(out, ini) {
  rbind(data.frame(time = out$time, value = out$x, species = "x", ini = ini), 
        data.frame(time = out$time, value = out$y, species = "y", ini = ini))
}


out2phase <- function(out, ini) {
  data.frame(x = out[,2], y = out[,3], ini = ini)
}

outConfig <- out2config(out1, "x0 < y0")
outPhase <- out2phase(out1, "x0 < y0")
outConfig <- rbind(outConfig, out2config(out2, "x0 > y0"))
outPhase <- rbind(outPhase, out2phase(out2, "x0 > y0"))

# configuration space
P1 <- ggplot(outConfig, aes(x = time, y = value, group = species, lty = species)) + geom_line() + facet_wrap(~ini)
print(P1)

# phase space
P2 <- ggplot(outPhase, aes(x = x, y = y, group = ini, lty = ini)) + geom_path()
print(P2)

• Add a noise term to the euler integration: $x_{n+1} = x_n + f(x_n)\cdot dt + \epsilon \cdot \sqrt{dt}$ with $\epsilon\sim \mathcal{N}(0, {\tt sd)}$. Use sd = c(sdx, sdy) as additional argument of your euler function and avoid $x_n < 0$. Compare the solutions obtained with sd = c(.3, .3) to the standard solutions without noise by plotting.

# Lotka-Volterra euler method with noise
eulerN_LV <- function(ini, pars, dt, tmax, sd) {
  with(as.list(pars),{
  N <- tmax/dt
  x <- y <- numeric(N+1)
  x[1] <- ini[1]
  y[1] <- ini[2]
  
  for(i in 1:N){
    dx <- a*x[i] - b*x[i]*y[i]
    dy <- c*x[i]*y[i] - d*y[i]
    noisex <- sqrt(dt)*rnorm(1, sd = sd[1])
    noisey <- sqrt(dt)*rnorm(1, sd = sd[2])
    x[i+1] <- x[i] + dt*dx + noisex
    y[i+1] <- y[i] + dt*dy + noisey
    if(x[i+1] < 0) x[i+1] <- 0
    if(y[i+1] < 0) y[i+1] <- 0
    }
  times <- seq(from = 0, by = dt, length.out = N+1)
  return(data.frame(time = times, x = x, y = y))
  })
}


pars <- c(a = 2/3, b = 4/30, c = 1/10 , d = 1)
ini <- c(x = 10, y = 10)
outN <- eulerN_LV(ini, pars, 1e-3, 100, sd=c(.3,.3))
outnoN <- eulerN_LV(ini, pars, 1e-3, 100, sd=c(0,0))
outConfig <- rbind(out2config(outN, "noise"),out2config(outnoN, "no noise"))
outPhase <- rbind(out2phase(outN, "noise"), out2phase(outnoN, "no noise"))
P1 <- ggplot(outConfig, aes(x = time, y = value, group = ini , color = ini)) + facet_wrap(~species) + geom_line() 
print(P1)

P2 <- ggplot(outPhase, aes(x = x, y = y, group = ini, lty = ini, color = ini)) + geom_path() + scale_color_manual(values = c("black", "red"))
print(P2)

Exercise 2: Lotka-Volterra extended

The LV system can be extended to: $$\begin{align} \dot{x}(t) & = a \cdot x(t) \cdot \left(1 - \frac{x(t)}{K} \right) - b \cdot \frac{x(t)}{x(t) + S} \cdot y(t), \\ \label{eq:ELV2} \dot{y}(t) & = c \cdot \frac{x(t)}{x(t) + S} \cdot y(t) - d \cdot y(t), \end{align}$$ where also $K$ and $S$ are positive parameters.

• Repeat Exercise 1 for the extended equations using $a = b = c = 1$; $d= \tfrac{1}{3}$; $K = 30$, $S=10$ and $(x(0),y(0)) = \{(1,1), (10,10)\}$.

# extended Lotka-Volterra euler method with noise
eulerN_LVE <- function(ini, pars, dt, tmax, sd) {
  with(as.list(pars),{
  N <- tmax/dt
  x <- y <- numeric(N+1)
  x[1] <- ini[1]
  y[1] <- ini[2]
  
  for(i in 1:N){
    #dx <- a*x[i]*(1-x[i]/K)- b*x[i]*y[i]
    #dy <- c*x[i]*y[i] - d*y[i]
    dx <- a*x[i]*(1-x[i]/K)- b*x[i]*y[i]/(x[i]+S)
    dy <- c*x[i]*y[i]/(x[i]+S) - d*y[i]
    #dx <- a*x[i] - b*x[i]*y[i]/(x[i]+S)
    #dy <- c*x[i]*y[i]/(x[i]+S) - d*y[i]
    noisex <- sqrt(dt)*rnorm(1, sd = sd[1])
    noisey <- sqrt(dt)*rnorm(1, sd = sd[2])
    x[i+1] <- x[i] + dt*dx + noisex
    y[i+1] <- y[i] + dt*dy + noisey
    if(x[i+1] < 0) x[i+1] <- 0
    if(y[i+1] < 0) y[i+1] <- 0
    }
  times <- seq(from = 0, by = dt, length.out = N+1)
  return(data.frame(time = times, x = x, y = y))
  })
}

pars <- c(a = 1, b = 1, c = 1 , d = 1/3, K = 30, S = 10)

ini <- c(x = 1, y = 1)
outN1 <- eulerN_LVE(ini, pars, 1e-3, 100, sd=c(0.5,0.5))
outnoN1 <- eulerN_LVE(ini, pars, 1e-3, 100, sd=c(0,0))

ini <- c(x = 10, y = 10)
outN2 <- eulerN_LVE(ini, pars, 1e-3, 100, sd=c(0.5,0.5))
outnoN2 <- eulerN_LVE(ini, pars, 1e-3, 100, sd=c(0,0))

outConfig <- cbind(rbind(out2config(outnoN1, "x=y=1"),out2config(outnoN2, "x=y=10")), noise = "no")
outConfig <- rbind(outConfig,cbind(rbind(out2config(outN1, "x=y=1"),out2config(outN2, "x=y=10")), noise = "yes"))
outPhase <- cbind(rbind(out2phase(outnoN1, "x=y=1"),out2phase(outnoN2, "x=y=10")), noise = "no")
outPhase <- rbind(outPhase,cbind(rbind(out2phase(outN1, "x=y=1"),out2phase(outN2, "x=y=10")), noise = "yes"))

P1 <- ggplot(outConfig, aes(x = time, y = value, group = ini , color = ini)) + facet_grid(noise~species, scales = "free_y") + geom_line() 
print(P1)

P2 <- ggplot(outPhase, aes(x = x, y = y, group = ini, lty = ini, color = ini)) + geom_path() + facet_wrap(~noise) + scale_color_manual(values = c("black", "red"))
print(P2)

• Interpret the extended system:
  • What is the meaning of the parameters?
  • What are the qualitative differences of both systems?
  • Compare how the noise affects both systems.

# -(a/K)x^2 = intraspecific competition
# x -> x/(x+S) = saturation of the production

Cathedral exercise

The choir of the cathedral is slightly inclined relative to the main nave. Why?