SIR models with Kermack and McKendrick

Posted byPierre JacobPosted inGeneral
kermit-bake-rolling-pin-egg-preview
Time to bake and to read Kermack and McKendrick!

Hi!

It seems about the right time to read Kermack & McKendrick, 1927, “A contribution to the Mathematical Theory of Epidemics”. It is an early article on the “Susceptible-Infected-Removed” or “SIR” model, a milestone in the mathematical modelling of infectious disease. In this blog post, I will go through the article, describe the model and the data considered by the authors (plague in Bombay in 1905-1906), which will turn out to be a questionable choice. Some references and R code are given at the end of the article. All of this comes with the disclaimer that I have no expertise in epidemiology.

The article starts by crediting Ronald Ross and Hilda Hudson for their articles in the 1910s on malaria; other related works are cited in Anderson (1991) and Bacaër (2011b) (full references are given below). The topic is the following. Some individuals in a closed population get infected by some new disease. Over time they go through various stages, might infect other people, and eventually recover, or die. People who get infected might, in turn, infect other people, and thus the disease can spread to a large part of the population. But it is also possible that the first infected individuals recover before infecting anyone, and the disease could rapidly disappear. The goal here is to try to understand what drives the spread of a disease: why some become large epidemics and some don’t, how many people get affected, etc. The paper is about a model, so “to understand” here means something like “to propose a convincing, simple and intuitive model that is still rich enough to describe faithfully some aspects of reality”.

The authors first explain a general model before delving into special cases including the celebrated SIR model. An infected person eventually recovers or dies, but never gets infected again. The population contains N individuals in a given area where people are in contact with one another, e.g. we can think of a (medieval walled) city. It is really helpful to think of N as indicating a population density in the context of the article, rather than just a population size. What does this mean? It means that the area under consideration remains fixed as we imagine variations in N, thus the density varies linearly with the number of individuals. Let’s assume that at time zero, a single person is infected: y_0 = 1 . The other N-1 individuals are “susceptible” of becoming infected, x_0 = N-1 , and initially no one has yet recovered, z_0 = 0 . Throughout the disease outbreak, the population size remains constant, x_t + y_t + z_t = N .

The general model of the article is first formulated in discrete time. It differentiates infected people according to the duration of their infection. At time t, the number of infected individuals can be written y_t = \sum_{\theta = 0}^t v_{t,\theta}, . Here v_{t,\theta} counts the people who have been infected for \theta units of time already. Over one time interval, for each \theta , we assume that v_{t,\theta} individuals generate \kappa(\theta) x_t v_{t,\theta} new infections and \ell(\theta) v_{t,\theta} removals.

  • Indeed \kappa(\theta) x_t v_{t,\theta} can be understood as the product of a rate of transmission \kappa(\theta) , which accounts for both infectiousness of the pathogen and contact rate between individuals, with x_t v_{t,\theta} counting all pairs of individuals with one susceptible and one infected for \theta time units, i.e. all possible contacts that can lead to a new infection.
  • Meanwhile \ell(\theta) v_{t,\theta} is the product of a rate of recovery/death \ell(\theta) and the number of people infected for \theta time units.

The rates of transmission and recovery, \kappa(\theta) and \ell(\theta) , are generally allowed to vary with the length of infection \theta . Indeed an individual infected for seven days might be more or less infectious than an individual infected for one day, for instance. Counting all the transfers of individuals from and to the different “compartments” (susceptible, infected for one unit of time, infected for two units of time, …, or recovered) the paper gives formulae that describe what happens to these numbers of individuals as time progresses.

From there the authors send the time period to zero. This means that they look at the time period in the eyes and they say: go to zero! Thus they go from discrete to continuous time. The numbers (x_t,y_t,z_t) , representing the number of susceptible, infected and removed individuals, are then shown to follow a system of differential equations. These equations do not have an analytical solution, so there are no explicit formulae giving  (x_t,y_t,z_t) as a function of t. The authors comment on various aspects of the equations, including connections with Volterra integral equations, and Fredholm integral equations. There are also some remarks on limits as time goes to infinity, how these limits depend on the parameters of the model, and how the equations behave for small time t and large population density N.

