Introduction

In this lesson, we delve into the effective reproduction number, Re, a key epidemiological measure that gauges the current potential of an infectious disease to spread within a population. Unlike R0, which is calculated under the assumption of a fully susceptible population, Re takes into account the real-world dynamics of immunity and intervention impacts. Importantly, Re is time-dependent, reflecting changes in the population and intervention measures over time.

Learning Objectives

  1. Define the effective reproduction number Re for infectious diseases.
  2. Calculate Re for a given infectious disease model.
  3. Derive the Re expression from SIR model ODEs.
  4. Understand the implications and limitations of Re.

1 Overview of Reproduction Numbers

1.1 Basic Reproduction Number (R0)

R0, the basic reproduction number, signifies the average number of secondary infections that an infected individual will produce in a completely susceptible population. It is a foundational concept in epidemiology for assessing the potential severity of infectious diseases at the outset of an outbreak.

For the SIR model with frequency-dependent transmission, R0 is calculated with the formula:

\[ R_0 = \frac{\beta}{\gamma} \]

where \(\beta\) is the transmission rate, and \(\gamma\) is the recovery rate. Initially, R0 offers a baseline measure of an epidemic’s potential to expand, but it becomes less relevant as individuals gain immunity through infection or vaccination.

1.2 Effective Reproduction Number (Re)

Re, or the effective reproduction number, adapts the concept of R0 to reflect the ongoing epidemic conditions. It estimates the number of secondary infections that occur at any given time, taking into account the proportion of the population that remains susceptible and the effects of control measures like vaccination.

Re is inherently time-dependent, changing as the epidemic progresses and as interventions are implemented.

1.3 Transitioning from R0 to Re

An outbreak initially might have an R0 of 3, suggesting rapid spread. However, as individuals recover and become immune, the pool of susceptibles diminishes, thereby reducing the reproduction number. For instance, once a third of the population have become infected and immune, an infectious person who would have infected three others at the beginning of an outbreak now ‘wastes’ one of their transmissions on an already immune person, and therefore only infects two others. Similarly, once 2/3 of the population have become infected and immune, this will have dropped to 1. The outbreak has reached its peak. From our definition, when R=1, the outbreak does not grow further nor does it decay (recall, we need R>1 for growth or R<1 for decay). However, this point of the outbreak, the peak with R=1, is somewhat instantaneous as, now, more than 2/3 of the population have become infected or immune which causes R to drop below 1. Thus, the outbreak starts to decline, with R continuing to go down as fewer and fewer susceptibles remain. This illustrates the time-dependent nature of Re.

This shows the impact of increasing immunity on the disease’s ability to spread over time.

Practice Exercise: SIR Model Simulation

Let’s simulate the example described in words above, using real R code. The first step is to choose a modelling framework. For this question, we will be using the classic SIR model with three compartments and a closed population (N stays constant).

As a reminder, the compartments and parameters of the SIR model are the following:

  • Susceptible (S): Number of individuals who can become infected (total population at risk).
  • Infected (I): Number of individuals who are infected and can transmit the disease.
  • Recovered (R): Number of individuals who have recovered and are no longer infectious.
  • Transmission coefficient (β): The rate at which susceptible individuals become infected.
  • Recovery rate (γ): The rate at which infected individuals recover. This is the inverse of the infectious period.

Step 0: Load the packages needed for SIR modelling in R

This code has been written for you:

# Load packages 
if(!require(pacman)) install.packages("pacman")
pacman::p_load(deSolve, # to solve differential equations
               tidyverse) # includes dplyr and ggplot2

Step 1: Setting up Initial Conditions and Parameters

The initial conditions for this question are:

  • Total population (N): 120
  • Initial infected (I): 1 (patient zero)
  • Initial susceptible (S): 119 (120 total - 1 infected)
  • Initial recovered (R): 0 (none recovered at the start)
# Define the initial conditions
initial_state <- "WRITE_YOUR_CODE_HERE"

initial_state

Next, let’s calculate the values of \(\beta\) and \(\gamma\)

  1. If the duration of the infectious period is 2 days. Calculate the value of \(\gamma\).

  2. The \(R_0\) for this disease is 3. If we assume frequency-dependent transmission for this disease, calculate the value of \(\beta\).

Make the duration of the simulation 1 month (30 days), with a time interval of 0.5 days.

# Time span for the simulation
times <- seq(0, 30, by = 0.5)

Step 2: Writing the SIR Model Functions

Now that we’ve set up our initial conditions and parameters, we’ll implement the SIR model using the {deSolve} package.

Start by defining a function called sir_model with differential equations for S, I, and R. Note that we are modelling frequency-dependent transmission in a closed population.

As a reminder, here are the three components of the function:

