diff --git a/tests/testthat/test-model_default.R b/tests/testthat/test-model_default.R
index 00c844ee..a10b97d1 100644
--- a/tests/testthat/test-model_default.R
+++ b/tests/testthat/test-model_default.R
@@ -243,6 +243,22 @@ test_that("Default model: vaccination and stats. correctness", {
   expect_true(
     all(epidemic_size(data_baseline) > epidemic_size(data))
   )
+
+  # expect that high vaccination rates (± 1% per day, e.g. Covid-19 vax)
+  # give statistically correct results (no negative susceptibles at any time)
+  high_rate_vax <- vaccination(
+    time_begin = matrix(200, nrow(contact_matrix)),
+    time_end = matrix(200 + 150, nrow(contact_matrix)),
+    nu = matrix(0.01, nrow = nrow(contact_matrix))
+  )
+  data <- model_default(
+    uk_population,
+    vaccination = high_rate_vax, time_end = 600
+  )
+  checkmate::expect_numeric(
+    data[grepl("susceptible", data$compartment, fixed = TRUE), ]$value,
+    lower = 0
+  )
 })
 
 test_that("Default model: time dependence", {
diff --git a/tests/testthat/test-model_vacamole.R b/tests/testthat/test-model_vacamole.R
index c79237f8..70ad0f1f 100644
--- a/tests/testthat/test-model_vacamole.R
+++ b/tests/testthat/test-model_vacamole.R
@@ -343,6 +343,31 @@ test_that("Vacamole model: two dose vaccination and stats. correctness", {
     data[grepl("dose", data$compartment, fixed = TRUE) &
       time %in% seq(50, 55), ]
   )
+
+  # expect that high vaccination rates (± 1% per day, e.g. Covid-19 vax)
+  # give statistically correct results (no negative susceptibles at any time)
+  high_rate_vax <- vaccination(
+    nu = matrix(
+      c(1e-2, 1e-2), # 1% per day
+      nrow = 2, ncol = 2, byrow = TRUE
+    ),
+    time_begin = matrix(
+      c(0, 50),
+      nrow = 2, ncol = 2, byrow = TRUE # second dose given from t = 50 onwards
+    ),
+    time_end = matrix(
+      100,
+      nrow = 2, ncol = 2
+    )
+  )
+  data <- model_vacamole(
+    uk_population,
+    vaccination = high_rate_vax, time_end = 600
+  )
+  checkmate::expect_numeric(
+    data[grepl("susceptible", data$compartment, fixed = TRUE), ]$value,
+    lower = 0
+  )
 })
 
 test_that("Vacamole model: time dependence", {
diff --git a/vignettes/modelling_vaccination.Rmd b/vignettes/modelling_vaccination.Rmd
index 6975ddda..64e04454 100644
--- a/vignettes/modelling_vaccination.Rmd
+++ b/vignettes/modelling_vaccination.Rmd
@@ -22,8 +22,6 @@ knitr::opts_chunk$set(
 
 ```{r setup}
 library(epidemics)
-library(dplyr)
-library(ggplot2)
 ```
 
 ## Prepare population and initial conditions
@@ -96,7 +94,15 @@ vaccinate_elders <- vaccination(
 vaccinate_elders
 ```
 
-## Run epidemic model
+::: {.alert .alert-warning}
+**Note that** when the vaccination rate is high, the number of individuals who could _potentially_ transition from 'susceptible' to 'vaccinated' may be greater than the number of individuals in the 'susceptible' compartment: $dS_{i_{t}} = \nu_{i_t}S_{i_{0}} > S_{i_{t}}$ for each demographic group $i$ at time $t$.
+
+This could realistically occur when there are more doses available than individuals eligible to receive them, for instance towards the end of an epidemic.
+
+_epidemics_ automatically handles such situations by setting $dS_{i_{t}}$ to be the minimum of doses or eligible individuals: : $dS_{i_t} = \text{Min}(\nu_{i_t}S_{i_0}, S_{i_t})$, such that $S_{i_{t+1}}$ does not take a negative value.
+:::
+
+The `<vaccination>` object can be passed to the `vaccination` argument of a model function call.
 
 ```{r}
 # run an epidemic model using `epidemic`
@@ -107,42 +113,6 @@ output <- model_default(
 )
 ```
 
-## Prepare data and visualise infections
-
-Plot epidemic over time, showing only the number of individuals in the exposed, infected, and vaccinated compartments.
-
-```{r class.source = 'fold-hide'}
-# plot figure of epidemic curve
-filter(output, compartment %in% c("exposed", "infectious")) %>%
-  ggplot(
-    aes(
-      x = time,
-      y = value,
-      col = demography_group,
-      linetype = compartment
-    )
-  ) +
-  geom_line() +
-  scale_y_continuous(
-    labels = scales::comma
-  ) +
-  scale_colour_brewer(
-    palette = "Dark2",
-    name = "Age group"
-  ) +
-  expand_limits(
-    y = c(0, 500e3)
-  ) +
-  coord_cartesian(
-    expand = FALSE
-  ) +
-  theme_bw() +
-  theme(
-    legend.position = "top"
-  ) +
-  labs(
-    x = "Simulation time (days)",
-    linetype = "Compartment",
-    y = "Individuals"
-  )
-```
+::: {.alert .alert-info}
+**Note that** vaccination modelling is currently only supported for the 'default' (`model_default()`) and 'Vacamole' (`model_vacamole()`) models.
+:::