The celebrated SIR model is obtained as a special case, when the rates of transmission and removal are assumed constant. The system of differential equations becomes:

\begin{cases}\frac{dx}{dt} &= -\kappa x y  \\ \frac{dy}{dt} &= \kappa xy - \ell y \\ \frac{dz}{dt} &= \ell y \end{cases}

Often the letters S,I,R are used in place of x,y,z. Sometimes the model is written in terms of the proportions of individuals of each type, whereas here it is describing the numbers of individuals of each type (per unit area). In my experience, it is easy to get confused by this. Since a product x y appears on the right-hand side, replacing all numbers (x,y,z) by proportions (x/N,y/N,z/N) requires an extra “N” to multiply \kappa on the right-hand side. This is sometimes referred to as “density-dependent versus frequency-dependent”, see e.g. this blog post.

An interesting aspect of the equations is that the sign of the change in numbers of infected individuals, dy/dt , depends on \kappa x_t / \ell being larger or smaller than one. At the start of the epidemic this is very close to \kappa N / \ell . If this is less than one, the number of infected people will decrease; if it is larger than one, it will increase… until  \kappa x_t = \ell and then it will decrease.

Thus what drives the occurrence of an epidemic is 1) the parameters \kappa, \ell , that appear directly in the equations, but also 2) the population density N, in this model. An epidemic occurs or not according to how large N is relative to  \ell / \kappa . There is a “critical threshold” of population density for any \kappa, \ell , above which epidemics occur. In later works, \kappa N / \ell would be called the basic reproduction number “R0”, and epidemics occur when it is larger than one. To illustrate this, here are curves of (x_t/N, y_t/N, z_t/N) against time, all obtained with \kappa = 4\times 10^{-4} , \ell = 0.15 and varying N between 200 and 1200.

test

The graph serves to illustrate several crucial points:

  • the population density N plays a key role in the occurrence of an epidemic under the model,
  • the model is rich enough to generate widespread epidemics or non-epidemics depending on the parameters,
  • contrarily to Kermack & McKendrick, we might not be too concerned about the lack of analytical solutions for the differential equations; we’re used to it and our computers can compute accurate numerical solutions, e.g. using deSolve,
  • it’s fun to play with gganimate.

I was initially surprised about the emphasis on population density in the article. Kermack and McKendrick speculate (towards the end of the article) that epidemics might “regulate” population densities and that many cities in the world might have population densities around the critical threshold (for many pathogens?)… otherwise they would be liable to catastrophic epidemics. Some sentences are quite chilling: “The longer the epidemic is withheld the greater will be the catastrophe, provided that the population continues to increase, and the threshold density remains unchanged”. Does it apply to the current pandemic? Well, if you take the first letter of each of the first thirteen paragraphs of the article, you get “coronavirus”. Or maybe you don’t.

Anderson (1991) provides some useful context: “Two explanations for the termination of an epidemic were most in favour amongst medical circles at that time [circa 1927], namely: (1) that the supply of susceptible people had been exhausted and (2) that during the course of the epidemic the virulence of the infectious agent had gradually (or rapidly) decreased.”

In this debate, the model of Kermack & McKendrick describes an alternative hypothesis: the removal of susceptible people lowers the density below some critical threshold, leading to the termination of the epidemic, even when the number of susceptible individuals might remain large in absolute terms, and even if the virulence of the pathogen remains constant throughout the epidemic.

plaguebombay

