The Basic Reproduction Number

The Basic Reproduction Number

The basic reproduction number, \(R_0\), is defined as the expected number of secondary cases produced by a single (typical) infection in a completely susceptible population. It is important to note that \(R_0\) is a dimensionless number and not a rate, which would have units of \(\mathrm{time}^{-1}\). Some authors incorrectly call \(R_0\) the “basic reproductive rate.”"

We can use the fact that \(R_0\) is a dimensionless number to help us in calculating it.

\[ R_0 \propto \left( \frac{{\rm infection}}{{\rm contact}} \right) \cdot \left( \frac{{\rm contact}}{{\rm time}}\right) \cdot \left( \frac{{\rm time}}{{\rm infection}}\right) \]

More specifically:

\[ R_0 = \tau \cdot \bar{c} \cdot d, \]

where \(\tau\) is the transmissibility (i.e., probability of infection given contact between a susceptible and infected individual), \(\bar{c}\) is the average rate of contact between susceptible and infected individuals, and \(d\) is the duration of infectiousness.

SIR Model

\[ \begin{aligned} \dot{S} &= -\beta SI \\ \dot{I} &= \beta SI - \nu I \\ \dot{R} &= \nu I \end{aligned} \]

where \(\dot{S}\) indicates the time-derivative of \(S\) (and similarly for \(\dot{I}\) and \(\dot{R}\)), \(\beta = \tau \bar{c}\) and is known as the effective contact rate, and \(\nu\) is the removal rate.

By assumption all rates are constant. This means that the expected duration of infection is simply the inverse of the removal rate: \(d = \nu^{-1}\).

What are the conditions for an epidemic? An epidemic occurs if the number of infected individuals increases, i.e., \(\dot{I} > 0\)

\[ 0 > \beta SI - \nu I \] \[ \frac{\beta S}{\nu} > 1 \] At the outset of an epidemic, nearly everyone (except the index case) is susceptible. So we can say that \(s \approx 1\). Substituting \(s=1\), we arrive at the following inequality:

\[ \frac{\beta}{\nu} = R_0 > 1 \]

Since \(\beta = \tau \bar{c}\) and \(d = \nu^{-1}\), we see that we have derived our expression for \(R_0\) given in the first section. This little bit of mathematical trickery explains why we have that cumbersome phrase “in a completely susceptible population”" tacked onto our definition for \(R_0\).

Next Generation Matrix: Inutition

If \(R_0\) is the number of secondary infections produced by a single typical infection in a rarefied population, how do we define it when there are multiple types of infected individuals. For example, what is a typical infection in a vector-borne disease like malaria? What about a sexually transmitted infection where there are large asymmetries in transmissibility (like HIV)? Or what about a multi-host pathogen like influenza?

It turns out that there is a straightforward extension of the theory for structured epidemic models. The mathematics behind this theory is not especially difficult, but it does involve scary German terms that are not familiar to the non-engineers in our midst. The key concept is that we now need to average the expected number of new infections over all possible infected types.

Assume that we have a system in which there are multiple discrete types of infected individuals (e.g., mosquitoes and humans; women and men; or humans, dogs, and chickens). We define the next generation matrix as the square matrix \(\mathbf{G}\) in which the \(ij\)th element of \(\mathbf{G}\), \(g_{ij}\), is the expected number of secondary infections of type \(i\) caused by a single infected individual of type \(j\), again assuming that the population of type \(i\) is entirely susceptible. That is, each element of the matrix \(\mathbf{G}\) is a reproduction number, but one where who infects whom is accounted for.

Once we have \(\mathbf{G}\), we are one step away from \(R_0\). The basic reproduction number is given by the spectral radius of \(\mathbf{G}\). The spectral radius is the also known as the dominant eigenvalue of \(\mathbf{G}\). The next generation matrix has a number of desirable properties from a mathematical standpoint. In particular, it is a non-negative matrix and, as such, it is guaranteed that there will be a single, unique eigenvalue which is positive, real, and strictly greater than all the others. This is \(R_0\).

