diff --git a/NEWS.md b/NEWS.md
index 228aab65..0a4161fe 100644
--- a/NEWS.md
+++ b/NEWS.md
@@ -1,3 +1,8 @@
+# epichains 0.1.9999
+
+## Documentation
+* A vignette outlining how to simulate interventions has been added
+
 # epichains 0.1.0
 
 ## Package name change
diff --git a/R/intervention.R b/R/intervention.R
deleted file mode 100644
index 3ab7339c..00000000
--- a/R/intervention.R
+++ /dev/null
@@ -1,51 +0,0 @@
-#' Reduce the mean of the offspring distribution
-#'
-#' @description
-#' `intvn_reduce_mean()` is a helper for the \code{simulate_*()} functions. It
-#' reduces/scales the mean of the offspring distribution in order to
-#' mimic the impact of a population-level intervention. Currently, it can only
-#' handle the poisson and negative binomial distributions and errors when other
-#' offspring distributions are specified alongside the `intvn_mean_reduction`
-#' argument.
-#'
-#' @inheritParams simulate_tree
-#' @param intvn_mean_reduction A number between 0
-#' and 1 for scaling/reducing the mean of `offspring_dist`. Serves as
-#' population-level intervention. `intvn_mean_reduction` = 0
-#' implies no intervention impact and `intvn_mean_reduction` = 1 implies full
-#' impact.
-#' @param pars_list Parameter(s) for poisson or negative binomial offspring
-#' distribution.
-#' @return List of the offspring distribution parameter(s) with the mean
-#' scaled by \code{1 - intvn_mean_reduction}.
-#' @details
-#' `intvn_reduce_mean()` scales the mean of the offspring distribution
-#' by \eqn{1 - {\sf intvn\_mean\_reduction}} so that the new mean is given as:
-#' \deqn{(1 - {\sf intvn\_mean\_reduction}) \times {\sf mean,}} This
-#' scaling when applied to the poisson and negative binomial offspring
-#' distributions corresponds to the population-level reduction of R0 as
-#' described in Lloyd-Smith et al, (2005). `intvn_reduce_mean()` is therefore
-#' only implemented for the aforementioned distributions and errors when other
-#' offspring distributions are specified along with the `intvn_mean_reduction`
-#' argument in the \code{simulate_*()} functions.
-#'
-#' @author James M. Azam
-#' @keywords internal
-#' @references Lloyd-Smith, J., Schreiber, S., Kopp, P. et al. Superspreading
-#' and the effect of individual variation on disease emergence. Nature 438,
-#' 355–359 (2005). \doi{10.1038/nature04153}
-intvn_reduce_mean <- function(intvn_mean_reduction, offspring_dist, pars_list) {
-  # Intervention only works for pois and nbinom
-  if (!offspring_dist %in% c("pois", "nbinom")) {
-    stop(
-      "`offspring_dist` must be one of c(\"pois\", \"nbinom\"), ",
-      "if intvn_mean_reduction is specified."
-    )
-  }
-  if (offspring_dist == "pois") {
-    pars_list$lambda <- (1 - intvn_mean_reduction) * pars_list$lambda
-  } else {
-    pars_list$mu <- (1 - intvn_mean_reduction) * pars_list$mu
-  }
-  return(pars_list)
-}
diff --git a/R/simulate.r b/R/simulate.r
index 11f34afe..c5b81d1b 100644
--- a/R/simulate.r
+++ b/R/simulate.r
@@ -1,6 +1,5 @@
 #' Simulate transmission trees from an initial number of infections
 #'
-#' @inheritParams intvn_reduce_mean
 #' @param nchains Number of chains to simulate.
 #' @param offspring_dist Offspring distribution: a character string
 #' corresponding to the R distribution function (e.g., "pois" for Poisson,
@@ -98,19 +97,6 @@
 #'   serials_dist = function(x) 3,
 #'   lambda = 2
 #' )
-#'
-#' # Run model with intervention a 50% reduction in R0.
-#' chains_with_intvn <- simulate_tree(
-#'   nchains = 10,
-#'   statistic = "size",
-#'   offspring_dist = "pois",
-#'   intvn_mean_reduction = 0.5,
-#'   stat_max = 10,
-#'   serials_dist = function(x) 3,
-#'   lambda = 2
-#' )
-#'
-#' chains_with_intvn
 #' @references
 #' Lehtinen S, Ashcroft P, Bonhoeffer S. On the relationship
 #' between serial interval, infectiousness profile and generation time.
