Integrating ODEs with deSolve

Install and load the package deSolve. Take a look at the first example of the documentation of ode().

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.

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)
}
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)

# 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)

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)

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)

Cathedral exercise

Close to the outside staircase leading to the top panorama platform there are seven huge figures. During the renovation around 1900 one of the figures was replaced by a figure of the Canon. What about this fact was so funny for a church representative from Konstanz when he visited the inauguration ceremony after the renovation?