For illustrative purposes, we will limit our discussion to the case where there are two classes of infected individual. The next generation matrix is thus \(2 \times 2\). For the \(2 \times 2\) matrix, there is a fairly simple formula for the eigenvalues that everyone who took algebra in high school should know. Define

\[ \mathbf{G} = \left[ \begin{array}{cc} a & b \\ c & d \end{array} \right] \]

The eigenvalues of \(\mathbf{G}\) are:

\[ \lambda_{\pm} = \frac{T}{2} \pm \sqrt{(T/2)^2 - D} \]

where \(T=a+d\) is the trace and \(D=ad-bc\) is the determinant of matrix \({\bf G}\).

Say you have a sexually transmitted disease in a completely heterosexual population. Define \(f\) as the expected number of infected women and \(m\) as the expected number of infected men given contact with a single infected member of the opposite sex in a completely susceptible population. The next generation matrix is

\[ \mathbf{G} = \left[ \begin{array}{cc} 0 & f \\ m & 0 \end{array} \right] \]

\(R_0\) is thus \(\sqrt{mf}\). It is worth noting that this is the geometric mean of the expected number of female and male secondary cases.

When you have more than two discrete types of infectious classes, the formula for calculating the eigenvalues is not as simple, but the same ideas apply. In practice, we typically calculate eigenvalues numerically using software like R.

(G <- matrix(c(0,2,1,0), nr=2, nc=2, byrow=TRUE))

##      [,1] [,2]
## [1,]    0    2
## [2,]    1    0

ev <- eigen(G)
ev$values[1]

## [1] 1.414214

A More Formal Approach

A review paper by Heffernan et al. (2005) provides a nice readable introduction for calculating \(R_0\) in structured population models. The notation I use here follows their usage.

Consider the next generation matrix \(\mathbf{G}\). It is comprised of two parts: \(F\) and \(V^{-1}\), where

\[ F = \left[ \frac{\partial F_i(x_0)}{\partial x_j} \right] \] and

\[ V = \left[ \frac{\partial V_i(x_0)}{\partial x_j} \right] \]

These are square matrices of the partial derivatives of new infections (\(F_i\)) and transfers between different compartments (\(V_i\)). The rank of these matrices is the number of distinct classes of infections. For example, in a simple two-sex SIR model, the next generation matrix would be \(2 \times 2\) since there are two classes of infection (i.e., female and male). \(x_0\) is the disease-free equilibrium state. This matrix should be non-negative, irreducible, and primitive.

\(R_0\) is then the dominant eigenvalue of the matrix \(\mathbf{G}=FV^{-1}\).

Example: SEIR Epidemic Consider a Susceptible-Exposed-Infected-Removed (SEIR) Epidemic. This is an appropriate model for a disease where there is a considerable post-infection incubation period in which the exposed person is not yet infectious.

State diagram for the SEIR model. \(\beta\) is the effective contact rate, \(\lambda\) is the “birth” rate of susceptibles, \(\mu\) is the mortality rate, \(k\) is the progression rate from exposed (latent) to infected, \(\gamma\) is the removal rate.

The simple SEIR model consists of a set of four differential equations:

\[ \begin{aligned} \dot{S} &= -\beta S I + \lambda - \mu S \\ \dot{E} &= \beta SI - (\mu + k) E \\ \dot{I} &= kE - (\gamma + \mu) I \\ \dot{R} &= \gamma I - \mu R \end{aligned} \]

where \(\beta\) is the effective contact rate, \(\lambda\) is the “birth” rate of susceptibles, \(\mu\) is the mortality rate, \(k\) is the progression rate from exposed (latent) to infected, \(\gamma\) is the removal rate.

To calculate the next generation matrix for the SEIR model, we need to enumerate the number of ways that (1) new infections can arise and (2) the number of ways that individuals can move between compartments. There are two disease states but only one way to create new infections:

