Dynamical Systems in Biology

Tutorial 1

Exercise 1: Plotting Verhulst dynamics

• Write a function f <- function(t, p) for Eq. (2) returning a data.frame with columns t, x, x0, K, r.

f <- function(t, p) {
  with(as.list(p), {
    x <- K*x0/(x0+(K-x0)*exp(-r*t))
    out <- data.frame(t = t, x = x, x0 = x0, r = r, K = K)  
    return(out)
  })
}

• Plot \(x(t)\) on the interval \([0,1]\) for \(K=0.5, r=4\) and a set of values for \(x_0\) in [0,1]. Group the plot by values of \(x_0\).

library(ggplot2)
theme_set(theme_bw())
times <- seq(0, 1, len = 100)
x0Sequence <- seq(0, 1, by = 0.1)

xlist <- lapply(x0Sequence, function(x0){
  p <- c(x0 = x0, K = 0.5, r = 4)
  out <- f(times, p)
  return(out)
})

out <- do.call(rbind, xlist)

P1 <- ggplot(out, aes(x = t, y = x, group = x0, color = x0)) + geom_line()
print(P1)

In addition, vary \(K\) within [0.2,1] and plot the curves in panels separated by their \(K\) values.

times <- seq(0, 1, len = 100)
x0Sequence <- seq(0, 1, by = 0.1)
KSequence <- seq(0.2, 1, len = 5)

x0KComb <- expand.grid(x0 = x0Sequence, K = KSequence)
NComb <- nrow(x0KComb)

xlist <- lapply(1:NComb, function(n){
  p <- c(x0 = x0KComb$x0[n], K = x0KComb$K[n], r = 4)
  out <- f(times, p)
  return(out)
})
out <- do.call(rbind, xlist)

P2 <- ggplot(out, aes(x = t, y = x, group = x0, color = x0)) + facet_wrap(~K) + geom_line()
print(P2)

• Make a grid plot by varying both, parameters \(K\) and \(r\).

times <- seq(0, 1, len = 100)
x0Sequence <- seq(0, 1, by = 0.1)
KSequence <- seq(0.2, 1, len = 5)
rSequence <- seq(1, 10, by = 2)
parComb <- expand.grid(x0 = x0Sequence, K = KSequence, r = rSequence)
NComb <- nrow(parComb)

xlist <- lapply(1:NComb, function(n){
  p <- c(x0 = parComb$x0[n], K = parComb$K[n], r = parComb$r[n])
  out <- f(times, p)
  return(out)
})
out <- do.call(rbind, xlist)

P2 <- ggplot(out, aes(x = t, y = x, group = x0, color = x0)) + facet_grid(r~K) + geom_line()
print(P2)

• Write a function VH <- function(x, p) for Eq. (1) returning a data.frame with columns x, dx.

VH <- function(x, p) {
  with(as.list(p), {
    dx <- r*x*(1-x/K)  
    out <- data.frame(x = x, dx = dx)
    return(out)
  })
}

• Plot \(\dot{x}(t)\) on the interval \(x\in [0,1]\) for \(K=0.75, r=4, x_0 = 1\).

x <- seq(0, 1, len = 100)
p <- c(K = 0.75, r = 4, x0 = 1)
out <- VH(x, p)
P2 <- ggplot(out, aes(x = x, y = dx)) + geom_line()
print(P2)

• Add arrows to the plot along the \(x\)-axis pointing from \(x\) to \(x+\dot x {\rm d}t\) with \({\rm d}t=0.1\).

x2 <- seq(0.01, 1, len = 10) 
out2 <- VH(x2, p)
P2 <- P2 + geom_segment(data = out2, aes(x = x, y = 0, xend = x + 0.1*dx, yend = 0),  arrow = arrow(angle=20, length = unit(0.12, "inches")), color="red")
print(P2)

Exercise 2: Data simulation and error bars

• Simulate data for 10 timepoints in [0,1] generating a data.frame with columns t, x, sigma.

sigma <- 0.05
timesData <- seq(0, 1, len = 10)
p <- c(x0 = 0.1, K = 0.5, r = 4)
outData <- f(timesData, p)
outData$sigma <- sigma 
noise <- rnorm(length(timesData), 0, sigma)
outData$x <- outData$x + noise

• Generate a data.frame for 100 timepoints in [0,1] for \(x(t)\) without noise.

times <- seq(0, 1, len = 100)
outSmooth <- f(times, p)
outSmooth$sigma <- sigma

• Plot the datasets, including data points (geom_point()), error bars (geom_errorbar()), curves \(x(t)\) (geom_line()) and error bands around the curves (geom_ribbon()).

P3 <- ggplot(outSmooth, aes(x = t, y = x, ymin = x - sigma, ymax = x + sigma)) + 
  geom_line() + geom_ribbon(alpha = 0.2) 
P3 <- P3 + geom_point(data = outData) + geom_errorbar(data = outData, width = 0)
print(P3)