A simple model for vaccination (w4 SIR)

1 Etivity

In this etivity, we are going back to the simple SIR model from week 3, without births or deaths, to look at the effect of vaccination. The aim of this etivity is to represent vaccination in a very simple way - we are assuming it already happened before we run our model! By changing the initial conditions, we can prepare the population so that it has received a certain coverage of vaccination. Even though this is very simplistic, it will allow us to study some important effects of vaccination on the infection dynamics.

We are starting with the transmission and recovery parameters beta = 0.4 days\(^{-1}\) and gamma = 0.1 days \(^{-1}\). To incorporate immunity from vaccination in the model, we assume that a proportion p of the total population starts in the recovered compartment, representing the vaccine coverage and assuming the vaccine is perfectly effective. Again, we assume the epidemic starts with a single infected case introduced into the population.

Model this scenario for a duration of 2 years, assuming that the vaccine coverage is 50%, and plot the prevalence in each compartment over time. Confirm that you observe an epidemic peaking at around 125 days after introduction of the infectious person in the population. Also have a look at the proportion susceptible and recovered, to double-check this looks like what you would expect given your initial conditions.

1.1 Question: Does everyone in the population need to be vaccinated in order to prevent an epidemic? What do you observe if you model the infection dynamics with different values for p? Can you explain why?

Hopefully, you have seen that once a certain proportion of the population is immune, no epidemic occurs. This is called the critical vaccination threshold or herd immunity threshold.

Now, use your code above to investigate how the herd immunity threshold changes if we are modelling a disease with a different infection and recovery rate.

1.2 Question: What proportion of the population needs to be vaccinated in order to prevent an epidemic if beta = 0.4 and gamma = 0.2 days\(^{-1}\)? What if beta = 0.6 and gamma = 0.1 days\(^{-1}\)?

As you can see, the proportion of the population that needs to be vaccinated varies with different infection-related parameters. Think about why that is so in the context of herd immunity, and what is different between the scenarios you have just modelled.

1.3 Question: Remember that vaccination changes the effective reproduction number, by reducing the number of people who are susceptible. Based on your answers to the previous questions, can you use the formula for the effective reproduction number R_eff to derive a formula for calculating the critical vaccination threshold?

2 Solutions - A Simple Model for Vaccination

Modelling a disease where \(\beta\) = 0.4 days\(^{-1}\), \(\gamma\) = 0.1 days \(^{-1}\) and the vaccine coverage p = 0.5

#<img src="../../IDM1/Graphics and Data/w4_nb1_model_diagram.png">

# LOAD THE PACKAGES:
library(deSolve)
library(reshape2)
library(ggplot2)

# MODEL INPUTS:

# Vaccine coverage
p <- 0.5

# Total population size
N <- 10^6

# Vector storing the initial number of people in each compartment (at timestep 0)
initial_state_values <- c(S = (1-p)*(N-1),   # a proportion 1-p of the total population is susceptible
                          I = 1,             # the epidemic starts with a single infected person
                          R = p*(N-1))       # a proportion p of the total population is vaccinated/immune

# Vector storing the parameters describing the transition rates in units of days^-1
parameters <- c(beta = 0.4,      # the infection rate, which acts on susceptibles
                gamma = 0.1)     # the rate of recovery, which acts on those infected

# TIMESTEPS:

# Vector storing the sequence of timesteps to solve the model at
times <- seq(from = 0, to = 730, by = 1)   # from 0 to 730 days in daily intervals

# SIR MODEL FUNCTION: 

# The model function takes as input arguments (in the following order): time, state and parameters
sir_model <- function(time, state, parameters) {  

    with(as.list(c(state, parameters)), {  # tell R to look for variable names within the state and parameters objects    
        
    # Calculating the total population size N (the sum of the number of people in each compartment)
      N <- S+I+R
      
    # Defining lambda as a function of beta and I:
      lambda <- beta * I/N
        
    # The differential equations
      dS <- -lambda * S               # people move out of (-) the S compartment at a rate lambda (force of infection)
      dI <- lambda * S - gamma * I    # people move into (+) the I compartment from S at a rate lambda, 
                                      # and move out of (-) the I compartment at a rate gamma (recovery)
      dR <- gamma * I                 # people move into (+) the R compartment from I at a rate gamma
      
    # Return the number of people in the S, I and R compartments at each timestep 
    # (in the same order as the input state variables)
    return(list(c(dS, dI, dR))) 
    })
  
}