\[ F = \left( \begin{array}{ll} \frac{\beta \lambda }{\mu } & 0 \\ 0 & 0 \end{array} \right) \label{eq:seir-F} \]

In contrast, there are various ways to move between the states:

\[ V = \left( \begin{array}{ll} 0 & k+\mu \\ \gamma +\mu & -k \end{array} \right) \]

\(R_0\) is the leading eigenvalue of the matrix \(FV^{-1}\). This is reasonably straightforward to calculate since \(FV^{-1}\) is simply a \(2\times2\) matrix.

\[ R_0 = \frac{k \beta \lambda }{\mu (k+\mu ) (\gamma +\mu )} \]

It is interesting to note that \(R_0\) is also the product of the rate of production of (1) new exposures and (2) new infections, as it should be.

What Is A Generation?

In demography, \(R_0\) represents the ratio of total population size from the start to the end of a generation, which is, roughly, the mean age of childbearing. \(R_0 = e^{rT}\), where \(r\) is the instantaneous rate of increase of the population. So what is a generation in an epidemic model? Generations in epidemic models are the waves of secondary infection that flow from each previous infection. So, the first generation of an epidemic is all the secondary infections that result from infectious contact with the index case, who is of generation zero. If \(R_i\) denotes the reproduction number of the \(i\)th generation, then \(R_0\) is simply the number of infections generated by the index case, i.e., generation zero. Now, these numbers are typically small and are therefore susceptible to random sampling error. Consequently, we talk about expected (i.e., averaged over many epidemics) numbers of secondary cases produced by generation zero.

The following figure shows a schematic representation of an epidemic. The zeroth generation of the epidemic is the index case – patient zero – indicated in red. The number of secondary infections generated by the case in generation zero is \(R_0=3\). In the first generation (blue), \(R_1=6/3=2\). The second generation (cyan) has \(R_1=12/6=2\).

Graphical depiction of generations in an epidemic.

Will An Epidemic Infect Everyone?

Will an epidemic, once it has taken off in a population, eventually infect everyone? In order to answer this question, we want to know how \(I\) changes with respect to the “fuel” for the epidemic, \(S\).

We thus divide equation for \(\dot{I}\) by equation for \(\dot{S}\).

\[ \frac{dI}{dS} = -1 + \frac{\nu}{\beta S} \]

We solve this equation by first multiplying both sides by \(dS\)

\[dI = (-1 + \frac{\nu}{\beta S})dS \]

We then integrate and do a little algebra, yielding

\[ \log(s_{\infty}) = R_0 (s_{\infty}-1) \]

This is an implicit equation for \(s_{\infty}\), the number of susceptibles at the end of the epidemic. When \(R_0 > 1\), this equation has exactly two roots, only one of which lies in the interval \((0,1)\). Subtract \(R_0 (s_{\infty}-1)\) from both sides and we get \(\log(s_{\infty}) - R_0 (s_{\infty}-1)=0\). Call the whole left-hand side \(y\). \(y\) will have different values for different values of \(\log(s_{\infty})\). Only a couple of those will satisfy the final-size equation. \(s_{\infty}=1\) will always satisfy the requirement of \(y=0\) (plug it in and see!). When \(R_0>1\), the other solution to \(y=0\) is the actual value of the final size. This is the one we really care about. If \(R_0<1\), the only value that satisfies the final-size equation is \(\log(s_{\infty})=0\). In words, at the end of the epidemic, everyone will still be susceptible (i.e., no one gets infected).

fs <- function(ro,si) log(si) - ro*(si-1) 
si <- seq(0,1,length=500)
ro <- c(3,2,1.1,0.9)
plot(si,fs(ro[1],si), type="l", xlab="Fraction Susceptible", ylab = "y")
for(i in 2:3) lines(si,fs(ro[i],si))
lines(si,fs(ro[4],si), col="red")
abline(h=0,lty=2)