# Define function and arguments
sir_model <- function(times, state, parameters) {
  with(as.list(c(state, parameters)), {
    N <- "WRITE_YOUR_CODE_HERE"
    # Rate of change
    dS <- "WRITE_YOUR_CODE_HERE"
    dI <- "WRITE_YOUR_CODE_HERE"
    dR <- "WRITE_YOUR_CODE_HERE"
    
    # Return the rate of change
    list(c("WRITE_YOUR_CODE_HERE"))
  })   # end with (as.list...) 
}

Step 3: Solving the Differential Equations

Use the ode() function from the {deSolve} package to solve the differential equations of your SIR model.

# Solve the ordinary differential equations
sir_out <- ode(
  # inital state
  "WRITE_YOUR_CODE_HERE",
  # sequence of time points
  "WRITE_YOUR_CODE_HERE",
  # model function
  "WRITE_YOUR_CODE_HERE", 
  # parameters in the model
  "WRITE_YOUR_CODE_HERE"
)

Convert the output to a data frame for easier manipulation and visualization

# Convert to data frame
sir_out <- "WRITE_YOUR_CODE_HERE"

sir_out

Step 4: Visualize the results

Create a plot using {ggplot2} to graph all three predicted lines (S, I, R).

# Create a plot to visualize the dynamics of the SIR model
ggplot(sir_out, aes(x = time)) +
  # Graph Susceptible line
  geom_line(aes(y = S, color = "Susceptible"), size = 1) +
  # Graph Infected line
  geom_line(aes(y = I, color = "Infectious"), size = 1) +
  # Graph Recovered line
  geom_line(aes(y = R, color = "Recovered"), size = 1) +
  # Add title and labels
  labs(title = "SIR Model Simulation, R0 of 3", x = "Time (days)") +
  # Specify color for each line
  scale_color_manual(values = c("Susceptible" = "#2F5597", "Infectious" = "#C00000", "Recovered" = "#548235")) +
  # Select theme 
  theme_light() +
  scale_x_continuous(expand = c(0,0)) +
  scale_y_continuous(expand = c(0,0))

2 Defining Re

The effective reproductive number (Re) is the average number of secondary cases per infectious case in a population.

Re Terminology and Notation

The effective reproduction number (Re) is sometimes denoted Reff or simply R. As with R0, the symbol R stands for reproduction.

It can also be referred to as Rt, defined as the average number of secondary cases per infectious case at time \(t\) in the epidemic.

2.1 Expression for Re

Re can be expressed as a modified version of R0, adjusted for the current number of susceptible individuals:

\[ R_e = R_0 \times \frac{S}{N} \]

  • \(R_0\) is the basic reproduction number.
  • \(S\) is the current number of susceptible individuals.
  • \(N\) is the total population size.

This formula underscores that Re decreases as the number of susceptible individuals, \(S\), declines relative to the total population, \(N\). Since \(S\) changes over time, Re is also time-dependent.

3 Derivation of \(R_e\) from SIR ODEs

3.1 Calculating Re for Density-Dependent Transmission

Theoretical Framework

Given the SIR model dynamics where:

  • S = Susceptible individuals
  • I = Infected individuals
  • R = Recovered individuals

The rates of change in these compartments are governed by the following ordinary differential equations (ODEs):

\[ \frac{dS}{dt} = -\beta S I, \quad \frac{dI}{dt} = \beta S I - \gamma I, \quad \frac{dR}{dt} = \gamma I \]

Here, \(\beta\) represents the transmission rate per contact and \(\gamma\) represents the recovery rate.

Let’s consider a scenario where the disease is at a point where it has already infected some of the population, and thus not everyone is susceptible. To determine Re:

  1. Start with the expression for \(\frac{dI}{dt}\):

Given \(\frac{dI}{dt} = \beta S I - \gamma I\), for the disease to persist in the population, \(\frac{dI}{dt}\) must be greater than zero.

  1. Rearrange to highlight the infection terms:

\[ \beta S I - \gamma I > 0 \] \[ \beta S I > \gamma I \] Dividing all terms by \(I\) (assuming \(I > 0\)): \[ \beta S > \gamma \]

  1. Isolate the susceptible ratio:

Divide by \(\gamma\) to find the condition under which the infection grows:

\[ \frac{\beta S_t}{\gamma} > 1 \]

NOTE HERE: THIS IS Re.

  1. Relate to R0:

Recognizing that \(R_0 = \frac{\beta N}{\gamma}\) for density-dependent transmission, substitute to express in terms of \(S\) and \(N\):

\[ \frac{\beta S}{\gamma} = \frac{\beta N}{\gamma} \times \frac{S}{N} = R_0 \times \frac{S}{N} \]

Thus, we have:

\[ R_e = R_0 \times \frac{S}{N} \]

3.2 Calculating Re for Frequency-Dependent Transmission

For frequency-dependent transmission, the contact rate is adjusted for the population size, thereby making the transmission rate depend not just on the number of contacts an infected individual has, but also on the proportion of those contacts that are with susceptible individuals. This model is particularly relevant in large, mixed populations where the probability of contact between any two individuals is low.

