SIR models of epidemics
A big thank you goes to the Theoretical Biology Group at the ETH Zurich, who openly provide their code from their course module on SIR models. A special thank you goes to Florence Débarre and Sebastian Bonhoeffer who developed this module.
General approach
The SIR model tracks the numbers of susceptible, infected and recovered individuals during an epidemic with the help of ordinary differential equations (ODE). The model can be coded in a few lines in R. We will learn how to simulate the model and how to plot and interpret the results. We will use simulation to verify some analytical results. We will plot time courses, phase diagrams and contour plots.
Basic problem
Infectious diseases are a major cause of death worldwide, and have in the past killed many more people than all the wars (think, for instance, of the Spanish flu). Mathematical modelling of infectious diseases was initiated by Bernoulli in 1760. The work of Kermack and McKendrick, published in 1927, had a major influence on the modelling framework. Their SIR model is still used to model epidemics of infectious diseases. We will study this basic model, and some of its extensions.
Course module handout
Simulation of SIR model
The aim of the following script is to implement the basic SIR model, and plot the simulation results.
Some basic things to know first:
# To make comments
<- Assigning a value to a symbol
c To create a vector c(v1, v2, …, vn)
e is *10^
seq(a,b,c) To create a sequence from a to b with intervals c
Function definition
Input:
t: time (not used here, because there is no explicit time dependence)
x: vector of the current values of all variables (S, I, R)
parms: vector of all parameter values
Output: der: vector of derivatives
To use the lsoda function, the model function has to be a function of t (time), x (the vector of the variables) and parms (the parameters of the model).
# the SIR function calculates the derivatives for the SIR model
SIR <- function(t, x, parms){
# Beta and r are not global variables. This means that if you type beta in R, the output will be 'beta',
# and not its value. You have to specify that you want to use the value of beta from 'parms' to solve the ODEs.
# Similarly, the variables of the model are taken from the vector x. This is done by the 'with' function.
with(as.list(c(parms,x)),{
dS <- - beta*S*I
dI <- + beta*S*I - r*I
dR <- r*I # Note: because S+I+R=constant, this equation could actually be omitted,
# and R at any time point could simply be calculated as N-S-I.
der <- c(dS, dI,dR)
list(der) # the output must be returned
}) # end of 'with'
} # end of function definition
# MAIN PROGRAM
### LOAD LIBRARIES
#load R library for ordinary differential equation solvers
library(deSolve)
### INITIALIZE PARAMETER SETTINGS
parms <- c(beta=1e-3, r=1e-1) # set the parameters of the model
inits <- c(S=499, I=1, R=0) # set the initial values
dt <- seq(0,100,0.1) # set the time points for evaluation
# Calculate and print R_0 on the screen
N <- sum(inits)
R_0 <- with(as.list(parms),{beta*N/r})
print(paste("R_0 =",R_0),quote=FALSE)
### SIMULATE THE MODEL
## Use lsoda to solve the differential equations numerically. The syntax should be
## lsoda(initial values, time points, function, parameters)
simulation <- as.data.frame(lsoda(inits, dt, SIR, parms=parms)) # this way our set 'parms' will be used as default
### PLOT THE OUTPUT
# Plot S according to time, in blue, and add the graph I and R according to time,
# in red and dark green respectively. Call help(plot) for further details.
attach(simulation) # this command allows you to refer to the columns of the data frame directly.
plot(dt, S, type="l", col="blue", ylim=c(0,sum(inits)), xlab="time", ylab="number of individuals",lwd=3)
lines(dt, I, type="l", col="red",lwd=3)
lines(dt, R, type="l", col="darkgreen",lwd=3)
# Add a legend to the graph
legend(70,400, legend=c("S","I","R"), col=c("blue", "red", "darkgreen"), lty=1,lwd=2)
#dev.off()
detach(simulation) # clean up the search path
Interesting questions that you can investigate now
What conditions are necessary for the outbreak of an epidemic?
What fraction of a population is going to be infected in a transient epidemic?
Can partial vaccination in a population protect against the outbreak of an epidemic?
Advanced questions
Modify the model to
allow for the loss of immunity
model treatment of the disease
model the emergence of drug resistance and find the optimal rate of treatment
model longer time scales that allow for the birth and death of individuals.
Glossary
Compartment models: models where the population is divided between several spatial compartments or classes, which are connected by the flow (migration or transformation) of individuals.
SIR models: models where the population is divided into 3 classes – susceptible individuals are uninfected and susceptible to the disease; infected individuals are infected and can infect susceptibles; recovered individuals have recovered from the infection and are immune to re-infection.
Basic reproduction number: the key parameter of epidemiology. It represents the number of secondary cases initiated by the introduction of a single infected individual into a ‘naive’ population.
Phase portrait: also called a phase diagram, it shows the temporal evolution of two or three variables of a system. It consists of trajectories plotted in a coordinate system that has axes corresponding to the variables of the system.
Herd immunity: immunity of a population to the outbreak of an epidemic, provided by the immunity of only a fraction of the individuals.
Literature and weblinks
Kermack, W. and McKendrick, A., 1927. A contribution to the mathematical theory of epidemics. Proc. R. Soc. London A 115, 700-721.
Anderson, R. M. and May, R. M. 1991. Infectious Diseases of Humans. Oxford. Oxford University Press
Earn D. J. et al. 2000. A Simple Model for Complex Dynamical Transitions in Epidemics. Science 287, 667-670.
Matt Keeling’s article in Plus: The mathematics of diseases open access.
Compartment models of epidemiology at Wikipedia.
Reproductive immunity, microbes and insect health 