diff --git a/.Rbuildignore b/.Rbuildignore
index 14fb78e2..d969fb5a 100644
--- a/.Rbuildignore
+++ b/.Rbuildignore
@@ -9,4 +9,3 @@ docker
 LICENSE.md
 codecov.yml
 .github
-vignettes/*
diff --git a/DESCRIPTION b/DESCRIPTION
index ebf727d7..1affbf8b 100644
--- a/DESCRIPTION
+++ b/DESCRIPTION
@@ -1,6 +1,6 @@
 Package: malariasimulation
 Title: An individual based model for malaria
-Version: 1.1.0
+Version: 1.2.0
 Authors@R: c(
   person(
     given = "Giovanni",
@@ -19,9 +19,9 @@ License: MIT + file LICENSE
 Encoding: UTF-8
 LazyData: true
 Remotes:
-  mrc-ide/malariaEquilibrium
+  mrc-ide/malariaEquilibrium,
 Imports:
-  individual,
+  individual (>= 0.1.7),
   malariaEquilibrium,
   Rcpp,
   statmod,
@@ -39,8 +39,9 @@ Suggests:
   DiagrammeR,
   cowplot,
   ggplot2,
-  covr
-RoxygenNote: 7.1.1
+  covr,
+  mgcv 
+RoxygenNote: 7.1.2
 Roxygen: list(markdown = TRUE)
 LinkingTo:
   Rcpp,
diff --git a/NEWS.md b/NEWS.md
new file mode 100644
index 00000000..0a813cf8
--- /dev/null
+++ b/NEWS.md
@@ -0,0 +1,7 @@
+# malariasimulation 1.2.0
+
+  * added a `news.md` file to track changes to the package.
+  * n_inc_clinical includes treated humans
+  * disaggregate EIR and state counts by mosquito species
+  * remove redundant Total_M and and global EIR outputs
+  * new vignette on how to calibrate EIR to prevalence
diff --git a/R/RcppExports.R b/R/RcppExports.R
index 249955b9..ccf76a71 100644
--- a/R/RcppExports.R
+++ b/R/RcppExports.R
@@ -65,6 +65,10 @@ bernoulli_multi_p_cpp <- function(p) {
     .Call(`_malariasimulation_bernoulli_multi_p_cpp`, p)
 }
 
+bitset_index_cpp <- function(a, b) {
+    .Call(`_malariasimulation_bitset_index_cpp`, a, b)
+}
+
 fast_weighted_sample <- function(size, probs) {
     .Call(`_malariasimulation_fast_weighted_sample`, size, probs)
 }
diff --git a/R/biting_process.R b/R/biting_process.R
index 15f65e2c..97d801f0 100644
--- a/R/biting_process.R
+++ b/R/biting_process.R
@@ -85,12 +85,13 @@ simulate_bites <- function(
   if (parameters$individual_mosquitoes) {
     infectious_index <- variables$mosquito_state$get_index_of('Im')
     susceptible_index <- variables$mosquito_state$get_index_of('Sm')
-    adult_index <- variables$mosquito_state$get_index_of('NonExistent')$not()
+    adult_index <- variables$mosquito_state$get_index_of('NonExistent')$not(TRUE)
   }
 
   EIR <- 0
 
   for (s_i in seq_along(parameters$species)) {
+    species_name <- parameters$species[[s_i]]
     solver_states <- solver_get_states(solvers[[s_i]])
     p_bitten <- prob_bitten(timestep, variables, s_i, parameters)
     Q0 <- parameters$Q0[[s_i]]
@@ -118,20 +119,24 @@ simulate_bites <- function(
 
     lagged_eir[[s_i]]$save(n_infectious * a, timestep)
     species_eir <- lagged_eir[[s_i]]$get(timestep - parameters$de)
+    renderer$render(paste0('EIR_', species_name), species_eir, timestep)
     EIR <- EIR + species_eir
-    n_bites <- rpois(1, species_eir * mean(psi))
-    if (n_bites > 0) {
-      bitten_humans$insert(
-        fast_weighted_sample(n_bites, lambda)
-      )
+    expected_bites <- species_eir * mean(psi)
+    if (expected_bites > 0) {
+      n_bites <- rpois(1, expected_bites)
+      if (n_bites > 0) {
+        bitten_humans$insert(
+          fast_weighted_sample(n_bites, lambda)
+        )
+      }
     }
 
     infectivity <- lagged_infectivity$get(timestep - parameters$delay_gam)
     lagged_infectivity$save(sum(human_infectivity * .pi), timestep)
     foim <- calculate_foim(a, infectivity)
-    renderer$render(paste0('FOIM_', s_i), foim, timestep)
+    renderer$render(paste0('FOIM_', species_name), foim, timestep)
     mu <- death_rate(f, W, Z, s_i, parameters)
-    renderer$render(paste0('mu_', s_i), mu, timestep)
+    renderer$render(paste0('mu_', species_name), mu, timestep)
 
     if (parameters$individual_mosquitoes) {
       # update the ODE with stats for ovoposition calculations
@@ -167,7 +172,6 @@ simulate_bites <- function(
     }
   }
 
-  renderer$render('EIR', EIR, timestep)
   renderer$render('n_bitten', bitten_humans$size(), timestep)
   bitten_humans
 }
@@ -191,7 +195,7 @@ effective_biting_rates <- function(a, .pi, p_bitten) {
 
 calculate_infectious <- function(species, solvers, variables, parameters) {
   if (parameters$individual_mosquitoes) {
-    adult_index <- variables$mosquito_state$get_index_of('NonExistent')$not()
+    adult_index <- variables$mosquito_state$get_index_of('NonExistent')$not(TRUE)
     return(
       calculate_infectious_individual(
         species,
@@ -221,7 +225,7 @@ calculate_infectious_individual <- function(
 }
 
 calculate_infectious_compartmental <- function(solver_states) {
-  solver_states[[ADULT_ODE_INDICES['Im']]]
+  max(solver_states[[ADULT_ODE_INDICES['Im']]], 0)
 }
 
 intervention_coefficient <- function(p_bitten) {
@@ -234,7 +238,7 @@ human_pi <- function(zeta, psi) {
 
 blood_meal_rate <- function(v, z, parameters) {
   gonotrophic_cycle <- get_gonotrophic_cycle(v, parameters)
-  interrupted_foraging_time <- parameters$foraging_time / (1 - z)
+  interrupted_foraging_time <- parameters$foraging_time[[v]] / (1 - z)
   1 / (interrupted_foraging_time + gonotrophic_cycle)
 }
 
diff --git a/R/compartmental.R b/R/compartmental.R
index 79d52d76..1e712b49 100644
--- a/R/compartmental.R
+++ b/R/compartmental.R
@@ -77,7 +77,7 @@ parameterise_solvers <- function(models, parameters) {
   )
 }
 
-create_ode_rendering_process <- function(renderer, solvers, parameters) {
+create_compartmental_rendering_process <- function(renderer, solvers, parameters) {
   if (parameters$individual_mosquitoes) {
     indices <- ODE_INDICES
   } else {
@@ -86,18 +86,17 @@ create_ode_rendering_process <- function(renderer, solvers, parameters) {
 
   function(timestep) {
     counts <- rep(0, length(indices))
-    for (i in seq_along(solvers)) {
-      row <- solver_get_states(solvers[[i]])
+    for (s_i in seq_along(solvers)) {
+      row <- solver_get_states(solvers[[s_i]])
+      for (i in seq_along(indices)) {
+        renderer$render(
+          paste0(names(indices)[[i]], '_', parameters$species[[s_i]], '_count'),
+          row[[i]],
+          timestep
+        )
+      }
       counts <- counts + row
     }
-
-    for (i in seq_along(indices)) {
-      renderer$render(
-        paste0(names(indices)[[i]], '_count'),
-        counts[[i]],
-        timestep
-      )
-    }
   }
 }
 
diff --git a/R/events.R b/R/events.R
index 66905042..941de05d 100644
--- a/R/events.R
+++ b/R/events.R
@@ -9,9 +9,6 @@ create_events <- function(parameters) {
     clinical_infection = individual::TargetedEvent$new(parameters$human_population),
     asymptomatic_infection = individual::TargetedEvent$new(parameters$human_population),
 
-    # whether the infection is detected
-    detection = individual::TargetedEvent$new(parameters$human_population), 
-
     # MDA events
     mda_administer = individual::Event$new(),
     smc_administer = individual::Event$new(),
@@ -135,6 +132,11 @@ attach_event_listeners <- function(
       parameters
     )
   )
+  events$asymptomatic_progression$add_listener(
+    function(timestep, target) {
+      variables$is_severe$queue_update('no', target)
+    }
+  )
 
   # Recovery events
   events$subpatent_progression$add_listener(
@@ -165,24 +167,6 @@ attach_event_listeners <- function(
   # Progression
   # ===========
   # When infection events fire, schedule the next stages of infection
-
-  events$clinical_infection$add_listener(
-    create_clinical_incidence_renderer(
-      variables$birth,
-      parameters,
-      renderer
-    )
-  )
-
-  events$detection$add_listener(
-    create_incidence_renderer(
-      variables$birth,
-      variables$is_severe,
-      parameters,
-      renderer
-    )
-  )
-
   if (parameters$individual_mosquitoes) {
     events$mosquito_infection$add_listener(
       individual::update_category_listener(variables$mosquito_state, 'Im')
diff --git a/R/human_infection.R b/R/human_infection.R
index 625a96ad..c0b246c3 100644
--- a/R/human_infection.R
+++ b/R/human_infection.R
@@ -33,6 +33,7 @@ simulate_infection <- function(
     variables,
     bitten_humans,
     parameters,
+    renderer,
     timestep
   )
 
@@ -56,7 +57,9 @@ simulate_infection <- function(
   clinical_infections <- calculate_clinical_infections(
     variables,
     infected_humans,
-    parameters
+    parameters,
+    renderer,
+    timestep
   )
 
   if (parameters$severe_enabled) {
@@ -65,16 +68,15 @@ simulate_infection <- function(
       clinical_infections,
       variables,
       infected_humans,
-      parameters
+      parameters,
+      renderer
     )
   }
 
-
   treated <- calculate_treated(
     variables,
     clinical_infections,
     events$recovery,
-    events$detection,
     parameters,
     timestep,
     renderer
@@ -98,12 +100,14 @@ simulate_infection <- function(
 #' @param variables a list of all of the model variables
 #' @param bitten_humans bitset of bitten humans
 #' @param parameters model parameters
+#' @param renderer model render object
 #' @param timestep current timestep
 #' @noRd
 calculate_infections <- function(
   variables,
   bitten_humans,
   parameters,
+  renderer,
   timestep
   ) {
   source_humans <- variables$state$get_index_of(
@@ -150,10 +154,22 @@ calculate_infections <- function(
     )
   }
 
-  bitset_at(
+  prob <- b * (1 - prophylaxis) * (1 - vaccine_efficacy)
+  infected <- bitset_at(source_humans, bernoulli_multi_p(prob))
+
+  incidence_renderer(
+    variables$birth,
+    renderer,
+    infected,
     source_humans,
-    bernoulli_multi_p(b * (1 - prophylaxis) * (1 - vaccine_efficacy))
+    prob,
+    'inc_',
+    parameters$incidence_rendering_min_ages,
+    parameters$incidence_rendering_max_ages,
+    timestep
   )
+
+  infected
 }
 
 #' @title Calculate clinical infections
@@ -162,12 +178,32 @@ calculate_infections <- function(
 #' @param variables a list of all of the model variables
 #' @param infections bitset of infected humans
 #' @param parameters model parameters
+#' @param renderer model render
+#' @param timestep current timestep
 #' @noRd
-calculate_clinical_infections <- function(variables, infections, parameters) {
+calculate_clinical_infections <- function(
+  variables,
+  infections,
+  parameters,
+  renderer,
+  timestep
+  ) {
   ica <- variables$ica$get_values(infections)
   icm <- variables$icm$get_values(infections)
   phi <- clinical_immunity(ica, icm, parameters)
-  bitset_at(infections, bernoulli_multi_p(phi))
+  clinical_infections <- bitset_at(infections, bernoulli_multi_p(phi))
+  incidence_renderer(
+    variables$birth,
+    renderer,
+    clinical_infections,
+    infections,
+    phi,
+    'inc_clinical_',
+    parameters$clinical_incidence_rendering_min_ages,
+    parameters$clinical_incidence_rendering_max_ages,
+    timestep
+  )
+  clinical_infections
 }
 
 #' @title Calculate severe infections
@@ -178,13 +214,15 @@ calculate_clinical_infections <- function(variables, infections, parameters) {
 #' @param variables a list of all of the model variables
 #' @param infections indices of all infected humans (for immunity boosting)
 #' @param parameters model parameters
+#' @param renderer model outputs
 #' @noRd
 update_severe_disease <- function(
   timestep,
   clinical_infections,
   variables,
   infections,
-  parameters
+  parameters,
+  renderer
   ) {
   if (clinical_infections$size() > 0) {
     age <- get_age(variables$birth$get_values(clinical_infections), timestep)
@@ -196,9 +234,22 @@ update_severe_disease <- function(
       parameters
     )
     develop_severe <- bernoulli_multi_p(theta)
+    severe_infections <- bitset_at(clinical_infections, develop_severe)
+    variables$is_severe$queue_update('yes', severe_infections)
     variables$is_severe$queue_update(
-      'yes',
-      bitset_at(clinical_infections, develop_severe)
+      'no',
+      severe_infections$copy()$not(TRUE)$and(infections)
+    ) # Prevent any unwanted severe reinfections
+    incidence_renderer(
+      variables$birth,
+      renderer,
+      severe_infections,
+      clinical_infections,
+      theta,
+      'inc_severe_',
+      parameters$severe_incidence_rendering_min_ages,
+      parameters$severe_incidence_rendering_max_ages,
+      timestep
     )
     boost_immunity(
       variables$iva,
@@ -224,7 +275,6 @@ calculate_treated <- function(
   variables,
   clinical_infections,
   recovery,
-  detection,
   parameters,
   timestep,
   renderer
@@ -266,7 +316,6 @@ calculate_treated <- function(
       timestep,
       treated_index
     )
-    detection$schedule(treated_index, 0)
   }
   treated_index
 }
@@ -288,24 +337,21 @@ schedule_infections <- function(
   parameters,
   asymptomatics
   ) {
-  included <- treated$not()
+  included <- treated$not(TRUE)
 
   to_infect <- clinical_infections$and(included)
-  to_infect_asym <- clinical_infections$not()$and(infections)$and(
+  to_infect_asym <- clinical_infections$copy()$not(TRUE)$and(infections)$and(
     included
   )
 
-  # change to symptomatic
   if(to_infect$size() > 0) {
     infection_times <- log_uniform(to_infect$size(), parameters$de)
     events$clinical_infection$schedule(to_infect, 0)
-    events$detection$schedule(to_infect, 0)
   }
 
   if(to_infect_asym$size() > 0) {
     infection_times <- log_uniform(to_infect_asym$size(), parameters$de)
     events$asymptomatic_infection$schedule(to_infect_asym, 0)
-    events$detection$schedule(to_infect_asym, 0)
   }
 }
 
diff --git a/R/mda_processes.R b/R/mda_processes.R
index 002df03c..3833a152 100644
--- a/R/mda_processes.R
+++ b/R/mda_processes.R
@@ -50,7 +50,7 @@ create_mda_listeners <- function(
       }
 
       # Move everyone else
-      other <- to_move$copy()$and(diseased$not())
+      other <- to_move$copy()$and(diseased$not(TRUE))
       if (other$size() > 0) {
         variables$state$queue_update('S', other)
       }
diff --git a/R/model.R b/R/model.R
index ce8c0d86..2b636991 100644
--- a/R/model.R
+++ b/R/model.R
@@ -8,14 +8,10 @@
 #'
 #'  * timestep: the timestep for the row
 #'  * infectivity: the infectivity from humans towards mosquitoes
-#'  * lambda: the effective biting rate on humans (per timestep) (per
-#'species). This defaults to lambda_All
-#'  * normal_lambda: the effective biting rate on adult humans (per
-#' timestep) (per species). This defaults to normal_lambda_All
 #'  * FOIM: the force of infection towards mosquitoes (per species)
 #'  * mu: the death rate of adult mosquitoes (per species)
-#'  * EIR: the Entomological Inoculation Rate (per timestep, over the whole
-#'  population)
+#'  * EIR: the Entomological Inoculation Rate (per timestep, per species, over 
+#'  the whole population)
 #'  * n_bitten: number of humans bitten by an infectious mosquito
 #'  * n_treated: number of humans treated for clinical or severe malaria this timestep
 #'  * n_infections: number of humans who get an asymptomatic, clinical or severe malaria this timestep
@@ -32,31 +28,45 @@
 #'  * n: number of humans between an inclusive age range at this timestep. This
 #' defaults to n_730_3650. Other age ranges can be set with
 #' prevalence_rendering_min_ages and prevalence_rendering_max_ages parameters.
-#' * n_detect: number of humans with an infection detectable by microscopy between an inclusive age range at this timestep. This
+#'  * n_detect: number of humans with an infection detectable by microscopy between an inclusive age range at this timestep. This
 #' defaults to n_detect_730_3650. Other age ranges can be set with
 #' prevalence_rendering_min_ages and prevalence_rendering_max_ages parameters.
-#'  * n_severe: number of humans with a severe infection detectable by microscopy 
-#'  between an inclusive age range at this timestep. Age ranges can be set with
+#'  * p_detect: the sum of probabilities of detection by microscopy between an
+#' inclusive age range at this timestep. This
+#' defaults to p_detect_730_3650. Other age ranges can be set with
+#' prevalence_rendering_min_ages and prevalence_rendering_max_ages parameters.
+#'  * n_severe: number of humans with a severe infection between an inclusive
+#' age range at this timestep. Age ranges can be set with
 #' severe_prevalence_rendering_min_ages and severe_prevalence_rendering_max_ages parameters.
 #'  * n_inc: number of new infections for humans between an inclusive age range at this timestep.
 #' incidence columns can be set with
 #' incidence_rendering_min_ages and incidence_rendering_max_ages parameters.
+#'  * p_inc: sum of probabilities of infection for humans between an inclusive age range at this timestep.
+#' incidence columns can be set with
+#' incidence_rendering_min_ages and incidence_rendering_max_ages parameters.
 #'  * n_inc_clinical: number of new clinical infections for humans between an inclusive age range at this timestep. 
 #' clinical incidence columns can be set with
 #' clinical_incidence_rendering_min_ages and clinical_incidence_rendering_max_ages parameters.
+#'  * p_inc_clinical: sub of probabilities of clinical infection for humans between an inclusive age range at this timestep. 
+#' clinical incidence columns can be set with
+#' clinical_incidence_rendering_min_ages and clinical_incidence_rendering_max_ages parameters.
 #'  * n_inc_severe: number of new severe infections for humans between an inclusive age range at this timestep.
 #' severe incidence columns can be set with
 #' severe_incidence_rendering_min_ages and severe_incidence_rendering_max_ages parameters.
+#'  * p_inc_severe: the sum of probabilities of severe infection for humans between an inclusive age range at this timestep.
+#' severe incidence columns can be set with
+#' severe_incidence_rendering_min_ages and severe_incidence_rendering_max_ages parameters.
 #'  * severe_deaths: number of deaths due to severe malaria. severe_enabled must be
 #' set to TRUE
-#'  * E_count: number of mosquitoes in the early larval stage
-#'  * L_count: number of mosquitoes in the late larval stage
-#'  * P_count: number of mosquitoes in the pupal stage
-#'  * Sm_count: number of adult female mosquitoes who are Susceptible
-#'  * Pm_count: number of adult female mosquitoes who are incubating
-#'  * Im_count: number of adult female mosquitoes who are infectious
-#'  * total_M: number of adult female mosquitoes. Variables with suffixes total_M_# refer to the number of adult 
-#' female mosquitoes from the different indicies of the set_species() vector
+#'  * E_count: number of mosquitoes in the early larval stage (per species)
+#'  * L_count: number of mosquitoes in the late larval stage (per species)
+#'  * P_count: number of mosquitoes in the pupal stage (per species)
+#'  * Sm_count: number of adult female mosquitoes who are Susceptible (per
+#'  species)
+#'  * Pm_count: number of adult female mosquitoes who are incubating (per
+#'  species)
+#'  * Im_count: number of adult female mosquitoes who are infectious (per
+#'  species)
 #'  * rate_D_A: rate that humans transition from clinical disease to
 #' asymptomatic
 #'  * rate_A_U: rate that humans transition from asymptomatic to
diff --git a/R/mortality_processes.R b/R/mortality_processes.R
index 7c07b911..2bcb1253 100644
--- a/R/mortality_processes.R
+++ b/R/mortality_processes.R
@@ -16,7 +16,7 @@ create_mortality_process <- function(variables, events, renderer, parameters) {
     ) / 365)
 
     died <- individual::Bitset$new(parameters$human_population)
-    died <- sample_bitset(died$not(), 1 / parameters$average_age)
+    died <- sample_bitset(died$not(TRUE), 1 / parameters$average_age)
 
     renderer$render('natural_deaths', died$size(), timestep)
 
@@ -72,7 +72,6 @@ create_mortality_process <- function(variables, events, renderer, parameters) {
         'recovery',
         'clinical_infection',
         'asymptomatic_infection',
-        'detection',
         'throw_away_net',
         'rtss_mass_doses',
         'rtss_mass_booster',
@@ -123,5 +122,5 @@ create_mortality_process <- function(variables, events, renderer, parameters) {
 died_from_severe <- function(severe, diseased, v, treated, fvt) {
   at_risk <- severe$copy()$and(diseased)
   treated_at_risk <- severe$copy()$and(treated)
-  sample_bitset(at_risk, v)$or(sample_bitset(treated_at_risk, fvt))
+  sample_bitset(at_risk$or(sample_bitset(treated_at_risk, fvt)), v)
 }
diff --git a/R/mosquito_biology.R b/R/mosquito_biology.R
index d30f3714..0b1076cc 100644
--- a/R/mosquito_biology.R
+++ b/R/mosquito_biology.R
@@ -148,7 +148,7 @@ peak_season_offset <- function(parameters) {
 #' @noRd
 death_rate <- function(f, W, Z, species, parameters) {
   mum <- parameters$mum[[species]]
-  p1_0 <- exp(-mum * parameters$foraging_time)
+  p1_0 <- exp(-mum * parameters$foraging_time[[species]])
   gonotrophic_cycle <- get_gonotrophic_cycle(species, parameters)
   p2 <- exp(-mum * gonotrophic_cycle)
   p1 <- p1_0 * W / (1 - Z * p1_0)
@@ -157,7 +157,7 @@ death_rate <- function(f, W, Z, species, parameters) {
 
 get_gonotrophic_cycle <- function(v, parameters) {
   f <- parameters$blood_meal_rates[[v]]
-  gonotrophic_cycle <- 1 / f - parameters$foraging_time
+  gonotrophic_cycle <- 1 / f - parameters$foraging_time[[v]]
 }
 
 #' @title Update the individual mosquito model after biting
diff --git a/R/parameters.R b/R/parameters.R
index d3ed8d6a..847e2b5f 100644
--- a/R/parameters.R
+++ b/R/parameters.R
@@ -346,8 +346,8 @@ get_parameters <- function(overrides = list()) {
     # treatment
     drug_efficacy          = numeric(0),
     drug_rel_c             = numeric(0),
-    drug_prophilaxis_shape = numeric(0),
-    drug_prophilaxis_scale = numeric(0),
+    drug_prophylaxis_shape = numeric(0),
+    drug_prophylaxis_scale = numeric(0),
     clinical_treatment_drugs     = list(),
     clinical_treatment_timesteps = list(),
     clinical_treatment_coverages = list(),
diff --git a/R/processes.R b/R/processes.R
index 6a7621a8..58b5a3b9 100644
--- a/R/processes.R
+++ b/R/processes.R
@@ -164,20 +164,16 @@ create_processes <- function(
       parameters,
       renderer
     ),
-    create_ode_rendering_process(renderer, solvers, parameters)
+    create_compartmental_rendering_process(renderer, solvers, parameters)
   )
 
   if (parameters$individual_mosquitoes) {
     processes <- c(
       processes,
-      individual::categorical_count_renderer_process(
-        renderer,
-        variables$mosquito_state,
-        c('Sm', 'Pm', 'Im')
-      ),
-      create_total_M_renderer_individual(
+      create_vector_count_renderer_individual(
         variables$mosquito_state,
         variables$species,
+        variables$mosquito_state,
         renderer,
         parameters
       )
@@ -187,7 +183,8 @@ create_processes <- function(
       processes,
       create_total_M_renderer_compartmental(
         renderer,
-        solvers
+        solvers,
+        parameters
       )
     )
   }
diff --git a/R/render.R b/R/render.R
index 4fd2cb15..90755891 100644
--- a/R/render.R
+++ b/R/render.R
@@ -24,19 +24,15 @@ create_prevelance_renderer <- function(
   renderer
   ) {
   function(timestep) {
-    age <- get_age(birth$get_values(), timestep)
     asymptomatic <- state$get_index_of('A')
-    asymptomatic_vector <- asymptomatic$to_vector()
-    asymptomatic_detected <- bernoulli_multi_p(
-      probability_of_detection(
-        age[asymptomatic_vector],
-        immunity$get_values(asymptomatic),
-        parameters
-      )
-    )
-    detected <- state$get_index_of(c('Tr', 'D'))$or(
-      bitset_at(asymptomatic, asymptomatic_detected)
+    prob <- probability_of_detection(
+      get_age(birth$get_values(asymptomatic), timestep),
+      immunity$get_values(asymptomatic),
+      parameters
     )
+    asymptomatic_detected <- bitset_at(asymptomatic, bernoulli_multi_p(prob))
+    clinically_detected <- state$get_index_of(c('Tr', 'D'))
+    detected <- clinically_detected$copy()$or(asymptomatic_detected)
 
     severe <- is_severe$get_index_of('yes')$and(detected)
     for (i in seq_along(parameters$prevalence_rendering_min_ages)) {
@@ -50,7 +46,14 @@ create_prevelance_renderer <- function(
       ) 
       renderer$render(
         paste0('n_detect_', lower, '_', upper),
-        in_age$and(detected)$size(),
+        in_age$copy()$and(detected)$size(),
+        timestep
+      )
+      renderer$render(
+        paste0('p_detect_', lower, '_', upper),
+        in_age$copy()$and(clinically_detected)$size() + sum(
+          prob[bitset_index(asymptomatic, in_age)]
+        ),
         timestep
       )
     }
@@ -62,7 +65,7 @@ create_prevelance_renderer <- function(
         paste0('n_', lower, '_', upper),
         in_age$size(),
         timestep
-      ) 
+      )
       renderer$render(
         paste0('n_severe_', lower, '_', upper),
         in_age$and(severe)$size(),
@@ -75,64 +78,45 @@ create_prevelance_renderer <- function(
 #' @title Render incidence statistics
 #' 
 #' @description renders incidence (new for this timestep) for indivduals
-#' detected by microscopy and with severe malaria
 #' 
 #' @param birth variable for birth of the individual
-#' @param is_severe variable for if individual has severe malaria
-#' @param parameters model parameters
-#' @param renderer model renderer
+#' @param renderer object for model outputs
+#' @param target incidence population
+#' @param source_pop the population which is sampled for infection
+#' @param prob probability of infection
+#' @param prefix for model outputs
+#' @param lowers age bounds
+#' @param uppers age bounds
+#' @param timestep current target
 #' 
 #' @noRd
-create_incidence_renderer <- function(birth, is_severe, parameters, renderer) {
-  function(timestep, target) {
-    severe <- is_severe$get_index_of('yes')$and(target)
-    for (i in seq_along(parameters$incidence_rendering_min_ages)) {
-      lower <- parameters$incidence_rendering_min_ages[[i]]
-      upper <- parameters$incidence_rendering_max_ages[[i]]
-      in_age <- in_age_range(birth, timestep, lower, upper)
-      renderer$render(paste0('n_', lower, '_', upper), in_age$size(), timestep)
-      renderer$render(
-        paste0('n_inc_', lower, '_', upper),
-        in_age$and(target)$size(),
-        timestep
-      )
-    }
-    for (i in seq_along(parameters$severe_incidence_rendering_min_ages)) {
-      lower <- parameters$severe_incidence_rendering_min_ages[[i]]
-      upper <- parameters$severe_incidence_rendering_max_ages[[i]]
-      in_age <- in_age_range(birth, timestep, lower, upper)
-      renderer$render(paste0('n_', lower, '_', upper), in_age$size(), timestep)
-      renderer$render(
-        paste0('n_inc_severe_', lower, '_', upper),
-        in_age$and(target)$and(severe)$size(),
-        timestep
-      )
-    }
-  }
-}
+incidence_renderer <- function(
+  birth,
+  renderer,
+  target,
+  source_pop,
+  prob,
+  prefix,
+  lowers,
+  uppers,
+  timestep
+  ) {
+  for (i in seq_along(lowers)) {
+    lower <- lowers[[i]]
+    upper <- uppers[[i]]
+    in_age <- in_age_range(birth, timestep, lower, upper)
+    renderer$render(paste0('n_', lower, '_', upper), in_age$size(), timestep)
+    renderer$render(
+      paste0('n_', prefix, lower, '_', upper),
+      in_age$copy()$and(target)$size(),
+      timestep
+    )
 
-#' @title Render incidence statistics
-#' 
-#' @description renders clinical incidence (new for this timestep)
-#' 
-#' @param birth variable for birth of the individual
-#' @param parameters model parameters
-#' @param renderer model renderer
-#' 
-#' @noRd
-create_clinical_incidence_renderer <- function(birth, parameters, renderer) {
-  function(timestep, target) {
-    for (i in seq_along(parameters$clinical_incidence_rendering_min_ages)) {
-      lower <- parameters$clinical_incidence_rendering_min_ages[[i]]
-      upper <- parameters$clinical_incidence_rendering_max_ages[[i]]
-      in_age <- in_age_range(birth, timestep, lower, upper)
-      renderer$render(paste0('n_', lower, '_', upper), in_age$size(), timestep)
-      renderer$render(
-        paste0('n_inc_clinical_', lower, '_', upper),
-        in_age$and(target)$size(),
-        timestep
-      )
-    }
+    renderer$render(
+      paste0('p_', prefix, lower, '_', upper),
+      sum(prob[bitset_index(source_pop, in_age)]),
+      timestep
+    )
   }
 }
 
@@ -152,36 +136,37 @@ create_variable_mean_renderer_process <- function(
   }
 }
 
-create_total_M_renderer_individual <- function(
+create_vector_count_renderer_individual <- function(
   mosquito_state,
   species,
+  state,
   renderer,
   parameters
   ) {
   function(timestep) {
-    adult <- mosquito_state$get_index_of('NonExistent')$not()
+    adult <- mosquito_state$get_index_of('NonExistent')$not(TRUE)
     for (i in seq_along(parameters$species)) {
-      renderer$render(
-        paste0('total_M_', i),
-        species$get_index_of(parameters$species[[i]])$and(adult)$size(),
-        timestep
-      )
+      species_name <- parameters$species[[i]]
+      species_index <- species$get_index_of(species_name)
+      for (s in state$get_categories()) {
+        renderer$render(
+          paste0(s, '_', species_name, '_count'),
+          state$get_index_of(s)$and(species_index)$size(),
+          timestep
+        )
+      }
     }
-
-    renderer$render('total_M', adult$size(), timestep)
   }
 }
 
-create_total_M_renderer_compartmental <- function(renderer, solvers) {
+create_total_M_renderer_compartmental <- function(renderer, solvers, parameters) {
   function(timestep) {
     total_M <- 0
     for (i in seq_along(solvers)) {
       row <- solver_get_states(solvers[[i]])
       species_M <- sum(row[ADULT_ODE_INDICES])
       total_M <- total_M + species_M
-      renderer$render(paste0('total_M_', i), species_M, timestep)
+      renderer$render(paste0('total_M_', parameters$species[[i]]), species_M, timestep)
     }
-
-    renderer$render('total_M', total_M, timestep)
   }
 }
diff --git a/R/rtss.R b/R/rtss.R
index 638deb7a..6d603777 100644
--- a/R/rtss.R
+++ b/R/rtss.R
@@ -28,7 +28,7 @@ create_rtss_epi_process <- function(
       not_recently_vaccinated <- variables$rtss_vaccinated$get_index_of(
         a = max(timestep - parameters$rtss_epi_min_wait, 0),
         b = timestep
-      )$not()
+      )$not(TRUE)
       target <- to_vaccinate$and(not_recently_vaccinated)$to_vector()
     }
 
@@ -78,7 +78,7 @@ create_rtss_mass_listener <- function(
       not_recently_vaccinated <- variables$rtss_vaccinated$get_index_of(
         a = max(timestep - parameters$rtss_mass_min_wait, 0),
         b = timestep
-      )$not()
+      )$not(TRUE)
       target <- in_age_group$and(not_recently_vaccinated)$to_vector()
     }
     
diff --git a/R/utils.R b/R/utils.R
index f9d046f9..b946c7ec 100644
--- a/R/utils.R
+++ b/R/utils.R
@@ -15,6 +15,15 @@ bitset_at <- function(b, i) {
 
 bernoulli_multi_p <- function(p) bernoulli_multi_p_cpp(p)
 
+
+#' @title find the indices of a where it intersects with b
+#' @description synonymous with \code{which(a$to_vector() %in%
+#' b$to_vector())} but faster
+#' @param a the bitset to index
+#' @param b the bitset to check
+#' @noRd
+bitset_index <- function(a, b) bitset_index_cpp(a$.bitset, b$.bitset)
+
 #' @importFrom stats runif
 log_uniform <- function(size, rate) -rate * log(runif(size))
 
diff --git a/R/vector_control.R b/R/vector_control.R
index dc39dd1c..f197a7d1 100644
--- a/R/vector_control.R
+++ b/R/vector_control.R
@@ -39,8 +39,7 @@ prob_bitten <- function(
 
   if (parameters$spraying) {
     phi_indoors <- parameters$phi_indoors[[species]]
-    unprotected <- variables$spray_time$get_index_of(set=-1)
-    protected <- unprotected$not()
+    protected <- variables$spray_time$get_index_of(set=-1)$not(TRUE)
     spray_time <- variables$spray_time$get_values(protected)
     matches <- match(spray_time, parameters$spraying_timesteps)
     ls_theta <- parameters$spraying_ls_theta[matches, species]
diff --git a/R/vector_parameters.R b/R/vector_parameters.R
index 888477e1..8ff25c54 100644
--- a/R/vector_parameters.R
+++ b/R/vector_parameters.R
@@ -2,6 +2,7 @@
 #' @details Default parameters:
 #' species: "gamb"
 #' blood_meal_rates: 0.3333333
+#' foraging_time: .69
 #' Q0: 0.92
 #' phi_bednets: 0.85
 #' phi_indoors: 0.90
@@ -13,6 +14,7 @@
 gamb_params <- list(
   species = 'gamb',
   blood_meal_rates = 1/3,
+  foraging_time = .69,
   Q0 = .92,
   phi_bednets = .85,
   phi_indoors = .90,
@@ -23,6 +25,7 @@ gamb_params <- list(
 #' @details Default parameters:
 #' species: "arab"
 #' blood_meal_rates: 0.3333333
+#' foraging_time: .69
 #' Q0: 0.71
 #' phi_bednets: 0.8
 #' phi_indoors: 0.86
@@ -34,6 +37,7 @@ gamb_params <- list(
 arab_params <- list(
   species = 'arab',
   blood_meal_rates = 1/3,
+  foraging_time = .69,
   Q0 = .71,
   phi_bednets = .8,
   phi_indoors = .86,
@@ -44,6 +48,7 @@ arab_params <- list(
 #' @details Default parameters:
 #' species: "fun"
 #' blood_meal_rates: 0.3333333
+#' foraging_time: .69
 #' Q0: 0.94
 #' phi_bednets: 0.78
 #' phi_indoors: 0.87
@@ -55,6 +60,7 @@ arab_params <- list(
 fun_params <- list(
   species = 'fun',
   blood_meal_rates = 1/3,
+  foraging_time = .69,
   Q0 = .94,
   phi_bednets = .78,
   phi_indoors = .87,
@@ -77,6 +83,7 @@ set_species <- function(parameters, species, proportions) {
   keys <- c(
     'species',
     'blood_meal_rates',
+    'foraging_time',
     'Q0',
     'phi_bednets',
     'phi_indoors',
@@ -86,6 +93,15 @@ set_species <- function(parameters, species, proportions) {
     parameters[[key]] <- rep(NA, length(species))
   }
   for (v in seq_along(species)) {
+    if (species[[v]]$foraging_time > 1 / species[[v]]$blood_meal_rates) {
+      stop(
+        sprintf(
+          "blood meal time (%f) must be >= foraging time (%f)",
+          1 / species[[v]]$blood_meal_rates,
+          species[[v]]$foraging_time
+        )
+      )
+    }
     for (key in keys) {
       parameters[[key]][[v]] <- species[[v]][[key]]
     }
diff --git a/man/arab_params.Rd b/man/arab_params.Rd
index b5655b25..144bbf70 100644
--- a/man/arab_params.Rd
+++ b/man/arab_params.Rd
@@ -5,7 +5,7 @@
 \alias{arab_params}
 \title{Preset parameters for the An. arabiensis vector}
 \format{
-An object of class \code{list} of length 6.
+An object of class \code{list} of length 7.
 }
 \usage{
 arab_params
@@ -17,6 +17,7 @@ Preset parameters for the An. arabiensis vector
 Default parameters:
 species: "arab"
 blood_meal_rates: 0.3333333
+foraging_time: .69
 Q0: 0.71
 phi_bednets: 0.8
 phi_indoors: 0.86
diff --git a/man/fun_params.Rd b/man/fun_params.Rd
index b5d0eaff..6063bfbc 100644
--- a/man/fun_params.Rd
+++ b/man/fun_params.Rd
@@ -5,7 +5,7 @@
 \alias{fun_params}
 \title{Preset parameters for the An. funestus vector}
 \format{
-An object of class \code{list} of length 6.
+An object of class \code{list} of length 7.
 }
 \usage{
 fun_params
@@ -17,6 +17,7 @@ Preset parameters for the An. funestus vector
 Default parameters:
 species: "fun"
 blood_meal_rates: 0.3333333
+foraging_time: .69
 Q0: 0.94
 phi_bednets: 0.78
 phi_indoors: 0.87
diff --git a/man/gamb_params.Rd b/man/gamb_params.Rd
index 434c46fa..fef4b29f 100644
--- a/man/gamb_params.Rd
+++ b/man/gamb_params.Rd
@@ -5,7 +5,7 @@
 \alias{gamb_params}
 \title{Preset parameters for the An. gambiae s.s vector}
 \format{
-An object of class \code{list} of length 6.
+An object of class \code{list} of length 7.
 }
 \usage{
 gamb_params
@@ -17,6 +17,7 @@ Preset parameters for the An. gambiae s.s vector
 Default parameters:
 species: "gamb"
 blood_meal_rates: 0.3333333
+foraging_time: .69
 Q0: 0.92
 phi_bednets: 0.85
 phi_indoors: 0.90
diff --git a/man/get_parameters.Rd b/man/get_parameters.Rd
index 0f6f8d11..5fea735c 100644
--- a/man/get_parameters.Rd
+++ b/man/get_parameters.Rd
@@ -233,8 +233,8 @@ outputs; default = 730
 \item incidence_rendering_min_ages - the minimum ages for incidence
 outputs (includes asymptomatic microscopy +); default = turned off
 \item incidence_rendering_max_ages - the corresponding max ages; default = turned off
-clinical_incidence_rendering_min_ages - the minimum ages for clinical incidence outputs (symptomatic); default = 0
-clinical_incidence_rendering_max_ages - the corresponding max ages; default = 1825
+\item clinical_incidence_rendering_min_ages - the minimum ages for clinical incidence outputs (symptomatic); default = 0
+\item clinical_incidence_rendering_max_ages - the corresponding max ages; default = 1825
 \item severe_prevalence_rendering_min_ages - the minimum ages for severe
 prevalence outputs; default = turned off
 \item severe_prevalence_rendering_max_ages - the corresponding max ages; default = turned off
@@ -253,7 +253,6 @@ individually or compartmentally; default = TRUE
 \item r_tol - the relative tolerance for the ode solver; default = 1e-4
 \item a_tol - the absolute tolerance for the ode solver; default = 1e-4
 \item ode_max_steps - the max number of steps for the solver; default = 1e6
-\item ode_debug - whether to print debug output for ode parameter updates; default = FALSE
 \item enable_heterogeneity - boolean whether to include heterogeneity in biting
 rates; default = TRUE
 }}
diff --git a/man/run_simulation.Rd b/man/run_simulation.Rd
index d25a6526..5dd1c8ec 100644
--- a/man/run_simulation.Rd
+++ b/man/run_simulation.Rd
@@ -25,14 +25,10 @@ The resulting dataframe contains the following columns:
 \itemize{
 \item timestep: the timestep for the row
 \item infectivity: the infectivity from humans towards mosquitoes
-\item lambda: the effective biting rate on humans (per timestep) (per
-species). This defaults to lambda_All
-\item normal_lambda: the effective biting rate on adult humans (per
-timestep) (per species). This defaults to normal_lambda_All
 \item FOIM: the force of infection towards mosquitoes (per species)
 \item mu: the death rate of adult mosquitoes (per species)
-\item EIR: the Entomological Inoculation Rate (per timestep, over the whole
-population)
+\item EIR: the Entomological Inoculation Rate (per timestep, per species, over
+the whole population)
 \item n_bitten: number of humans bitten by an infectious mosquito
 \item n_treated: number of humans treated for clinical or severe malaria this timestep
 \item n_infections: number of humans who get an asymptomatic, clinical or severe malaria this timestep
@@ -46,32 +42,48 @@ population)
 \item icm_mean: the mean maternal immunity to clinical infection over the population of humans
 \item ib_mean: the mean blood immunity to all infection over the population of humans
 \item id_mean: the mean immunity from detection through microscopy over the population of humans
-\item pv: prevalence for humans between an inclusive age range (in timesteps). This
-defaults to pv_730_3650. Other prevalence columns can be set with
+\item n: number of humans between an inclusive age range at this timestep. This
+defaults to n_730_3650. Other age ranges can be set with
 prevalence_rendering_min_ages and prevalence_rendering_max_ages parameters.
-\item pv_severe: prevalence for severe malaria in humans between an inclusive age range (in timesteps).
-These columns can be set with the
+\item n_detect: number of humans with an infection detectable by microscopy between an inclusive age range at this timestep. This
+defaults to n_detect_730_3650. Other age ranges can be set with
+prevalence_rendering_min_ages and prevalence_rendering_max_ages parameters.
+\item p_detect: the sum of probabilities of detection by microscopy between an
+inclusive age range at this timestep. This
+defaults to p_detect_730_3650. Other age ranges can be set with
+prevalence_rendering_min_ages and prevalence_rendering_max_ages parameters.
+\item n_severe: number of humans with a severe infection between an inclusive
+age range at this timestep. Age ranges can be set with
 severe_prevalence_rendering_min_ages and severe_prevalence_rendering_max_ages parameters.
-\item inc: incidence for humans between an inclusive age range (in timesteps).
+\item n_inc: number of new infections for humans between an inclusive age range at this timestep.
+incidence columns can be set with
+incidence_rendering_min_ages and incidence_rendering_max_ages parameters.
+\item p_inc: sum of probabilities of infection for humans between an inclusive age range at this timestep.
 incidence columns can be set with
 incidence_rendering_min_ages and incidence_rendering_max_ages parameters.
-\item clin_inc: clinical incidence for humans between an inclusive age range (in timesteps)
+\item n_inc_clinical: number of new clinical infections for humans between an inclusive age range at this timestep.
 clinical incidence columns can be set with
 clinical_incidence_rendering_min_ages and clinical_incidence_rendering_max_ages parameters.
-\item inc_severe: severe incidence for humans between an inclusive age range (in
-timesteps).
+\item p_inc_clinical: sub of probabilities of clinical infection for humans between an inclusive age range at this timestep.
+clinical incidence columns can be set with
+clinical_incidence_rendering_min_ages and clinical_incidence_rendering_max_ages parameters.
+\item n_inc_severe: number of new severe infections for humans between an inclusive age range at this timestep.
+severe incidence columns can be set with
+severe_incidence_rendering_min_ages and severe_incidence_rendering_max_ages parameters.
+\item p_inc_severe: the sum of probabilities of severe infection for humans between an inclusive age range at this timestep.
 severe incidence columns can be set with
 severe_incidence_rendering_min_ages and severe_incidence_rendering_max_ages parameters.
 \item severe_deaths: number of deaths due to severe malaria. severe_enabled must be
 set to TRUE
-\item E_count: number of mosquitoes in the early larval stage
-\item L_count: number of mosquitoes in the late larval stage
-\item P_count: number of mosquitoes in the pupal stage
-\item Sm_count: number of adult female mosquitoes who are Susceptible
-\item Pm_count: number of adult female mosquitoes who are incubating
-\item Im_count: number of adult female mosquitoes who are infectious
-\item total_M: number of adult female mosquitoes. Variables with suffixes total_M_# refer to the number of adult
-female mosquitoes from the different indicies of the set_species() vector
+\item E_count: number of mosquitoes in the early larval stage (per species)
+\item L_count: number of mosquitoes in the late larval stage (per species)
+\item P_count: number of mosquitoes in the pupal stage (per species)
+\item Sm_count: number of adult female mosquitoes who are Susceptible (per
+species)
+\item Pm_count: number of adult female mosquitoes who are incubating (per
+species)
+\item Im_count: number of adult female mosquitoes who are infectious (per
+species)
 \item rate_D_A: rate that humans transition from clinical disease to
 asymptomatic
 \item rate_A_U: rate that humans transition from asymptomatic to
diff --git a/man/set_mass_rtss.Rd b/man/set_mass_rtss.Rd
index 8622e912..aa1a479d 100644
--- a/man/set_mass_rtss.Rd
+++ b/man/set_mass_rtss.Rd
@@ -26,7 +26,9 @@ set_mass_rtss(
 
 \item{max_ages}{for the target population, inclusive (in timesteps)}
 
-\item{min_wait}{the minimum acceptable time since the last vaccination (in timesteps);}
+\item{min_wait}{the minimum acceptable time since the last vaccination (in timesteps);
+When using both set_mass_rtss and set_rtss_epi, this represents the minimum
+time between an individual being vaccinated under one scheme and vaccinated under another.}
 
 \item{boosters}{the timesteps (following the initial vaccination) at which booster vaccinations are administered}
 
diff --git a/man/set_rtss_epi.Rd b/man/set_rtss_epi.Rd
index 77b82309..3d759ef0 100644
--- a/man/set_rtss_epi.Rd
+++ b/man/set_rtss_epi.Rd
@@ -28,8 +28,10 @@ set_rtss_epi(
 \item{age}{for the target population, (in timesteps)}
 
 \item{min_wait}{the minimum acceptable time since the last vaccination (in
-timesteps); This includes seasonal boosters and vaccinations between EPI and mass
-strategies}
+timesteps); When seasonal_boosters = TRUE, this represents the minimum time
+between an individual receiving the final dose and the first booster. When using
+both set_mass_rtss and set_rtss_epi, this represents the minimum time between
+an individual being vaccinated under one scheme and vaccinated under another.}
 
 \item{boosters}{the timesteps (following the final dose) at which booster vaccinations are administered}
 
diff --git a/src/RcppExports.cpp b/src/RcppExports.cpp
index 0d778d51..a73a67bf 100644
--- a/src/RcppExports.cpp
+++ b/src/RcppExports.cpp
@@ -232,6 +232,18 @@ BEGIN_RCPP
     return rcpp_result_gen;
 END_RCPP
 }
+// bitset_index_cpp
+std::vector<size_t> bitset_index_cpp(Rcpp::XPtr<individual_index_t> a, Rcpp::XPtr<individual_index_t> b);
+RcppExport SEXP _malariasimulation_bitset_index_cpp(SEXP aSEXP, SEXP bSEXP) {
+BEGIN_RCPP
+    Rcpp::RObject rcpp_result_gen;
+    Rcpp::RNGScope rcpp_rngScope_gen;
+    Rcpp::traits::input_parameter< Rcpp::XPtr<individual_index_t> >::type a(aSEXP);
+    Rcpp::traits::input_parameter< Rcpp::XPtr<individual_index_t> >::type b(bSEXP);
+    rcpp_result_gen = Rcpp::wrap(bitset_index_cpp(a, b));
+    return rcpp_result_gen;
+END_RCPP
+}
 // fast_weighted_sample
 Rcpp::IntegerVector fast_weighted_sample(size_t size, std::vector<double> probs);
 RcppExport SEXP _malariasimulation_fast_weighted_sample(SEXP sizeSEXP, SEXP probsSEXP) {
@@ -263,6 +275,7 @@ static const R_CallMethodDef CallEntries[] = {
     {"_malariasimulation_solver_step", (DL_FUNC) &_malariasimulation_solver_step, 1},
     {"_malariasimulation_random_seed", (DL_FUNC) &_malariasimulation_random_seed, 1},
     {"_malariasimulation_bernoulli_multi_p_cpp", (DL_FUNC) &_malariasimulation_bernoulli_multi_p_cpp, 1},
+    {"_malariasimulation_bitset_index_cpp", (DL_FUNC) &_malariasimulation_bitset_index_cpp, 2},
     {"_malariasimulation_fast_weighted_sample", (DL_FUNC) &_malariasimulation_fast_weighted_sample, 2},
     {"run_testthat_tests", (DL_FUNC) &run_testthat_tests, 0},
     {NULL, NULL, 0}
diff --git a/src/test-helper.cpp b/src/test-helper.cpp
deleted file mode 100644
index 898ae2bb..00000000
--- a/src/test-helper.cpp
+++ /dev/null
@@ -1,12 +0,0 @@
-/*
- * test-helper.cpp
- *
- *  Created on: 20 Aug 2020
- *      Author: gc1610
- */
-
-#include "test-helper.h"
-
-individual_index_t individual_index(size_t size, std::vector<size_t> i) {
-    return individual_index_t(size, i.cbegin(), i.cend());
-}
diff --git a/src/test-helper.h b/src/test-helper.h
deleted file mode 100644
index b1065d8a..00000000
--- a/src/test-helper.h
+++ /dev/null
@@ -1,15 +0,0 @@
-/*
- * test-helper.h
- *
- *  Created on: 20 Aug 2020
- *      Author: gc1610
- */
-
-#ifndef SRC_TEST_HELPER_H_
-#define SRC_TEST_HELPER_H_
-
-#include <individual.h>
-
-individual_index_t individual_index(size_t size, std::vector<size_t> i);
-
-#endif /* SRC_TEST_HELPER_H_ */
diff --git a/src/test-mock.h b/src/test-mock.h
deleted file mode 100644
index dcc5ffad..00000000
--- a/src/test-mock.h
+++ /dev/null
@@ -1,38 +0,0 @@
-// -*- mode: C++; c-indent-level: 4; c-basic-offset: 4; indent-tabs-mode: nil; -*-
-/*
- * test-emergence.h
- *
- *  Created on: 5 Aug 2020
- *      Author: gc1610
- */
-
-#ifndef SRC_TEST_MOCK_H_
-#define SRC_TEST_MOCK_H_
-
-#include "solver.h"
-#include <testthat/trompeloeil.hpp>
-#include <individual.h>
-#include "Random.h"
-
-struct MockCategory : public CategoricalVariable {
-    MockCategory(
-        const std::vector<std::string> categories,
-        const std::vector<std::string> values)
-        : CategoricalVariable(categories, values) {}
-    MAKE_MOCK2(queue_update, void(const std::string, const individual_index_t&), override);
-};
-
-class MockRandom : public RandomInterface {
-public:
-    MAKE_MOCK2(bernoulli, std::vector<size_t>(size_t, double), override);
-    MAKE_MOCK3(sample, std::vector<size_t>(size_t, size_t, bool), override);
-};
-
-
-integration_function_t mock_integration = [](
-    const state_t&,
-    state_t&,
-    double t
-) {};
-
-#endif /* SRC_TEST_MOCK_H_ */
diff --git a/src/utils.cpp b/src/utils.cpp
index 585a75b7..838a9530 100644
--- a/src/utils.cpp
+++ b/src/utils.cpp
@@ -1,6 +1,7 @@
 
 #include <Rcpp.h>
 #include "Random.h"