In the case of frequency-dependent transmission, the transmission term is modified to \(\frac{\beta S I}{N}\), reflecting that each infectious individual makes contact with a fraction of the population that is proportional to the total number of susceptible individuals divided by the total population size. The ODEs are:

\[ \frac{dS}{dt} = -\frac{\beta S I}{N}, \quad \frac{dI}{dt} = \frac{\beta S I}{N} - \gamma I, \quad \frac{dR}{dt} = \gamma I \]

Given the modified infection term, Re can be derived similarly by analyzing the growth condition of \(\frac{dI}{dt}\):

  1. From the equation for \(\frac{dI}{dt}\):

\[ \frac{\beta S I}{N} - \gamma I > 0 \]

  1. Isolate the transmission terms:

\[ \frac{\beta S I}{N} > \gamma I \]

Dividing by \(I\) and rearranging:

\[ \frac{\beta S}{N} > \gamma \]

  1. Express in terms of R0:

For frequency-dependent transmission, if we consider \(R_0\) under initial conditions where \(S \approx N\), then:

\[ R_0 = \frac{\beta}{\gamma} \]

Thus, substituting back:

\[ R_e = R_0 \times \frac{S}{N} \]

This calculation reaffirms that Re, like R0, diminishes as the proportion of susceptible individuals in the population decreases, but does so in a manner that is adjusted for the total population, which is crucial in densely populated or highly mixed environments.

3.3 Implications of Re

How does Re reflect the progression of an outbreak?

Re > 1

Each infected patient will go on to infect on average more than one person. Number of cases may increase exponentially over time.

Re = 1

Each infected patient will go on to infect one other person. Statistically, this would suggest that the outbreak is stabilizing, with no changes in number of infected people over time.

Re < 1

Each infected patient will go on to infect on average less than one person, meaning fewer people are being infected which suggests the outbreak is slowing down over time.

Practice Exercise: Estimate \(R_0\) from SIR Curves

Hint: Look at the proportion of S at the epidemic peak.

4 Conclusion

Re is time-dependent and effectively captures the current dynamics of disease transmission, reflecting changes in population susceptibility. It is a critical measure for understanding the control measures’ impact on an outbreak and determining strategies for intervention.


5 Answer Key

5.1 SIR Model Simulation

# Load packages 
if(!require(pacman)) install.packages("pacman")
pacman::p_load(deSolve, # package to solve differential equations
               tidyverse, # package that includes ggplot2
               here)
# Define the parameters
parameters <- c(beta = 1.5, # Transmission rate
                gamma = 0.5)   # Recovery rate

initial_state <- c(S = 119, I = 1, R = 0)


# Time span for the simulation
times <- seq(0, 30, by = 0.5)

sir_model <- function(times, state, parameters) {
  with(as.list(c(state, parameters)), {
    N <- S+I+R
    # Rate of change
    dS <- -beta * S * I / N
    dI <- (beta * S * I / N) - gamma * I
    dR <- gamma * I
    
    # Return the rate of change
    list(c(dS, dI, dR))
  })   # end with (as.list...) 
}

sir_out <- ode(
  # inital state
  y = initial_state,
  # sequence of time points
  times = times,
  # model function
  func = sir_model, 
  # parameters in the model
  parms = parameters
)

sir_out <- as.data.frame(sir_out)

sir_out
# Create a plot to visualize the dynamics of the SIR model
ggplot(sir_out, aes(x = time)) +
  # Graph Susceptible line
  geom_line(aes(y = S, color = "Susceptible"), size = 1) +
  # Graph Infected line
  geom_line(aes(y = I, color = "Infectious"), size = 1) +
  # Graph Recovered line
  geom_line(aes(y = R, color = "Recovered"), size = 1) +
  # Add title and labels
  labs(title = "SIR Model Simulation, R0 of 3", x = "Time (days)") +
  # Specify color for each line
  scale_color_manual(values = c("Susceptible" = "#2F5597", "Infectious" = "#C00000", "Recovered" = "#548235")) +
  # Select theme 
  theme_light() +
  scale_x_continuous(expand = c(0,0)) +
  scale_y_continuous(expand = c(0,0))
sir_out <- sir_out %>% mutate(sir_out,
                   Prop = S/(S+I+R))

max_I_prop <- sir_out$Prop[which.max(sir_out$I)]
max_I_prop

1/max_I_prop

5.2 Estimate \(R_0\) from SIR Curves

\(R_0 = 2\) because \(S/N = 0.5\) when \(R_e = 1\). The calculation steps can be done as follows:

\[ R_e = R_0 \times \frac{S}{N} \]

\[ 1 = R_0 \times 0.5 \]

\[ R_0 = \frac{1}{0.5} = 2 \]