Let’s now look at the data set considered by Kermack & McKendrick, shown above. It shows the weekly deaths from the plague during thirty weeks over 1905-1906 in Bombay. Some context is provided by Bacaër (2011), who mentions some interesting concerns about the use of a simple SIR model for this particular outbreak. The plague appeared in Bombay in 1896 and reappeared with “strong seasonal character” in the following years. This seasonal aspect is not accounted for by a simple SIR model, but could play a big role in the decrease of the infections after week ~20.  With the simple SIR model, it is possible to obtain a good fit for the curve dz/dt to the data points, but the associated parameter values are unrealistic. For example you will find a nice fit with \kappa = 8.025\times 10^{-5}, \ell = 5.68, N=75,500 , but the population of Bombay was around one million individuals at that time. Bacaër (2011) proposes a fix: a modified SIR model with seasonal components that provides more satisfactory results. The SIR model seems to provide a useful template for more sophisticated models that account for various other factors, specific to each outbreak, but might not be an adequate model out of the box.

The basic SIR model can be an OK model for certain data sets. An example is the classical “boarding school” data set. This data set was reported in the British Medical Journal in 1978, and concerns an influenza outbreak in a boarding school in the north of England. The data include the number of children confined in bed day after day. There were N=763 boys in that school. The curve y_t of infected individuals can be made to match the data quite closely with \kappa = 2.2\times 10^{-3}, \ell = 0.44 .

fluboarding

Some final thoughts:

  • As illustrated by the two data examples above, measurements about disease outbreaks can be in the form of case counts per time unit, or numbers of individuals “removed”; we might know the size of the susceptible population exactly or not; we might know the exact times of infection, or not, etc. There seem to be as many scenarios as disease outbreaks.
  • In the early works, models are either in discrete or continuous time and they are often deterministic. We can interpret some quantities probabilistically if we want (e.g. the chance that some individual gets infected in some small time interval), but there are no random variables in the description of the model. We can fit curves to data by minimizing least squares, but there are no stochastic processes, likelihood functions or probabilistic models for measurement errors.
  • The choice of example made by Kermack and McKendrick is questionable: their data set would have been better modelled with the consideration of seasonal effects. Clearly that did not prevent the article from becoming extremely influential, and the SIR model from being widely used to this day.
  • Brauer (2005) mentions that “One of the products of the SARS epidemic of 2002-2003 was a variety of epidemic models including general contact rates, quarantine, and isolation.” Hopefully, these developments are proving useful now? What modelling developments will follow from the current epidemic?
  • According to Breda et al. (2012) the original article of Kermack & McKendrick is unfortunately hardly ever read. I have certainly found it useful to read it; certain parts of it, about differential equations, are a bit technical, but it is overall a very well-written article. It’s always interesting to see how pioneers explain their works themselves, and whether the writing style has aged. I was also glad to have, as reading companions, Anderson (1991) and Bacaër (2011).

Here’s a link to an R script producing the above figures and performing quick least squares fit: https://github.com/pierrejacob/statisfaction-code/blob/master/2020-04-sir.R

To read more on the topic:

  • Ronald Ross and Hilda P. Hudson (1917) An Application of the Theory of Probabilities to the Study of a priori Pathometry. Part I, II and III.
  • Roy Anderson (1991) Discussion: The Kermack-McKendrick epidemic threshold theorem. [link]
  • Fred Brauer (2005) The Kermack–McKendrick epidemic model revisited. [link]
  • Nicolas Bacaër (2011a) The model of Kermack and McKendrick for the plague epidemic in Bombay and the type reproduction number with seasonality. [link]
  • Nicolas Bacaër (2011b) A Short History of Mathematical Population Dynamics. [Chapter 16 is on McKendrick and Kermack, link]
  • D. Breda , O. Diekmann , W. F. de Graaf , A. Pugliese & R. Vermiglio (2012) On the formulation of epidemic models (an appraisal of Kermack and McKendrick). [link]
  • H. Hesterbeek et al (2013) Modeling infectious disease dynamics in the complex landscape of global health [a fairly recent review on the topic by some of the world experts https://science.sciencemag.org/content/347/6227/aaa4339]
  • Textbooks on the topic include
    • Bailey (1975) The mathematical theory of infectious diseases and its applications.
    • Anderson and May (1991) Infectious diseases of humans: dynamics and control.
    • Britton and Pardoux (2019) Stochastic Epidemic Models with Inference.

On the Internet:

Published by Pierre Jacob

Professor of statistics, ESSEC Business School

Leave a comment