A deeper dive into the SIR model

Semi-organised lecture notes

The SIR model of disease dynamics is foundational for a good reason. It is simple (three equations and two parameters in its simplest formulation), it describes the spread of some diseases well enough (the flu, the plague, for example), and it holds a lot of insights about the fate of a population in which a disease is introduced. In this entry, we will go through a few interesting information that are accessible after a bit of simple manipulation of the differential equations that make up the model.

The SIR model, and the conditions for an outbreak

The SIR model represents the changes in the population size of Susceptible, Infectious, and Recovered (or Removed) individuals. The model itself is written as a system of three differential equations (where $\dot x$ is $dx/dt$):

$$ \dot S = -\beta SI \\ \dot I = \beta SI - \gamma I \\ \dot R = \gamma I $$

The transmission rate is $\beta$, and the recovery rate is $\gamma$ (which can be approximated by the inverse of the contagious period).

In a previous post, I illustrated a number of ways to calculate $\mathcal{R}_0$, the basic reproductive number, which gives the number of expected new infectious cases created by a single infectious individual introduced in an otherwise entirely susceptible population. For this expression of the SIR model,

$$\mathcal{R}_0 = \frac{\beta}{\gamma}S(0) ,$$

where $N = S+I+R$. When $\mathcal{R}_0$ is positive, $\dot I > 0$ (in an otherwise susceptible population), and the disease prevalence increases.

There are a few things we can notice before we start.

One is a little thing that is going to help us simplify some of the expressions down the line: when is $\dot I$ positive? All that is required is $\beta SI - \gamma I > 0$, which assuming $I > 0$ (or thinking about the growht per capita), ends up being $\beta/\gamma \times S > 0$. This is remarkably similar to the expression for $\mathcal{R}_0$, whith $S(0)$ replaced by $S(t)$. We will call this quantity

$$\phi = \frac{\beta}{\gamma}S(t) ,$$

which represents the number of new cases created by a single infectious individual in a population defined by its number of suspectible individuals at a given time. Is there something relevant here? Yes! We can find the value of $S(t)$ for which $\phi < 1$, which gives us the threshold number of susceptible individuals below which the disease prevalence will start to drop, and it is going to be a function of $\beta$ and $\gamma$. We call this threshold $\hat S$, and its value is $\gamma/\beta$.

The second thing to notice is that $\dot S + \dot I + \dot R = 0$ – in other words, if we write $N(t) = S(t)+I(t)+R(t)$, then $\dot N = 0$, and the total population size remains constant. This will help. Later. You’ll see.

How high is the epidemic peak?

A very important question we can answer with the SIR model is: how many people will be infected at once? Having a way to express $I(t)$ directly would be really nice, but there is something far easier and equally useful we can do, which is to find the expression for $I(S)$.

How?

Let’s get back to the system of differential equations for a moment. In $\dot S$ and $\dot I$, the value of $R$ does not appear. This is good news, because we can assume (for now) that the SIR model only has two equations. On one hand, we have the differential w.r.t. $t$ of $I$, and on the other hand, the same thing for $S$. Whatever happens to $R$ is entirely irrelevant here – we can think of $R$ as a compartment in which individuals accumulate after contributing to the epidemic.

In order to figure out the maximum number of simultaneously infectious cases, we need to find the value of $S$ for which $\dot I$ (assumed to be initially positive as a condition for the epidemic to happen!) becomes 0 – this is easy, we have previously established this value as $\hat S = \gamma/\beta$.

The value $I_\text{max}$ is therefore going to be $I(\hat S)$, as soon as we have found a way to express $I(S)$. And this is very straightforward! We know the derivatives of $S$ and $I$, so we can use implicit differentiation to treat $I$ as a variable of $S$:

$$\frac{dI}{dS} = \frac{dI}{dt}/\frac{dt}{dS} .$$

This end up being

$$\frac{dI}{dS} = \frac{I\times(\beta S - \gamma)}{-\beta S I} ,$$

which simplifies really neatly by dropping the $I$, and simplifying the $\beta S$, to

$$\frac{dI}{dS} = \frac{\gamma}{\beta S}-1 .$$

Wait! There is something important here: $\gamma/(\beta S)$ is $\mathcal{R}_0^{-1}$ at $S(0)$. So, if we do the substitution, we have $\mathcal{R}_0^{-1}-1$, which is also $-(1-\mathcal{R}_0^{-1})$, and the term in the parentheses is the herd immunity threshold for this model.

It’s nice when things rhyme.

But we are still no closer to answering our question, which requires an expression of $I(S)$. All we have at the moment is the derivative, which means that we need to take the integral of $dI/dS$, and see where that leads us. The indefinite integral will do just fine (even though we actually know the biological bounds of the system for $S$, namely $0$ and $N$, but who cares?).

So let’s have a little go at this:

$$I(S) = \int \frac{\gamma}{\beta S}-1 \quad\text{d}S .$$

After checking the calculus cheatsheet that is on everyone’s desk at all times (right?), we get the solution to our problem:

$$I(S) = \frac{\gamma}{\beta}\text{ln}S - S + c .$$

There is now an integration constant $c$, and it is going to help us a lot. Let’s define a function $f(S,I)$ that is

$$f(S,I) = I + S - \frac{\gamma}{\beta}\text{ln}S .$$

