diff --git a/R/compatibility.R b/R/compatibility.R
index 9324f5db..ec71d3df 100644
--- a/R/compatibility.R
+++ b/R/compatibility.R
@@ -184,7 +184,7 @@ remove_unused_equilibrium <- function(params) {
 #' equilibrium parameters and set up the initial human and mosquito population
 #' to acheive init_EIR
 #' @param parameters model parameters to update
-#' @param init_EIR the desired initial EIR (infectious bites per day over the entire human
+#' @param init_EIR the desired initial EIR (infectious bites per person per day over the entire human
 #' population)
 #' @param eq_params parameters from the malariaEquilibrium package, if null.
 #' The default malariaEquilibrium parameters will be used
diff --git a/R/events.R b/R/events.R
index 2c17e30d..c01646ad 100644
--- a/R/events.R
+++ b/R/events.R
@@ -27,7 +27,7 @@ create_events <- function(parameters) {
   }
 
   # EPI vaccination events
-  if (!is.null(parameters$rtss_epi_start)) {
+  if (!is.null(parameters$rtss_epi_coverage)) {
     rtss_epi_doses <- lapply(
       seq_along(parameters$rtss_doses),
       function(.) individual::TargetedEvent$new(parameters$human_population)
diff --git a/R/processes.R b/R/processes.R
index 133ebb49..44402861 100644
--- a/R/processes.R
+++ b/R/processes.R
@@ -124,14 +124,16 @@ create_processes <- function(
   # =========
   # RTS,S EPI
   # =========
-  if (!is.null(parameters$rtss_epi_start)) {
+  if (!is.null(parameters$rtss_epi_coverages)) {
     processes <- c(
       processes,
       create_rtss_epi_process(
         variables,
         events,
         parameters,
-        correlations
+        correlations,
+        parameters$rtss_epi_coverages,
+        parameters$rtss_epi_timesteps
       )
     )
   }
diff --git a/R/rtss.R b/R/rtss.R
index 6d603777..5a407596 100644
--- a/R/rtss.R
+++ b/R/rtss.R
@@ -12,16 +12,23 @@ create_rtss_epi_process <- function(
   variables,
   events,
   parameters,
-  correlations
+  correlations,
+  coverages,
+  timesteps
   ) {
   function(timestep) {
-    if (!between(timestep, parameters$rtss_epi_start, parameters$rtss_epi_end)) {
+    timestep_index <- match_timestep(ts = timesteps, t = timestep)
+    if(timestep_index == 0){
+      return()
+    }
+    coverage <- coverages[timestep_index]
+    if(coverage == 0){
       return()
     }
+    
     to_vaccinate <- variables$birth$get_index_of(
       set = timestep - parameters$rtss_epi_age
     )
-
     if (parameters$rtss_epi_min_wait == 0) {
       target <- to_vaccinate$to_vector()
     } else {
@@ -36,10 +43,11 @@ create_rtss_epi_process <- function(
       sample_intervention(
         target,
        'rtss',
-        parameters$rtss_epi_coverage,
+        coverage,
         correlations
       )
     ]
+
     schedule_vaccination(
       target,
       events,
diff --git a/R/utils.R b/R/utils.R
index 343e6309..f25244bd 100644
--- a/R/utils.R
+++ b/R/utils.R
@@ -64,4 +64,3 @@ rtexp <- function(n, m, t) { itexp(runif(n), m, t) }
 match_timestep <- function(ts, t) {
   min(sum(ts <= t), length(ts))
 }
-
diff --git a/R/vaccine_parameters.R b/R/vaccine_parameters.R
index 164be3f9..2f5fe810 100644
--- a/R/vaccine_parameters.R
+++ b/R/vaccine_parameters.R
@@ -4,9 +4,8 @@
 #' age. Efficacy will take effect after the last dose
 #'
 #' @param parameters a list of parameters to modify
-#' @param start timestep to turn on epi vaccination
-#' @param end timestep to turn off epi vaccination
-#' @param coverage the coverage for the starter doses
+#' @param coverages a vector of coverages for the starter doses
+#' @param timesteps a vector of timesteps associated with coverages
 #' @param age for the target population, (in timesteps)
 #' @param min_wait the minimum acceptable time since the last vaccination (in
 #' timesteps); When seasonal_boosters = TRUE, this represents the minimum time
@@ -21,19 +20,18 @@
 #' @export
 set_rtss_epi <- function(
   parameters,
-  start,
-  end,
-  coverage,
+  coverages,
+  timesteps,
   age,
   min_wait,
   boosters,
   booster_coverage,
   seasonal_boosters = FALSE
   ) {
-  stopifnot(length(start) == 1 && start > 1)
-  stopifnot(length(end) == 1 && end > start)
+  if(length(coverages) != length(timesteps)){
+    stop("coverages and timesteps must align")
+  }
   stopifnot(min_wait >= 0)
-  stopifnot(coverage >= 0 & coverage <= 1)
   stopifnot(age >= 0)
   stopifnot(is.logical(seasonal_boosters))
   if (seasonal_boosters) {
@@ -47,9 +45,8 @@ set_rtss_epi <- function(
     stop('booster and booster_coverage does not align')
   }
   parameters$rtss <- TRUE
-  parameters$rtss_epi_start <- start
-  parameters$rtss_epi_end <- end
-  parameters$rtss_epi_coverage <- coverage
+  parameters$rtss_epi_coverages <- coverages
+  parameters$rtss_epi_timesteps <- timesteps
   parameters$rtss_epi_age <- age
   parameters$rtss_epi_boosters <- boosters
   parameters$rtss_epi_min_wait <- min_wait
diff --git a/man/set_equilibrium.Rd b/man/set_equilibrium.Rd
index f1c83c54..c1b8eb7b 100644
--- a/man/set_equilibrium.Rd
+++ b/man/set_equilibrium.Rd
@@ -9,8 +9,7 @@ set_equilibrium(parameters, init_EIR, eq_params = NULL)
 \arguments{
 \item{parameters}{model parameters to update}
 
-\item{init_EIR}{the desired initial EIR (infectious bites per day over the entire human
-population)}
+\item{init_EIR}{the desired initial EIR (infectious bites per person per day over the entire human population)}
 
 \item{eq_params}{parameters from the malariaEquilibrium package, if null.
 The default malariaEquilibrium parameters will be used}
diff --git a/man/set_rtss_epi.Rd b/man/set_rtss_epi.Rd
index 3d759ef0..78d6a2cb 100644
--- a/man/set_rtss_epi.Rd
+++ b/man/set_rtss_epi.Rd
@@ -6,9 +6,8 @@
 \usage{
 set_rtss_epi(
   parameters,
-  start,
-  end,
-  coverage,
+  coverages,
+  timesteps,
   age,
   min_wait,
   boosters,
@@ -19,11 +18,9 @@ set_rtss_epi(
 \arguments{
 \item{parameters}{a list of parameters to modify}
 
-\item{start}{timestep to turn on epi vaccination}
+\item{coverages}{a vector of coverages for the starter doses}
 
-\item{end}{timestep to turn off epi vaccination}
-
-\item{coverage}{the coverage for the starter doses}
+\item{timesteps}{a vector of timesteps associated with coverages}
 
 \item{age}{for the target population, (in timesteps)}
 
diff --git a/tests/testthat/test-rtss-epi.R b/tests/testthat/test-rtss-epi.R
index 54a593f9..0d4eef98 100644
--- a/tests/testthat/test-rtss-epi.R
+++ b/tests/testthat/test-rtss-epi.R
@@ -2,18 +2,16 @@ test_that('RTS,S epi strategy parameterisation works', {
   parameters <- get_parameters()
   parameters <- set_rtss_epi(
     parameters,
-    start = 10,
-    end = 100,
-    coverage = 0.8,
+    coverages = c(0.1, 0.8),
+    timesteps = c(10, 100),
     min_wait = 0,
     age = 5 * 30,
     boosters = c(18, 36) * 30,
     booster_coverage = c(.9, .8)
   )
   expect_equal(parameters$rtss, TRUE)
-  expect_equal(parameters$rtss_epi_start, 10)
-  expect_equal(parameters$rtss_epi_end, 100)
-  expect_equal(parameters$rtss_epi_coverage, .8)
+  expect_equal(parameters$rtss_epi_coverages, c(0.1, 0.8))
+  expect_equal(parameters$rtss_epi_timesteps, c(10, 100))
   expect_equal(parameters$rtss_epi_age, 5 * 30)
   expect_equal(parameters$rtss_epi_min_wait, 0)
   expect_equal(parameters$rtss_epi_boosters, c(18, 36) * 30)
@@ -24,9 +22,8 @@ test_that('RTS,S epi fails pre-emptively', {
   expect_error(
     set_rtss_epi(
       parameters,
-      start = 10,
-      end = 100,
-      coverages = 0.8,
+      coverages = c(0.1, 0.8),
+      timesteps = c(10, 100),
       min_wait = 0,
       min_ages = 5 * 30,
       max_ages = 17 * 30,
@@ -42,9 +39,8 @@ test_that('RTS,S epi targets correct age and respects min_wait', {
   parameters <- get_parameters(list(human_population = 5))
   parameters <- set_rtss_epi(
     parameters,
-    start = 10,
-    end = timestep,
-    coverage = 0.8,
+    timesteps = 10,
+    coverages = 0.8,
     min_wait = 2*365,
     age = 18 * 365,
     boosters = c(18, 36) * 30,
@@ -65,7 +61,9 @@ test_that('RTS,S epi targets correct age and respects min_wait', {
     variables,
     events,
     parameters,
-    get_correlation_parameters(parameters)
+    get_correlation_parameters(parameters),
+    parameters$rtss_epi_coverages,
+    parameters$rtss_epi_timesteps
   )
 
   mockery::stub(
@@ -103,9 +101,8 @@ test_that('RTS,S EPI respects min_wait when scheduling seasonal boosters', {
   parameters <- get_parameters(list(human_population = 5))
   parameters <- set_rtss_epi(
     parameters,
-    start = 10,
-    end = timestep,
-    coverage = 0.8,
+    timesteps = 10,
+    coverages = 0.8,
     min_wait = 6 * 30,
     age = 18 * 365,
     boosters = c(3, 12) * 30,
@@ -138,9 +135,8 @@ test_that('RTS,S EPI schedules for the following year with seasonal boosters', {
   parameters <- get_parameters(list(human_population = 5))
   parameters <- set_rtss_epi(
     parameters,
-    start = 10,
-    end = timestep,
-    coverage = 0.8,
+    timesteps = 10,
+    coverages = 0.8,
     min_wait = 6 * 30,
     age = 18 * 365,
     boosters = c(3, 12) * 30,
diff --git a/vignettes/Vaccines.Rmd b/vignettes/Vaccines.Rmd
index fb434ed0..6e75b950 100644
--- a/vignettes/Vaccines.Rmd
+++ b/vignettes/Vaccines.Rmd
@@ -25,171 +25,123 @@ We are going to set the default parameters to run the simulation from an equilib
 
 ```{r}
 year <- 365
+month <- 30
 sim_length <- 3 * year
-human_population <- 1000
-starting_EIR <- 50
+human_population <- 10000
+starting_EIR <- 20
 
 simparams <- get_parameters(list(
     human_population = human_population,
-    model_seasonality = TRUE, # Let's try a bi-modal model
-    g0 = 0.28605,
-    g = c(0.20636, -0.0740318, -0.0009293),
-    h = c(0.173743, -0.0730962, -0.116019),
-    prevalence_rendering_min_ages = 0,
-    prevalence_rendering_max_ages = 18 * year,
     incidence_rendering_min_ages = 0,
-    incidence_rendering_max_ages = 18 * year,
+    incidence_rendering_max_ages = 5 * year,
     individual_mosquitoes = FALSE
   )
 )
 
 simparams <- set_equilibrium(simparams, starting_EIR)
 
-# Plotting functions
-plot_prevalence <- function(output) {
-  ggplot(output) + geom_line(
-    aes(x = timestep, y = (n_inc_0_6570 / n_0_6570))) + 
-    labs(x = "timestep", y = "incidence 0-18")
-}
-
-add_intervention_lines <- function(plot, events) {
-  plot + geom_vline(
-    data = events,
-    mapping = aes(xintercept=timestep),
-    color="blue"
-  ) + geom_text(
-    data = events,
-    mapping = aes(x = timestep, y = 0, label = name),
-    size = 4,
-    angle = 90,
-    vjust = -0.4,
-    hjust = 0
-  )
-}
 ```
 
 Then we can run the simulation for a variety of Vaccination strategies:
 
 ## Mass RTS,S
 
-This is a round of RTS,S vaccine for individuals between 5 - 17 months and a booster after 18 months; coverage of 80%; for 10 years:
+This is a round of RTS,S vaccine for individuals between 5 months and 10 years followed by a booster after 18 months:
 
 ```{r}
 rtssparams <- simparams
 
-peak <- peak_season_offset(rtssparams)
-month <- 30
-
-# Add RTS,S strategy
-rtssevents = data.frame(
-  timestep = c(1, 2) * year + peak - month, #vaccine efficacy kicks off a month before the peak
-  name=c("RTS,S 1", "RTS,S 2")
-)
-
 rtssparams <- set_mass_rtss(
   rtssparams,
-  timesteps = rtssevents$timestep,
-  coverages = rep(0.8, 2),
-  min_wait = 2 * 365,
+  timesteps = 1 * year,
+  coverages = 1,
+  min_wait = 0,
   min_ages = 5 * month,
-  max_ages = 17 * year,
+  max_ages = 10 * year,
   boosters = 18 * month,
-  booster_coverage = .7
+  booster_coverage = 0.95
 )
 
 output <- run_simulation(sim_length, rtssparams)
+output$clinical_incidence <- 1000 * output$n_inc_0_1825 / output$n_0_1825
 
-add_intervention_lines(plot_prevalence(output), rtssevents)
-```
-
+ggplot(data = output, aes(x = timestep / 365, y = clinical_incidence)) + 
+  geom_line(col = "grey80") +
+  stat_smooth(col = "darkblue", se = FALSE) +
+  geom_vline(xintercept = 1, col = "red") +
+  xlab("year") +
+  ylab("Clinical incidence \n (per 1000 children aged 0-5)") +
+  theme_bw()
 
-You can plot the distribution of doses using the n_rtss_mass_dose_* outputs:
+```
 
+You can look at the distribution of doses using the n_rtss_mass_dose\_\* outputs:
 
 ```{r}
-plot_dosage <- function(output, strategy, doses, boosters) {
-  distributed <- NULL
-  label <- NULL
-  timestep <- NULL
-  for (dose in seq(doses)) {
-    distributed <- c(
-      distributed,
-      cumsum(output[[paste0('n_rtss_', strategy, '_dose_', dose)]])
-    )
-    label <- c(label, rep(paste0('dose ', dose), nrow(output)))
-    timestep <- c(timestep, output$timestep)
-  }
-  for (booster in seq(boosters)) {
-    distributed <- c(
-      distributed,
-      cumsum(output[[paste0('n_rtss_', strategy, '_booster_', booster)]])
-    )
-    label <- c(label, rep(paste0('booster ', booster), nrow(output)))
-    timestep <- c(timestep, output$timestep)
-  }
-  long_output <- data.frame(
-    distributed = distributed,
-    label = label,
-    timestep = timestep
-  )
-  ggplot(long_output) + geom_line(aes(x = timestep, y = distributed, group=label, color = label))
-}
+dose_data <- data.frame(timestep = output$timestep,
+                        dose = rep(c("mass 1", "mass 2", "mass 3", "mass booster"), each = length(output$timestep)),
+                        n = c(output$n_rtss_mass_dose_1, output$n_rtss_mass_dose_2, output$n_rtss_mass_dose_3,
+                              output$n_rtss_mass_booster_1))
+dose_data[dose_data$n > 0, ]
 
-add_intervention_lines(plot_dosage(output, 'mass', 3, 1), rtssevents)
 ```
 
 ## RTS,S EPI
 
-You can opt for a more gradual dosing using the EPI strategy. Individuals will be vaccinated once they reach a certain age...
+You can opt for a more gradual dosing using the EPI strategy. 
+Individuals will be vaccinated once they reach 5 months. For this intervention
+we see a much more gradual increase in impact following implementation
 
 ```{r}
 rtssepiparams <- simparams
 
 # Add RTS,S strategy
-rtssepievents = data.frame(
-  timestep = c(1, 2) * year,
-  name=c("RTS,S EPI start", "RTS,S EPI end")
-)
-
 rtssepiparams <- set_rtss_epi(
   rtssepiparams,
-  start = rtssepievents$timestep[[1]],
-  end = rtssepievents$timestep[[2]],
-  age = 5 * month,
-  coverage = 0.8,
+  timesteps = 1 * year,
+  coverages = 1,
   min_wait = 0,
+  age = 5 * month,
   boosters = 18 * month,
-  booster_coverage = .7
+  booster_coverage = 0.95
 )
 
-output <- run_simulation(sim_length, rtssepiparams)
+output <- run_simulation(sim_length * 2, rtssepiparams)
+output$clinical_incidence <- 1000 * output$n_inc_0_1825 / output$n_0_1825
 
-add_intervention_lines(plot_dosage(output, 'epi', 3, 1), rtssepievents)
+ggplot(data = output, aes(x = timestep / 365, y = clinical_incidence)) + 
+  geom_line(col = "grey80") +
+  stat_smooth(col = "darkblue", se = FALSE) +
+  geom_vline(xintercept = 1, col = "red") +
+  xlab("year") +
+  ylab("Clinical incidence \n (per 1000 children aged 0-5)") +
+  theme_bw()
 ```
 
 ## RTS,S seasonal boosters
 
-We can set booster timesteps relative to the start of the year. This allows us to target seasonal dynamics
+In a seasonal setting, we can set booster timesteps relative to the start of the year.
+This allows us to target seasonal dynamics
 
 ```{r}
 rtssepiseasonalparams <- simparams
+rtssepiseasonalparams$model_seasonality = TRUE
+rtssepiseasonalparams$g0 = 0.28605
+rtssepiseasonalparams$g = c(0.20636, -0.0740318, -0.0009293)
+rtssepiseasonalparams$h = c(0.173743, -0.0730962, -0.116019)
+peak <- peak_season_offset(rtssepiseasonalparams)
 
 # Add RTS,S seasonal strategy
 rtssepiseasonalparams <- set_rtss_epi(
   rtssepiseasonalparams,
-  start = rtssepievents$timestep[[1]],
-  end = rtssepievents$timestep[[2]],
-  age = 5 * month,
-  coverage = 0.8,
+  timesteps = 1 * year,
+  coverages = 1,
   min_wait = 6 * month,
+  age = 5 * month,
   boosters = (peak - 3 * month) + c(0, year),
   booster_coverage = rep(.7, 2),
   seasonal_boosters = TRUE
 )
-
-output <- run_simulation(5 * year, rtssepiseasonalparams)
-
-add_intervention_lines(plot_dosage(output, 'epi', 3, 2), rtssepievents)
 ```
 
 ## RTS,S dosing
@@ -201,40 +153,40 @@ rtssepiparams2 <- rtssepiparams
 rtssepiparams2$rtss_doses <- c(0, 30, 60)
 rtssepiparams2 <- set_rtss_epi(
   rtssepiparams2,
-  start = rtssepievents$timestep[[1]],
-  end = rtssepievents$timestep[[2]],
+  timesteps = 1 * year,
+  coverages = 1,
   age = 5 * month,
-  coverage = 0.8,
   min_wait = 0,
   boosters = c(12 * month, 18 * month, 24 * month),
   booster_coverage = c(1, 1, 1)
 )
-
-output <- run_simulation(5 * year, rtssepiparams2)
-add_intervention_lines(plot_dosage(output, 'epi', 3, 3), rtssepievents)
 ```
 
 ## TBV
 
-This is a round of the TBV (transmission blocking) vaccine at ages 1, 2, 3 and 18; coverage of 80%; for 10 years:
+We can also model vaccines with completely different modes of actions. For example
+a transmission blocking vaccines (TBV). This example shows 5 rounds of the TBV
+to everyone aged over 5
 
 ```{r}
 tbvparams <- simparams
 
-# Add TBV strategy
-tbvevents = data.frame(
-  timestep = c(1, 2) * 365,
-  name=c("TBV 1", "TBV 2")
-)
-
 tbvparams <- set_tbv(
   tbvparams,
-  timesteps = tbvevents$timestep,
-  coverages = rep(0.9, 2),
-  ages = seq(18)
+  timesteps = round(c(1, 1.25, 1.5, 1.75, 2) * 365),
+  coverages = rep(0.99, 5),
+  ages = 5:60
 )
 
 output <- run_simulation(sim_length, tbvparams)
+output$clinical_incidence <- 1000 * output$n_inc_0_1825 / output$n_0_1825
+
+ggplot(data = output, aes(x = timestep / 365, y = clinical_incidence)) + 
+  geom_line(col = "grey80") +
+  stat_smooth(col = "darkblue", se = FALSE) +
+  geom_vline(xintercept = 1, col = "red") +
+  xlab("year") +
+  ylab("Clinical incidence \n (per 1000 children aged 0-5)") +
+  theme_bw()
 
-add_intervention_lines(plot_prevalence(output), tbvevents)
 ```