The many ways of finding R₀

(all this ecological modeling had a point, after all)

This morning, I gave a 90 minutes introduction to compartment models in epidemiology, and discussed the notion of R₀. Specifically, I talked about how we can reach the same result using very different intuitions, and how all paths lead to the same result.

Let’s take a very simple SIR model:

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

This is your run-of-the-mill model, with an infection term ($\beta$), and a recovery term ($\gamma$). It’s nice, the perfect model when it’s 10:30 am on a Thursday 3 years into a pandemic and you’ll have to calculate determinants live in front of a class.

So how do we usually figure out R₀? Well, we start from the definition, which is to say, the number of infectious cases that are supposed to emerge from the introduction of a single infectious individual in an otherwise fully susceptible population. In a previous post, I talked about the fact that this verbal definition boils down to the product of transmission ($\beta$), contacts ($SI$), and the duration of the infectious period ($\gamma^{-1}$); if we get a fully susceptible population, then we can write $S=N$, and the value of $R_0$ becomes

$$R_0 = \frac{\beta N}{\gamma} \,.$$

OK, but what if I’m an ecologist? First of all, sorry. Been there. Let’s reframe what $R_0$ is: essentially, we have a resident population, which is susceptible. The question we may want to ask is: can it be invaded? Of course our trick of choice for this is to calculate the invasion fitness of our new population (infectious) in the demographic equilibrium of the susceptible population. We will refer to this equilibrium as the “Disease-Free Equilibrium”, or $\text{DFE}$ for short.

What is the invasion fitness of $I$? Quite simply, it is the per capita growth rate of $I$ when initially rare, which we can write as

$$\frac{\dot I}{I} = \beta S - \gamma \,.$$

We want the invasion fitness to be positive for the new population to invade, so basically

$$\left(\beta S - \gamma\right) \biggr |_{\text{DFE}} > 0 \,.$$

Because we are working at the DFE, $S = N$, so this amounts to $\beta N > \gamma$, which you may recognize as

$$\frac{\beta N}{\gamma} > 1 \,.$$

By working out the condition for an infectious population to invade a susceptible equilibrium, we get back to what we wanted: $R_0 > 1$.

OK, but what if I’m an ecologist and want to flex a little? Now we’re talking. Let’s think about $R_0$ some more. That the population of infectious individuals is able to invade means that the DFE is locally unstable, and so, that’s right, it’s Jacobian matrix time. Let’s write the Jacobian matrix for this system - we have three equations, so it’s a $3\times 3$ matrix.

$$\mathbf{J} = \begin{pmatrix} -\beta I & -\beta S & 0\\ \beta I & \beta S - \gamma & 0\\ 0 & \gamma & 0 \end{pmatrix}$$

We want to evaluate $\mathbf{J}$ at the DFE, and because in this model $\dot S + \dot I + \dot R = 0$, we know that $S + I + R$ is a constant, which we call $N$. Problem solved, the DFE is therefore $(N, 0, 0)$.

This gives us the Jacobian matrix for this equilibrium:

$$\begin{pmatrix} 0 & -\beta N & 0\\ 0 & \beta N - \gamma & 0\\ 0 & \gamma & 0 \end{pmatrix}$$

In order for the population of infectious to invade, this equilibrium must be unstable; or in other words, its eigenvalues must have real parts that are positive. How do we get the eigenvalues? But by looking at the characteristic polynomial of this matrix, of course. Let’s rewrite $\textbf{M} = \textbf{J}|_{\text{DFE}}-\lambda \mathbf{I}$ first:

$$\begin{pmatrix} -{\color{red}\lambda} & -\beta N & 0\\ 0 & \beta N - \gamma - {\color{red}\lambda} & 0\\ 0 & \gamma & -{\color{red}\lambda} \end{pmatrix}$$

The $\lambda$ are highlighted; the characteristic polynomical is given by the determinant of this matrix, which at a glance looks bad because we have three rows, so this will probabily involve something cubed, but it actually won’t. Bear with me for a minute.

The determinant of this matrix is

$$(-\lambda)\times \text{Det}{\color{teal}\begin{pmatrix}\beta N - \gamma - \lambda & 0 \\ \gamma & -\lambda\end{pmatrix}} - (-\beta N)\times \text{Det}{\color{purple}\begin{pmatrix}0 & 0 \\ 0 & -\lambda\end{pmatrix}}.$$

Note that the purple matrix determinant is ${\color{purple}-\lambda\times 0 - 0\times 0} = 0$, which very neatly removes the entire right of the above formula. We won’t have no such luck with the left side, but because there is an off-diagonal 0 in the left matrix, its determinant is ${\color{teal}(-\lambda)(\beta N - \gamma - \lambda)}$.

And at this point, I reminded students that hard work (specifically, the type of hard work that’s really not hard but would consist in expandin this expression) is over-rated. Instead, let’s write the characteristic polynomial as (again, $\lambda$ highlighted because it should now be clearer):

$$(-{\color{red}\lambda})\times(-{\color{red}\lambda})\times(\beta N - \gamma - {\color{red}\lambda})\,.$$

Would you look at that? We have three terms that we multiply, and each of them has one root. So when does our equilibrium shifts from stable to unstable? Well, we have $\lambda_3 = 0$ (we read this from the first factor), $\lambda_2 = 0$ (we read that from the second factor), and finally… $\lambda_1 = \beta N - \gamma$. When is this positive? When

$$ \frac{\beta N}{\gamma} > 1 \,. $$

Isn’t that nice? The leading eigenvalue gives has the expression of $R_0$ hidden inside it! So does the invasibility condition! This is how you know something is a big deal: it shows up in multiple ways in your system. I also like this, as a teachable moment: often concepts that don’t appear immediately related are, in fact, the same thing.