Getting started with {ringbp}

Introduction

Understanding how an infectious disease might spread, and under what conditions it can be controlled, is crucial for public health planning. The {ringbp} R package provides tools for simulating infectious disease transmission using a branching process model with explicit representation of case detection, isolation, and contact tracing. The package was developed in the early phase of the COVID-19 pandemic to support epidemiological analyses of outbreak dynamics and control feasibility.

This vignette provides a brief introduction to the package and demonstrates how to:

• load {ringbp}
• specify epidemiological parameters
• run outbreak simulations
• plot and summarise simulation results

It concludes with a simplified version of the analysis of COVID-19 control feasibility from Hellewell et al. (2020).

Firstly, we load the {ringbp} package, as well as the {data.table} and {tinyplot} packages which are used in this vignette.

library(ringbp)
library(data.table)
library(tinyplot)
#> Warning: package 'tinyplot' was built under R version 4.4.3

Model overview

At its core, {ringbp} simulates outbreaks using a branching process model. Each infected individual generates secondary infections according to an offspring distribution, with modifications arising from:

• delays between infection, symptom onset, and isolation
• probabilities governing asymptomatic infection and case detection
• optional intervention mechanisms such as quarantine

Simulations are constructed by combining a set of options that describe these components and then running one or more stochastic outbreak realisations.

For a more detailed exposition of the {ringbp} epidemiological model and how non-pharmaceutical interventions influence disease transmission see the ringbp-model.Rmd vignette.

Specifying model components

offspring <- offspring_opts(
  community = \(n) rnbinom(n = n, mu = 2.5, size = 0.16),
  isolated = \(n) rgeom(n = n, prob = 0.8), 
  asymptomatic = \(n) rpois(n = n, lambda = 1.5)
)

delays <- delay_opts(
  incubation_period = \(n) rlnorm(n, meanlog = 0.9, sdlog = 0.5),
  onset_to_isolation = \(n) rgamma(n, shape = 3, scale = 2)
)

event_probs <- event_prob_opts(
  asymptomatic = 0.1,
  presymptomatic_transmission = 0.5,
  symptomatic_traced = 0.2
)

interventions <- intervention_opts(quarantine = FALSE)

sim <- sim_opts(
  cap_max_days = 350,
  cap_cases = 4500
)

Running outbreak simulations

With all components defined, simulations can be run using scenario_sim(). The example below runs 100 independent outbreak simulations starting from a single initial case.

We set the seed to ensure we have the same output each time the vignette is rendered. When using {ringbp}, setting the seed is not required unless you need to simulate the same outbreak multiple times.

set.seed(1)

outbreak <- scenario_sim(
  n = 100,
  initial_cases = 1,
  offspring = offspring,
  delays = delays,
  event_probs = event_probs,
  interventions = interventions,
  sim = sim
)

The returned object is a data.table containing the simulated outbreak trajectories. The columns are:

• sim: simulation replicate ID
• week: outbreak week, zero indexed
• weekly_cases: number of new cases that week
• cumulative: cumulative weekly cases
• effective_r0: the effective reproduction number, the same for each simulation replicate
• cases_per_gen: the number of cases per generation of the outbreak

Visualising results

A simple way to explore results is to plot cumulative cases over time.

tinyplot(
  cumulative ~ week | as.factor(sim), 
  data = outbreak, 
  type = "l", 
  lwd = 3,
  ylab = "Cumulative number of cases", 
  xlab = "Week",
  legend = FALSE,
  theme = "clean"
)

We can also plot the weekly incidence of cases. Note: the maximum number of cases in the simulation is set to 4,500, so the decline in cases to zero is due to reaching that upper bound.

tinyplot(
  weekly_cases ~ week | as.factor(sim), 
  data = outbreak, 
  type = "l", 
  lwd = 3,
  ylab = "Weekly number of cases", 
  xlab = "Week",
  legend = FALSE,
  theme = "clean"
)

Summarising outcomes

A common quantity of interest is the probability that an outbreak goes extinct (i.e. dies out without sustained transmission). This can be calculated using extinct_prob().

extinct_prob(outbreak)
#> Calculating extinction using the extinction status from the simulation.
#> [1] 0.8

By default, extinction is defined as all infectious individuals have had the opportunity to transmit but no new infections are generated. The extinction_week argument in extinct_prob() can also be used to specify whether extinction has happened by a certain week, e.g. extinct_prob(..., extinction_week = 5) checks whether extinction has occurred by week 5 of the outbreak. See ?extinct_prob for other methods of checking for extinction.