# MODEL OUTPUT (solving the differential equations):

# Solving the differential equations using the ode integration algorithm
output <- as.data.frame(ode(y = initial_state_values, 
                            times = times, 
                            func = sir_model,
                            parms = parameters))

output_long <- melt(as.data.frame(output), id = "time")                  # turn output dataset into long format

# Adding a column for the prevalence proportion to the long-format output
output_long$prevalence <- output_long$value/sum(initial_state_values)

# Plot the prevalence proportion
ggplot(data = output_long,                                               # specify object containing data to plot
       aes(x = time, y = prevalence, colour = variable, group = variable)) +  # assign columns to axes and groups
  geom_line() +                                                          # represent data as lines
  xlab("Time (days)")+                                                   # add label for x axis
  ylab("Prevalence (proportion)") +                                      # add label for y axis
  labs(colour = "Compartment",                                           # add legend title
       title = "Prevalence of infection, susceptibility and recovery over time")   # add plot title    

Increasing the vaccine coverage to 75%

# MODEL INPUTS:

# Vaccine coverage
p <- 0.75

# Vector storing the initial number of people in each compartment (at timestep 0)
initial_state_values <- c(S = (1-p)*(N-1),   # a proportion 1-p of the total population is susceptible
                          I = 1,             # the epidemic starts with a single infected person
                          R = p*(N-1))       # a proportion p of the total population is vaccinated/immune

# MODEL OUTPUT (solving the differential equations):

# Solving the differential equations using the ode integration algorithm
output <- as.data.frame(ode(y = initial_state_values, 
                            times = times, 
                            func = sir_model,
                            parms = parameters))

output_long <- melt(as.data.frame(output), id = "time")                  # turn output dataset into long format

# Adding a column for the prevalence proportion to the long-format output
output_long$prevalence <- output_long$value/sum(initial_state_values)

# Plot the prevalence proportion
ggplot(data = output_long,                                               # specify object containing data to plot
       aes(x = time, y = prevalence, colour = variable, group = variable)) +  # assign columns to axes and groups
  geom_line() +                                                          # represent data as lines
  xlab("Time (days)")+                                                   # add label for x axis
  ylab("Prevalence (proportion)") +                                      # add label for y axis
  labs(colour = "Compartment",                                           # add legend title
       title = "Prevalence of infection, susceptibility and recovery over time")   # add plot title    

2.0.1 Does everyone in the population need to be vaccinated in order to prevent an epidemic? What do you observe if you model the infection dynamics with different values for p? Can you explain why?

No, not everyone in the population needs to be vaccinated in order to prevent an epidemic. In this scenario, if p equals 0.75 or higher, no epidemic occurs - 75% is the critical vaccination/herd immunity threshold. Remember, as you heard in week 2, herd immunity describes the phenomenon in which there is sufficient immunity in a population to interrupt transmission. Because of this, not everyone needs to be vaccinated to prevent an outbreak.

2.0.2 What proportion of the population needs to be vaccinated in order to prevent an epidemic if \(\beta\) = 0.4 and \(\gamma\) = 0.2 days\(^{-1}\)? What if \(\beta\) = 0.6 and \(\gamma\) = 0.1 days\(^{-1}\)?

If \(\beta\) = 0.4 and \(\gamma\) = 0.2 days\(^{-1}\), the herd immunity threshold is 50%. If \(\beta\) = 0.6 and \(\gamma\) = 0.1 days\(^{-1}\), the required vaccination coverage is around 83%.

2.0.3 Remember that vaccination changes the effective reproduction number, by reducing the number of people who are susceptible. Based on your answers to the previous questions, can you use the formula for the effective reproduction number R_eff to derive a formula for calculating the critical vaccination threshold?

In mathematical modelling terms, herd immunity is just the same as saying that R_eff < 1. We can derive the herd immunity threshold by solving the formula for R_eff for p when R_eff = 1:

\[\begin{align} R_{eff} & = R_{0} \frac{S}{N} \\ R_{eff} & = R_{0} (1-p) \\ p = 1-\frac{1}{R_{0}} \end{align}\]

Remember, we can calculate R₀ by dividing \(\beta\) by \(\gamma\).