Modelling infectious diseases

Over the next few weeks, we will learn about how to model the spread of an infectious diseases through the population. We will begin with the simplest model, and slowly add more features to our model to make it more expressive of the real world.

The SIR model

The simplest model of infection is known as the SIR model. The three letters stand for:

  • Susceptible: people who can get infected
  • Infectives: people who are infected, and who can spread the disease (also labelled Infected or Infectious)
  • Removed: people who have either recovered or died from the disease, and cannot catch it again (also labelled Recovered)

The SIR model does not apply to all diseases: for example, you can catch a cold or a flu again after catching it once. In such a case, we would use an SIS model (the last S is Susceptible) instead.

However, for diseases like chicken pox and measles, where getting the disease once offers immunity against getting it again, the SIR model is suitable.

Before reading on, watch the following video (11 minutes) to get a better understanding of the SIR model. Pay particular attention to how they determine the rates of infection ($b$ and $a$ in the video):

In [8]:
from IPython.display import YouTubeVideo
YouTubeVideo('XWXqXzAYe4E', width=800, height=300)

Modelling SIR dynamics

The parameters of the SIR model may be summarized in the following diagram:

$$ S \xrightarrow{\quad \beta\,SI\quad} I \xrightarrow{\quad \gamma\,I \quad} R $$

This diagram includes:

  • state variables: $S,I,R$
  • parameters: $\beta$ and $\gamma$ (the video uses $b$ and $a$)
  • transition rates: $\beta\,SI$ and $\gamma\,I$
  • arrows to indicate the direction of each transition

The state variables $S,I,R$ depend on time $t$, so they will sometimes be written $S_t, I_t, R_t$ (or $S(t), I(t), R(t)$) if we wish to emphasize this. So the diagram should really look like:

In contrast, the parameters $\beta$ and $\gamma$ are constant and do not change with time.

For our purposes, we will assume that time is discrete and takes values $0,1,2,\dots$, rather than being continuous.

Let's look at each of the arrows in the diagram one at a time.

Recovery equations

Let's start with the second arrow: $I \xrightarrow{\gamma I} R$. This tells us that Infectives become Removed at a rate proportional to the number of Infectives.

Question 1:

Which of the following equations describes how $I$ changes over one time step? (ignore any effects from the first arrow)

  1. $I_{t+1} - I_t = \gamma I_t $
  2. $I_{t+1} - I_t = -\gamma I_t $

It might be clearer if you move $I_t$ to the RHS of the equations

Question 2:

Which of the following equations describes how $R$ changes over one time step?

  1. $R_{t+1} - R_t = \gamma I_t $
  2. $R_{t+1} - R_t = -\gamma I_t $

Infection equations

Now let's look at the first arrow $S \xrightarrow{\beta SI} I$. This is a bit more complicated.

Question 3:

Which of the following equations describes how $I$ changes over one time step? (ignore any effects from the second arrow)

  1. $I_{t+1} - I_t = \beta S_t I_t $
  2. $I_{t+1} - I_t = -\beta S_t I_t $

Question 4:

Which of the following equations describes how $S$ changes over one time step?

  1. $S_{t+1} - S_t = \beta S_t I_t $
  2. $S_{t+1} - S_t = -\beta S_t I_t $

Question 5:

Putting it all together, complete the following equations:

  1. $S_{t+1} - S_t = ?$
  2. $I_{t+1} - I_t = ?$ (this is the tricky one)
  3. $R_{t+1} - R_t = ?$

Question 6:

Using your answers to question 5, replace the 3 question marks in the function update(S,I,R) with the appropriate code. You may just write out those 3 lines when submitting your answer (instead of the whole code).

In [59]:
# Parameters
beta = 0.6
gamma = 0.2

# Initial state (represented as fraction of the total population)
S = 0.99
I = 1-S
R = 0

# Total population
N = 10000

# Number of iterations
iters = 50

# Function to update state
def update(S,I,R):
    newS = ?
    newI = ?
    newR = ?
    return (newS, newI, newR)

# Run the simulation and print the numbers at every time step

print("t", "\t", "S", "\t", "I", "\t", "R")

for t in range(iters):
    print(t, "\t", int(S*N), "\t", int(I*N), "\t", int(R*N))
    S,I,R = update(S,I,R)
    if S < 0: # If number of susceptibles goes below zero, stop the simulation

To think about

  • Play around with the constants beta and gamma (and also N, iters and the initial S). How do these values affect what happens in the long run?
  • This model assumes that everyone comes into contact with everyone else with equal probability, which is not very realistic. What might you want to include into the model to account for this?
  • How might we model the use of quarantines in curbing the spread of a disease? What about the use of vaccines?

Next time, we will try to answer some of these questions. We will also learn how to make nice graphs showing the progress of infection over time - this makes it easier to quickly see what's happening, instead of having to look at a table of numbers.