Dynamical Systems in Biology

Exercise 1: SIR model

The ODEs of the SIR model, \(\dot{\vec x} = f(\vec x)\), are given by: \[ \begin{align} \dot{S}(t) & = - r \cdot S(t) \cdot I(t), \\ \dot{I}(t) & = r \cdot S(t) \cdot I(t) - a \cdot I(t), \\ \dot{R}(t) & = a \cdot I(t), \end{align} \] where \(\mathbb{R} \ni t \mapsto S(t),I(t),R(t) \in \mathbb{R}\) and the parameters \((a, r) \in \mathbb{R}_{+}^2\) are positive.

Write a function SIR <- function(time, ini, pars) that can be passed as func argument to ode(). Add two functions out2config <- function(out, parameter) and out2phase <- function(out, parameter) that arrange the output of ode() in a config or phase space data.frame with columns:
- time, value, species, parameter, where species \(\in \{S,I,R\}\) and parameter, labelling the parameter choice, are character variables.
- S, I, R, parameter.

library(deSolve)

# SIR ODE
SIR <- function(time, ini, pars) {
  with(as.list(c(ini, pars)), {
    dS <- -r*S*I
    dI <- r*S*I - a*I
    dR <- a*I
    return(list(c(dS, dI, dR)))
  })  
}


## Write result of ode() in config space data.frame
out2config <- function(out, parameter) {
  rbind(data.frame(time = out[,"time"], value = out[,"S"], species = "S", parameter = parameter),
        data.frame(time = out[,"time"], value = out[,"I"], species = "I", parameter = parameter),
        data.frame(time = out[,"time"], value = out[,"R"], species = "R", parameter = parameter)
  )
}

## Write result of ode() in phase space data.frame 
out2phase <- function(out, parameter) {
  data.frame(S = out[,2], I= out[,3], R = out[,4], parameter=parameter)
}

Use ode() to integrate the system for \(r = 0.1\); \(a=50\) and the initial values \(R(t=0) = 0\), \(S(t=0) = 450\), \(N = R(t=0) + S(t=0) + I(t=0) = 1000\).
- Plot solutions for \(S(t=0) = \{450,950\}\) in configuration space.
- Plot solutions for \(r\in [0.1,0.5]\) in phase space.
- Plot solutions for \(r\in [0.1,0.5]\) but \(I(t=0)=1\) in phase space and compare it to the previous plot.

library(ggplot2)
theme_set(theme_bw())

# Parameters and initial values:
N <- 1000
r <- 0.1
a <- 50
R <- 0
S <- 950
I <- N-R-S

# getting the config space data.frames:
yini <- c(S=S, I=I, R=R)
pars <- c(r=r, a=a)
times <- seq(0, 0.3, len=100)
out <- out2config(ode(yini, times, SIR, pars), "S=950")

S <- 450
I <- N-R-S # need to override initial value for I
yini <- c(S=S, I=I, R=R)
pars <- c(r=r, a=a)
out <- rbind(out, out2config(ode(yini, times, SIR, pars), "S=450"))

# config space plot:
P <- ggplot(out, aes(x=time, y=value, group=parameter, color=parameter)) + facet_wrap(~species, scales="free") +  geom_line()
print(P)

# phase space:
times <- seq(0, 0.3, len=1000)
S <- 450
I <- N-R-S
myr <- seq(0.1, 0.5, by=0.1)
yini <- c(S=S, I=I, R=R)

out <- do.call(rbind, lapply(myr, function(ri){
  pars <- c(r=ri, a=a)
  out2phase(ode(yini, times, SIR, pars), paste0("r=", ri))
}))

P <- ggplot(out, aes(x=S, y=I, group=parameter, color=parameter)) + geom_path()
print(P)

# add S(t=0)=1 solution:
I <- 1
S <- N-R-I
myr <- seq(0.1, 0.5, by=0.1)
yini <- c(S=S, I=I, R=R)

outI <- do.call(rbind, lapply(myr, function(ri){
  pars <- c(r=ri, a=a)
  out2phase(ode(yini, times, SIR, pars), paste0("r=", ri))
}))

out$ini <- "550"
outI$ini <- "1"

P <- ggplot(rbind(out,outI), aes(x=S, y=I, group=parameter, color=parameter)) + facet_wrap(~ini, scales = "free_x") + geom_path()
print(P)

How many individuals need to be infected at \(t=0\) such that an epidemic occurs?
- Plot the derivative dI as a function of \(I(t=0)\).
- Find the root of dI using uniroot.

# function to get dI for given I:
deriv <- function(I){
  yini <- c(R=0, I=I, S=N-I)
  pars <- c(a=a, r=r)
  out <- SIR(0, yini, pars)[[1]]
  return(out[2])
}

# plot dI as a function of I:
myI <- seq(0, 1000, by=10)
mydI <- sapply(myI, deriv)
P <- ggplot(data.frame(I = myI, dI = mydI), aes(x = I, y=dI)) + geom_line() + ylab("dI/dt")
print(P)

# finding the root of dI and adding it to the plot:
myroot <- uniroot(deriv, c(20, 1e3))$root
P <- P + geom_point(data = data.frame(I = myroot, dI = 0), color = 2, size = 3) + 
  geom_vline(xintercept = myroot, lty=2) + geom_hline(yintercept = 0, lty=2)
print(P)

Evaluate the severity of the disease as a function of the parameters a and r for \(N = 1000\), \(S(t=0)=950\), \(R(t=0)=0\).
- plot Imax, the maximal value of \(I(t)\), in a 2d grid plot via geom_tile()
- mark the non-epidemic region by adding a line via annotate("line",...)

times <- seq(0, 0.3, len=500)
S <- 950
R <- 0
I <- N-R-S
r <- seq(0.005, 0.5, len=20)
a <- seq(1, 200, len=20)
yini <- c(S=S, I=I, R=R)

pars <- expand.grid(a=a, r=r)

Imax <- unlist(lapply(1:nrow(pars), function(i) max(ode(yini, times, SIR, pars[i,])[,"I"])))

out <- cbind(pars, Imax=Imax)

P <- ggplot(out, aes(x=a, y=r, z=Imax, fill=Imax)) + geom_tile() + scale_fill_gradientn(colours=rainbow(7))
print(P)

# add a line to the plot seperating the epidemic region
P <- P + annotate("line", x=a, y=a/S, lty=2) + 
  annotate("text", x=1.5*mean(a), y=0.5*mean(a/S), label="not epidemic")
print(P)

Simulate the SIR model in configuration space for \(I(t=0) = 50\), \(S(t=0) \in [100,1000]\) and \(r=0.1, a=10\). What difference would you expect for an epidemic in Manhattan (27.500 citizens per km\(^2\)) and Mongolia (1,7 citizens per km\(^2\))?

I <- 50

myS <- seq(100,1000,length.out = 10)
pars <- c(r=0.1, a=10)
times <- seq(0, 0.3, len=1000)

out <- do.call(rbind, lapply(myS, function(Si){
  yini <- c(S=Si, I=I, R=R)
  out2config(ode(yini, times, SIR, pars), paste0("Sini=", round(Si,digits = 2)))
}))

# config space plot:
P <- ggplot(out, aes(x=time, y=value, group=parameter, color=parameter)) + facet_wrap(~species, scales="free") +  geom_line()
print(P)

Dynamical Systems in Biology

Tutorial 3

Daniel Kaschek, Mirjam Fehling-Kaschek

Summer term 2017

Integrating ODEs with deSolve

Exercise 1: SIR model

Cathedral exercise