Simplified COVID-19 contact tracing effectiveness analysis

Next we demonstrate how to use {ringbp} to reproduce a simplified version of the outbreak control analysis from Hellewell et al. (2020). This study, published in the first months of the COVID-19 pandemic, modelled the probability that an introduced COVID-19 outbreak is contained by isolation and contact tracing. Containment was defined as no disease transmission between weeks 12-16 of the outbreak and the outbreak not reaching 5,000 total cases.

Using the functions described above we define several parameter sets to analyse how varying the basic reproduction number in the community, the proportion of contacts successfully traced, and the number of initial infectors that seed independent outbreaks influence the likelihood of outbreak control.

This is a simplified version of the analysis from Hellewell et al. (2020), for the full analysis see the GitHub repository with the analysis scripts.

Define epidemiological parameters

scenarios <- data.table(
  expand.grid(
    initial_cases = c(5, 20),
    r0_community = c(1.5, 2.5),
    symptomatic_traced = c(0, 0.5, 1)
  )
)
scenarios[, scenario :=  1:.N]
scenarios <- scenarios[, list(data = list(.SD)), by = scenario]

Run the simulation

n <- 10

scenarios[, sims := lapply(data, \(x, n) {
  scenario_sim(
    n = n,
    initial_cases = x$initial_cases,
    offspring = offspring_opts(
      community = \(n) rnbinom(n = n, mu = x$r0_community, size = 0.16),
      isolated = \(n) rnbinom(n = n, mu = 0, size = 1)
    ),
    delays = delay_opts(
      incubation_period = \(n) rweibull(n, shape = 1.65, scale = 4.28),
      onset_to_isolation = \(n) rweibull(n, shape = 1.65, scale = 2.31)
    ),
    event_probs = event_prob_opts(
       asymptomatic = 0.1,
       presymptomatic_transmission = 0.1,
       symptomatic_traced = x$symptomatic_traced
    ),
    interventions = intervention_opts(quarantine = FALSE),
    sim = sim_opts(cap_max_days = 365, cap_cases = 5000)
  )
  },
  n = n
)]

For a more detailed example of running {ringbp} across multiple parameter sets see the parameter-sweep.Rmd vignette, which also includes information on how to parallelise the simulation.

scenarios[, 
  pext := extinct_prob(sims[[1]], extinction_week = 12:16), 
  by = scenario
]
#> Calculating extinction as no new cases within weeks: 12 to 16 (inclusive).
#> Calculating extinction as no new cases within weeks: 12 to 16 (inclusive).
#> Calculating extinction as no new cases within weeks: 12 to 16 (inclusive).
#> Calculating extinction as no new cases within weeks: 12 to 16 (inclusive).
#> Calculating extinction as no new cases within weeks: 12 to 16 (inclusive).
#> Calculating extinction as no new cases within weeks: 12 to 16 (inclusive).
#> Calculating extinction as no new cases within weeks: 12 to 16 (inclusive).
#> Calculating extinction as no new cases within weeks: 12 to 16 (inclusive).
#> Calculating extinction as no new cases within weeks: 12 to 16 (inclusive).
#> Calculating extinction as no new cases within weeks: 12 to 16 (inclusive).
#> Calculating extinction as no new cases within weeks: 12 to 16 (inclusive).
#> Calculating extinction as no new cases within weeks: 12 to 16 (inclusive).

# probability of extinction for each scenario
scenarios$pext
#>  [1] 1.0 0.7 0.7 0.0 1.0 1.0 0.8 0.7 1.0 1.0 1.0 1.0

Hellewell et al. (2020) used {ringbp} to show that:

• Outbreaks are controllable if contact tracing is fast and covers a high proportion of contacts
• Delays in isolating cases or many undetected asymptomatic cases make control unlikely
• Small initial clusters are easier to contain than larger ones

{ringbp} Use Cases

This vignette has focused on the high-level workflow for running simulations with {ringbp}. Further uses of the package include, but are not limited to:

• Simulate outbreaks under different assumptions about \(R_0\), incubation period, and proportion of asymptomatic infections.
• Explore the impact of delays in isolation and the effectiveness of contact tracing.
• Estimate the probability of outbreak extinction under different intervention strategies.