Projecting re-emergence risk after waning or new births

Modelling the initial epidemic

We follow the approach described in the “Modelling uncertainty in R₀” vignette, which includes uncertainty in final size calculations by running final_size() for values drawn from a distribution, and summarising the outcomes.

First we define the future window of time we are considering, as well as the rate of waning (assumed to be 5% per year) and influx of new susceptibles (assuming 14 births per 1,000 per year).

We also define a distribution of \(R_0\) values for the initial epidemic and the second later re-emerging outbreak, both following a normal distribution \(N(\mu = 2, \sigma = 0.3)\).

# future time period considered
years_range <- seq(0, 20)

# proportion losing immunity per year
annual_wane <- 0.05

# proportion newly susceptible per year
birth_rate <- 14 / 1000

# define r0 values
r0_mean <- 2.0
r0_samples <- rnorm(n = 100, mean = r0_mean, sd = 0.3) # for first outbreak
r0_samples2 <- rnorm(n = 100, mean = r0_mean, sd = 0.3) # for second outbreak

Next we simulate initial epidemics using final_size(), iterating over \(R_0\) values sampled above.

# run final size for each value and collate estimates
# use `Map()` to iterate over values and an index to add to the final data
final_size_est <- Map(
  function(r0, i) {
    fs <- final_size(r0)$p_infected # run final size model, get estimate only
    data.frame(fs = fs, i = i, r0 = r0) # data.frame for each value considered
  },
  r0 = r0_samples,
  i = seq_along(r0_samples)
)

# use Reduce to combine the list of data.frames into a single data.frame
final_size_df <- Reduce(rbind, final_size_est)

Modelling waning immunity and influx of susceptibles

We calculate the relative population immunity, accounting for the waning of immunity from the initial epidemic, and accounting for the influx (due to new births) of new susceptibles over a range of 20 years.

This is calculated as \(I_y = (1 - \text{annual waning})^y \times (1 - b y)\), where \(b\) is the birth rate; this is vectorised over the range of years in the code.

# define changes in susceptibility over time from waning
# and population turnover (assuming) rectangular age distribution
relative_immunity <- (1 - annual_wane)^(years_range) *
  pmax(1 - birth_rate * years_range, 0) # use pmax in case of longer year range

Projected estimates of effective \(R\)

Having simulated the distribution of initial epidemic sizes based on the distribution of \(R_0\), we next implement waning and an influx of new susceptibles to calculate the level of immunity over coming years, relative to the level of immunity at the end of the epidemic. This is a step-wise process:

• The initial final size calculations have given us estimates of the population’s ‘relative immunity’ in year 0 (i.e., the final size; assuming all recovered individuals are immune);

• This allows us to calculate the relative susceptibility from the relative immunity, as \(\text{susceptibility} = 1 - I_y\); this can be calculated for any year \(y\) using the predicted immunity for that year \(I_y\);

• The effective reproduction number (i.e. \(R\) for a partially susceptible population) can then be calculated as \(R = R_0 \times \text{susceptibility}\) for each year;

• We can then predict \(R\) over the coming years as population susceptibility increases.

When \(R\) rises above the critical value of 1, there is potential for sustained transmission in the population and hence another large epidemic.

We calculate the level of susceptibility - and hence \(R\) - over the distribution of \(R_0\) values we have assumed for a re-emerging epidemic.

# calculate effective R for each value and collate estimates
r_eff_est <- Map(
  function(r0, fs) {
    r_eff <- (1 - fs * relative_immunity) * r0

    # data.frame for each value considered
    data.frame(r_eff = r_eff, yr = years_range)
  },
  r0 = r0_samples2,
  fs = final_size_df$fs
)

# use Reduce to combine the list of data.frames into a single data.frame
r_eff_df <- Reduce(rbind, r_eff_est)

In this scenario, the majority of new susceptibility comes from waning, rather than influx of new births (i.e. 5% vs 1.4% per year). As a result, the predicted value of \(R\) rises substantially in the years following the initial epidemic, suggesting potential for re-emergence within 5-10 years.

# plot results
ggplot(r_eff_df) +
  geom_hline(yintercept = 1, linetype = "dashed", color = "gray") +
  geom_vline(xintercept = 0, colour = "red", linetype = "dashed") +
  annotate(
    geom = "text",
    x = 0, y = 2, angle = 90, vjust = "outward",
    colour = "red",
    label = "Initial epidemic at year 0"
  ) +
  stat_summary(
    aes(
      yr, r_eff
    ),
    fun = mean,
    fun.min = function(x) {
      quantile(x, 0.025)
    },
    fun.max = function(x) {
      quantile(x, 0.975)
    }
  ) +
  scale_y_continuous(
    limits = c(0, 3)
  ) +
  theme_classic() +
  coord_cartesian(
    expand = TRUE
  ) +
  labs(
    x = "Years after initial epidemic",
    y = "Effective R (mean and 95% CI)"
  )

Figure 1: Predicted value of the effective reproduction number \(R\) following an initial epidemic of an emerging infection in year 0. We assume 5% of the immune population becomes susceptible again each year following the epidemic, and there is an influx of new susceptibles via a birth rate of 14/1000 per year. We also assume \(R_0\) has a mean of 2, with a standard deviation around this estimate of 0.3.