+#include <individual.h>
 
 //[[Rcpp::export]]
 void random_seed(size_t seed) {
@@ -16,6 +17,22 @@ std::vector<size_t> bernoulli_multi_p_cpp(const std::vector<double> p) {
     return values;
 }
 
+//[[Rcpp::export]]
+std::vector<size_t> bitset_index_cpp(
+    Rcpp::XPtr<individual_index_t> a,
+    Rcpp::XPtr<individual_index_t> b
+    ) {
+    auto values = std::vector<size_t>();
+    auto i = 1u;
+    for (const auto& v : *a) {
+        if (b->find(v) != b->cend()) {
+            values.push_back(i);
+        }
+        ++i;
+    }
+    return values;
+}
+
 //[[Rcpp::export(rng = false)]]
 Rcpp::IntegerVector fast_weighted_sample(
     size_t size,
diff --git a/tests/testthat/helper-integration.R b/tests/testthat/helper-integration.R
index 8a3d6f4e..b69b7dcb 100644
--- a/tests/testthat/helper-integration.R
+++ b/tests/testthat/helper-integration.R
@@ -23,7 +23,23 @@ mock_method <- function(class, method, mock) {
     paste0('Mock', class$classname),
     inherit = class
   )
-  MockClass$set('public', method, function(...) mock(...))
+  MockClass$set(
+    'public',
+    method,
+    function(...) {
+      args <- lapply(
+        list(...),
+        function(arg) {
+          # copy a any bitsets because these will be lazily evaluated
+          if (inherits(arg, 'Bitset')) {
+            return(arg$copy())
+          }
+          arg
+        }
+      )
+      do.call(mock, args)
+    }
+  )
   MockClass$set('public', paste0(method, '_mock'), function() mock)
   MockClass
 }
diff --git a/tests/testthat/test-biting.R b/tests/testthat/test-biting.R
new file mode 100644
index 00000000..e6c19fb9
--- /dev/null
+++ b/tests/testthat/test-biting.R
@@ -0,0 +1,12 @@
+test_that('compartmental always gives positive infectious', {
+  solver_states <- rep(0, length(ODE_INDICES) + length(ADULT_ODE_INDICES))
+  solver_states[[ADULT_ODE_INDICES['Im']]] <- -1e-10
+  expect_gte(calculate_infectious_compartmental(solver_states), 0)
+})
+
+test_that('gonotrophic_cycle cannot be negative', {
+  params <- get_parameters()
+  vparams <- gamb_params
+  vparams$blood_meal_rates <- 5
+  expect_error(set_species(params, list(vparams), 1))
+})
diff --git a/tests/testthat/test-infection-integration.R b/tests/testthat/test-infection-integration.R
index f55e484c..2dbd2ad5 100644
--- a/tests/testthat/test-infection-integration.R
+++ b/tests/testthat/test-infection-integration.R
@@ -38,6 +38,8 @@ test_that('simulate_infection integrates different types of infection and schedu
   mockery::stub(simulate_infection, 'update_severe_disease', severe_infection_mock)
   mockery::stub(simulate_infection, 'calculate_treated', treated_mock)
   mockery::stub(simulate_infection, 'schedule_infections', schedule_mock)
+  mockery::stub(simulate_infection, 'incidence_renderer', mockery::mock())
+  mockery::stub(simulate_infection, 'clinical_incidence_renderer', mockery::mock())
   simulate_infection(
     variables,
     events,
@@ -64,6 +66,7 @@ test_that('simulate_infection integrates different types of infection and schedu
     variables,
     bitten,
     parameters,
+    renderer,
     timestep
   )
 
@@ -72,7 +75,9 @@ test_that('simulate_infection integrates different types of infection and schedu
     1,
     variables,
     infected,
-    parameters
+    parameters,
+    renderer,
+    timestep
   )
 
   mockery::expect_args(
@@ -82,7 +87,8 @@ test_that('simulate_infection integrates different types of infection and schedu
     clinical,
     variables,
     infected,
-    parameters
+    parameters,
+    renderer
   )
 
   mockery::expect_args(
@@ -91,7 +97,6 @@ test_that('simulate_infection integrates different types of infection and schedu
     variables,
     clinical,
     events$recovery,
-    events$detection,
     parameters,
     timestep,
     renderer
@@ -147,6 +152,7 @@ test_that('calculate_infections works various combinations of drug and vaccinati
     variables,
     bitten_humans, 
     parameters,
+    mock_render(timestep),
     timestep
   )
 
@@ -245,9 +251,7 @@ test_that('calculate_treated correctly samples treated and updates the drug stat
   )
 
   recovery_mock <- mockery::mock()
-  detection_mock <- mockery::mock()
   mockery::stub(calculate_treated, 'recovery$schedule', recovery_mock)
-  mockery::stub(calculate_treated, 'detection$schedule', detection_mock)
 
   seek_treatment <- individual::Bitset$new(4)$insert(c(1, 2, 4))
   mockery::stub(
@@ -267,7 +271,6 @@ test_that('calculate_treated correctly samples treated and updates the drug stat
     variables,
     clinical_infections,
     events$recovery,
-    events$detection,
     parameters,
     timestep,
     mock_render(timestep)
@@ -296,12 +299,6 @@ test_that('calculate_treated correctly samples treated and updates the drug stat
   )
   expect_bitset_update(variables$drug$queue_update, c(2, 1), c(1, 4))
   expect_bitset_update(variables$drug_time$queue_update, 5, c(1, 4))
-
-  expect_bitset_schedule(
-    detection_mock,
-    c(1, 4),
-    0
-  )
 })
 
 test_that('schedule_infections correctly schedules new infections', {
@@ -309,10 +306,8 @@ test_that('schedule_infections correctly schedules new infections', {
 
   clinical_mock <- mockery::mock()
   asym_mock <- mockery::mock()
-  detection_mock <- mockery::mock()
 
   events <- list(
-    detection = list(schedule = detection_mock),
     clinical_infection = list(
       schedule = clinical_mock
     ),
@@ -341,24 +336,11 @@ test_that('schedule_infections correctly schedules new infections', {
     0
   )
 
-  expect_bitset_schedule(
-    detection_mock,
-    c(5, 6, 13, 14, 15),
-    0
-  )
-
   expect_bitset_schedule(
     asym_mock,
     c(1, 2, 3, 4, 16, 17, 18, 19, 20),
     0,
   )
-
-  expect_bitset_schedule(
-    detection_mock,
-    c(1, 2, 3, 4, 16, 17, 18, 19, 20),
-    0,
-    call = 2
-  )
 })
 
 test_that('prophylaxis is considered for medicated humans', {
@@ -383,7 +365,13 @@ test_that('prophylaxis is considered for medicated humans', {
   m <- mockery::mock(seq(3))
   mockery::stub(calculate_infections, 'bernoulli_multi_p', m)
 
-  calculate_infections(variables, bitten, parameters, timestep)
+  calculate_infections(
+    variables,
+    bitten,
+    parameters,
+    mock_render(timestep),
+    timestep
+  )
 
   expect_equal(
     mockery::mock_args(m)[[1]][[1]],
diff --git a/tests/testthat/test-mda.R b/tests/testthat/test-mda.R
index 689cad3c..d9d65ada 100644
--- a/tests/testthat/test-mda.R
+++ b/tests/testthat/test-mda.R
@@ -50,7 +50,6 @@ test_that('MDA moves the diseased and non-diseased population correctly', {
     mockery::mock_args(mock_correlation)[[1]][[1]],
     c(3, 4)
   )
-
   expect_bitset_update(
     variables$state$queue_update_mock(),
     'Tr',
diff --git a/tests/testthat/test-render.R b/tests/testthat/test-render.R
index ea980230..a1e68f43 100644
--- a/tests/testthat/test-render.R
+++ b/tests/testthat/test-render.R
@@ -14,7 +14,6 @@ test_that('that default rendering works', {
     c('no', 'no', 'yes', 'no')
   )
 
-
   renderer <- mock_render(1)
   process <- create_prevelance_renderer(
     state,
@@ -25,6 +24,7 @@ test_that('that default rendering works', {
     renderer
   )
 
+  mockery::stub(process, 'probability_of_detection', mockery::mock(.5))
   mockery::stub(process, 'bernoulli_multi_p', mockery::mock(1))
   process(timestep)
 
@@ -44,6 +44,14 @@ test_that('that default rendering works', {
     timestep
   )
 
+  mockery::expect_args(
+    renderer$render_mock(),
+    3,
+    'p_detect_730_3650',
+    1.5,
+    timestep
+  )
+
 })
 
 test_that('that default rendering works when no one is in the age range', {
@@ -151,79 +159,22 @@ test_that('that severe rendering works', {
 test_that('that clinical incidence rendering works', {
   timestep <- 0
   year <- 365
-  parameters <- get_parameters(list(
-    clinical_incidence_rendering_min_ages = c(0, 2) * year,
-    clinical_incidence_rendering_max_ages = c(5, 10) * year,
-    prevalence_rendering_min_ages = NULL,
-    prevalence_rendering_max_ages = NULL
-  ))
-
   birth <- individual::IntegerVariable$new(
     -c(2, 5, 10, 11) * year
   )
 
   renderer <- mock_render(1)
-  process <- create_clinical_incidence_renderer(birth, parameters, renderer)
-
-  process(timestep, individual::Bitset$new(4)$insert(c(1, 2, 4)))
-
-  mockery::expect_args( 
-    renderer$render_mock(),
-    1,
-    'n_0_1825',
-    2,
-    timestep
-  )
-
-  mockery::expect_args( 
-    renderer$render_mock(),
-    2,
-    'n_inc_clinical_0_1825',
-    2,
-    timestep
-  )
-
-  mockery::expect_args( 
-    renderer$render_mock(),
-    3,
-    'n_730_3650',
-    3,
-    timestep
-  )
-
-  mockery::expect_args( 
-    renderer$render_mock(),
-    4,
-    'n_inc_clinical_730_3650',
-    2,
+  incidence_renderer(
+    birth,
+    renderer,
+    individual::Bitset$new(4)$insert(c(1, 2, 4)),
+    individual::Bitset$new(4)$insert(seq(4)),
+    c(.1, .2, .3, .4),
+    'inc_clinical_',
+    c(0, 2) * year,
+    c(5, 10) * year,
     timestep
   )
-})
-
-
-test_that('that incidence rendering works', {
-  timestep <- 0
-  year <- 365
-  parameters <- get_parameters(list(
-    incidence_rendering_min_ages = c(0, 2) * year,
-    incidence_rendering_max_ages = c(5, 10) * year,
-    prevalence_rendering_min_ages = NULL,
-    prevalence_rendering_max_ages = NULL
-  ))
-
-  birth <- individual::IntegerVariable$new(
-    -c(2, 5, 10, 11) * 365
-  )
-
-  is_severe <- individual::CategoricalVariable$new(
-    c('yes', 'no'),
-    c('no', 'yes', 'no', 'no')
-  )
-
-  renderer <- mock_render(1)
-  process <- create_incidence_renderer(birth, is_severe, parameters, renderer)
-
-  process(timestep, individual::Bitset$new(4)$insert(c(1, 2, 4)))
 
   mockery::expect_args(
     renderer$render_mock(),
@@ -236,7 +187,7 @@ test_that('that incidence rendering works', {
   mockery::expect_args(
     renderer$render_mock(),
     2,
-    'n_inc_0_1825',
+    'n_inc_clinical_0_1825',
     2,
     timestep
   )
@@ -244,73 +195,32 @@ test_that('that incidence rendering works', {
   mockery::expect_args(
     renderer$render_mock(),
     3,
-    'n_730_3650',
-    3,
+    'p_inc_clinical_0_1825',
+    .3,
     timestep
   )
 
   mockery::expect_args(
     renderer$render_mock(),
     4,
-    'n_inc_730_3650',
-    2,
-    timestep
-  )
-})
-
-test_that('that severe incidence rendering works', {
-  timestep <- 0
-  year <- 365
-  parameters <- get_parameters(list(
-    severe_incidence_rendering_min_ages = c(0, 2) * year,
-    severe_incidence_rendering_max_ages = c(5, 10) * year,
-    prevalence_rendering_min_ages = NULL,
-    prevalence_rendering_max_ages = NULL
-  ))
-
-  birth <- individual::IntegerVariable$new(
-    -c(2, 6, 10, 11) * 365
-  )
-
-  is_severe <- individual::CategoricalVariable$new(
-    c('yes', 'no'),
-    c('no', 'yes', 'no', 'no')
-  )
-
-  renderer <- mock_render(1)
-  process <- create_incidence_renderer(birth, is_severe, parameters, renderer)
-
-  process(timestep, individual::Bitset$new(4)$insert(c(1, 2, 4)))
-
-  mockery::expect_args(
-    renderer$render_mock(),
-    1,
-    'n_0_1825',
-    1,
+    'n_730_3650',
+    3,
     timestep
   )
 
   mockery::expect_args(
     renderer$render_mock(),
+    5,
+    'n_inc_clinical_730_3650',
     2,
-    'n_inc_severe_0_1825',
-    0,
     timestep
   )
 
   mockery::expect_args(
     renderer$render_mock(),
-    3,
-    'n_730_3650',
-    3,
-    timestep
-  )
-
-  mockery::expect_args(
-    renderer$render_mock(),
-    4,
-    'n_inc_severe_730_3650',
-    1,
+    6,
+    'p_inc_clinical_730_3650',
+    .6,
     timestep
   )
 })
diff --git a/tests/testthat/test-rtss.R b/tests/testthat/test-rtss.R
index 060a8a21..3ee0619a 100644
--- a/tests/testthat/test-rtss.R
+++ b/tests/testthat/test-rtss.R
@@ -90,6 +90,7 @@ test_that('Infection considers vaccine efficacy', {
     variables,
     bitten_humans = individual::Bitset$new(4)$insert(seq(4)),
     parameters,
+    mock_render(timestep),
     timestep
   )
 
diff --git a/vignettes/EIRprevmatch.Rmd b/vignettes/EIRprevmatch.Rmd
new file mode 100644
index 00000000..db608a14
--- /dev/null
+++ b/vignettes/EIRprevmatch.Rmd
@@ -0,0 +1,176 @@
+---
+title: "Matching PfPR2-10 to EIR"
+output: rmarkdown::html_vignette
+vignette: >
+  %\VignetteIndexEntry{EIRprevmatch}
+  %\VignetteEngine{knitr::rmarkdown}
+  %\VignetteEncoding{UTF-8}
+---
+
+```{r, include = FALSE}
+knitr::opts_chunk$set(
+  collapse = TRUE,
+  comment = "#>"
+)
+```
+
+```{r setup, message=F, warning=F}
+library(malariasimulation)
+library(mgcv)
+```
+
+Begin by parameterizing the model you wish to run. Be sure to include any seasonality factors and interventions you want to include at baseline. 
+
+```{r parameters}
+year <- 365
+human_population <- 5000
+
+params <- get_parameters(list(
+  human_population = human_population,
+  average_age = 8453.323,     # to match flat_demog
+  model_seasonality = TRUE,   # assign seasonality
+  g0 = 0.284596,                              
+  g = c(-0.317878, -0.0017527, 0.116455),     
+  h = c(-0.331361, 0.293128, -0.0617547),     
+  prevalence_rendering_min_ages = 2 * year,    # set prev age range
+  prevalence_rendering_max_ages = 10 * year,
+  fvt = 0,
+  v = 0,
+  severe_enabled = T,
+  individual_mosquitoes = FALSE))
+
+  # set species / drugs / treatment parameters
+params <- set_species(params, species = list(arab_params, fun_params, gamb_params), 
+                      proportions = c(0.25, 0.25, 0.5))
+params <- set_drugs(params, list(AL_params))
+params <- set_clinical_treatment(params, 1, c(1), c(0.45))
+
+```
+
+Now run `malariasimulation` over a range of EIRs. The goal is to run enough points to generate a curve to which you can match PfPR2-10 to EIR effectively.
+
+```{r model}
+# loop over malariasimulation runs
+init_EIR <- c(0.01, 0.1, 1, 5, 10, 25, 50) # set EIR values
+
+# run model
+outputs <- lapply(
+  init_EIR,
+  function(init) {
+    p_i <- set_equilibrium(params, init)
+    run_simulation(5 * year, p_i) # sim time = 5 years
+  }
+)
+
+# output EIR values
+EIR <- lapply(
+  outputs,
+  function(output) {
+    mean(
+      rowSums(
+        output[
+          output$timestep %in% seq(4 * 365, 5 * 365), # just use data from the last year
+          grepl('EIR_', names(output))
+        ] / human_population * year
+      )
+    )
+  }
+)
+
+# output prev 2-10 values
+prev <- lapply(
+  outputs,
+  function(output) {
+    mean(
+      output[
+        output$timestep %in% seq(4 * 365, 5 * 365),
+        'n_detect_730_3650'
+      ] / output[
+        output$timestep %in% seq(4 * 365, 5 * 365),
+        'n_730_3650'
+      ]
+    )
+  }
+)
+
+# create dataframe of initial EIR, output EIR, and prev 2-10 results
+EIR_prev <- cbind.data.frame(init_EIR, output_EIR = unlist(EIR), prev = unlist(prev))
+
+```
+
+Now plot your results! Code is included to compare the results of matching PfPR2-10 to EIR based on `malariaEquilibrium` (blue line) versus matching based on parameterized `malariasimulation` runs (red line). Notice that the generated points do not form a smooth curve. We ran `malariasimulation` using a population of just 5,000. Increasing the population to 10,000 or even 100,000 will generate more accurate estimates, but will take longer to run. 
+
+```{r plot}
+# calculate EIR / prev 2-10 relationship from malariaEquilibrium
+eir <- seq(from = 0.1, to = 50, by=.5)
+eq_params <- malariaEquilibrium::load_parameter_set("Jamie_parameters.rds")
+
+prev <- vapply( # calculate prevalence between 2:10 for a bunch of EIRs
+  eir,
+  function(eir) {
+    eq <- malariaEquilibrium::human_equilibrium(
+      eir,
+      ft=0,
+      p=eq_params,
+      age=0:100
+    )
+    sum(eq$states[2:10, 'pos_M']) / sum(eq$states[2:10, 'prop'])
+  },
+  numeric(1)
+)
+
+prevmatch <- cbind.data.frame(eir, prev)
+
+# calculate best fit line through malariasimulation data
+fit <- predict(gam(prev~s(init_EIR, k=5), data=EIR_prev), 
+               newdata = data.frame(init_EIR=c(0,seq(0.1,50,0.1))),  type="response")
+fit <- cbind(fit, data.frame(init_EIR=c(0,seq(0.1,50,0.1))))
+
+# plot
+plot(x=1, type = "n", 
+     frame = F, 
+     xlab = "initial EIR", 
+     ylab = expression(paste(italic(Pf),"PR"[2-10])),
+     xlim = c(0,50), 
+     ylim = c(0,.7))
+
+points(x=EIR_prev$init_EIR, 
+       y=EIR_prev$prev,
+       pch = 19,
+       col = 'black')
+
+lines(x=fit$init_EIR, 
+      y=fit$fit, 
+      col = "red", 
+      type = "l", 
+      lty = 1)
+
+lines(x=prevmatch$eir, 
+      y=prevmatch$prev, 
+      col = "blue", 
+      type = "l", 
+      lty = 1)
+
+legend("bottomright", legend=c("malariasimulation", "malariaEquilibrium"), col=c("red", "blue"), lty = c(1,1), cex=0.8)
+
+```
+
+Now extract values from the fitted line to generate EIR estimates at your desired PfPR2-10 values.
+
+```{r table, warning=F}
+# Pre-intervention baseline PfPR2-10 starting at values 10, 25, 35, 55 
+PfPR <- c(.10, .25, .35, .45)
+
+# match via stat_smooth predictions
+match <- function(x){
+  
+  m <- which.min(abs(fit$fit-x)) # index of closest prev match
+  fit[m,2]                       # print match
+  
+}
+
+eir <- unlist(lapply(PfPR, match))
+
+cbind.data.frame(PfPR, eir)
+
+```
diff --git a/vignettes/Treatment.Rmd b/vignettes/Treatment.Rmd
index d2be6bfd..3986c202 100644
--- a/vignettes/Treatment.Rmd
+++ b/vignettes/Treatment.Rmd
@@ -24,17 +24,11 @@ Set up some baseline parameters:
 
 ```{r}
 year <- 365
-sim_length <- 1 * year
+sim_length <- 2 * year
 human_population <- 1000
 
 simparams <- get_parameters(list(human_population = human_population))
 simparams <- set_drugs(simparams, list(AL_params, DHA_PQP_params))
-
-# NOTE: this equilibrium is calculated using a different prophylaxis model
-# it is just used as a starting point
-get_equilibrium <- function(EIR, ft) {
-  human_equilibrium(EIR = EIR, ft = ft, p = jamie_params, age = 0:99)
-}
 ```
 
 Run the model for some starting EIRs and fts:
@@ -45,11 +39,12 @@ starting_EIR <- 10
 fts <- c(0, .25, .5, .75, 1.)
 
 for (ft in fts) {
-  eq_params <- load_parameter_set()
   simparams <- set_clinical_treatment(simparams, 1, 1, ft / 2)
   simparams <- set_clinical_treatment(simparams, 2, 1, ft / 2)
-  simparams <- set_equilibrium(simparams, starting_EIR, eq_params)
+  simparams <- set_equilibrium(simparams, starting_EIR)
   output <- run_simulation(sim_length, simparams)
+  # cut out the warmup
+  output <- output[output$timestep > year,]
   outputs[[length(outputs) + 1]] <- output
 }
 ```
diff --git a/vignettes/Vaccines.Rmd b/vignettes/Vaccines.Rmd
index 994af5c0..fb434ed0 100644
--- a/vignettes/Vaccines.Rmd
+++ b/vignettes/Vaccines.Rmd
@@ -17,7 +17,6 @@ knitr::opts_chunk$set(
 ```{r setup}
 suppressPackageStartupMessages(library(ggplot2))
 library(malariasimulation)
-library(malariaEquilibrium)
 ```
 
 # Parameterisation
diff --git a/vignettes/Variation.Rmd b/vignettes/Variation.Rmd
new file mode 100644
index 00000000..e8f39746
--- /dev/null
+++ b/vignettes/Variation.Rmd
@@ -0,0 +1,91 @@
+---
+title: "Variation"
+output: html_document
+vignette: >
+  %\VignetteIndexEntry{Variation}
+  %\VignetteEngine{knitr::rmarkdown}
+  %\VignetteEncoding{UTF-8}
+---
+
+```{r setup, include=FALSE}
+library(malariasimulation)
+library(cowplot)
+library(ggplot2)
+```
+
+## Variation in outputs
+
+Malariasimulation is a stochastic model, so there will be be variation in model outputs. For example, we could look at detected malaria cases over a year...
+
+```{r}
+# A small population
+p <- get_parameters(list(
+  human_population = 1000,
+  individual_mosquitoes = FALSE
+))
+p <- set_equilibrium(p, 2)
+small_o <- run_simulation(365, p)
+p <- get_parameters(list(
+  human_population = 10000,
+  individual_mosquitoes = FALSE
+))
+p <- set_equilibrium(p, 2)
+big_o <- run_simulation(365, p)
+
+plot_grid(
+  ggplot(small_o) + geom_line(aes(timestep, n_detect_730_3650)) + ggtitle('n_detect at n = 1,000'),
+  ggplot(small_o) + geom_line(aes(timestep, p_detect_730_3650)) + ggtitle('p_detect at n = 1,000'),
+  ggplot(big_o) + geom_line(aes(timestep, n_detect_730_3650)) + ggtitle('n_detect at n = 10,000'),
+  ggplot(big_o) + geom_line(aes(timestep, p_detect_730_3650)) + ggtitle('p_detect at n = 10,000')
+)
+```
+
+The n_detect output shows the result of sampling individuals who would be detected by microscopy.
+
+While the p_detect output shows the sum of probabilities of detection in the population.
+
+Notice that the p_detect output is slightly smoother. That's because it forgoes the sampling step. At low population sizes, p_detect will be smoother than the n_detect counterpart.
+
+## Estimating variation 
+
+We can estimate the variation in the number of detectable cases by repeating the simulation several times...
+
+```{r}
+outputs <- run_simulation_with_repetitions(
+  timesteps = 365,
+  repetitions = 10,
+  overrides = p,
+  parallel=TRUE
+)
+
+df <- aggregate(
+  outputs$n_detect_730_3650,
+  by=list(outputs$timestep),
+  FUN = function(x) {
+    c(
+      median = median(x),
+      lowq = unname(quantile(x, probs = .25)),
+      highq = unname(quantile(x, probs = .75)),
+      mean = mean(x),
+      lowci = mean(x) + 1.96*sd(x),
+      highci = mean(x) - 1.96*sd(x)
+    )
+  }
+)
+
+df <- data.frame(cbind(t = df$Group.1, df$x))
+
+plot_grid(
+  ggplot(df)
+  + geom_line(aes(x = t, y = median, group = 1))
+  + geom_ribbon(aes(x = t, ymin = lowq, ymax = highq), alpha = .2)
+  + xlab('timestep') + ylab('n_detect') + ggtitle('IQR spread'),
+  ggplot(df)
+  + geom_line(aes(x = t, y = mean, group = 1))
+  + geom_ribbon(aes(x = t, ymin = lowci, ymax = highci), alpha = .2)
+  + xlab('timestep') + ylab('n_detect') + ggtitle('95% confidence interval')
+)
+
+```
+
+These are useful techniques for comparing the expected variance from model outputs to outcomes from trials with lower population sizes.
\ No newline at end of file