The figure shows the solutions of the final-size equation for various values of \(R_0 > 1\) in black. The point where the curve crosses the horizontal axis is the value for \(s_{\infty}\), the total fraction of the population infected at the end of the epidemic. As \(R_0\) gets larger, the final size of the epidemic gets larger as well. The figure also shows the solution when \(R_0 < 1\) in red. The curve never crosses the horizontal axis, meaning that essentially none of the total population becomes infected when an infection is sub-critical.

Optimal Virulence: Pathogen Life History Evolution

But enough about you, let’s talk about me for a while. It’s instructive to think about epidemics from the pathogen’s perspective. Pathogens bear biological information in their nucleic acids. This information varies from one copy of a pathogen to another, and the ability of a pathogen to persist and multiply can be a function of this variability. We therefore have fulfilled the necessary and sufficient conditions for natural selection. Pathogens evolve.

We will consider a model in which transmissibility and disease-induced mortality trade-off introduced by van Baalen and Sabelis (1995). This interaction is mediated by virulence, which we will take to be proportional to parasite burden or parasitemia. More copies of a virus (say) means that conditional on contact with an infected individual, the pathogen is more likely to be transmitted. However, more viral copies means the host is sicker – and potentially dead. Dead hosts do not transmit and very sick hosts are less likely to be up and interacting with susceptibles.

Consider a directly-transmitted infection from which there is no recovery (e.g., HSV-2). The population experiences a baseline mortality rate, \(\mu\), and a disease-induced mortality \(\delta\).

\[ \begin{aligned} \dot{S} & = -\beta SI - \mu S \\ \dot{I} & = \beta SI - (\mu + \delta) I, \end{aligned} \]

It is straightforward to show that the basic reproduction number is given by:

\[ R_0 = \frac{\beta}{\mu+\delta} > 1, \]

The parameter \(\mu\) is independent of the epidemic, but the parameters, \(\beta\) and \(\delta\) can conceivably be functions of virulence, which we denote by \(x\).

An Evolutionary Stable Strategy (ESS) is a phenotype that can not be invaded by a rare mutant. Loosely speaking, it represents the optimal phenotype. The ESS virulence occurs where \(dR_0/dx = 0\).

Differentiate the equation for \(R_0\) with respect to \(x\) using the quotient rule:

\[ \frac{dR_0}{dx} = \frac{\beta'(\mu + \delta) - \delta' \beta}{(\mu + \delta)^2} = 0 \]

Rearranging and evaluating \(\beta(x)\) and \(\delta(x)\) at the ESS values of \(x\) (denoted \(x^*\)), we have:

\[ \frac{d\beta(x)}{d\delta(x)} = \frac{\beta(x^*)}{\mu+\delta(x^*)} \]

This result has a nice geometric interpretation. The ESS virulence occurs where a line rooted at the origin is tangent to the curve that relates \(\beta\) to \(\delta\). This result is known as the Marginal Value Theorem and has applications in economics and ecology as well as epidemiology. The MVT model for optimal virulence is plotted in the following figure. In the lower curve, the tangent line hits further out on the horizontal axis and mortality is higher.

Marginal value theorem for optimal virulence. The ESS virulence occurs where a line rooted at the origin is tangent to the curve that relates \(\beta\) to \(\delta\). Two curves are depicted. The first curve shows a pathogen in which transmissibility increases relatively rapidly with mortality. Point \(A\) indicates the optimal balance between \(\beta(x)\) and \(\delta(x)\) under this case, and the optimal virulence is indicated \(x^*\). For the second curve, relative transmission is less efficient. Therefore, the tangent line from the origin to the curve hits further out (\(B\)) along the mortality axis and the optimal virulence is higher (\(x^{**}\)).

Frank (1996) has a very thorough review of models of the evolution of virulence.

Can a More Virulent Strain Invade?

Following Diekman (2004), we use adaptive dynamics to explore this problem.