This is an implicit function, $f(S,I) = c$, and it defines the orbits of the system in the space defined by $S$ and $I$. What is the value of $c$? This is where our little observation ($N(t)$ is a constant of motion) helps, as we know that $f(S,I)$ has the same solution for any points on the orbit. This includes the initial values of our system, and we can write

$$f(S,I) = f(S(0), I(0)) ,$$

which amounts to

$$I + S - \frac{\gamma}{\beta}\text{ln}S = I(0) + S(0) - \frac{\gamma}{\beta}\text{ln}S(0) .$$

Problem solved!

Almost – we need to do a few substitution. We know that $I_\text{max}$ is reached at $\hat S$, and further that $f(\hat S,I_\text{max}) = f(S(0), I(0))$. Oh boy this notation is a mess. Regardless, we have the following relationship:

$$I_\text{max} + \hat S - \frac{\gamma}{\beta}\text{ln}\hat S = I(0) + S(0) - \frac{\gamma}{\beta}\text{ln}S(0) .$$

Here, it is important to note that $\hat S = \gamma/\beta$, so there is a lot we can simplify. Let’s write $\psi = \hat S = \gamma/\beta$, and the expression above becomes

$$I_\text{max} + \psi - \psi\text{ln}\psi = I(0) + S(0) - \psi\text{ln}S(0) ,$$

and now we can start moving terms around, to end up with

$$I_\text{max} = I(0) + S(0) - \psi\text{ln}S(0) - \psi + \psi\text{ln}\psi ,$$

which obviously means that we should gather the factors of $\psi$, to finally end up with

$$I_\text{max} = I(0) + S(0) + \psi \left( \text{ln}\frac{\psi}{S(0)} - 1\right) .$$

There is something really neat hidden here, namely the fact that $\psi/S(0)$ is $\mathcal{R}_0^{-1}$, so we get to write the maximum number of simultaneously infectious cases as

$$I_\text{max} = I(0) + S(0) + \frac{\gamma}{\beta} \left( \text{ln}\frac{1}{\mathcal{R}_0} - 1\right) .$$

How many people will be infected by the end?

Epidemics, like epidemiologists, experience burnout.

Epidemic burnout is essentially the fact that the new cases are going to drop to 0 at some point, and an interesting question is “why”? – specifically, is this because of lack of susceptibles to infect, or lack of infectious?

Believe it or not, we have already solved this problem, because we know that the orbit of out system is given by $f(S,I)$, which has a term in $\text{ln} S$, and so $S > 0$. We can conclude from this that epidemics that follow an SIR model end because of a lack of infectious cases high enough to sustain the infection process. But this is not really a satisfying answer. So let’s get back to where we were before:

$$I + S - \frac{\gamma}{\beta}\text{ln}S = I(0) + S(0) - \frac{\gamma}{\beta}\text{ln}S(0) .$$

Because we have established that $S > 0$ (and because $\dot R$ is always positive or null, so every infectious individual eventually recovers), we know that the point $(S(\infty), 0)$ is somewhere on the orbit (at the end, actually). We can replace these values in the expression above:

$$S(\infty) - \frac{\gamma}{\beta}\text{ln}S(\infty) = I(0) + S(0) - \frac{\gamma}{\beta}\text{ln}S(0) .$$

Let’s make a reasonable assumption: $I(0)$ is going to be so small that it might as well be 0, and $S(0)$ is going to be so not-small that it might as well be called anything else that strikes your fancy, but I don’t want to introduce another variable, so let’s just keep it as is:

$$S(\infty) - \frac{\gamma}{\beta}\text{ln}S(\infty) = S(0) - \frac{\gamma}{\beta}\text{ln}S(0) .$$

The next big-brain move here is to put the $S(…)$ on one side, and the logs on the other:

$$S(0) - S(\infty) = \frac{\gamma}{\beta}\text{ln}S(0) - \frac{\gamma}{\beta}\text{ln}S(\infty) ,$$

which becomes

$$\frac{\beta}{\gamma}\left(S(0) - S(\infty)\right) = \text{ln}S(0) - \text{ln}S(\infty) .$$

This is actually useful because we can expand the l.h.s. like so:

$$\frac{\beta}{\gamma}S(0) - \frac{\beta}{\gamma}\frac{S(0)}{S(0)}S(\infty) ,$$

and we have manifested $\mathcal{R}_0$ out of thin air, so that the expression above becomes

$$\mathcal{R}_0 \left( 1 - \frac{S(\infty)}{S(0)}\right) = \text{ln}S(0) - \text{ln}S(\infty) .$$

This ties the severity of the epidemic (how many individuals will remain uninfected by the end) to the value of $\mathcal{R}_0$, which is neat. But maybe more importantly, we know that $S(\infty) > 0$, and as $\dot S$ is always negative or null, $S(\infty) < S(0)$. This means that the left hand side is always (i) finite and (ii) positive, and is another way of getting at the fact that by the end of the epidemic, some individuals will remain uninfected.

What have we learned?

There is plenty that we can tell about the behavior of the model by not jumping to simulations straight away. The same approach rapidly breaks down for more complex variations of the SIR model (but also, can be applied to look at $R(S)$ and get another persepctive on epidemic burnout), because the complexity increases rapidly; nevertheless, it’s an interesting exercise to do, and can be used as a benchmark to see whether the model implementation is correct!