diff --git a/R/correlation.R b/R/correlation.R
index 41ea6a82..68119369 100644
--- a/R/correlation.R
+++ b/R/correlation.R
@@ -139,7 +139,7 @@ CorrelationParameters <- R6::R6Class(
 #'   min_wait = 0,
 #'   min_ages = 100,
 #'   max_ages = 1000,
-#'   booster_timestep = numeric(0),
+#'   booster_spacing = numeric(0),
 #'   booster_coverage = numeric(0),
 #'   booster_profile = NULL
 #' )
diff --git a/R/events.R b/R/events.R
index 879d72bd..d7ad3416 100644
--- a/R/events.R
+++ b/R/events.R
@@ -18,7 +18,7 @@ create_events <- function(parameters) {
       function(.) individual::TargetedEvent$new(parameters$human_population)
     )
     mass_pev_boosters <- lapply(
-      seq_along(parameters$mass_pev_booster_timestep),
+      seq_along(parameters$mass_pev_booster_spacing),
       function(.) individual::TargetedEvent$new(parameters$human_population)
     )
     events$mass_pev <- individual::Event$new()
@@ -33,7 +33,7 @@ create_events <- function(parameters) {
       function(.) individual::TargetedEvent$new(parameters$human_population)
     )
     pev_epi_boosters <- lapply(
-      seq_along(parameters$pev_epi_booster_timestep),
+      seq_along(parameters$pev_epi_booster_spacing),
       function(.) individual::TargetedEvent$new(parameters$human_population)
     )
     events$pev_epi_doses <- pev_epi_doses
@@ -129,9 +129,10 @@ attach_event_listeners <- function(
     attach_pev_dose_listeners(
       variables,
       parameters,
+      parameters$mass_pev_timesteps,
       events$mass_pev_doses,
       events$mass_pev_boosters,
-      parameters$mass_pev_booster_timestep,
+      parameters$mass_pev_booster_spacing,
       parameters$mass_pev_booster_coverage,
       parameters$mass_pev_profile_indices,
       'mass',
@@ -143,9 +144,10 @@ attach_event_listeners <- function(
     attach_pev_dose_listeners(
       variables,
       parameters,
+      parameters$pev_epi_timesteps,
       events$pev_epi_doses,
       events$pev_epi_boosters,
-      parameters$pev_epi_booster_timestep,
+      parameters$pev_epi_booster_spacing,
       parameters$pev_epi_booster_coverage,
       parameters$pev_epi_profile_indices,
       'epi',
diff --git a/R/pev.R b/R/pev.R
index 039b6044..bc27fc9f 100644
--- a/R/pev.R
+++ b/R/pev.R
@@ -179,6 +179,7 @@ create_pev_efficacy_listener <- function(variables, pev_profile_index) {
 create_pev_booster_listener <- function(
   variables,
   coverage,
+  pev_distribution_timesteps,
   booster_number,
   pev_profile_index,
   next_booster_event,
@@ -192,7 +193,11 @@ create_pev_booster_listener <- function(
   force(next_booster_delay)
   force(coverage)
   function(timestep, target) {
-    target <- sample_bitset(target, coverage)
+    cov_t <- coverage[
+      match_timestep(pev_distribution_timesteps, timestep),
+      booster_number
+    ]
+    target <- sample_bitset(target, cov_t)
     variables$last_pev_timestep$queue_update(timestep, target)
     variables$last_eff_pev_timestep$queue_update(timestep, target)
     variables$pev_profile$queue_update(pev_profile_index, target)
@@ -236,6 +241,7 @@ create_dosage_renderer <- function(renderer, strategy, dose) {
 attach_pev_dose_listeners <- function(
   variables,
   parameters,
+  pev_distribution_timesteps,
   dose_events,
   booster_events,
   booster_delays,
@@ -302,7 +308,8 @@ attach_pev_dose_listeners <- function(
     booster_events[[b]]$add_listener(
       create_pev_booster_listener(
         variables = variables,
-        coverage = booster_coverages[[b]],
+        coverage = booster_coverages,
+        pev_distribution_timesteps = pev_distribution_timesteps,
         booster_number = b,
         pev_profile_index = pev_profile_indices[[b + 1]],
         next_booster_event = next_booster_event,
diff --git a/R/pev_parameters.R b/R/pev_parameters.R
index aeb6bd5b..7edcb39c 100644
--- a/R/pev_parameters.R
+++ b/R/pev_parameters.R
@@ -62,7 +62,7 @@ rtss_booster_profile <- create_pev_profile(
 #' age. Efficacy will take effect after the last dose
 #'
 #' @param parameters a list of parameters to modify
-#' @param profile primary vaccine profile of type PEVProfile
+#' @param profile a list of details for the vaccine profile, create with `create_pev_profile`
 #' @param coverages a vector of coverages for the primary doses
 #' @param timesteps a vector of timesteps associated with coverages
 #' @param age the age when an individual will receive the first dose of the
@@ -72,12 +72,11 @@ rtss_booster_profile <- create_pev_profile(
 #' between an individual receiving the final dose and the first booster. When using
 #' both set_mass_pev and set_pev_epi, this represents the minimum time between
 #' an individual being vaccinated under one scheme and vaccinated under another.
-#' @param booster_timestep the timesteps (following the final dose) at which booster vaccinations are administered
-#' @param booster_coverage the proportion of the vaccinated population relative to the last vaccination (whether a previous booster or the primary series)
-#' @param booster_profile list of booster vaccine profiles, of type
-#' PEVProfile, for each timestep in booster_timeteps
+#' @param booster_spacing the timesteps (following the final primary dose) at which booster vaccinations are administered
+#' @param booster_coverage a matrix of coverages (timesteps x boosters) specifying the proportion the previously vaccinated population to continue receiving booster doses. The rows of the matrix must be the same size as `timesteps`. The columns of the matrix must be the same size as `booster_spacing`.
+#' @param booster_profile list of lists representing each booster profile, the outer list must be the same length as `booster_spacing`. Create vaccine profiles with `create_pev_profile`
 #' @param seasonal_boosters logical, if TRUE the first booster timestep is
-#' relative to the start of the year, otherwise they are relative to the last dose
+#' relative to the start of the year, otherwise they are relative to the last primary dose
 #' @export
 set_pev_epi <- function(
   parameters,
@@ -86,33 +85,47 @@ set_pev_epi <- function(
   timesteps,
   age,
   min_wait,
-  booster_timestep,
+  booster_spacing,
   booster_coverage,
   booster_profile,
   seasonal_boosters = FALSE
   ) {
   stopifnot(all(coverages >= 0) && all(coverages <= 1))
+  stopifnot(is.matrix(booster_coverage))
 
   # Check that the primary timing parameters make sense
   if(length(coverages) != length(timesteps)){
     stop("coverages and timesteps must align")
   }
 
+  # Check that booster_spacing are monotonically increasing
+  if (length(booster_spacing) > 1) {
+    if (!all(diff(booster_spacing) > 0)) {
+      stop('booster_spacing must be monotonically increasing')
+    }
+  }
+
   # Check that seasonal booster parameters make sense
   stopifnot(min_wait >= 0)
   stopifnot(age >= 0)
   stopifnot(is.logical(seasonal_boosters))
   if (seasonal_boosters) {
-    if(booster_timestep[[1]] < 0) {
-      booster_timestep <- booster_timestep + 365
+    if(booster_spacing[[1]] < 0) {
+      booster_spacing <- booster_spacing + 365
     }
   }
 
   # Check that the booster timing parameters make sense
-  stopifnot((length(booster_timestep) == 0) || all(booster_timestep > 0))
+  stopifnot((length(booster_spacing) == 0) || all(booster_spacing > 0))
   stopifnot((length(booster_coverage)) == 0 || all(booster_coverage >= 0 & booster_coverage <= 1))
-  if (!all(c(length(booster_coverage), length(booster_timestep), length(booster_profile)) == length(booster_timestep))) {
-    stop('booster_timestep and booster_coverage and booster_profile does not align')
+  if (!all(c(ncol(booster_coverage), length(booster_profile)) == length(booster_spacing))) {
+    stop('booster_spacing, booster_coverage and booster_profile do not align')
+  }
+  # Check that booster_coverage and timesteps align
+  if (length(booster_coverage) > 0) {
+    if (nrow(booster_coverage) != length(timesteps)) {
+      stop('booster_coverage and timesteps do not align')
+    }
   }
 
   # Index the new vaccine profiles
@@ -125,7 +138,7 @@ set_pev_epi <- function(
   parameters$pev_epi_coverages <- coverages
   parameters$pev_epi_timesteps <- timesteps
   parameters$pev_epi_age <- age
-  parameters$pev_epi_booster_timestep <- booster_timestep
+  parameters$pev_epi_booster_spacing <- booster_spacing
   parameters$pev_epi_min_wait <- min_wait
   parameters$pev_epi_booster_coverage <- booster_coverage
   parameters$pev_epi_profile_indices <- profile_indices
@@ -139,7 +152,7 @@ set_pev_epi <- function(
 #' Efficacy will take effect after the last dose
 #'
 #' @param parameters a list of parameters to modify
-#' @param profile primary vaccine profile of type PEVProfile
+#' @param profile a list of details for the vaccine profile, create with `create_pev_profile`
 #' @param timesteps a vector of timesteps for each round of vaccinations
 #' @param coverages the coverage for each round of vaccinations
 #' @param min_wait the minimum acceptable time since the last vaccination (in timesteps);
@@ -147,10 +160,9 @@ set_pev_epi <- function(
 #' time between an individual being vaccinated under one scheme and vaccinated under another.
 #' @param min_ages for the target population, inclusive (in timesteps)
 #' @param max_ages for the target population, inclusive (in timesteps)
-#' @param booster_timestep the timesteps (following the initial vaccination) at which booster vaccinations are administered
-#' @param booster_coverage the proportion of the vaccinated population relative to the last vaccination (whether a previous booster or the primary series)
-#' @param booster_profile list of booster vaccine profiles, of type
-#' PEVProfile, for each timestep in booster_timeteps
+#' @param booster_spacing the timesteps (following the final primary dose) at which booster vaccinations are administered
+#' @param booster_coverage a matrix of coverages (timesteps x boosters) specifying the proportion the previously vaccinated population to continue receiving booster doses. The rows of the matrix must be the same size as `timesteps`. The columns of the matrix must be the same size as `booster_spacing`.
+#' @param booster_profile list of lists representing each booster profile, the outer list must be the same length as `booster_spacing`. Create vaccine profiles with `create_pev_profile`
 #' @export
 set_mass_pev <- function(
   parameters,
@@ -160,7 +172,7 @@ set_mass_pev <- function(
   min_ages,
   max_ages,
   min_wait,
-  booster_timestep,
+  booster_spacing,
   booster_coverage,
   booster_profile
   ) {
@@ -168,13 +180,28 @@ set_mass_pev <- function(
   stopifnot(min_wait >= 0)
   stopifnot(all(coverages >= 0) && all(coverages <= 1))
   stopifnot(all(min_ages >= 0 & max_ages >= 0))
-  stopifnot(all(booster_timestep > 0))
+  stopifnot(all(booster_spacing > 0))
   stopifnot(all(booster_coverage >= 0 & booster_coverage <= 1))
   if (length(min_ages) != length(max_ages)) {
     stop('min and max ages do not align')
   }
-  if (!all(c(length(booster_coverage), length(booster_timestep), length(booster_profile)) == length(booster_timestep))) {
-    stop('booster_timestep, booster_coverage and booster_profile does not align')
+
+  # Check that booster_spacing are monotonically increasing
+  if (length(booster_spacing) > 1) {
+    if (!all(diff(booster_spacing) > 0)) {
+      stop('booster_spacing must be monotonically increasing')
+    }
+  }
+
+  stopifnot((length(booster_coverage)) == 0 || all(booster_coverage >= 0 & booster_coverage <= 1))
+  if (!all(c(ncol(booster_coverage), length(booster_profile)) == length(booster_spacing))) {
+    stop('booster_spacing, booster_coverage and booster_profile do not align')
+  }
+  # Check that booster_coverage and timesteps align
+  if (length(booster_coverage) > 0) {
+    if (nrow(booster_coverage) != length(timesteps)) {
+      stop('booster_coverage and timesteps do not align')
+    }
   }
 
   # Index the new vaccine profiles
@@ -189,7 +216,7 @@ set_mass_pev <- function(
   parameters$mass_pev_min_ages <- min_ages
   parameters$mass_pev_max_ages <- max_ages
   parameters$mass_pev_min_wait <- min_wait
-  parameters$mass_pev_booster_timestep <- booster_timestep
+  parameters$mass_pev_booster_spacing <- booster_spacing
   parameters$mass_pev_booster_coverage <- booster_coverage
   parameters$mass_pev_profile_indices <- profile_indices
   parameters
diff --git a/man/CorrelationParameters.Rd b/man/CorrelationParameters.Rd
index c2e6ada7..b1d22578 100644
--- a/man/CorrelationParameters.Rd
+++ b/man/CorrelationParameters.Rd
@@ -14,17 +14,17 @@ Describes an event in the simulation
 \section{Methods}{
 \subsection{Public methods}{
 \itemize{
-\item \href{#method-CorrelationParameters-new}{\code{CorrelationParameters$new()}}
-\item \href{#method-CorrelationParameters-inter_round_rho}{\code{CorrelationParameters$inter_round_rho()}}
-\item \href{#method-CorrelationParameters-inter_intervention_rho}{\code{CorrelationParameters$inter_intervention_rho()}}
-\item \href{#method-CorrelationParameters-sigma}{\code{CorrelationParameters$sigma()}}
-\item \href{#method-CorrelationParameters-mvnorm}{\code{CorrelationParameters$mvnorm()}}
-\item \href{#method-CorrelationParameters-clone}{\code{CorrelationParameters$clone()}}
+\item \href{#method-new}{\code{CorrelationParameters$new()}}
+\item \href{#method-inter_round_rho}{\code{CorrelationParameters$inter_round_rho()}}
+\item \href{#method-inter_intervention_rho}{\code{CorrelationParameters$inter_intervention_rho()}}
+\item \href{#method-sigma}{\code{CorrelationParameters$sigma()}}
+\item \href{#method-mvnorm}{\code{CorrelationParameters$mvnorm()}}
+\item \href{#method-clone}{\code{CorrelationParameters$clone()}}
 }
 }
 \if{html}{\out{<hr>}}
-\if{html}{\out{<a id="method-CorrelationParameters-new"></a>}}
-\if{latex}{\out{\hypertarget{method-CorrelationParameters-new}{}}}
+\if{html}{\out{<a id="method-new"></a>}}
+\if{latex}{\out{\hypertarget{method-new}{}}}
 \subsection{Method \code{new()}}{
 initialise correlation parameters
 \subsection{Usage}{
@@ -40,8 +40,8 @@ initialise correlation parameters
 }
 }
 \if{html}{\out{<hr>}}
-\if{html}{\out{<a id="method-CorrelationParameters-inter_round_rho"></a>}}
-\if{latex}{\out{\hypertarget{method-CorrelationParameters-inter_round_rho}{}}}
+\if{html}{\out{<a id="method-inter_round_rho"></a>}}
+\if{latex}{\out{\hypertarget{method-inter_round_rho}{}}}
 \subsection{Method \code{inter_round_rho()}}{
 Add rho between rounds
 \subsection{Usage}{
@@ -60,8 +60,8 @@ the intervention}
 }
 }
 \if{html}{\out{<hr>}}
-\if{html}{\out{<a id="method-CorrelationParameters-inter_intervention_rho"></a>}}
-\if{latex}{\out{\hypertarget{method-CorrelationParameters-inter_intervention_rho}{}}}
+\if{html}{\out{<a id="method-inter_intervention_rho"></a>}}
+\if{latex}{\out{\hypertarget{method-inter_intervention_rho}{}}}
 \subsection{Method \code{inter_intervention_rho()}}{
 Add rho between interventions
 \subsection{Usage}{
@@ -83,8 +83,8 @@ the intervention}
 }
 }
 \if{html}{\out{<hr>}}
-\if{html}{\out{<a id="method-CorrelationParameters-sigma"></a>}}
-\if{latex}{\out{\hypertarget{method-CorrelationParameters-sigma}{}}}
+\if{html}{\out{<a id="method-sigma"></a>}}
+\if{latex}{\out{\hypertarget{method-sigma}{}}}
 \subsection{Method \code{sigma()}}{
 Standard deviation of each intervention between rounds
 \subsection{Usage}{
@@ -93,8 +93,8 @@ Standard deviation of each intervention between rounds
 
 }
 \if{html}{\out{<hr>}}
-\if{html}{\out{<a id="method-CorrelationParameters-mvnorm"></a>}}
-\if{latex}{\out{\hypertarget{method-CorrelationParameters-mvnorm}{}}}
+\if{html}{\out{<a id="method-mvnorm"></a>}}
+\if{latex}{\out{\hypertarget{method-mvnorm}{}}}
 \subsection{Method \code{mvnorm()}}{
 multivariate norm draws for these parameters
 \subsection{Usage}{
@@ -103,8 +103,8 @@ multivariate norm draws for these parameters
 
 }
 \if{html}{\out{<hr>}}
-\if{html}{\out{<a id="method-CorrelationParameters-clone"></a>}}
-\if{latex}{\out{\hypertarget{method-CorrelationParameters-clone}{}}}
+\if{html}{\out{<a id="method-clone"></a>}}
+\if{latex}{\out{\hypertarget{method-clone}{}}}
 \subsection{Method \code{clone()}}{
 The objects of this class are cloneable with this method.
 \subsection{Usage}{
diff --git a/man/get_correlation_parameters.Rd b/man/get_correlation_parameters.Rd
index 21995b91..068aed64 100644
--- a/man/get_correlation_parameters.Rd
+++ b/man/get_correlation_parameters.Rd
@@ -27,7 +27,7 @@ parameters <- set_mass_pev(
   min_wait = 0,
   min_ages = 100,
   max_ages = 1000,
-  booster_timestep = numeric(0),
+  booster_spacing = numeric(0),
   booster_coverage = numeric(0),
   booster_profile = NULL
 )
diff --git a/man/set_mass_pev.Rd b/man/set_mass_pev.Rd
index 7c28547f..ef3263e5 100644
--- a/man/set_mass_pev.Rd
+++ b/man/set_mass_pev.Rd
@@ -12,7 +12,7 @@ set_mass_pev(
   min_ages,
   max_ages,
   min_wait,
-  booster_timestep,
+  booster_spacing,
   booster_coverage,
   booster_profile
 )
@@ -20,7 +20,7 @@ set_mass_pev(
 \arguments{
 \item{parameters}{a list of parameters to modify}
 
-\item{profile}{primary vaccine profile of type PEVProfile}
+\item{profile}{a list of details for the vaccine profile, create with \code{create_pev_profile}}
 
 \item{timesteps}{a vector of timesteps for each round of vaccinations}
 
@@ -34,12 +34,11 @@ set_mass_pev(
 When using both set_mass_pev and set_pev_epi, this represents the minimum
 time between an individual being vaccinated under one scheme and vaccinated under another.}
 
-\item{booster_timestep}{the timesteps (following the initial vaccination) at which booster vaccinations are administered}
+\item{booster_spacing}{the timesteps (following the final primary dose) at which booster vaccinations are administered}
 
-\item{booster_coverage}{the proportion of the vaccinated population relative to the last vaccination (whether a previous booster or the primary series)}
+\item{booster_coverage}{a matrix of coverages (timesteps x boosters) specifying the proportion the previously vaccinated population to continue receiving booster doses. The rows of the matrix must be the same size as \code{timesteps}. The columns of the matrix must be the same size as \code{booster_spacing}.}
 
-\item{booster_profile}{list of booster vaccine profiles, of type
-PEVProfile, for each timestep in booster_timeteps}
+\item{booster_profile}{list of lists representing each booster profile, the outer list must be the same length as \code{booster_spacing}. Create vaccine profiles with \code{create_pev_profile}}
 }
 \description{
 distribute pre-erythrocytic vaccine to a population in an age range.
diff --git a/man/set_pev_epi.Rd b/man/set_pev_epi.Rd
index 4dea14ed..a3a4d195 100644
--- a/man/set_pev_epi.Rd
+++ b/man/set_pev_epi.Rd
@@ -11,7 +11,7 @@ set_pev_epi(
   timesteps,
   age,
   min_wait,
-  booster_timestep,
+  booster_spacing,
   booster_coverage,
   booster_profile,
   seasonal_boosters = FALSE
@@ -20,7 +20,7 @@ set_pev_epi(
 \arguments{
 \item{parameters}{a list of parameters to modify}
 
-\item{profile}{primary vaccine profile of type PEVProfile}
+\item{profile}{a list of details for the vaccine profile, create with \code{create_pev_profile}}
 
 \item{coverages}{a vector of coverages for the primary doses}
 
@@ -35,15 +35,14 @@ between an individual receiving the final dose and the first booster. When using
 both set_mass_pev and set_pev_epi, this represents the minimum time between
 an individual being vaccinated under one scheme and vaccinated under another.}
 
-\item{booster_timestep}{the timesteps (following the final dose) at which booster vaccinations are administered}
+\item{booster_spacing}{the timesteps (following the final primary dose) at which booster vaccinations are administered}
 
-\item{booster_coverage}{the proportion of the vaccinated population relative to the last vaccination (whether a previous booster or the primary series)}
+\item{booster_coverage}{a matrix of coverages (timesteps x boosters) specifying the proportion the previously vaccinated population to continue receiving booster doses. The rows of the matrix must be the same size as \code{timesteps}. The columns of the matrix must be the same size as \code{booster_spacing}.}
 
-\item{booster_profile}{list of booster vaccine profiles, of type
-PEVProfile, for each timestep in booster_timeteps}
+\item{booster_profile}{list of lists representing each booster profile, the outer list must be the same length as \code{booster_spacing}. Create vaccine profiles with \code{create_pev_profile}}
 
 \item{seasonal_boosters}{logical, if TRUE the first booster timestep is
-relative to the start of the year, otherwise they are relative to the last dose}
+relative to the start of the year, otherwise they are relative to the last primary dose}
 }
 \description{
 distribute vaccine when an individual becomes a certain
diff --git a/tests/testthat/test-biology.R b/tests/testthat/test-biology.R
index 5895d45b..c739f7e7 100644
--- a/tests/testthat/test-biology.R
+++ b/tests/testthat/test-biology.R
@@ -92,13 +92,13 @@ test_that('FOIM is consistent with equilibrium', {
       psi <- unique_biting_rate(age, parameters)
       zeta <- variables$zeta$get_values()
       .pi <- human_pi(psi, zeta)
-      calculate_foim(a, sum(.pi * variables$infectivity$get_values()))
+      calculate_foim(a, sum(.pi * variables$infectivity$get_values()), 1.)
     }
   )
   expect_equal(
     expected_foim,
     actual_foim,
-    tolerance = 1e-4
+    tolerance = 1e-3
   )
 })
 
diff --git a/tests/testthat/test-infection-integration.R b/tests/testthat/test-infection-integration.R
index 238a2ae6..0400a361 100644
--- a/tests/testthat/test-infection-integration.R
+++ b/tests/testthat/test-infection-integration.R
@@ -121,8 +121,8 @@ test_that('calculate_infections works various combinations of drug and vaccinati
     min_ages = 0,
     max_ages = 100 * 365,
     min_wait = 0,
-    booster_timestep = 365,
-    booster_coverage = 1,
+    booster_spacing = 365,
+    booster_coverage = matrix(1),
     booster_profile = list(rtss_booster_profile)
   )
 
diff --git a/tests/testthat/test-pev-epi.R b/tests/testthat/test-pev-epi.R
index ef9148a3..6555700c 100644
--- a/tests/testthat/test-pev-epi.R
+++ b/tests/testthat/test-pev-epi.R
@@ -7,8 +7,8 @@ test_that('pev epi strategy parameterisation works', {
     timesteps = c(10, 100),
     min_wait = 0,
     age = 5 * 30,
-    booster_timestep = c(18, 36) * 30,
-    booster_coverage = c(.9, .8),
+    booster_spacing = c(18, 36) * 30,
+    booster_coverage = matrix(c(.9, .8, .9, .8), nrow=2, ncol=2),
     booster_profile = list(rtss_booster_profile, rtss_booster_profile)
   )
   expect_equal(parameters$pev, TRUE)
@@ -16,7 +16,7 @@ test_that('pev epi strategy parameterisation works', {
   expect_equal(parameters$pev_epi_timesteps, c(10, 100))
   expect_equal(parameters$pev_epi_age, 5 * 30)
   expect_equal(parameters$pev_epi_min_wait, 0)
-  expect_equal(parameters$pev_epi_booster_timestep, c(18, 36) * 30)
+  expect_equal(parameters$pev_epi_booster_spacing, c(18, 36) * 30)
   expect_equal(parameters$pev_profiles, list(rtss_profile, rtss_booster_profile, rtss_booster_profile))
   expect_equal(parameters$pev_epi_profile_indices, seq(3))
 
@@ -28,8 +28,8 @@ test_that('pev epi strategy parameterisation works', {
       timesteps = 10,
       min_wait = 0,
       age = 5 * 30,
-      booster_timestep = c(18, 36) * 30,
-      booster_coverage = c(.9, .8),
+      booster_spacing = c(18, 36) * 30,
+      booster_coverage = matrix(c(.9, .8, .9, .8), nrow=2, ncol=2),
       booster_profile = list(rtss_booster_profile, rtss_booster_profile)
     ), "all(coverages >= 0) && all(coverages <= 1) is not TRUE",
     fixed = TRUE
@@ -43,30 +43,53 @@ test_that('pev epi strategy parameterisation works', {
       timesteps = 10,
       min_wait = 0,
       age = 5 * 30,
-      booster_timestep = c(18, 36) * 30,
-      booster_coverage = c(.9, .8),
+      booster_spacing = c(18, 36) * 30,
+      booster_coverage = matrix(c(.9, .8, .9, .8), nrow=2, ncol=2),
       booster_profile = list(rtss_booster_profile, rtss_booster_profile)
     ), "all(coverages >= 0) && all(coverages <= 1) is not TRUE",
     fixed = TRUE
   )
 })
 
-test_that('pev epi fails pre-emptively with unaligned booster parameters', {
+test_that('set_pev_epi checks booster coverage matrix shape', {
   parameters <- get_parameters()
   expect_error(
-    set_pev_epi(
+    parameters <- set_pev_epi(
+      parameters,
+      profile = rtss_profile,
+      coverages = c(0.1, 0.8),
+      timesteps = c(10, 100),
+      min_wait = 0,
+      age = 5 * 30,
+      booster_spacing = c(18, 36) * 30,
+      booster_coverage = matrix(c(.9, .8), nrow=2, ncol=1),
+      booster_profile = list(rtss_booster_profile, rtss_booster_profile)
+    ),
+    'booster_spacing, booster_coverage and booster_profile do not align',
+    fixed = TRUE
+  )
+})
+
+test_that('set_pev_epi checks that booster_spacing are increasing', {
+  parameters <- get_parameters()
+  expect_error(
+    parameters <- set_pev_epi(
+      parameters,
       profile = rtss_profile,
       coverages = c(0.1, 0.8),
       timesteps = c(10, 100),
       min_wait = 0,
       age = 5 * 30,
-      booster_timestep = c(18, 36) * 30,
-      booster_coverage = .9,
+      booster_spacing = c(5, 5) * 30,
+      booster_coverage = matrix(c(.9, .8), nrow=2, ncol=1),
       booster_profile = list(rtss_booster_profile, rtss_booster_profile)
-    )
+    ),
+    'booster_spacing must be monotonically increasing',
+    fixed = TRUE
   )
 })
 
+
 test_that('pev epi targets correct age and respects min_wait', {
   timestep <- 5*365 
   parameters <- get_parameters(list(human_population = 5))
@@ -77,8 +100,8 @@ test_that('pev epi targets correct age and respects min_wait', {
     coverages = 0.8,
     min_wait = 2*365,
     age = 18 * 365,
-    booster_timestep = c(18, 36) * 30,
-    booster_coverage = c(.9, .8),
+    booster_spacing = c(18, 36) * 30,
+    booster_coverage = matrix(c(.9, .8), nrow=1, ncol=2),
     booster_profile = list(rtss_booster_profile, rtss_booster_profile)
   )
   events <- create_events(parameters)
@@ -141,8 +164,8 @@ test_that('EPI ignores individuals scheduled for mass vaccination', {
     min_wait = 0,
     min_ages = c(1, 2, 3, 18) * 365,
     max_ages = (c(1, 2, 3, 18) + 1) * 365 - 1,
-    booster_timestep = c(18, 36) * 30,
-    booster_coverage = c(.9, .8),
+    booster_spacing = c(18, 36) * 30,
+    booster_coverage = matrix(c(.9, .8, .9, .8), nrow=2, ncol=2),
     booster_profile = list(rtss_booster_profile, rtss_booster_profile)
   )
   parameters <- set_pev_epi(
@@ -152,8 +175,8 @@ test_that('EPI ignores individuals scheduled for mass vaccination', {
     coverages = 0.8,
     min_wait = 0,
     age = 18 * 365,
-    booster_timestep = c(18, 36) * 30,
-    booster_coverage = c(.9, .8),
+    booster_spacing = c(18, 36) * 30,
+    booster_coverage = matrix(c(.9, .8), nrow=1, ncol=2),
     booster_profile = list(rtss_booster_profile, rtss_booster_profile)
   )
   events <- create_events(parameters)
@@ -210,8 +233,8 @@ test_that('pev EPI respects min_wait when scheduling seasonal boosters', {
     coverages = 0.8,
     min_wait = 6 * 30,
     age = 18 * 365,
-    booster_timestep = c(3, 12) * 30,
-    booster_coverage = c(.9, .8),
+    booster_spacing = c(3, 12) * 30,
+    booster_coverage = matrix(c(.9, .8), nrow=1, ncol=2),
     booster_profile = list(rtss_booster_profile, rtss_booster_profile),
     seasonal_boosters = TRUE
   )
@@ -246,8 +269,8 @@ test_that('pev EPI schedules for the following year with seasonal boosters', {
     coverages = 0.8,
     min_wait = 6 * 30,
     age = 18 * 365,
-    booster_timestep = c(3, 12) * 30,
-    booster_coverage = c(.9, .8),
+    booster_spacing = c(3, 12) * 30,
+    booster_coverage = matrix(c(.9, .8), nrow=1, ncol=2),
     booster_profile = list(rtss_booster_profile, rtss_booster_profile),
     seasonal_boosters = TRUE
   )
diff --git a/tests/testthat/test-pev.R b/tests/testthat/test-pev.R
index 30007cd3..88b95b35 100644
--- a/tests/testthat/test-pev.R
+++ b/tests/testthat/test-pev.R
@@ -8,8 +8,8 @@ test_that('Mass vaccination strategy parameterisation works', {
     min_wait = 0,
     min_ages = 5 * 30,
     max_ages = 17 * 30,
-    booster_timestep = c(18, 36) * 30,
-    booster_coverage = c(.9, .8),
+    booster_spacing = c(18, 36) * 30,
+    booster_coverage = matrix(c(.9, .8), nrow=1, ncol=2),
     booster_profile = list(rtss_booster_profile, rtss_booster_profile)
   )
   expect_equal(parameters$pev, TRUE)
@@ -17,7 +17,7 @@ test_that('Mass vaccination strategy parameterisation works', {
   expect_equal(parameters$mass_pev_coverages, .8)
   expect_equal(parameters$mass_pev_min_ages, 5 * 30)
   expect_equal(parameters$mass_pev_max_ages, 17 * 30)
-  expect_equal(parameters$mass_pev_booster_timestep, c(18, 36) * 30)
+  expect_equal(parameters$mass_pev_booster_spacing, c(18, 36) * 30)
   expect_equal(parameters$pev_profiles, list(rtss_profile, rtss_booster_profile, rtss_booster_profile))
   expect_equal(parameters$mass_pev_profile_indices, seq(3))
 
@@ -30,8 +30,8 @@ test_that('Mass vaccination strategy parameterisation works', {
       min_wait = 0,
       min_ages = 5 * 30,
       max_ages = 17 * 30,
-      booster_timestep = c(18, 36) * 30,
-      booster_coverage = c(.9, .8),
+      booster_spacing = c(18, 36) * 30,
+      booster_coverage = matrix(c(.9, .8), nrow=1, ncol=2),
       booster_profile = list(rtss_booster_profile, rtss_booster_profile)
     ), "all(coverages >= 0) && all(coverages <= 1) is not TRUE",
     fixed = TRUE
@@ -46,14 +46,54 @@ test_that('Mass vaccination strategy parameterisation works', {
       min_wait = 0,
       min_ages = 5 * 30,
       max_ages = 17 * 30,
-      booster_timestep = c(18, 36) * 30,
-      booster_coverage = c(.9, .8),
+      booster_spacing = c(18, 36) * 30,
+      booster_coverage = matrix(c(.9, .8), nrow=1, ncol=2),
       booster_profile = list(rtss_booster_profile, rtss_booster_profile)
     ), "all(coverages >= 0) && all(coverages <= 1) is not TRUE",
     fixed = TRUE
   )
 })
 
+test_that('set_mass_pev checks booster coverage matrix shape', {
+  parameters <- get_parameters()
+  expect_error(
+    parameters <- set_mass_pev(
+      parameters,
+      profile = rtss_profile,
+      coverages = c(0.1),
+      timesteps = c(10),
+      min_wait = 0,
+      min_ages = 5 * 30,
+      max_ages = 17 * 30,
+      booster_spacing = c(18, 36) * 30,
+      booster_coverage = matrix(c(.9, .8), nrow=2, ncol=1),
+      booster_profile = list(rtss_booster_profile, rtss_booster_profile)
+    ),
+    'booster_spacing, booster_coverage and booster_profile do not align',
+    fixed = TRUE
+  )
+})
+
+test_that('set_mass_pev checks booster_spacing is increasing', {
+  parameters <- get_parameters()
+  expect_error(
+    parameters <- set_mass_pev(
+      parameters,
+      profile = rtss_profile,
+      coverages = c(0.1),
+      timesteps = c(10),
+      min_wait = 0,
+      min_ages = 5 * 30,
+      max_ages = 17 * 30,
+      booster_spacing = c(5, 5) * 30,
+      booster_coverage = matrix(c(.9, .8), nrow=2, ncol=1),
+      booster_profile = list(rtss_booster_profile, rtss_booster_profile)
+    ),
+    'booster_spacing must be monotonically increasing',
+    fixed = TRUE
+  )
+})
+
 test_that('Mass vaccination fails pre-emptively for unaligned booster parameters', {
   parameters <- get_parameters()
   expect_error(
@@ -65,8 +105,8 @@ test_that('Mass vaccination fails pre-emptively for unaligned booster parameters
       min_wait = 0,
       min_ages = 5 * 30,
       max_ages = 17 * 30,
-      booster_timestep = c(18, 36) * 30,
-      booster_coverage = c(.9),
+      booster_spacing = c(18, 36) * 30,
+      booster_coverage = matrix(.9, nrow=1, ncol=1),
       booster_profile = list(rtss_booster_profile, rtss_booster_profile)
     )
   )
@@ -83,8 +123,8 @@ test_that('Infection considers pev efficacy', {
     min_wait = 0,
     min_ages = 5 * 30,
     max_ages = 17 * 30,
-    booster_timestep = c(18, 36) * 30,
-    booster_coverage = c(.9, .8),
+    booster_spacing = c(18, 36) * 30,
+    booster_coverage = matrix(c(.9, .8), nrow=1, ncol=2),
     booster_profile = list(rtss_booster_profile, rtss_booster_profile)
   )
   events <- create_events(parameters)
@@ -150,8 +190,8 @@ test_that('Mass vaccinations update vaccination time', {
     min_wait = 0,
     min_ages = c(1, 2, 3, 18) * 365,
     max_ages = (c(1, 2, 3, 18) + 1) * 365 - 1,
-    booster_timestep = c(18, 36) * 30,
-    booster_coverage = c(.9, .8),
+    booster_spacing = c(18, 36) * 30,
+    booster_coverage = matrix(c(.9, .8, .9, .8), nrow=2, ncol=2),
     booster_profile = list(rtss_booster_profile, rtss_booster_profile)
   )
   events <- create_events(parameters)
@@ -220,8 +260,8 @@ test_that('Mass vaccinations ignore EPI individuals', {
     min_wait = 0,
     min_ages = c(1, 2, 3, 18) * 365,
     max_ages = (c(1, 2, 3, 18) + 1) * 365 - 1,
-    booster_timestep = c(18, 36) * 30,
-    booster_coverage = c(.9, .8),
+    booster_spacing = c(18, 36) * 30,
+    booster_coverage = matrix(c(.9, .8, .9, .8), nrow=2, ncol=2),
     booster_profile = list(rtss_booster_profile, rtss_booster_profile)
   )
   parameters <- set_pev_epi(
@@ -231,8 +271,8 @@ test_that('Mass vaccinations ignore EPI individuals', {
     coverages = 0.8,
     min_wait = 2*365,
     age = 18 * 365,
-    booster_timestep = c(18, 36) * 30,
-    booster_coverage = c(.9, .8),
+    booster_spacing = c(18, 36) * 30,
+    booster_coverage = matrix(c(.9, .8), nrow=1, ncol=2),
     booster_profile = list(rtss_booster_profile, rtss_booster_profile)
   )
   events <- create_events(parameters)
@@ -287,8 +327,8 @@ test_that('Mass boosters update profile params and reschedule correctly', {
     min_wait = 0,
     min_ages = c(1, 2, 3, 18) * 365,
     max_ages = (c(1, 2, 3, 18) + 1) * 365 - 1,
-    booster_timestep = c(1, 6) * 30,
-    booster_coverage = c(1, 1),
+    booster_spacing = c(1, 6) * 30,
+    booster_coverage = matrix(1, nrow=2, ncol=2),
     booster_profile = list(rtss_booster_profile, rtss_booster_profile)
   )
   events <- create_events(parameters)
@@ -312,7 +352,8 @@ test_that('Mass boosters update profile params and reschedule correctly', {
 
   listener <- create_pev_booster_listener(
     variables = variables,
-    coverage = 1,
+    coverage = parameters$mass_pev_booster_coverage,
+    parameters$mass_pev_timesteps,
     booster_number = 1,
     pev_profile_index = 2,
     next_booster_event = events$mass_pev_boosters[[2]],
@@ -363,8 +404,8 @@ test_that('Mass booster coverages sample subpopulations correctly', {
     min_ages = c(1, 2, 3, 18) * 365,
     max_ages = (c(1, 2, 3, 18) + 1) * 365 - 1,
     min_wait = 0,
-    booster_timestep = c(1, 6) * 30,
-    booster_coverage = c(.9, .8),
+    booster_spacing = c(1, 6) * 30,
+    booster_coverage = matrix(c(.9, .8), nrow=1, ncol=2),
     booster_profile = list(rtss_booster_profile, rtss_booster_profile)
   )
 
@@ -389,7 +430,8 @@ test_that('Mass booster coverages sample subpopulations correctly', {
 
   listener <- create_pev_booster_listener(
     variables = variables,
-    coverage = .9,
+    coverage = parameters$mass_pev_booster_coverage,
+    pev_distribution_timesteps = parameters$mass_pev_timesteps,
     booster_number = 1,
     pev_profile_index = 2,
     next_booster_event = events$mass_pev_boosters[[2]],
@@ -442,8 +484,8 @@ test_that('mass pev targets correct age and respects min_wait', {
     min_ages = 0,
     max_ages = 19 * 365,
     min_wait = 2*365,
-    booster_timestep = c(1, 6) * 30,
-    booster_coverage = c(.9, .8),
+    booster_spacing = c(1, 6) * 30,
+    booster_coverage = matrix(c(.9, .8, .9, .8), nrow=2, ncol=2),
     booster_profile = list(rtss_booster_profile, rtss_booster_profile)
   )
   events <- create_events(parameters)
@@ -504,8 +546,8 @@ test_that('Mass efficacy listener works correctly', {
     min_ages = c(1, 2, 3, 18) * 365,
     max_ages = (c(1, 2, 3, 18) + 1) * 365 - 1,
     min_wait = 0,
-    booster_timestep = c(1, 6) * 30,
-    booster_coverage = c(.9, .8),
+    booster_spacing = c(1, 6) * 30,
+    booster_coverage = matrix(c(.9, .8), nrow=1, ncol=2),
     booster_profile = list(rtss_booster_profile, rtss_booster_profile)
   )
 
@@ -542,8 +584,8 @@ test_that('Mass dose events are not ruined by lazy evaluation', {
     min_ages = c(1, 2, 3, 18) * 365,
     max_ages = (c(1, 2, 3, 18) + 1) * 365 - 1,
     min_wait = 0,
-    booster_timestep = c(1, 6, 12) * 30,
-    booster_coverage = c(.9, .8, .7),
+    booster_spacing = c(1, 6, 12) * 30,
+    booster_coverage = matrix(c(.9, .8, .7), nrow=1, ncol=3),
     booster_profile = list(rtss_booster_profile, rtss_booster_profile, rtss_booster_profile)
   )
 
@@ -553,9 +595,10 @@ test_that('Mass dose events are not ruined by lazy evaluation', {
   attach_pev_dose_listeners(
     variables = variables,
     parameters = parameters,
+    pev_distribution_timesteps = parameters$mass_pev_timesteps,
     dose_events = events$mass_pev_doses,
     booster_events = events$mass_pev_boosters,
-    booster_delays = parameters$mass_pev_booster_timestep,
+    booster_delays = parameters$mass_pev_booster_spacing,
     booster_coverages = parameters$mass_pev_booster_coverage,
     pev_profile_indices = parameters$mass_pev_profile_indices,
     strategy = 'mass',
@@ -577,8 +620,8 @@ test_that('Mass dose events are not ruined by lazy evaluation', {
   expect_equal(
     as.list(environment(
       events$mass_pev_boosters[[1]]$.listeners[[1]]
-    ))$coverage,
-    .9
+    ))$booster_number,
+    1
   )
 })
 
@@ -612,4 +655,101 @@ test_that('Efficacies are calculated correctly', {
   )
 })
 
+test_that('pev timed booster coverage can select the first coverage for the first booster', {
+  timestep <- 5 * 365
+  parameters <- get_parameters(list(human_population = 5))
+  parameters <- set_pev_epi(
+    parameters,
+    profile = rtss_profile,
+    timesteps = 10,
+    coverages = 0.8,
+    min_wait = 6 * 30,
+    age = 18 * 365,
+    booster_spacing = c(3, 12) * 30,
+    booster_coverage = matrix(c(.9, .8), nrow=1, ncol=2),
+    booster_profile = list(rtss_booster_profile, rtss_booster_profile),
+  )
+  events <- create_events(parameters)
+
+  booster_event <- mock_event(events$pev_epi_boosters[[1]])
 
+  listener <- create_pev_booster_listener(
+    variables = create_variables(parameters),
+    coverage = parameters$pev_epi_booster_coverage,
+    pev_distribution_timesteps = parameters$pev_epi_timesteps,
+    booster_number = 1,
+    pev_profile_index = 2,
+    next_booster_event = booster_event,
+    next_booster_delay = 9 * 30,
+    renderer = mock_render(timestep),
+    strategy = 'epi'
+  )
+
+  target <- individual::Bitset$new(5)$insert(seq(5))
+
+  mock_sample_bitset = mockery::mock(individual::Bitset$new(5)$insert(c(1, 2)))
+  mockery::stub(
+    listener,
+    'sample_bitset',
+    mock_sample_bitset
+  )
+
+  listener(timestep, target)
+
+  mockery::expect_args(
+    mock_sample_bitset,
+    1,
+    target,
+    .9
+  )
+})
+
+
+test_that('pev boosters can select the second coverage for the first booster', {
+  timestep <- 5 * 365
+  parameters <- get_parameters(list(human_population = 5))
+  parameters <- set_pev_epi(
+    parameters,
+    profile = rtss_profile,
+    timesteps = c(10, 30),
+    coverages = c(0.8, 0.4),
+    min_wait = 6 * 30,
+    age = 18 * 365,
+    booster_spacing = c(3, 12) * 30,
+    booster_coverage = matrix(c(.9, .45, .8, .8), nrow=2, ncol=2),
+    booster_profile = list(rtss_booster_profile, rtss_booster_profile),
+  )
+  events <- create_events(parameters)
+
+  booster_event <- mock_event(events$pev_epi_boosters[[1]])
+
+  listener <- create_pev_booster_listener(
+    variables = create_variables(parameters),
+    coverage = parameters$pev_epi_booster_coverage,
+    pev_distribution_timesteps = parameters$pev_epi_timesteps,
+    booster_number = 1,
+    pev_profile_index = 2,
+    next_booster_event = booster_event,
+    next_booster_delay = 9 * 30,
+    renderer = mock_render(timestep),
+    strategy = 'epi'
+  )
+
+  target <- individual::Bitset$new(5)$insert(seq(5))
+
+  mock_sample_bitset = mockery::mock(individual::Bitset$new(5)$insert(c(1, 2)))
+  mockery::stub(
+    listener,
+    'sample_bitset',
+    mock_sample_bitset
+  )
+
+  listener(timestep, target)
+
+  mockery::expect_args(
+    mock_sample_bitset,
+    1,
+    target,
+    .45
+  )
+})
diff --git a/vignettes/Vaccines.Rmd b/vignettes/Vaccines.Rmd
index 51bffd0e..02eb9358 100644
--- a/vignettes/Vaccines.Rmd
+++ b/vignettes/Vaccines.Rmd
@@ -158,8 +158,8 @@ rtssmassparams <- set_mass_pev(
   min_wait = 0, # The minimum acceptable time since the last vaccination is 0 because in our case we are only implementing one round of vaccination.
   min_ages = 5 * month, # The minimum age for the target population to be vaccinated. 
   max_ages = 50 * year, # The maximum age for the target population to be vaccinated.
-  booster_timestep = 12 * month, # The booster is given at 12 months after the primary series. 
-  booster_coverage = 0.95, # Coverage of the booster dose is 95%.
+  booster_spacing = 12 * month, # The booster is given at 12 months after the primary series. 
+  booster_coverage = matrix(0.95), # Coverage of the booster dose is 95%.
   booster_profile = list(rtss_booster_profile) # We will model implementation of the RTSS booster.
 )
 
@@ -228,8 +228,8 @@ seasmass_simparams <- set_mass_pev(
   min_ages = 5 * month, # The minimum age for the target population to be vaccinated. 
   max_ages = 50 * year, # The maximum age for the target population to be vaccinated.
   min_wait = 0, # There is no minimum wait between the last vaccination.
-  booster_timestep = round(c(12 * month + 2 * month)), # The booster is given 14 months after the first dose. 
-  booster_coverage = 1, # 100% of the vaccinated population is boosted.
+  booster_spacing = round(c(12 * month + 2 * month)), # The booster is given 14 months after the first dose. 
+  booster_coverage = matrix(1), # 100% of the vaccinated population is boosted.
   booster_profile = list(rtss_booster_profile) # We will model implementation of the RTSS booster.
   ) 
 
@@ -269,8 +269,8 @@ rtssepiparams <- set_pev_epi(
   coverages = 1, # Vaccine coverage is 100%.
   min_wait = 0, # There is no minimum wait since the last vaccination.
   age = 5 * month, # Individuals will be vaccinated once they reach 5 months of age.
-  booster_timestep = 12 * month, # The booster is administered 12 months following the third dose. 
-  booster_coverage = 0.95, # 95% of those vaccinated with the primary series will be boosted.
+  booster_spacing = 12 * month, # The booster is administered 12 months following the third dose. 
+  booster_coverage = matrix(0.95), # 95% of those vaccinated with the primary series will be boosted.
   booster_profile = list(rtss_booster_profile) # We will model implementation of the RTSS booster.
 )
 
@@ -322,8 +322,8 @@ rtssepiseasonalparams <- set_pev_epi(
   coverages = 1, # Vaccine coverage is 100%.
   min_wait = 6 * month, # When seasonal_boosters = TRUE, this is the minimum time between an individual receiving the final dose and the first booster.
   age = 5 * month, # Individuals will be vaccinated once they reach 5 months of age.
-  booster_timestep = peak - month * 3.5 , # Because seasonal_boosters = TRUE, the timestep here is relative to the start of the year. Here, we will give a booster at 3.5 months prior to peak transmission. 
-  booster_coverage = 0.95, # 95% of the vaccinated population is boosted.
+  booster_spacing = peak - month * 3.5 , # Because seasonal_boosters = TRUE, the timestep here is relative to the start of the year. Here, we will give a booster at 3.5 months prior to peak transmission. 
+  booster_coverage = matrix(0.95), # 95% of the vaccinated population is boosted.
   seasonal_boosters = TRUE, # Boosters will be given based on a seasonal schedule, so the timing in the boosters= argument above will be relative to the start of the year instead of relative to the 3rd dose.
   booster_profile = list(rtss_booster_profile) # We will model implementation of the RTSS booster.
 )
@@ -362,8 +362,8 @@ rtssepiparams2 <- set_pev_epi(
   coverages = 1, # Vaccine coverage is 100%.
   age = 5 * month, # Individuals will be vaccinated once they reach 5 months of age.
   min_wait = 0, # When seasonal_boosters = FALSE, this is the minimum time between doses.
-  booster_timestep = c(12 * month, 24 * month), # Here, we are testing a strategy with 2 boosters, one at 1 year after the 3rd dose and the second 2 years after the 3rd dose.
-  booster_coverage = c(1, 1), # For each of the two boosters, coverage is 100%.
+  booster_spacing = c(12 * month, 24 * month), # Here, we are testing a strategy with 2 boosters, one at 1 year after the 3rd dose and the second 2 years after the 3rd dose.
+  booster_coverage = matrix(c(1, 1), nrow=1, ncol=2), # For each of the two boosters, coverage is 100%.
   booster_profile = list(rtss_booster_profile, rtss_booster_profile) # We will model implementation of the RTSS booster.
 )