@@ -127,7 +113,6 @@
 #' 1186–1204. \doi{https://doi.org/10.3390/ijerph7031204}
 simulate_tree <- function(nchains, statistic = c("size", "length"),
                           offspring_dist, stat_max = Inf,
-                          intvn_mean_reduction = 0,
                           serials_dist, t0 = 0,
                           tf = Inf, ...) {
   statistic <- match.arg(statistic)
@@ -137,13 +122,6 @@ simulate_tree <- function(nchains, statistic = c("size", "length"),
   # check that offspring is properly specified
   check_offspring_valid(offspring_dist)
 
-  # Check that the intvn_mean_reduction is well specified
-  checkmate::assert_number(
-    intvn_mean_reduction,
-    lower = 0,
-    upper = 1
-  )
-
   # check that offspring function exists in base R
   roffspring_name <- paste0("r", offspring_dist)
   check_offspring_func_valid(roffspring_name)
@@ -151,15 +129,6 @@ simulate_tree <- function(nchains, statistic = c("size", "length"),
   # Gather offspring distribution parameters
   pars <- list(...)
 
-  # Prepare interventions if specified
-  if (intvn_mean_reduction > 0) {
-    pars <- intvn_reduce_mean(
-      intvn_mean_reduction = intvn_mean_reduction,
-      offspring_dist = offspring_dist,
-      pars_list = pars
-    )
-  }
-
   if (!missing(serials_dist)) {
     check_serial_valid(serials_dist)
   } else if (!missing(tf)) {
@@ -284,7 +253,6 @@ simulate_tree <- function(nchains, statistic = c("size", "length"),
 #' Simulate transmission chains sizes/lengths
 #'
 #' @inheritParams simulate_tree
-#' @inheritParams intvn_reduce_mean
 #' @param stat_max A cut off for the chain statistic (size/length) being
 #' computed. Results above the specified value, are set to `Inf`.
 #' @inheritSection simulate_tree Calculating chain sizes and lengths
@@ -303,22 +271,9 @@ simulate_tree <- function(nchains, statistic = c("size", "length"),
 #'   stat_max = 10,
 #'   lambda = 2
 #' )
-#'
-#' # Run model with intervention a 50% reduction in R0.
-#' chain_summary_with_intvn <- simulate_summary(
-#'   nchains = 10,
-#'   statistic = "size",
-#'   offspring_dist = "pois",
-#'   intvn_mean_reduction = 0.5,
-#'   stat_max = 10,
-#'   lambda = 2
-#' )
-#'
-#' chain_summary_with_intvn
 #' @export
 simulate_summary <- function(nchains, statistic = c("size", "length"),
                              offspring_dist,
-                             intvn_mean_reduction = 0,
                              stat_max = Inf, ...) {
   statistic <- match.arg(statistic)
 
@@ -327,13 +282,6 @@ simulate_summary <- function(nchains, statistic = c("size", "length"),
   # check that offspring is properly specified
   check_offspring_valid(offspring_dist)
 
-  # Check that the intvn_mean_reduction is well specified
-  checkmate::assert_number(
-    intvn_mean_reduction,
-    lower = 0,
-    upper = 1
-  )
-
   # check that offspring function exists in base R
   roffspring_name <- paste0("r", offspring_dist)
   check_offspring_func_valid(roffspring_name)
@@ -341,15 +289,6 @@ simulate_summary <- function(nchains, statistic = c("size", "length"),
   # Gather offspring distribution parameters
   pars <- list(...)
 
-  # Prepare interventions if specified
-  if (intvn_mean_reduction > 0) {
-    pars <- intvn_reduce_mean(
-      intvn_mean_reduction = intvn_mean_reduction,
-      offspring_dist = offspring_dist,
-      pars_list = pars
-    )
-  }
-
   # Initialisations
   stat_track <- rep(1, nchains) ## track length or size (depending on `stat`)
   n_offspring <- rep(1, nchains) ## current number of offspring
@@ -404,7 +343,6 @@ simulate_summary <- function(nchains, statistic = c("size", "length"),
 #' population
 #'
 #' @inheritParams simulate_tree
-#' @inheritParams intvn_mean_reduction
 #' @param pop The susceptible population size.
 #' @param offspring_dist Offspring distribution: a character string
 #' corresponding to the R distribution function (e.g., "pois" for Poisson,
@@ -465,21 +403,9 @@ simulate_summary <- function(nchains, statistic = c("size", "length"),
 #' size = 1.1,
 #' serials_dist = function(x) 3
 #' )
-#'
-#' # Simulate with negative binomial offspring with intervention (50%
-#' # reduction in R0)
-#' simulate_tree_from_pop(
-#'   pop = 100,
-#'   offspring_dist = "nbinom",
-#'   intvn_mean_reduction = 0.5,
-#'   mu = 0.5,
-#'   size = 1.1,
-#'   serials_dist = function(x) 3
-#' )
 #' @export
 simulate_tree_from_pop <- function(pop,
                                    offspring_dist = c("pois", "nbinom"),
-                                   intvn_mean_reduction = 0,
                                    serials_dist,
                                    initial_immune = 0,
                                    t0 = 0,
@@ -487,25 +413,9 @@ simulate_tree_from_pop <- function(pop,
                                    ...) {
   offspring_dist <- match.arg(offspring_dist)
 
-  # Check that the intvn_mean_reduction is well specified
-  checkmate::assert_number(
-    intvn_mean_reduction,
-    lower = 0,
-    upper = 1
-  )
-
   # Gather offspring distribution parameters
   pars <- list(...)
 
-  # Prepare interventions if specified
-  if (intvn_mean_reduction > 0) {
-    pars <- intvn_reduce_mean(
-      intvn_mean_reduction = intvn_mean_reduction,
-      offspring_dist = offspring_dist,
-      pars_list = pars
-    )
-  }
-
   if (offspring_dist == "pois") {
     ## Use a right truncated poisson distribution
     ## to avoid more cases than susceptibles
diff --git a/_pkgdown.yml b/_pkgdown.yml
index cc8cc70d..a84a41a0 100644
--- a/_pkgdown.yml
+++ b/_pkgdown.yml
@@ -10,6 +10,7 @@ articles:
   navbar: Package vignettes
   contents:
   - projecting_incidence
+  - interventions
 - title: Modelling guides and background
   navbar: Modelling guides and background
   contents:
diff --git a/man/intvn_reduce_mean.Rd b/man/intvn_reduce_mean.Rd
deleted file mode 100644
index 080ec2a9..00000000
--- a/man/intvn_reduce_mean.Rd
+++ /dev/null
@@ -1,55 +0,0 @@
-% Generated by roxygen2: do not edit by hand
-% Please edit documentation in R/intervention.R
-\name{intvn_reduce_mean}
-\alias{intvn_reduce_mean}
-\title{Reduce the mean of the offspring distribution}
-\usage{
-intvn_reduce_mean(intvn_mean_reduction, offspring_dist, pars_list)
-}
-\arguments{
-\item{intvn_mean_reduction}{A number between 0
-and 1 for scaling/reducing the mean of \code{offspring_dist}. Serves as
-population-level intervention. \code{intvn_mean_reduction} = 0
-implies no intervention impact and \code{intvn_mean_reduction} = 1 implies full
-impact.}
-
-\item{offspring_dist}{Offspring distribution: a character string
-corresponding to the R distribution function (e.g., "pois" for Poisson,
-where \code{\link{rpois}} is the R function to generate Poisson random
-numbers).}
-
-\item{pars_list}{Parameter(s) for poisson or negative binomial offspring
-distribution.}
-}
-\value{
-List of the offspring distribution parameter(s) with the mean
-scaled by \code{1 - intvn_mean_reduction}.
-}
-\description{
-\code{intvn_reduce_mean()} is a helper for the \code{simulate_*()} functions. It
-reduces/scales the mean of the offspring distribution in order to
-mimic the impact of a population-level intervention. Currently, it can only
-handle the poisson and negative binomial distributions and errors when other
-offspring distributions are specified alongside the \code{intvn_mean_reduction}
-argument.
-}
-\details{
-\code{intvn_reduce_mean()} scales the mean of the offspring distribution
-by \eqn{1 - {\sf intvn\_mean\_reduction}} so that the new mean is given as:
-\deqn{(1 - {\sf intvn\_mean\_reduction}) \times {\sf mean,}} This
-scaling when applied to the poisson and negative binomial offspring
-distributions corresponds to the population-level reduction of R0 as
-described in Lloyd-Smith et al, (2005). \code{intvn_reduce_mean()} is therefore
-only implemented for the aforementioned distributions and errors when other
-offspring distributions are specified along with the \code{intvn_mean_reduction}
-argument in the \code{simulate_*()} functions.
-}
-\references{
-Lloyd-Smith, J., Schreiber, S., Kopp, P. et al. Superspreading
-and the effect of individual variation on disease emergence. Nature 438,
-355–359 (2005). \doi{10.1038/nature04153}
-}
-\author{
-James M. Azam
-}
-\keyword{internal}
diff --git a/man/simulate_summary.Rd b/man/simulate_summary.Rd
index 9c1283af..a453d680 100644
--- a/man/simulate_summary.Rd
+++ b/man/simulate_summary.Rd
@@ -8,7 +8,6 @@ simulate_summary(
   nchains,
   statistic = c("size", "length"),
   offspring_dist,
-  intvn_mean_reduction = 0,
   stat_max = Inf,
   ...
 )
@@ -29,12 +28,6 @@ corresponding to the R distribution function (e.g., "pois" for Poisson,
 where \code{\link{rpois}} is the R function to generate Poisson random
 numbers).}
 
-\item{intvn_mean_reduction}{A number between 0
-and 1 for scaling/reducing the mean of \code{offspring_dist}. Serves as
-population-level intervention. \code{intvn_mean_reduction} = 0
-implies no intervention impact and \code{intvn_mean_reduction} = 1 implies full
-impact.}
-
 \item{stat_max}{A cut off for the chain statistic (size/length) being
 computed. Results above the specified value, are set to \code{Inf}.}
 
@@ -105,18 +98,6 @@ simulate_summary(
   stat_max = 10,
   lambda = 2
 )
-
-# Run model with intervention a 50\% reduction in R0.
-chain_summary_with_intvn <- simulate_summary(
-  nchains = 10,
-  statistic = "size",
-  offspring_dist = "pois",
-  intvn_mean_reduction = 0.5,
-  stat_max = 10,
-  lambda = 2
-)
-
-chain_summary_with_intvn
 }
 \seealso{
 \itemize{
diff --git a/man/simulate_tree.Rd b/man/simulate_tree.Rd
index cb3d0629..1d6722dc 100644
--- a/man/simulate_tree.Rd
+++ b/man/simulate_tree.Rd
@@ -9,7 +9,6 @@ simulate_tree(
   statistic = c("size", "length"),
   offspring_dist,
   stat_max = Inf,
-  intvn_mean_reduction = 0,
   serials_dist,
   t0 = 0,
   tf = Inf,
@@ -36,12 +35,6 @@ numbers).}
 computed. Results above the specified value, are set to this value.
 Defaults to \code{Inf}.}
 
-\item{intvn_mean_reduction}{A number between 0
-and 1 for scaling/reducing the mean of \code{offspring_dist}. Serves as
-population-level intervention. \code{intvn_mean_reduction} = 0
-implies no intervention impact and \code{intvn_mean_reduction} = 1 implies full
-impact.}
-
 \item{serials_dist}{The serial interval distribution function; the name
 of a user-defined named or anonymous function with only one argument \code{n},
 representing the number of serial intervals to generate. See details.}
@@ -128,19 +121,6 @@ chains <- simulate_tree(
   serials_dist = function(x) 3,
   lambda = 2
 )
-
-# Run model with intervention a 50\% reduction in R0.
-chains_with_intvn <- simulate_tree(
-  nchains = 10,
-  statistic = "size",
-  offspring_dist = "pois",
-  intvn_mean_reduction = 0.5,
-  stat_max = 10,
-  serials_dist = function(x) 3,
-  lambda = 2
-)
-
-chains_with_intvn
 }
 \references{
 Lehtinen S, Ashcroft P, Bonhoeffer S. On the relationship
diff --git a/man/simulate_tree_from_pop.Rd b/man/simulate_tree_from_pop.Rd
index 030f9fe8..4f42d429 100644
--- a/man/simulate_tree_from_pop.Rd
+++ b/man/simulate_tree_from_pop.Rd
@@ -8,7 +8,6 @@ population}
 simulate_tree_from_pop(
   pop,
   offspring_dist = c("pois", "nbinom"),
-  intvn_mean_reduction = 0,
   serials_dist,
   initial_immune = 0,
   t0 = 0,
@@ -24,12 +23,6 @@ corresponding to the R distribution function (e.g., "pois" for Poisson,
 where \code{\link{rpois}} is the R function to generate Poisson random
 numbers). Only supports "pois" and "nbinom".}
 
-\item{intvn_mean_reduction}{A number between 0
-and 1 for scaling/reducing the mean of \code{offspring_dist}. Serves as
-population-level intervention. \code{intvn_mean_reduction} = 0
-implies no intervention impact and \code{intvn_mean_reduction} = 1 implies full
-impact.}
-
 \item{serials_dist}{The serial interval distribution function; the name
 of a user-defined named or anonymous function with only one argument \code{n},
 representing the number of serial intervals to generate. See details.}
@@ -137,17 +130,6 @@ mu = 0.5,
 size = 1.1,
 serials_dist = function(x) 3
 )
-
-# Simulate with negative binomial offspring with intervention (50\%
-# reduction in R0)
-simulate_tree_from_pop(
-  pop = 100,
-  offspring_dist = "nbinom",
-  intvn_mean_reduction = 0.5,
-  mu = 0.5,
-  size = 1.1,
-  serials_dist = function(x) 3
-)
 }
 \seealso{
 \itemize{
diff --git a/tests/testthat/test-simulate.R b/tests/testthat/test-simulate.R
index b6fd787e..4852fd9c 100644
--- a/tests/testthat/test-simulate.R
+++ b/tests/testthat/test-simulate.R
@@ -20,25 +20,6 @@ test_that("Simulators work", {
     size = 1.1,
     serials_dist = serial_func
   )
-  #' Simulate an outbreak from a susceptible population (pois) with
-  #' 50% R0 reduction
-  susc_outbreak_raw_intvn <- simulate_tree_from_pop(
-    pop = 100,
-    offspring_dist = "pois",
-    lambda = 1.5,
-    serials_dist = serial_func,
-    intvn_mean_reduction = 0.5
-  )
-  #' Simulate an outbreak from a susceptible population (nbinom) with
-  #' 50% R0 reduction
-  susc_outbreak_raw_intvn2 <- simulate_tree_from_pop(
-    pop = 100,
-    offspring_dist = "nbinom",
-    mu = 1.5,
-    size = 1.1,
-    serials_dist = serial_func,
-    intvn_mean_reduction = 0.5
-  )
   #' Simulate a tree of infections without serials
   tree_sim_raw <- simulate_tree(
     nchains = 2,
@@ -55,25 +36,6 @@ test_that("Simulators work", {
     serials_dist = function(x) 3,
     lambda = 2
   )
-  #' Simulate a tree of infections without serials and with 50% reduction
-  #' in R0
-  tree_sim_raw_intvn <- simulate_tree(
-    nchains = 2,
-    offspring_dist = "pois",
-    statistic = "length",
-    lambda = 0.9,
-    intvn_mean_reduction = 0.5
-  )
-  #' Simulate a tree of infections with nbinom offspring and with 50% reduction
-  #' in R0
-  tree_sim_raw_intvn2 <- simulate_tree(
-    nchains = 2,
-    offspring_dist = "nbinom",
-    statistic = "length",
-    mu = 0.9,
-    size = 1.1,
-    intvn_mean_reduction = 0.5
-  )
   #' Simulate chain statistics
   chain_summary_raw <- simulate_summary(
     nchains = 2,
@@ -81,37 +43,11 @@ test_that("Simulators work", {
     statistic = "length",
     lambda = 0.9
   )
-  #' Simulate chain statistics and with a 50% reduction in R0
-  chain_summary_raw_intvn <- simulate_summary(
-    nchains = 2,
-    offspring_dist = "pois",
-    statistic = "length",
-    lambda = 0.9,
-    intvn_mean_reduction = 0.5
-  )
-  #' Simulate chain statistics with nbinom offspring and with a 50% reduction
-  #' in R0
-  chain_summary_raw_intvn2 <- simulate_summary(
-    nchains = 2,
-    offspring_dist = "nbinom",
-    statistic = "length",
-    mu = 1.9,
-    size = 1.1,
-    intvn_mean_reduction = 0.5
-  )
   #' Expectations
   expect_length(
     chain_summary_raw,
     2
   )
-  expect_length(
-    chain_summary_raw_intvn,
-    2
-  )
-  expect_length(
-    chain_summary_raw_intvn2,
-    2
-  )
   expect_gte(
     nrow(tree_sim_raw),
     2
@@ -120,14 +56,6 @@ test_that("Simulators work", {
     nrow(tree_sim_raw2),
     2
   )
-  expect_identical(
-    nrow(tree_sim_raw_intvn),
-    3L
-  )
-  expect_identical(
-    nrow(tree_sim_raw_intvn2),
-    2L
-  )
   expect_gte(
     nrow(susc_outbreak_raw),
     1
@@ -136,14 +64,6 @@ test_that("Simulators work", {
     nrow(susc_outbreak_raw2),
     1
   )
-  expect_identical(
-    nrow(susc_outbreak_raw_intvn),
-    2L
-  )
-  expect_identical(
-    nrow(susc_outbreak_raw_intvn2),
-    1L
-  )
   expect_true(
     all(
       simulate_tree(
@@ -218,17 +138,6 @@ test_that("simulate_tree throws errors", {
     ),
     "must be specified"
   )
-  expect_error(
-    simulate_tree(
-      nchains = 2,
-      offspring_dist = "binom",
-      statistic = "length",
-      size = 1,
-      prob = 0.5,
-      intvn_mean_reduction = 0.5
-    ),
-    "must be one of"
-  )
 })
 
 test_that("simulate_summary throws errors", {
@@ -270,17 +179,6 @@ test_that("simulate_summary throws errors", {
     ),
     "character string"
   )
-  expect_error(
-    simulate_summary(
-      nchains = 2,
-      offspring_dist = "binom",
-      statistic = "length",
-      size = 1,
-      prob = 0.5,
-      intvn_mean_reduction = 0.5
-    ),
-    "must be one of"
-  )
 })
 
 test_that("simulate_tree_from_pop throws errors", {
@@ -333,18 +231,8 @@ test_that("simulate_tree is numerically correct", {
     statistic = "length",
     lambda = 0.9
   )
-  #' Simulate a tree of infections without serials and with 50% reduction
-  #' in R0
-  tree_sim_raw_intvn <- simulate_tree(
-    nchains = 2,
-    offspring_dist = "pois",
-    statistic = "length",
-    lambda = 0.9,
-    intvn_mean_reduction = 0.5
-  )
   #' summarise the results
   tree_sim_summary <- summary(tree_sim_raw)
-  tree_sim_intvn_summary <- summary(tree_sim_raw_intvn)
   #' Expectations
   expect_identical(
     tree_sim_summary$chains_run,
@@ -414,17 +302,8 @@ test_that("simulate_summary is numerically correct", {
     statistic = "length",
     lambda = 0.9
   )
-  #' Simulate chain statistics and with a 50% reduction in R0
-  chain_summary_raw_intvn <- simulate_summary(
-    nchains = 2,
-    offspring_dist = "pois",
-    statistic = "length",
-    lambda = 0.9,
-    intvn_mean_reduction = 0.5
-  )
   #' Summarise the results
   chain_summary_summaries <- summary(chain_summary_raw)
-  chain_summary_intvn_summaries <- summary(chain_summary_raw_intvn)
   #' Expectations
   expect_identical(
     chain_summary_summaries$chains_run,
@@ -442,22 +321,6 @@ test_that("simulate_summary is numerically correct", {
     as.vector(chain_summary_raw),
     c(1.00, 3.00)
   )
-  expect_identical(
-    chain_summary_intvn_summaries$chains_run,
-    2.00
-  )
-  expect_identical(
-    chain_summary_intvn_summaries$max_chain_stat,
-    2.00
-  )
-  expect_identical(
-    chain_summary_intvn_summaries$min_chain_stat,
-    1.00
-  )
-  expect_identical(
-    as.vector(chain_summary_raw_intvn),
-    c(2.00, 1.00)
-  )
 })
 
 test_that("simulate_tree_from_pop is numerically correct", {
@@ -469,19 +332,8 @@ test_that("simulate_tree_from_pop is numerically correct", {
     lambda = 0.9,
     serials_dist = serial_func
   )
-  #' Simulate an outbreak from a susceptible population (pois) with
-  #' 50% R0 reduction
-  set.seed(7)
-  susc_outbreak_raw_intvn <- simulate_tree_from_pop(
-    pop = 100,
-    offspring_dist = "pois",
-    lambda = 1.5,
-    serials_dist = serial_func,
-    intvn_mean_reduction = 0.5
-  )
   #' Summarise the results
   susc_outbreak_summary <- summary(susc_outbreak_raw)
-  susc_outbreak_summary_intvn <- summary(susc_outbreak_raw_intvn)
   #' Expectations
   expect_identical(
     susc_outbreak_summary$unique_ancestors,
@@ -512,25 +364,4 @@ test_that("simulate_tree_from_pop is numerically correct", {
     susc_outbreak_raw$time,
     0.00
   )
-  #' Expectations for intervention simulation
-  expect_identical(
-    susc_outbreak_summary_intvn$unique_ancestors,
-    12L
-  )
-  expect_identical(
-    round(
-      susc_outbreak_summary_intvn$max_time,
-      1
-    ),
-    72.1
-  )
-  expect_identical(
-    susc_outbreak_summary_intvn$max_generation,
-    10L
-  )
-  expect_null(susc_outbreak_summary_intvn$chains_run)
-  expect_identical(
-    sum(aggregate(susc_outbreak_raw_intvn, "time")$cases),
-    20L
-  )
 })
diff --git a/vignettes/epichains.Rmd b/vignettes/epichains.Rmd
index 3f8c0b90..7d361795 100644
--- a/vignettes/epichains.Rmd
+++ b/vignettes/epichains.Rmd
@@ -255,72 +255,6 @@ sim_tree_from_pop_eg <- simulate_tree_from_pop(
 head(sim_tree_from_pop_eg)
 ```
 
-#### Simulating chains with interventions
-
-All the `simulate_*()` functions can model interventions that reduce the $R_0$,
-using the `intvn_mean_reduction` argument. In general, these can be
-interpreted as population-level interventions.
-
-To illustrate this, we will use the previous examples for each function and specify
-a population-level intervention that reduces $R_0$ by $50\%$.
-
-Using `simulate_tree()`, we can specify an initial number of cases
-and a population level intervention, `intvn_mean_reduction`, that reduces $R_0$ by $50\%$.
-
-```{r}
-set.seed(123)
-# Define serial distribution
-serial_func <- function(x) {
-  return(3)
-}
-
-sim_tree_intvn_eg <- simulate_tree(
-  nchains = 10,
-  statistic = "size",
-  offspring_dist = "pois",
-  intvn_mean_reduction = 0.5,
-  stat_max = 10,
-  serials_dist = serial_func,
-  lambda = 0.9
-)
-
-head(sim_tree_intvn_eg)
-```
-
-Here is an example with `simulate_summary()`, modelling an intervention that reduces $R_0$ by $50\%$.
-```{r}
-simulate_summary_intvn_eg <- simulate_summary(
-  nchains = 10,
-  statistic = "size",
-  offspring_dist = "pois",
-  intvn_mean_reduction = 0.5,
-  stat_max = 10,
-  lambda = 0.9
-)
-
-# Print the results
-simulate_summary_intvn_eg
-```
-
-Finally, let's use `simulate_tree_from_pop()`.
-```{r}
-set.seed(7)
-# Define serial distribution
-serial_func <- function(x) {
-  return(3)
-}
-
-sim_tree_from_pop_intvn_eg <- simulate_tree_from_pop(
-  pop = 1000,
-  offspring_dist = "pois",
-  intvn_mean_reduction = 0.5,
-  lambda = 1,
-  serials_dist = serial_func
-)
-
-head(sim_tree_from_pop_intvn_eg)
-```
-
 ## Other functionalities
 
 ### Summarising
diff --git a/vignettes/interventions.Rmd b/vignettes/interventions.Rmd
new file mode 100644
index 00000000..2176283f
--- /dev/null
+++ b/vignettes/interventions.Rmd
@@ -0,0 +1,167 @@
+---
+title: "Modelling disease control interventions"
+author: "Sebastian Funk and James M. Azam"
+output:
+  bookdown::html_vignette2:
+    fig_caption: yes
+    code_folding: show
+pkgdown:
+  as_is: true
+bibliography: references.json
+link-citations: true
+vignette: >
+  %\VignetteIndexEntry{Modelling disease control interventions}
+  %\VignetteEncoding{UTF-8}
+  %\VignetteEngine{knitr::rmarkdown}
+editor_options: 
+  chunk_output_type: console
+---
+
+```{r include=FALSE}
+knitr::opts_chunk$set(
+  echo = TRUE,
+  message = FALSE,
+  warning = FALSE,
+  collapse = TRUE,
+  comment = "#>"
+)
+```
+
+_epichains_ does not provide any direct functionality for studying public health interventions.
+However, the flexible simulation functionality that it includes can be used to consider some specific changes to the parameters that can be interpreted as the result of control measures.
+Here we investigate the effect on outbreak sizes, but the same approaches could be used for investigating chain lengths (using the `statistic` argument to `simulate_summary`) or the time progression of outbreaks (using the `simulate_tree` function).
+
+```{r load_libraries}
+## main package
+library("epichains")
+## for plotting
+library("ggplot2")
+## for truncating the offspring distribution later
+library("truncdist")
+```
+
+As a base case we consider the spread of an infection with a negative binomial offspring distribution with mean 1.2 and overdispersion parameter 0.5.
+We simulate 200 chains tracking up to 99 infections:
+
+```{r simulate_chains}
+sims <- simulate_summary(
+  nchains = 200, offspring_dist = "nbinom", stat_max = 99, mu = 1.2, size = 0.5
+)
+```
+
+We then plot the resulting distribution of chain sizes
+```{r uncontrolled_chains_plot}
+sims[is.infinite(sims)] <- 100 # Replace infections > 99 with 100 for plotting.
+ggplot(data.frame(x = sims), aes(x = x)) +
+  geom_histogram(breaks = seq(0, 100, by = 5), closed = "left") +
+  scale_x_continuous(breaks = c(0, 25, 50, 75, 100),
+                     labels = c(0, 25, 50, 75, ">99")) +
+  theme_bw()
+```
+
+# Reducing the strength of transmission
+
+Following [@lloyd-smith2005] we consider two ways in which disease control interventions can reduce the reproduction number: _population-wide_ and _individual-specific_ control.
+
+## Population-wide control
+
+By population-level control we mean an intervention that reduces the mean number of offspring (i.e. the reproduction number) by a fixed proportion.
+For example, to reduce R by 25% at the population level we scale the `mu` parameter from 1.2 to 0.9:
+
+```{r simulate_chains_pop_control}
+sims <- simulate_summary(
+  nchains = 200, offspring_dist = "nbinom", stat_max = 99, mu = 0.9, size = 0.5
+)
+sims[is.infinite(sims)] <- 100 # Replace infections > 99 with 100 for plotting.
+ggplot(data.frame(x = sims), aes(x = x)) +
+  geom_histogram(breaks = seq(0, 100, by = 5), closed = "left") +
+  scale_x_continuous(breaks = c(0, 25, 50, 75, 100),
+                     labels = c(0, 25, 50, 75, ">99")) +
+  theme_bw()
+```
+
+## Individual-level control.
+
+In simulating population-level control we now apply the same reduction as before (25%) but instead of assuming that the mean is reduced we apply this such that 25% of individuals do not transmit further at all, whereas the remaining 75% generate offspring as in the uncontrolled case.
+
+To do this, we can no longer use the standard negative binomial distribution that comes with R.
+Instead, we define a random generator from a modified negative binomial distribution that includes our individual-level control as a `control` argument indicating the level of individual-level control (0: no control; 1: full control):
+
+```{r nbinom_ind_control}
+rnbinom_ind <- function(n, ..., control = 0) {
+  ## initialise number of offspring to 0
+  offspring <- rep(0L, n)
+  ## for each individual, decide whether they transmit further
+  transmits <- rbinom(n = n, prob = 1 - control, size = 1)
+  ## check if anyone transmits further
+  if (any(transmits == 1L)) {
+    ## for those that transmit, sample from negative binomial with given
+    ## parameters
+    offspring[which(transmits == 1L)] <- rnbinom(n = n, ...)
+  }
+  return(offspring)
+}
+```
+
+Having defined this, we can generate simulations as before:
+
+```{r simulate_chains_ind_control}
+sims <- simulate_summary(
+  nchains = 200, offspring_dist = "nbinom_ind", stat_max = 99, mu = 1.2,
+  size = 0.5, control = 0.25
+)
+sims[is.infinite(sims)] <- 100 # Replace infections > 99 with 100 for plotting.
+ggplot(data.frame(x = sims), aes(x = x)) +
+  geom_histogram(breaks = seq(0, 100, by = 5), closed = "left") +
+  scale_x_continuous(breaks = c(0, 25, 50, 75, 100),
+                     labels = c(0, 25, 50, 75, ">99")) +
+  theme_bw()
+```
+
+# Preventing superspreading events
+
+Another way of controlling a disease would be to prevent individuals from spreading to a large number of others, for example by preventing mass gatherings or, more generally, settings where superspreading events can occur.
+
+We can model this by truncating the offspring distribution at a certain size.
+This can be done, for example, using the [truncdist](https://cran.r-project.org/package=truncdist) R package.
+We use this to define a truncated negative binomial offspring distribution:
+
+```{r negbin_truncated}
+rnbinom_truncated <- function(n, ..., max = Inf) {
+  return(rtrunc(n = n, spec = "nbinom", b = max, ...))
+}
+```
+
+We use this to simulate chains in a situation where the maximum of secondary cases that each infected person can generate is 10.
+This can be likened to a disease control strategy where gatherings are limited to 10 people.
+
+```{r simulate_chains_truncated}
+sims <- simulate_summary(
+  nchains = 200, offspring_dist = "nbinom_truncated", stat_max = 99, mu = 1.2,
+  size = 0.5, max =  10
+)
+sims[is.infinite(sims)] <- 100 # Replace infections > 99 with 100 for plotting.
+ggplot(data.frame(x = sims), aes(x = x)) +
+  geom_histogram(breaks = seq(0, 100, by = 5), closed = "left") +
+  scale_x_continuous(breaks = c(0, 25, 50, 75, 100),
+                     labels = c(0, 25, 50, 75, ">99")) +
+  theme_bw()
+```
+
+# Truncating the generation interval
+
+Lastly, we consider a situation where the generation interval is shortened.
+We do not model this explicitly but instead consider the effect on the offspring distribution.
+
+For example, if our generation interval is from a gamma distribution with shape = 25 and rate = 5 (corresponding to a mean of 5 and standard deviation of 1), and we stop all transmission that would normally occur more than 6 days after infection, we can calculate the proportion of transmissions that are prevented as
+
+```{r truncate_gen_int}
+control <- 1 - pgamma(6, shape = 25, rate = 5)
+signif(control, 2)
+```
+
+
+In other words, this would prevent `r round(100 * control)`% of infections in this example.
+The value of `control` can be used in the examples above to study the effect on outbreak sizes.
+
+# References