diff --git a/DESCRIPTION b/DESCRIPTION
index 74445497..ab68ccb0 100644
--- a/DESCRIPTION
+++ b/DESCRIPTION
@@ -69,15 +69,13 @@ LazyData: true
 Remotes:
   mrc-ide/malariaEquilibrium,
   mrc-ide/individual
-Additional_repositories:
-  https://mrc-ide.r-universe.dev
 Imports:
     individual (>= 0.1.16),
     malariaEquilibrium (>= 1.0.1),
     Rcpp,
     statmod,
     MASS,
-    dqrng (>= 0.3.2.2),
+    dqrng (>= 0.4),
     sitmo,
     BH,
     R6,
diff --git a/R/biting_process.R b/R/biting_process.R
index 801aea20..279802ea 100644
--- a/R/biting_process.R
+++ b/R/biting_process.R
@@ -210,9 +210,7 @@ simulate_bites <- function(
 calculate_eir <- function(species, solvers, variables, parameters, timestep) {
   a <- human_blood_meal_rate(species, variables, parameters, timestep)
   infectious <- calculate_infectious(species, solvers, variables, parameters)
-  age <- get_age(variables$birth$get_values(), timestep)
-  psi <- unique_biting_rate(age, parameters)
-  infectious * a * mean(psi)
+  infectious * a
 }
 
 effective_biting_rates <- function(a, .pi, p_bitten) {
diff --git a/R/model.R b/R/model.R
index 67f79758..d59132f5 100644
--- a/R/model.R
+++ b/R/model.R
@@ -29,16 +29,13 @@
 #'  * 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_lm (or pcr): number of humans with an infection detectable by microscopy (or pcr) 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.
-#'  * p_detect: the sum of probabilities of detection by microscopy between an
+#'  * p_detect_lm (or pcr): the sum of probabilities of detection by microscopy (or pcr) 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.
diff --git a/R/parameters.R b/R/parameters.R
index 647b78b2..bb692873 100644
--- a/R/parameters.R
+++ b/R/parameters.R
@@ -8,17 +8,39 @@
 #' @param overrides a named list of parameter values to use instead of defaults
 #' The parameters are defined below.
 #'
-#' fixed state transitions:
+#' initial state proportions:
+#'
+#' * s_proportion - the proportion of `human_population` that begin as susceptible; default = 0.420433246
+#' * d_proportion - the proportion of `human_population` that begin with
+#' clinical disease; default = 0.007215064
+#' * a_proportion - the proportion of `human_population` that begin as
+#' asymptomatic; default = 0.439323667
+#' * u_proportion - the proportion of `human_population` that begin as
+#' subpatents; default = 0.133028023
+#' * t_proportion - the proportion of `human_population` that begin treated; default = 0
+#'
+#' human fixed state transitions:
 #'
 #' * dd - the delay for humans to move from state D to A; default = 5
 #' * dt - the delay for humans to move from state Tr to S; default = 5
 #' * da - the delay for humans to move from state A to U; default = 195
 #' * du - the delay for humans to move from state U to S; default = 110
-#' * del - the delay for mosquitoes to move from state E to L; default = 6.64
-#' * dl - the delay for mosquitoes to move from state L to P; default = 3.72
-#' * dpl - the delay mosquitoes to move from state P to Sm; default = 0.643
-#' * mup - the rate at which pupal mosquitoes die; default = 0.249
-#' * mum - the rate at which developed mosquitoes die; default (An. gambiae) = .132
+#' 
+#' human demography parameters:
+#'
+#' * human_population - the initial number of humans to model; default = 100
+#' * average_age - the average age of humans (in timesteps), this is only used
+#' if custom_demography is FALSE; default = 7663
+#' * custom_demography - population demography given; default = FALSE
+#'
+#' initial immunity values:
+#'
+#' * init_icm - the immunity from clinical disease at birth; default = 0
+#' * init_ivm - the immunity from severe disease at birth; default = 0
+#' * init_ib  - the initial pre-erythrocitic immunity; default = 0
+#' * init_ica - the initial acquired immunity from clinical disease; default = 0
+#' * init_iva - the initial acquired immunity from severe disease; default = 0
+#' * init_id  - the initial acquired immunity to detectability; default = 0
 #'
 #' immunity decay rates:
 #'
@@ -29,6 +51,26 @@
 #' * rva - decay rate for acquired immunity to severe disease; default = 10950
 #' * rid - decay rate for acquired immunity to detectability; default = 3650
 #'
+#' immunity boost grace periods:
+#'
+#' * ub - period in which pre-erythrocytic immunity cannot be boosted; default = 7.2
+#' * uc - period in which clinical immunity cannot be boosted; default = 6.06
+#' * uv - period in which severe immunity cannot be boosted; default = 11.4321
+#' * ud - period in which immunity to detectability cannot be boosted; default = 9.44512
+#'
+#' maternal immunity parameters:
+#'
+#' * pcm - new-born clinical immunity relative to mother's; default = 0.774368
+#' * pvm - new-born severe immunity relative to mother's; default = 0.195768
+#'
+#' unique biting rate:
+#'
+#' * a0 - age dependent biting parameter; default = 2920
+#' * rho - age dependent biting parameter; default = 0.85
+#' * sigma_squared - heterogeneity parameter; default = 1.67
+#' * n_heterogeneity_groups - number discretised groups for heterogeneity, used
+#' for sampling mothers; default = 5
+#'
 #' probability of pre-erythrocytic infection:
 #'
 #' * b0 - maximum probability due to no immunity; default = 0.59
@@ -36,6 +78,15 @@
 #' * ib0 - scale parameter; default = 43.9
 #' * kb - shape parameter; default = 2.16
 #'
+#' probability of detection by light-microscopy when asymptomatic:
+#'
+#' * fd0 - time-scale at which immunity changes with age; default = 0.007055
+#' * ad - scale parameter relating age to immunity; default = 7993.5
+#' * gammad - shape parameter relating age to immunity; default = 4.8183
+#' * d1 - minimum probability due to immunity; default = 0.160527
+#' * id0 - scale parameter; default = 1.577533
+#' * kd - shape parameter; default = 0.476614
+#'
 #' probability of clinical infection:
 #'
 #' * phi0 - maximum probability due to no immunity; default = 0.792
@@ -53,22 +104,6 @@
 #' * av - age dependent modifier; default = 2493.41
 #' * gammav - age dependent modifier; default = 2.91282
 #'
-#' immunity reducing probability of detection:
-#'
-#' * fd0 - time-scale at which immunity changes with age; default = 0.007055
-#' * ad - scale parameter relating age to immunity; default = 7993.5
-#' * gammad - shape parameter relating age to immunity; default = 4.8183
-#' * d1 - minimum probability due to immunity; default = 0.160527
-#' * id0 - scale parameter; default = 1.577533
-#' * kd - shape parameter; default = 0.476614
-#'
-#' immunity boost grace periods:
-#'
-#' * ub - period in which pre-erythrocytic immunity cannot be boosted; default = 7.2
-#' * uc - period in which clinical immunity cannot be boosted; default = 6.06
-#' * uv - period in which severe immunity cannot be boosted; default = 11.4321
-#' * ud - period in which immunity to detectability cannot be boosted; default = 9.44512
-#'
 #' infectivity towards mosquitoes:
 #'
 #' * cd - infectivity of clinically diseased humans towards mosquitoes; default = 0.068
@@ -76,62 +111,21 @@
 #' * cu - infectivity of sub-patent infection; default = 0.0062
 #' * ct - infectivity of treated infection; default = 0.021896
 #'
-#' unique biting rate:
-#'
-#' * a0 - age dependent biting parameter; default = 2920
-#' * rho - age dependent biting parameter; default = 0.85
-#' * sigma_squared - heterogeneity parameter; default = 1.67
-#' * n_heterogeneity_groups - number discretised groups for heterogeneity, used
-#' for sampling mothers; default = 5
-#'
-#' mortality parameters:
-#'
-#' * average_age - the average age of humans (in timesteps), this is only used
-#' if custom_demography is FALSE; default = 7663
-#' * pcm - new-born clinical immunity relative to mother's; default = 0.774368
-#' * pvm - new-born severe immunity relative to mother's; default = 0.195768
+#' mosquito fixed state transitions (including mortality):
+#' 
+#' * del - the delay for mosquitoes to move from state E to L; default = 6.64
+#' * dl - the delay for mosquitoes to move from state L to P; default = 3.72
+#' * dpl - the delay mosquitoes to move from state P to Sm; default = 0.643
 #' * me - early stage larval mortality rate; default = 0.0338
 #' * ml - late stage larval mortality rate; default = 0.0348
-#'
-#' carrying capacity parameters:
-#'
-#' * model_seasonality - boolean switch TRUE iff the simulation models seasonal rainfall; default = FALSE
-#' * g0 - rainfall fourier parameter; default = 2
-#' * g - rainfall fourier parameter; default = 0.3, 0.6, 0.9
-#' * h - rainfall fourier parameters; default = 0.1, 0.4, 0.7
-#' * gamma - effect of density dependence on late instars relative to early
-#' instars; default = 13.25
-#' * rainfall_floor - the minimum rainfall value (must be above 0); default 0.001
-#'
-#' initial state proportions:
-#'
-#' * s_proportion - the proportion of `human_population` that begin as susceptible; default = 0.420433246
-#' * d_proportion - the proportion of `human_population` that begin with
-#' clinical disease; default = 0.007215064
-#' * a_proportion - the proportion of `human_population` that begin as
-#' asymptomatic; default = 0.439323667
-#' * u_proportion - the proportion of `human_population` that begin as
-#' subpatents; default = 0.133028023
-#' * t_proportion - the proportion of `human_population` that begin treated; default = 0
-#'
-#' initial immunity values:
-#'
-#' * init_icm - the immunity from clinical disease at birth; default = 0
-#' * init_ivm - the immunity from severe disease at birth; default = 0
-#' * init_ib  - the initial pre-erythrocitic immunity; default = 0
-#' * init_ica - the initial acquired immunity from clinical disease; default = 0
-#' * init_iva - the initial acquired immunity from severe disease; default = 0
-#' * init_id  - the initial acquired immunity to detectability; default = 0
-#'
-#' incubation periods:
-#'
-#' * de - Duration of the human latent period of infection; default = 12
-#' * delay_gam - Lag from parasites to infectious gametocytes; default = 12.5
-#' * dem - Extrinsic incubation period in mosquito population model; default = 10
+#' * mup - the rate at which pupal mosquitoes die; default = 0.249
+#' * mum - the rate at which developed mosquitoes die; default (An. gambiae) = .132
 #'
 #' vector biology:
 #' species specific values are vectors
-#'
+#' please set species parameters using the convenience function
+#' \code{\link{set_species}}
+#' 
 #' * beta - the average number of eggs laid per female mosquito per day; default = 21.2
 #' * total_M - the initial number of adult mosquitos in the simulation; default = 1000
 #' * init_foim - the FOIM used to calculate the equilibrium state for mosquitoes; default = 0
@@ -141,18 +135,30 @@
 #' * Q0 - proportion of blood meals taken on humans; default = 0.92
 #' * foraging_time - time spent taking blood meals; default = 0.69
 #'
-#' feeding cycle:
-#' please set vector control strategies using `set_betnets` and `set_spraying`
+#' seasonality and carrying capacity parameters:
+#' please set flexible carrying capacity using the convenience function
+#' \code{\link{set_carrying_capacity}}
+#' 
+#' * model_seasonality - boolean switch TRUE iff the simulation models seasonal rainfall; default = FALSE
+#' * g0 - rainfall fourier parameter; default = 2
+#' * g - rainfall fourier parameter; default = 0.3, 0.6, 0.9
+#' * h - rainfall fourier parameters; default = 0.1, 0.4, 0.7
+#' * gamma - effect of density dependence on late instars relative to early
+#' instars; default = 13.25
+#' * rainfall_floor - the minimum rainfall value (must be above 0); default 0.001
+#' * carrying_capacity; default = FALSE
+#' * carrying_capacity_timesteps; default = NULL
+#' * carrying_capacity_values; default = NULL#'
 #'
-#' * bednets - boolean for if bednets are enabled; default = FALSE
-#' * phi_bednets - proportion of bites taken in bed; default = 0.85
-#' * k0 - proportion of females bloodfed with no net; default = 0.699
-#' * spraying - boolean for if indoor spraying is enabled; default = FALSE
-#' * phi_indoors - proportion of bites taken indoors; default = 0.90
+#' parasite incubation periods:
 #'
+#' * de - duration of the human latent period of infection; default = 12
+#' * delay_gam - lag from parasites to infectious gametocytes; default = 12.5
+#' * dem - extrinsic incubation period in mosquito population model; default = 10
+#' 
 #' treatment parameters:
-#' please set treatment parameters with the convenience functions in
-#' `drug_parameters.R`
+#' please set treatment parameters with the convenience functions
+#' \code{\link{set_drugs}} and \code{\link{set_clinical_treatment}}
 #'
 #' * drug_efficacy - a vector of efficacies for available drugs; default = turned off
 #' * drug_rel_c - a vector of relative onward infectiousness values for drugs; default = turned off
@@ -164,21 +170,32 @@
 #' clinically diseased (these values refer to the index in drug_* parameters); default = NULL, NULL, NULL
 #' * clinical_treatment_coverage - a vector of coverage values for each drug; default = NULL, NULL, NULL
 #'
-#' PEV parameters: 
-#' please set vaccine strategies with the convenience functions in
-#' `pev_parameters.R:set_pev_epi`
-#' `pev_parameters.R:set_mass_pev`
+#' MDA, SMC and PMC parameters:
+#' please set these parameters with the convenience functions
+#' \code{\link{set_mda}}, \code{\link{set_smc}} and \code{\link{set_pmc}},
+#' with \code{\link{peak_season_offset}}
+#'
+#' bednet, irs and mosquito feeding cycle parameters:
+#' please set vector control strategies using \code{\link{set_bednets}} and \code{\link{set_spraying}}
+#'
+#' * bednets - boolean for if bednets are enabled; default = FALSE
+#' * phi_bednets - proportion of bites taken in bed; default = 0.85
+#' * k0 - proportion of females bloodfed with no net; default = 0.699
+#' * spraying - boolean for if indoor spraying is enabled; default = FALSE
+#' * phi_indoors - proportion of bites taken indoors; default = 0.90
+#'
+#'
+#' PEV parameters:
+#' please set vaccine strategies with the convenience functions
+#' \code{\link{set_pev_epi}} and \code{\link{set_mass_pev}}
 #'
 #' * pev_doses - the dosing schedule before the vaccine takes effect; default =
 #' c(0, 1.5 * 30, 3 * 30)
 #' default = 365
 #'
-#' MDA, SMC and PMC parameters:
-#' please set these parameters with the convenience functions in `mda_parameters.R`
-#'
 #' TBV parameters:
 #' please set TBV parameters with the convenience functions in
-#' `vaccine_parameters.R:set_tbv`
+#' \code{\link{set_tbv}}
 #'
 #' * tbv_mt - effect on treated infectiousness; default = 35
 #' * tbv_md - effect on diseased infectiousness; default = 46.7
@@ -195,7 +212,7 @@
 #' 
 #' Antimalarial resistance parameters:
 #' please set antimalarial resistance parameters with the convenience functions in
-#' `antimalarial_resistance.R:set_antimalarial_resistance`
+#' \code{\link{set_antimalarial_resistance}}
 #' 
 #' * antimalarial_resistance - boolean for if antimalarial resistance is enabled; default = FALSE
 #' * antimalarial_resistance_drug - vector of drugs for which resistance can be parameterised; default = NULL
@@ -210,11 +227,11 @@
 #' * dt_slow_parasite_clearance - the delay for humans experiencing slow parasite clearance to move from state Tr to S; default = NULL
 #'
 #' rendering:
-#' All values are in timesteps and all ranges are inclusive
-#'
-#' * prevalence_rendering_min_ages - the minimum ages for clinical prevalence
-#' outputs; default = 730
-#' * prevalence_rendering_max_ages - the corresponding max ages; default = 3650
+#' All values are in timesteps and all ranges are inclusive.
+#' Please set rendered age groups using the convenience function
+#' 
+#' * age_group_rendering_min_ages - the minimum ages for population size outputs; default = turned off
+#' * age_group_rendering_max_ages - the corresponding max ages; default = turned off
 #' * incidence_rendering_min_ages - the minimum ages for incidence
 #' outputs (includes asymptomatic microscopy +); default = turned off
 #' * incidence_rendering_max_ages - the corresponding max ages; default = turned off
@@ -223,16 +240,18 @@
 #' * severe_incidence_rendering_min_ages - the minimum ages for severe incidence
 #' outputs; default = turned off
 #' * severe_incidence_rendering_max_ages - the corresponding max ages; default = turned off
+#' * prevalence_rendering_min_ages - the minimum ages for clinical prevalence
+#' outputs; default = 730
+#' * prevalence_rendering_max_ages - the corresponding max ages; default = 3650
 #'
 #' miscellaneous:
 #'
-#' * human_population - the initial number of humans to model; default = 100
-#' * human_population_timesteps - the timesteps at which the population should
-#' change; default = 0
 #' * mosquito_limit - the maximum number of mosquitoes to allow for in the
 #' simulation; default = 1.00E+05
 #' * individual_mosquitoes - boolean whether adult mosquitoes are modeled
 #' individually or compartmentally; default = TRUE
+#' * human_population_timesteps - the timesteps at which the population should
+#' change; default = 0
 #' * r_tol - the relative tolerance for the ode solver; default = 1e-4
 #' * a_tol - the absolute tolerance for the ode solver; default = 1e-4
 #' * ode_max_steps - the max number of steps for the solver; default = 1e6
@@ -242,17 +261,28 @@
 #' @export
 get_parameters <- function(overrides = list()) {
   parameters <- list(
+    # initial state proportions
+    s_proportion = 0.420433246,
+    d_proportion = 0.007215064,
+    a_proportion = 0.439323667,
+    u_proportion = 0.133028023,
+    t_proportion = 0,
+    # human fixed state transitions
     dd    = 5,
     dt    = 5,
     da    = 195,
     du    = 110,
-    del   = 6.64,
-    dl    = 3.72,
-    dpl   = .643,
-    mup   = .249,
-    mum   = .132,
-    sigma_squared   = 1.67,
-    n_heterogeneity_groups = 5,
+    # human demography parameters
+    human_population = 100,
+    average_age = 7663,
+    custom_demography = FALSE,
+    # initial immunities
+    init_ica = 0,
+    init_iva = 0,
+    init_icm = 0,
+    init_ivm = 0,
+    init_id  = 0,
+    init_ib  = 0,
     # immunity decay rates
     rm    = 67.6952,
     rvm   = 76.8365,
@@ -260,24 +290,32 @@ get_parameters <- function(overrides = list()) {
     rc    = 30 * 365,
     rva   = 30 * 365,
     rid   = 10 * 365,
-    # blood immunity parameters
-    b0    = 0.59,
-    b1    = 0.5,
-    ib0   = 43.9,
-    kb    = 2.16,
     # immunity boost grace periods
     ub    = 7.2,
     uc    = 6.06,
     uv    = 11.4321,
     ud    = 9.44512,
-    # infectivity towards mosquitos
-    cd    = 0.068,
-    gamma1= 1.82425,
-    cu    = 0.0062,
-    ct    = 0.021896,
+    # maternal immunity parameters
+    pcm   = .774368,
+    pvm   = .195768,
     # unique biting rate
     a0    = 8 * 365,
     rho   = .85,
+    sigma_squared   = 1.67,
+    n_heterogeneity_groups = 5,
+    enable_heterogeneity = TRUE,
+    # blood immunity parameters
+    b0    = 0.59,
+    b1    = 0.5,
+    ib0   = 43.9,
+    kb    = 2.16,
+    # asymptomatic immunity to lm detection parameters
+    fd0   = 0.007055,
+    ad    = 21.9 * 365,
+    gammad= 4.8183,
+    d1    = 0.160527,
+    id0   = 1.577533,
+    kd    = .476614,
     # clinical immunity parameters
     phi0 = .792,
     phi1 = .00074,
@@ -291,62 +329,44 @@ get_parameters <- function(overrides = list()) {
     av      = 2493.41,
     gammav  = 2.91282,
     iv0     = 1.09629,
-    # delay for infection
-    de        = 12,
-    delay_gam = 12.5,
-    dem       = 10,
-    # asymptomatic immunity parameters
-    fd0   = 0.007055,
-    ad    = 21.9 * 365,
-    gammad= 4.8183,
-    d1    = 0.160527,
-    id0   = 1.577533,
-    kd    = .476614,
-    # mortality parameters
-    average_age = 7663,
-    pcm   = .774368,
-    pvm   = .195768,
-    # carrying capacity parameters
-    g0    = 2,
-    g     = c(.3, .6, .9),
-    h     = c(.1, .4, .7),
-    gamma = 13.25,
-    model_seasonality = FALSE,
-    rainfall_floor = 0.001,
-    # larval mortality rates
+    # infectivity towards mosquitos
+    cd    = 0.068,
+    gamma1= 1.82425,
+    cu    = 0.0062,
+    ct    = 0.021896,
+    # mosquito fixed state transitions (inc. mortality)
+    del   = 6.64,
+    dl    = 3.72,
+    dpl   = .643,
     me    = .0338,
     ml    = .0348,
-    # initial state proportions
-    s_proportion = 0.420433246,
-    d_proportion = 0.007215064,
-    a_proportion = 0.439323667,
-    u_proportion = 0.133028023,
-    t_proportion = 0,
-    # initial immunities
-    init_ica = 0,
-    init_iva = 0,
-    init_icm = 0,
-    init_ivm = 0,
-    init_id  = 0,
-    init_ib  = 0,
-    # vector biology
-    beta     = 21.2,
-    total_M  = 1000,
-    init_foim= 0,
+    mup   = .249,
+    mum   = .132,
     # species-specific vector biology (default is An. gambiae s.s)
     species             = 'gamb',
     species_proportions = 1,
     blood_meal_rates    = 1/3,
     Q0                  = .92,
     foraging_time       = .69,
-    # bed nets
-    bednets = FALSE,
-    phi_bednets = .85,
-    k0 = .699,
-    # indoor spraying
-    spraying = FALSE,
-    phi_indoors = .90,
-    # treatment
+    beta     = 21.2,
+    total_M  = 1000,
+    init_foim= 0,
+    # carrying capacity parameters
+    g0    = 2,
+    g     = c(.3, .6, .9),
+    h     = c(.1, .4, .7),
+    gamma = 13.25,
+    model_seasonality = FALSE,
+    rainfall_floor = 0.001,
+    # flexible carrying capacity
+    carrying_capacity = FALSE,
+    carrying_capacity_timesteps = NULL,
+    carrying_capacity_values = NULL,
+    # parasite incubation delays
+    de        = 12,
+    delay_gam = 12.5,
+    dem       = 10,
+    # treatment parameters
     drug_efficacy          = numeric(0),
     drug_rel_c             = numeric(0),
     drug_prophylaxis_shape = numeric(0),
@@ -354,9 +374,6 @@ get_parameters <- function(overrides = list()) {
     clinical_treatment_drugs     = list(),
     clinical_treatment_timesteps = list(),
     clinical_treatment_coverages = list(),
-    # rts,s
-    pev = FALSE,
-    pev_doses = c(0, 1.5 * 30, 3 * 30),
     # MDA
     mda = FALSE,
     mda_drug = 0,
@@ -376,6 +393,16 @@ get_parameters <- function(overrides = list()) {
     pmc_timesteps = NULL,
     pmc_coverages = NULL,
     pcs_ages = -1,
+    # bed nets
+    bednets = FALSE,
+    phi_bednets = .85,
+    k0 = .699,
+    # indoor spraying
+    spraying = FALSE,
+    phi_indoors = .90,
+    # pev
+    pev = FALSE,
+    pev_doses = c(0, 1.5 * 30, 3 * 30),
     # tbv
     tbv = FALSE,
     tbv_mt = 35,
@@ -405,10 +432,6 @@ get_parameters <- function(overrides = list()) {
     late_parasitological_failure_probability = NULL,
     reinfection_during_prophylaxis_probability = NULL,
     dt_slow_parasite_clearance = NULL,
-    # flexible carrying capacity
-    carrying_capacity = FALSE,
-    carrying_capacity_timesteps = NULL,
-    carrying_capacity_values = NULL,
     # rendering
     prevalence_rendering_min_ages = 2 * 365,
     prevalence_rendering_max_ages = 10 * 365,
@@ -423,12 +446,9 @@ get_parameters <- function(overrides = list()) {
     age_group_rendering_min_ages = numeric(0),
     age_group_rendering_max_ages = numeric(0),
     # misc
-    custom_demography = FALSE,
-    human_population = 100,
-    human_population_timesteps = 0,
     mosquito_limit   = 100 * 1000,
     individual_mosquitoes = FALSE,
-    enable_heterogeneity = TRUE,
+    human_population_timesteps = 0,
     r_tol = 1e-4,
     a_tol = 1e-4,
     ode_max_steps = 1e6,
diff --git a/R/render.R b/R/render.R
index 60cb5b73..4f5816f3 100644
--- a/R/render.R
+++ b/R/render.R
@@ -4,8 +4,11 @@ in_age_range <- function(birth, timestep, lower, upper) {
 
 #' @title Render prevalence statistics
 #' 
-#' @description renders prevalence numerators and denominators for indivduals
-#' detected by microscopy and with severe malaria
+#' @description renders prevalence numerators and denominators for individuals
+#' detected by light microscopy and pcr, where those infected asymptomatically by
+#' P. falciparum have reduced probability of infection due to detectability
+#' immunity (reported as an integer sample: n_detect_lm, and summing over
+#' detection probabilities: p_detect_lm)
 #' 
 #' @param state human infection state
 #' @param birth variable for birth of the individual
@@ -32,7 +35,8 @@ create_prevelance_renderer <- function(
 
     clinically_detected <- state$get_index_of(c('Tr', 'D'))
     detected <- clinically_detected$copy()$or(asymptomatic_detected)
-
+    pcr_detected <- state$get_index_of(c('Tr', 'D', 'A', 'U'))
+    
     for (i in seq_along(parameters$prevalence_rendering_min_ages)) {
       lower <- parameters$prevalence_rendering_min_ages[[i]]
       upper <- parameters$prevalence_rendering_max_ages[[i]]
@@ -43,17 +47,22 @@ create_prevelance_renderer <- function(
         timestep
       ) 
       renderer$render(
-        paste0('n_detect_', lower, '_', upper),
+        paste0('n_detect_lm_', lower, '_', upper),
         in_age$copy()$and(detected)$size(),
         timestep
       )
       renderer$render(
-        paste0('p_detect_', lower, '_', upper),
+        paste0('p_detect_lm_', lower, '_', upper),
         in_age$copy()$and(clinically_detected)$size() + sum(
           prob[bitset_index(asymptomatic, in_age)]
         ),
         timestep
       )
+      renderer$render(
+        paste0('n_detect_pcr_', lower, '_', upper),
+        pcr_detected$copy()$and(in_age)$size(),
+        timestep
+      )
     }
   }
 }
diff --git a/man/get_parameters.Rd b/man/get_parameters.Rd
index e43792d5..fe4fb259 100644
--- a/man/get_parameters.Rd
+++ b/man/get_parameters.Rd
@@ -10,17 +10,42 @@ get_parameters(overrides = list())
 \item{overrides}{a named list of parameter values to use instead of defaults
 The parameters are defined below.
 
-fixed state transitions:
+initial state proportions:
+\itemize{
+\item s_proportion - the proportion of \code{human_population} that begin as susceptible; default = 0.420433246
+\item d_proportion - the proportion of \code{human_population} that begin with
+clinical disease; default = 0.007215064
+\item a_proportion - the proportion of \code{human_population} that begin as
+asymptomatic; default = 0.439323667
+\item u_proportion - the proportion of \code{human_population} that begin as
+subpatents; default = 0.133028023
+\item t_proportion - the proportion of \code{human_population} that begin treated; default = 0
+}
+
+human fixed state transitions:
 \itemize{
 \item dd - the delay for humans to move from state D to A; default = 5
 \item dt - the delay for humans to move from state Tr to S; default = 5
 \item da - the delay for humans to move from state A to U; default = 195
 \item du - the delay for humans to move from state U to S; default = 110
-\item del - the delay for mosquitoes to move from state E to L; default = 6.64
-\item dl - the delay for mosquitoes to move from state L to P; default = 3.72
-\item dpl - the delay mosquitoes to move from state P to Sm; default = 0.643
-\item mup - the rate at which pupal mosquitoes die; default = 0.249
-\item mum - the rate at which developed mosquitoes die; default (An. gambiae) = .132
+}
+
+human demography parameters:
+\itemize{
+\item human_population - the initial number of humans to model; default = 100
+\item average_age - the average age of humans (in timesteps), this is only used
+if custom_demography is FALSE; default = 7663
+\item custom_demography - population demography given; default = FALSE
+}
+
+initial immunity values:
+\itemize{
+\item init_icm - the immunity from clinical disease at birth; default = 0
+\item init_ivm - the immunity from severe disease at birth; default = 0
+\item init_ib  - the initial pre-erythrocitic immunity; default = 0
+\item init_ica - the initial acquired immunity from clinical disease; default = 0
+\item init_iva - the initial acquired immunity from severe disease; default = 0
+\item init_id  - the initial acquired immunity to detectability; default = 0
 }
 
 immunity decay rates:
@@ -33,6 +58,29 @@ immunity decay rates:
 \item rid - decay rate for acquired immunity to detectability; default = 3650
 }
 
+immunity boost grace periods:
+\itemize{
+\item ub - period in which pre-erythrocytic immunity cannot be boosted; default = 7.2
+\item uc - period in which clinical immunity cannot be boosted; default = 6.06
+\item uv - period in which severe immunity cannot be boosted; default = 11.4321
+\item ud - period in which immunity to detectability cannot be boosted; default = 9.44512
+}
+
+maternal immunity parameters:
+\itemize{
+\item pcm - new-born clinical immunity relative to mother's; default = 0.774368
+\item pvm - new-born severe immunity relative to mother's; default = 0.195768
+}
+
+unique biting rate:
+\itemize{
+\item a0 - age dependent biting parameter; default = 2920
+\item rho - age dependent biting parameter; default = 0.85
+\item sigma_squared - heterogeneity parameter; default = 1.67
+\item n_heterogeneity_groups - number discretised groups for heterogeneity, used
+for sampling mothers; default = 5
+}
+
 probability of pre-erythrocytic infection:
 \itemize{
 \item b0 - maximum probability due to no immunity; default = 0.59
@@ -41,6 +89,16 @@ probability of pre-erythrocytic infection:
 \item kb - shape parameter; default = 2.16
 }
 
+probability of detection by light-microscopy when asymptomatic:
+\itemize{
+\item fd0 - time-scale at which immunity changes with age; default = 0.007055
+\item ad - scale parameter relating age to immunity; default = 7993.5
+\item gammad - shape parameter relating age to immunity; default = 4.8183
+\item d1 - minimum probability due to immunity; default = 0.160527
+\item id0 - scale parameter; default = 1.577533
+\item kd - shape parameter; default = 0.476614
+}
+
 probability of clinical infection:
 \itemize{
 \item phi0 - maximum probability due to no immunity; default = 0.792
@@ -60,24 +118,6 @@ probability of severe infection:
 \item gammav - age dependent modifier; default = 2.91282
 }
 
-immunity reducing probability of detection:
-\itemize{
-\item fd0 - time-scale at which immunity changes with age; default = 0.007055
-\item ad - scale parameter relating age to immunity; default = 7993.5
-\item gammad - shape parameter relating age to immunity; default = 4.8183
-\item d1 - minimum probability due to immunity; default = 0.160527
-\item id0 - scale parameter; default = 1.577533
-\item kd - shape parameter; default = 0.476614
-}
-
-immunity boost grace periods:
-\itemize{
-\item ub - period in which pre-erythrocytic immunity cannot be boosted; default = 7.2
-\item uc - period in which clinical immunity cannot be boosted; default = 6.06
-\item uv - period in which severe immunity cannot be boosted; default = 11.4321
-\item ud - period in which immunity to detectability cannot be boosted; default = 9.44512
-}
-
 infectivity towards mosquitoes:
 \itemize{
 \item cd - infectivity of clinically diseased humans towards mosquitoes; default = 0.068
@@ -86,67 +126,21 @@ infectivity towards mosquitoes:
 \item ct - infectivity of treated infection; default = 0.021896
 }
 
-unique biting rate:
-\itemize{
-\item a0 - age dependent biting parameter; default = 2920
-\item rho - age dependent biting parameter; default = 0.85
-\item sigma_squared - heterogeneity parameter; default = 1.67
-\item n_heterogeneity_groups - number discretised groups for heterogeneity, used
-for sampling mothers; default = 5
-}
-
-mortality parameters:
+mosquito fixed state transitions (including mortality):
 \itemize{
-\item average_age - the average age of humans (in timesteps), this is only used
-if custom_demography is FALSE; default = 7663
-\item pcm - new-born clinical immunity relative to mother's; default = 0.774368
-\item pvm - new-born severe immunity relative to mother's; default = 0.195768
+\item del - the delay for mosquitoes to move from state E to L; default = 6.64
+\item dl - the delay for mosquitoes to move from state L to P; default = 3.72
+\item dpl - the delay mosquitoes to move from state P to Sm; default = 0.643
 \item me - early stage larval mortality rate; default = 0.0338
 \item ml - late stage larval mortality rate; default = 0.0348
-}
-
-carrying capacity parameters:
-\itemize{
-\item model_seasonality - boolean switch TRUE iff the simulation models seasonal rainfall; default = FALSE
-\item g0 - rainfall fourier parameter; default = 2
-\item g - rainfall fourier parameter; default = 0.3, 0.6, 0.9
-\item h - rainfall fourier parameters; default = 0.1, 0.4, 0.7
-\item gamma - effect of density dependence on late instars relative to early
-instars; default = 13.25
-\item rainfall_floor - the minimum rainfall value (must be above 0); default 0.001
-}
-
-initial state proportions:
-\itemize{
-\item s_proportion - the proportion of \code{human_population} that begin as susceptible; default = 0.420433246
-\item d_proportion - the proportion of \code{human_population} that begin with
-clinical disease; default = 0.007215064
-\item a_proportion - the proportion of \code{human_population} that begin as
-asymptomatic; default = 0.439323667
-\item u_proportion - the proportion of \code{human_population} that begin as
-subpatents; default = 0.133028023
-\item t_proportion - the proportion of \code{human_population} that begin treated; default = 0
-}
-
-initial immunity values:
-\itemize{
-\item init_icm - the immunity from clinical disease at birth; default = 0
-\item init_ivm - the immunity from severe disease at birth; default = 0
-\item init_ib  - the initial pre-erythrocitic immunity; default = 0
-\item init_ica - the initial acquired immunity from clinical disease; default = 0
-\item init_iva - the initial acquired immunity from severe disease; default = 0
-\item init_id  - the initial acquired immunity to detectability; default = 0
-}
-
-incubation periods:
-\itemize{
-\item de - Duration of the human latent period of infection; default = 12
-\item delay_gam - Lag from parasites to infectious gametocytes; default = 12.5
-\item dem - Extrinsic incubation period in mosquito population model; default = 10
+\item mup - the rate at which pupal mosquitoes die; default = 0.249
+\item mum - the rate at which developed mosquitoes die; default (An. gambiae) = .132
 }
 
 vector biology:
 species specific values are vectors
+please set species parameters using the convenience function
+\code{\link{set_species}}
 \itemize{
 \item beta - the average number of eggs laid per female mosquito per day; default = 21.2
 \item total_M - the initial number of adult mosquitos in the simulation; default = 1000
@@ -158,19 +152,32 @@ species specific values are vectors
 \item foraging_time - time spent taking blood meals; default = 0.69
 }
 
-feeding cycle:
-please set vector control strategies using \code{set_betnets} and \code{set_spraying}
+seasonality and carrying capacity parameters:
+please set flexible carrying capacity using the convenience function
+\code{\link{set_carrying_capacity}}
 \itemize{
-\item bednets - boolean for if bednets are enabled; default = FALSE
-\item phi_bednets - proportion of bites taken in bed; default = 0.85
-\item k0 - proportion of females bloodfed with no net; default = 0.699
-\item spraying - boolean for if indoor spraying is enabled; default = FALSE
-\item phi_indoors - proportion of bites taken indoors; default = 0.90
+\item model_seasonality - boolean switch TRUE iff the simulation models seasonal rainfall; default = FALSE
+\item g0 - rainfall fourier parameter; default = 2
+\item g - rainfall fourier parameter; default = 0.3, 0.6, 0.9
+\item h - rainfall fourier parameters; default = 0.1, 0.4, 0.7
+\item gamma - effect of density dependence on late instars relative to early
+instars; default = 13.25
+\item rainfall_floor - the minimum rainfall value (must be above 0); default 0.001
+\item carrying_capacity; default = FALSE
+\item carrying_capacity_timesteps; default = NULL
+\item carrying_capacity_values; default = NULL#'
+}
+
+parasite incubation periods:
+\itemize{
+\item de - duration of the human latent period of infection; default = 12
+\item delay_gam - lag from parasites to infectious gametocytes; default = 12.5
+\item dem - extrinsic incubation period in mosquito population model; default = 10
 }
 
 treatment parameters:
-please set treatment parameters with the convenience functions in
-\code{drug_parameters.R}
+please set treatment parameters with the convenience functions
+\code{\link{set_drugs}} and \code{\link{set_clinical_treatment}}
 \itemize{
 \item drug_efficacy - a vector of efficacies for available drugs; default = turned off
 \item drug_rel_c - a vector of relative onward infectiousness values for drugs; default = turned off
@@ -183,22 +190,33 @@ clinically diseased (these values refer to the index in drug_* parameters); defa
 \item clinical_treatment_coverage - a vector of coverage values for each drug; default = NULL, NULL, NULL
 }
 
+MDA, SMC and PMC parameters:
+please set these parameters with the convenience functions
+\code{\link{set_mda}}, \code{\link{set_smc}} and \code{\link{set_pmc}},
+with \code{\link{peak_season_offset}}
+
+bednet, irs and mosquito feeding cycle parameters:
+please set vector control strategies using \code{\link{set_bednets}} and \code{\link{set_spraying}}
+\itemize{
+\item bednets - boolean for if bednets are enabled; default = FALSE
+\item phi_bednets - proportion of bites taken in bed; default = 0.85
+\item k0 - proportion of females bloodfed with no net; default = 0.699
+\item spraying - boolean for if indoor spraying is enabled; default = FALSE
+\item phi_indoors - proportion of bites taken indoors; default = 0.90
+}
+
 PEV parameters:
-please set vaccine strategies with the convenience functions in
-\code{pev_parameters.R:set_pev_epi}
-\code{pev_parameters.R:set_mass_pev}
+please set vaccine strategies with the convenience functions
+\code{\link{set_pev_epi}} and \code{\link{set_mass_pev}}
 \itemize{
 \item pev_doses - the dosing schedule before the vaccine takes effect; default =
 c(0, 1.5 * 30, 3 * 30)
 default = 365
 }
 
-MDA, SMC and PMC parameters:
-please set these parameters with the convenience functions in \code{mda_parameters.R}
-
 TBV parameters:
 please set TBV parameters with the convenience functions in
-\code{vaccine_parameters.R:set_tbv}
+\code{\link{set_tbv}}
 \itemize{
 \item tbv_mt - effect on treated infectiousness; default = 35
 \item tbv_md - effect on diseased infectiousness; default = 46.7
@@ -216,7 +234,7 @@ please set TBV parameters with the convenience functions in
 
 Antimalarial resistance parameters:
 please set antimalarial resistance parameters with the convenience functions in
-\code{antimalarial_resistance.R:set_antimalarial_resistance}
+\code{\link{set_antimalarial_resistance}}
 \itemize{
 \item antimalarial_resistance - boolean for if antimalarial resistance is enabled; default = FALSE
 \item antimalarial_resistance_drug - vector of drugs for which resistance can be parameterised; default = NULL
@@ -232,11 +250,11 @@ please set antimalarial resistance parameters with the convenience functions in
 }
 
 rendering:
-All values are in timesteps and all ranges are inclusive
+All values are in timesteps and all ranges are inclusive.
+Please set rendered age groups using the convenience function
 \itemize{
-\item prevalence_rendering_min_ages - the minimum ages for clinical prevalence
-outputs; default = 730
-\item prevalence_rendering_max_ages - the corresponding max ages; default = 3650
+\item age_group_rendering_min_ages - the minimum ages for population size outputs; default = turned off
+\item age_group_rendering_max_ages - the corresponding max ages; default = turned off
 \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
@@ -245,17 +263,19 @@ outputs (includes asymptomatic microscopy +); default = turned off
 \item severe_incidence_rendering_min_ages - the minimum ages for severe incidence
 outputs; default = turned off
 \item severe_incidence_rendering_max_ages - the corresponding max ages; default = turned off
+\item prevalence_rendering_min_ages - the minimum ages for clinical prevalence
+outputs; default = 730
+\item prevalence_rendering_max_ages - the corresponding max ages; default = 3650
 }
 
 miscellaneous:
 \itemize{
-\item human_population - the initial number of humans to model; default = 100
-\item human_population_timesteps - the timesteps at which the population should
-change; default = 0
 \item mosquito_limit - the maximum number of mosquitoes to allow for in the
 simulation; default = 1.00E+05
 \item individual_mosquitoes - boolean whether adult mosquitoes are modeled
 individually or compartmentally; default = TRUE
+\item human_population_timesteps - the timesteps at which the population should
+change; default = 0
 \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
diff --git a/man/run_simulation.Rd b/man/run_simulation.Rd
index 383b38f6..54ee26ed 100644
--- a/man/run_simulation.Rd
+++ b/man/run_simulation.Rd
@@ -45,16 +45,13 @@ the whole population)
 \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 n_detect: number of humans with an infection detectable by microscopy between an inclusive age range at this timestep. This
+\item n_detect_lm (or pcr): number of humans with an infection detectable by microscopy (or pcr) 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
+\item p_detect_lm (or pcr): the sum of probabilities of detection by microscopy (or pcr) 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 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.
diff --git a/src/Random.cpp b/src/Random.cpp
index 2906ccc7..a4c55a7a 100644
--- a/src/Random.cpp
+++ b/src/Random.cpp
@@ -60,7 +60,7 @@ void Random::prop_sample_bucket(
 
     // all probabilities are the same
     if (heavy == n) {
-        for (auto i = 0; i < size; ++i) {
+        for (size_t i = 0; i < size; ++i) {
             *result = (*rng)((uint64_t)n);
             ++result;
         }
@@ -121,9 +121,9 @@ void Random::prop_sample_bucket(
     }
 
     // sample
-    for (auto i = 0; i < size; ++i) {
+    for (size_t i = 0; i < size; ++i) {
         size_t bucket = (*rng)((uint64_t)n);
-        double acceptance = dqrng::uniform01((*rng)());
+        double acceptance = rng->uniform01();
         *result = (acceptance < dividing_probs[bucket]) ? bucket :
             alternative_index[bucket];
         ++result;
diff --git a/tests/testthat/test-biology.R b/tests/testthat/test-biology.R
index c739f7e7..e429a48f 100644
--- a/tests/testthat/test-biology.R
+++ b/tests/testthat/test-biology.R
@@ -16,10 +16,7 @@ test_that('total_M and EIR functions are consistent with equilibrium EIR', {
     vector_models <- parameterise_mosquito_models(parameters, 1)
     solvers <- parameterise_solvers(vector_models, parameters)
     estimated_eir <- calculate_eir(1, solvers, variables, parameters, 0)
-    age <- get_age(variables$birth$get_values(), 0)
-    psi <- unique_biting_rate(age, parameters)
-    omega <- mean(psi)
-    mean(estimated_eir) / omega / population * 365
+    mean(estimated_eir) / population * 365
   })
 
   expect_equal(
@@ -44,10 +41,7 @@ test_that('total_M and EIR functions are consistent with equilibrium EIR (with h
     vector_models <- parameterise_mosquito_models(parameters, 1)
     solvers <- parameterise_solvers(vector_models, parameters)
     estimated_eir <- calculate_eir(1, solvers, variables, parameters, 0)
-    age <- get_age(variables$birth$get_values(), 0)
-    psi <- unique_biting_rate(age, parameters)
-    omega <- mean(psi)
-    mean(estimated_eir) / omega / population * 365
+    mean(estimated_eir) / population * 365
   })
 
   expect_equal(
diff --git a/tests/testthat/test-render.R b/tests/testthat/test-render.R
index 920f09bb..cc4e9f87 100644
--- a/tests/testthat/test-render.R
+++ b/tests/testthat/test-render.R
@@ -34,7 +34,7 @@ test_that('that default rendering works', {
   mockery::expect_args(
     renderer$render_mock(),
     2,
-    'n_detect_730_3650',
+    'n_detect_lm_730_3650',
     2,
     timestep
   )
@@ -42,11 +42,19 @@ test_that('that default rendering works', {
   mockery::expect_args(
     renderer$render_mock(),
     3,
-    'p_detect_730_3650',
+    'p_detect_lm_730_3650',
     1.5,
     timestep
   )
 
+  mockery::expect_args(
+    renderer$render_mock(),
+    4,
+    'n_detect_pcr_730_3650',
+    3,
+    timestep
+  )
+  
 })
 
 test_that('that default rendering works when no one is in the age range', {
diff --git a/vignettes/Antimalarial_Resistance.Rmd b/vignettes/Antimalarial_Resistance.Rmd
index 5b6acd40..1c59de93 100644
--- a/vignettes/Antimalarial_Resistance.Rmd
+++ b/vignettes/Antimalarial_Resistance.Rmd
@@ -107,9 +107,9 @@ With the outputs from the `run_simulation()` function, we can visualise the effe
 ```{r, fig.align = 'center', eval = TRUE}
 
 # In each timestep, calculate PfPR_2-10 and add it to as a new column for each simulation output:
-sim_out_baseline$pfpr210 = sim_out_baseline$n_detect_730_3650/sim_out_baseline$n_730_3650
-sim_out_clin_treatment$pfpr210 = sim_out_clin_treatment$n_detect_730_3650/sim_out_clin_treatment$n_730_3650
-sim_out_resistance$pfpr210 = sim_out_resistance$n_detect_730_3650/sim_out_resistance$n_730_3650
+sim_out_baseline$pfpr210 = sim_out_baseline$n_detect_lm_730_3650/sim_out_baseline$n_730_3650
+sim_out_clin_treatment$pfpr210 = sim_out_clin_treatment$n_detect_lm_730_3650/sim_out_clin_treatment$n_730_3650
+sim_out_resistance$pfpr210 = sim_out_resistance$n_detect_lm_730_3650/sim_out_resistance$n_730_3650
 
 # Plot the prevalence through time in each of the three simulated scenarios:
 plot.new(); par(mar = c(4, 4, 1, 1), new = TRUE)
@@ -225,7 +225,7 @@ We can visualise the effect of increasing resistance through time by plotting th
 ```{r, fig.align = 'center', eval = TRUE}
 
 # Calculate the prevalence (PfPR210) through time:
-dynamic_resistance_output$pfpr210 <- dynamic_resistance_output$n_detect_730_3650/dynamic_resistance_output$n_730_3650
+dynamic_resistance_output$pfpr210 <- dynamic_resistance_output$n_detect_lm_730_3650/dynamic_resistance_output$n_730_3650
 
 # Open a new plotting window and add a grid:
 plot.new(); par(mar = c(4, 4, 1, 1), new = TRUE)
@@ -365,7 +365,7 @@ plot.new(); par(mar = c(4, 4, 1, 1), new = TRUE)
 
 # Plot malaria prevalence in 2-10 years through time:
 plot(x = simulation_outputs$base$timestep,
-     y = simulation_outputs$base$n_detect_730_3650/simulation_outputs$base$n_730_3650,
+     y = simulation_outputs$base$n_detect_lm_730_3650/simulation_outputs$base$n_730_3650,
      xlab = "Time (days)",
      ylab = expression(paste(italic(Pf),"PR"[2-10])), cex = 0.8,
      ylim = c(0, 1), type = "l", lwd = 2, xaxs = "i", yaxs = "i",
@@ -373,19 +373,19 @@ plot(x = simulation_outputs$base$timestep,
 
 # Add the dynamics for no-intervention simulation
 lines(x = simulation_outputs$treatment$timestep,
-      y = simulation_outputs$treatment$n_detect_730_3650/simulation_outputs$treatment$n_730_3650,
+      y = simulation_outputs$treatment$n_detect_lm_730_3650/simulation_outputs$treatment$n_730_3650,
       col = cols[4])
 
 lines(x = simulation_outputs$etf$timestep,
-      y = simulation_outputs$etf$n_detect_730_3650/simulation_outputs$etf$n_730_3650,
+      y = simulation_outputs$etf$n_detect_lm_730_3650/simulation_outputs$etf$n_730_3650,
       col = cols[5])
 
 lines(x = simulation_outputs$spc$timestep,
-      y = simulation_outputs$spc$n_detect_730_3650/simulation_outputs$spc$n_730_3650,
+      y = simulation_outputs$spc$n_detect_lm_730_3650/simulation_outputs$spc$n_730_3650,
       col = cols[6])
 
 lines(x = simulation_outputs$etf_spc$timestep,
-      y = simulation_outputs$etf_spc$n_detect_730_3650/simulation_outputs$etf_spc$n_730_3650,
+      y = simulation_outputs$etf_spc$n_detect_lm_730_3650/simulation_outputs$etf_spc$n_730_3650,
       col = cols[7])
 
 # Add vlines to indicate when SP-AQ were administered:
@@ -407,4 +407,4 @@ legend(x = 3000, y = 0.99, legend = c("Baseline", "Treatment", "ETF-only", "SPC-
 ```
 
 ## References
-Slater, H.C., Griffin, J.T., Ghani, A.C. and Okell, L.C., 2016. Assessing the potential impact of artemisinin and partner drug resistance in sub-Saharan Africa. Malaria journal, 15(1), pp.1-11.
\ No newline at end of file
+Slater, H.C., Griffin, J.T., Ghani, A.C. and Okell, L.C., 2016. Assessing the potential impact of artemisinin and partner drug resistance in sub-Saharan Africa. Malaria journal, 15(1), pp.1-11.
diff --git a/vignettes/Carrying-capacity.Rmd b/vignettes/Carrying-capacity.Rmd
index 39558355..db38e54a 100644
--- a/vignettes/Carrying-capacity.Rmd
+++ b/vignettes/Carrying-capacity.Rmd
@@ -66,10 +66,10 @@ and can run the simulations and compare the outputs
 ```{r, fig.width = 7, fig.height = 4, out.width='100%'}
 set.seed(123)
 s <- run_simulation(timesteps = timesteps, parameters = p)
-s$pfpr <- s$n_detect_730_3650 / s$n_730_3650
+s$pfpr <- s$n_detect_lm_730_3650 / s$n_730_3650
 set.seed(123)
 s_seasonal <- run_simulation(timesteps = timesteps, parameters = p_seasonal)
-s_seasonal$pfpr <- s_seasonal$n_detect_730_3650 / s_seasonal$n_730_3650
+s_seasonal$pfpr <- s_seasonal$n_detect_lm_730_3650 / s_seasonal$n_730_3650
 
 par(mfrow = c(1, 2))
 plot(s$EIR_gamb ~ s$timestep, t = "l", ylim = c(0, 200), xlab = "Time", ylab = "EIR")
@@ -106,7 +106,7 @@ p_lsm <- p |>
   )
 set.seed(123)
 s_lsm <- run_simulation(timesteps = timesteps, parameters = p_lsm)
-s_lsm$pfpr <- s_lsm$n_detect_730_3650 / s_lsm$n_730_3650
+s_lsm$pfpr <- s_lsm$n_detect_lm_730_3650 / s_lsm$n_730_3650
 
 par(mfrow = c(1, 2))
 plot(s$EIR_gamb ~ s$timestep, t = "l", ylim = c(0, 150), xlab = "Time", ylab = "EIR")
@@ -144,7 +144,7 @@ p_stephensi <- p_stephensi |>
 
 set.seed(123)
 s_stephensi <- run_simulation(timesteps = timesteps, parameters = p_stephensi)
-s_stephensi$pfpr <- s_stephensi$n_detect_730_3650 / s_stephensi$n_730_3650
+s_stephensi$pfpr <- s_stephensi$n_detect_lm_730_3650 / s_stephensi$n_730_3650
 
 par(mfrow = c(1, 2))
 plot(s_stephensi$EIR_gamb ~ s_stephensi$timestep,
@@ -171,7 +171,7 @@ p_flexible <- p |>
 
 set.seed(123)
 s_flexible <- run_simulation(timesteps = timesteps, parameters = p_flexible)
-s_flexible$pfpr <- s_flexible$n_detect_730_3650 / s_flexible$n_730_3650
+s_flexible$pfpr <- s_flexible$n_detect_lm_730_3650 / s_flexible$n_730_3650
 
 par(mfrow = c(1, 2))
 plot(s$EIR_gamb ~ s$timestep, t = "l", ylim = c(0, 150), xlab = "Time", ylab = "EIR")
diff --git a/vignettes/EIRprevmatch.Rmd b/vignettes/EIRprevmatch.Rmd
index f519faa5..7040ed0c 100644
--- a/vignettes/EIRprevmatch.Rmd
+++ b/vignettes/EIRprevmatch.Rmd
@@ -117,7 +117,7 @@ simparams <- set_clinical_treatment(parameters = simparams,
 
 Having established a set of `malariasimulation` parameters, we are now ready to run simulations. In the following code chunk, we'll run the `run_simulation()` function across a range of initial EIR values to generate sufficient points to fit a curve matching *Pf*PR~2-10~ to the initial EIR. For each initial EIR, we first use the `set_equilibrium()` to update the model parameter list with the human and vector population parameter values required to achieve the specified EIR at equilibrium. This updated parameter list is then used to run the simulation.
 
-The `run_simulation()` outputs an EIR per time step, per species, across the entire human population. We first convert these to get the number of infectious bites experienced, on average, by each individual across the final year across all vector species. Next, the average *Pf*PR~2-10~ across the final year of the simulation is calculated by dividing the total number of individuals aged 2-10 by the number (`n_730_3650`) of detectable cases of malaria in individuals aged 2-10 (`n_detect_730_3650`) on each day and calculating the mean of these values. Finally, initial EIR, output EIR, and *Pf*PR~2-10~ are stored in a data frame.
+The `run_simulation()` outputs an EIR per time step, per species, across the entire human population. We first convert these to get the number of infectious bites experienced, on average, by each individual across the final year across all vector species. Next, the average *Pf*PR~2-10~ across the final year of the simulation is calculated by dividing the total number of individuals aged 2-10 by the number (`n_730_3650`) of detectable cases of malaria in individuals aged 2-10 (`n_detect_lm_730_3650`) on each day and calculating the mean of these values. Finally, initial EIR, output EIR, and *Pf*PR~2-10~ are stored in a data frame.
 
 ```{r}
 # Establish a vector of initial EIR values to simulate over and generate matching
@@ -158,7 +158,7 @@ malSim_prev <- lapply(
     mean(
       output[
         output$timestep %in% seq(4 * 365, 5 * 365),
-        'n_detect_730_3650'
+        'n_detect_lm_730_3650'
       ] / output[
         output$timestep %in% seq(4 * 365, 5 * 365),
         'n_730_3650'
@@ -276,7 +276,7 @@ library(cali)
 summary_mean_pfpr_2_10 <- function (x) {
  
  # Calculate the PfPR2-10:
- prev_2_10 <- mean(x$n_detect_730_3650/x$n_730_3650)
+ prev_2_10 <- mean(x$n_detect_lm_730_3650/x$n_730_3650)
  
  # Return the calculated PfPR2-10:
  return(prev_2_10)
@@ -325,8 +325,8 @@ malsim_sim <- run_simulation(timesteps = (simparams_malsim$timesteps),
                              parameters = simparams_malsim)
 
 # Extract the PfPR2-10 values for the cali and malsim simulation outputs:
-cali_pfpr2_10 <- cali_sim$n_detect_730_3650 / cali_sim$n_730_3650
-malsim_pfpr2_10 <- malsim_sim$n_detect_730_3650 / malsim_sim$n_730_3650
+cali_pfpr2_10 <- cali_sim$n_detect_lm_730_3650 / cali_sim$n_730_3650
+malsim_pfpr2_10 <- malsim_sim$n_detect_lm_730_3650 / malsim_sim$n_730_3650
 
 # Store the PfPR2-10 in each time step for the two methods:
 df <- data.frame(timestep = seq(1, length(cali_pfpr2_10)),
diff --git a/vignettes/MDA.Rmd b/vignettes/MDA.Rmd
index d464c533..1d71aeb2 100644
--- a/vignettes/MDA.Rmd
+++ b/vignettes/MDA.Rmd
@@ -143,7 +143,7 @@ plot.new(); par(mar = c(4, 4, 1, 1), new = TRUE)
 
 # Plot malaria prevalence in 2-10 years through time:
 plot(x = mda_output$timestep,
-     y = mda_output$n_detect_730_3650/mda_output$n_730_3650,
+     y = mda_output$n_detect_lm_730_3650/mda_output$n_730_3650,
      xlab = "Time (days)",
      ylab = expression(paste(italic(Pf),"PR"[2-10])), cex = 0.8,
      ylim = c(0, 1), type = "l", lwd = 2, xaxs = "i", yaxs = "i",
@@ -151,7 +151,7 @@ plot(x = mda_output$timestep,
 
 # Add the dynamics for no-intervention simulation
 lines(x = noint_output$timestep,
-     y = noint_output$n_detect_730_3650/noint_output$n_730_3650,
+     y = noint_output$n_detect_lm_730_3650/noint_output$n_730_3650,
      col = cols[4])
 
 # Add vlines to indicate when SP-AQ were administered:
@@ -282,7 +282,7 @@ legend("topleft",
 plot.new(); par(new = TRUE, mar = c(4, 4, 1, 1))
 
 # Plot malaria prevalence in 2-10 years through time:
-plot(x = smc_output$timestep, y = smc_output$n_detect_730_3650/smc_output$n_730_3650,
+plot(x = smc_output$timestep, y = smc_output$n_detect_lm_730_3650/smc_output$n_730_3650,
      xlab = "Time (days)", 
      ylab = expression(paste(italic(Pf),"PR"[2-10])), cex = 0.8,
      ylim = c(0, 1), type = "l", lwd = 2, xaxs = "i", yaxs = "i",
@@ -290,7 +290,7 @@ plot(x = smc_output$timestep, y = smc_output$n_detect_730_3650/smc_output$n_730_
 
 # Add the dynamics for no-intervention simulation
 lines(x = noint_output$timestep,
-     y = noint_output$n_detect_730_3650/noint_output$n_730_3650,
+     y = noint_output$n_detect_lm_730_3650/noint_output$n_730_3650,
      col = cols[4])
 
 # Add lines to indicate SMC events:
@@ -305,4 +305,4 @@ legend("topleft",
        legend = c("SMC", "No-Int"),
        col= c(cols[3:4]), box.col = "white",
        lwd = 1, lty = c(1, 1), cex = 0.8)
-```
\ No newline at end of file
+```
diff --git a/vignettes/Metapopulation.Rmd b/vignettes/Metapopulation.Rmd
index bece7840..880adfb7 100644
--- a/vignettes/Metapopulation.Rmd
+++ b/vignettes/Metapopulation.Rmd
@@ -107,7 +107,7 @@ output3$mix <- 'perfectly-mixed'
 output <- rbind(output1, output2, output3)
 
 # get mean PfPR 2-10 by year
-output$prev2to10 = output$p_detect_730_3650 / output$n_730_3650
+output$prev2to10 = output$p_detect_lm_730_3650 / output$n_730_3650
 output$year = ceiling(output$timestep / 365)
 output$mix = factor(output$mix, levels = c('isolated', 'semi-mixed', 'perfectly-mixed'))
 output <- aggregate(prev2to10 ~ mix + EIR + year, data = output, FUN = mean)
diff --git a/vignettes/Model.Rmd b/vignettes/Model.Rmd
index c09932c2..2eb55b75 100644
--- a/vignettes/Model.Rmd
+++ b/vignettes/Model.Rmd
@@ -154,8 +154,8 @@ The `run_simulation()` function then simulates malaria transmission dynamics and
 -   `iva_mean`: mean acquired immunity to severe infection
 -   `ivm_mean`: mean maternal immunity to severe infection
 -   `n_730_3650`: population size of an age group of interest (where the default is set to 730-3650 days old, or 2-10 years, but which may be adjusted (see [Demography](https://mrc-ide.github.io/malariasimulation/articles/Demography.html) vignette for more details)
--   `n_detect_730_3650`: number with possible detection through microscopy of a given age group
--   `p_detect_730_3650`: the sum of probabilities of detection through microscopy of a given age group
+-   `n_detect_lm_730_3650`: number with possible detection through microscopy of a given age group
+-   `p_detect_lm_730_3650`: the sum of probabilities of detection through microscopy of a given age group
 -   `E_gamb_count`, `L_gamb_count`, `P_gamb_count`, `Sm_gamb_count`, `Pm_gamb_count`, `Im_gamb_count`: species-specific mosquito population sizes in each state (default set to *An. gambiae*)
 -   `total_M_gamb`: species-specific number of adult mosquitoes (default set to *An. gambiae*)
 
@@ -173,7 +173,7 @@ Where **treatments** are specified, `n_treated` will report the number that have
 
 ### Output visualisation
 
-These outputs can then be visualised, such as the population changes in states. Another key output is the prevalence of detectable infections between the ages of 2-10 (*Pf*PR~2-10~), which can be obtained by dividing `n_detect_730_3650` by `n_730_3650`.
+These outputs can then be visualised, such as the population changes in states. Another key output is the prevalence of detectable infections between the ages of 2-10 (*Pf*PR~2-10~), which can be obtained by dividing `n_detect_lm_730_3650` by `n_730_3650`.
 
 ```{r, fig.align = 'center', out.width='100%', fig.asp=0.55, }
 # Define vector of column names to plot
@@ -202,7 +202,7 @@ par(mfrow = c(1,2))
 states_plot(test_sim)
 
 # Calculate Pf PR 2-10
-test_sim$PfPR2_10 <- test_sim$n_detect_730_3650/test_sim$n_730_3650
+test_sim$PfPR2_10 <- test_sim$n_detect_lm_730_3650/test_sim$n_730_3650
 
 # Plot Pf PR 2-10
 plot(x = test_sim$timestep, y = test_sim$PfPR2_10, type = "l",
@@ -226,7 +226,7 @@ par(mfrow = c(1,2))
 states_plot(test_sim_eq)
 
 # Calculate Pf PR 2-10
-test_sim_eq$PfPR2_10 <- test_sim_eq$n_detect_730_3650/test_sim_eq$n_730_3650
+test_sim_eq$PfPR2_10 <- test_sim_eq$n_detect_lm_730_3650/test_sim_eq$n_730_3650
 
 # Plot Pf PR 2-10
 plot(x = test_sim_eq$timestep, y = test_sim_eq$PfPR2_10, type = "l",
diff --git a/vignettes/Parameter_variation.Rmd b/vignettes/Parameter_variation.Rmd
index a0eef828..9b177b8a 100644
--- a/vignettes/Parameter_variation.Rmd
+++ b/vignettes/Parameter_variation.Rmd
@@ -40,7 +40,7 @@ sim <- run_simulation(timesteps = 1000, simparams)
 # plot the default median parameter
 plot(
   sim$timestep,
-  sim$n_detect_730_3650 / sim$n_730_3650,
+  sim$n_detect_lm_730_3650 / sim$n_730_3650,
   t = "l",
   ylim = c(0, 1),
   ylab = "PfPr",
@@ -60,7 +60,7 @@ Keep in mind that `set_parameter_draw` must be called prior to `set_equilibrium`
 # plot the default median parameter
 plot(
   sim$timestep[1:500],
-  sim$n_detect_730_3650[1:500] / sim$n_730_3650[1:500],
+  sim$n_detect_lm_730_3650[1:500] / sim$n_730_3650[1:500],
   t = "l",
   ylim = c(0, 1),
   ylab = "PfPr",
@@ -78,9 +78,9 @@ for (i in 1:7) {
   
   sim <- run_simulation(timesteps = 500, param_draw)
   
-  lines(sim$timestep, sim$n_detect_730_3650 / sim$n_730_3650, col = cols[i])
+  lines(sim$timestep, sim$n_detect_lm_730_3650 / sim$n_730_3650, col = cols[i])
 }
 ```
 
 
-For more information on uncertainty in parameters, please refer to [The US President's Malaria Initiative, Plasmodium falciparum transmission and mortality: A modelling study](https://journals.plos.org/plosmedicine/article?id=10.1371/journal.pmed.1002448), Supplemental material, Section 4. 
\ No newline at end of file
+For more information on uncertainty in parameters, please refer to [The US President's Malaria Initiative, Plasmodium falciparum transmission and mortality: A modelling study](https://journals.plos.org/plosmedicine/article?id=10.1371/journal.pmed.1002448), Supplemental material, Section 4. 
diff --git a/vignettes/SetSpecies.Rmd b/vignettes/SetSpecies.Rmd
index 10b1676d..d6a2e61f 100644
--- a/vignettes/SetSpecies.Rmd
+++ b/vignettes/SetSpecies.Rmd
@@ -40,11 +40,11 @@ We will create a plotting function to visualise the output.
 ```{r}
 # Plotting functions
 plot_prev <- function() {
-  plot(x = output_endophilic$timestep, y = output_endophilic$n_detect_730_3650 / output_endophilic$n_730_3650, 
+  plot(x = output_endophilic$timestep, y = output_endophilic$n_detect_lm_730_3650 / output_endophilic$n_730_3650, 
        type = "l", col = cols[3], lwd = 1,
        xlab = "Time (days)", ylab = expression(paste(italic(Pf),"PR"[2-10])),
        xaxs = "i", yaxs = "i", ylim = c(0,1))
-  lines(x = output_exophilic$timestep, y = output_exophilic$n_detect_730_3650 / output_exophilic$n_730_3650,
+  lines(x = output_exophilic$timestep, y = output_exophilic$n_detect_lm_730_3650 / output_exophilic$n_730_3650,
         col = cols[5], lwd = 1)
   abline(v = sprayingtimesteps, lty = 2, lwd = 2.5, col = "black")
   text(x = sprayingtimesteps + 10, y = 0.9, labels = "Spraying\nint.", adj = 0, cex = 0.8)
diff --git a/vignettes/Treatment.Rmd b/vignettes/Treatment.Rmd
index d325aa74..03efa301 100644
--- a/vignettes/Treatment.Rmd
+++ b/vignettes/Treatment.Rmd
@@ -124,20 +124,20 @@ output <- run_simulation(sim_length, treatment_params)
 
 ### Visualisation
 
-Following simulation of malaria transmission under this treatment regimen, we can now visualise the effect of the regimen on the number of detectable cases through time using the `n_detect_730_3650`.
+Following simulation of malaria transmission under this treatment regimen, we can now visualise the effect of the regimen on the number of detectable cases through time using the `n_detect_lm_730_3650`.
 
 ```{r, fig.align = 'center', out.width='100%'}
 # Plot results
-plot(x = output$timestep, y = output$n_detect_730_3650, type = "l",
+plot(x = output$timestep, y = output$n_detect_lm_730_3650, type = "l",
      xlab = "Days", ylab = "Detectable cases", col = cols[1],
-     ylim = c(min(output$n_detect_730_3650)-1, max(output$n_detect_730_3650)+7),
+     ylim = c(min(output$n_detect_lm_730_3650)-1, max(output$n_detect_lm_730_3650)+7),
      xaxs = "i", yaxs = "i")
 
 # Show treatment times
 abline(v = 300, lty = 2)
-text(x = 310, y = max(output$n_detect_730_3650), labels = "Treatment\nbegins", adj = 0, cex = 0.8)
+text(x = 310, y = max(output$n_detect_lm_730_3650), labels = "Treatment\nbegins", adj = 0, cex = 0.8)
 abline(v = 600, lty = 2)
-text(x = 610, y = max(output$n_detect_730_3650), labels = "Treatment\nends", adj = 0, cex = 0.8)
+text(x = 610, y = max(output$n_detect_lm_730_3650), labels = "Treatment\nends", adj = 0, cex = 0.8)
 
 # Add grid lines
 grid(lty = 2, col = "grey80", nx = 11, ny = 10, lwd = 0.5)
diff --git a/vignettes/Vaccines.Rmd b/vignettes/Vaccines.Rmd
index 02eb9358..4638b1a8 100644
--- a/vignettes/Vaccines.Rmd
+++ b/vignettes/Vaccines.Rmd
@@ -77,12 +77,12 @@ plot_prev <- function(type = "not seasonal"){
   output$time_year <- output$timestep / year
   comparison$time_year <- comparison$timestep / year
 
-  plot(x = output$time_year, y = output$n_detect_730_3650/output$n_730_3650, 
+  plot(x = output$time_year, y = output$n_detect_lm_730_3650/output$n_730_3650, 
        type = "l", col = cols[3], ylim=c(0,1), lwd = 3,
        xlab = "Time (years)", ylab = expression(paste(italic(Pf),"PR"[2-10])),
        xaxs = "i", yaxs = "i")
   grid(lty = 2, col = "grey80", lwd = 0.5)
-  lines(x = comparison$time_year, y = comparison$n_detect_730_3650/comparison$n_730_3650,
+  lines(x = comparison$time_year, y = comparison$n_detect_lm_730_3650/comparison$n_730_3650,
         col = cols[6], lwd = 3)
   abline(v = vaccinetime, col = "black", lty = 2, lwd = 2.5)
   text(x = vaccinetime + 0.05, y = 0.9, labels = "Start of\nvaccination", adj = 0, cex = 0.8)
@@ -388,7 +388,7 @@ plot_prev()
 ```{r, fig.align = 'center', out.width='100%', message=FALSE}
 # Plot doses
 output$month <- ceiling(output$timestep / month)
-doses <- output[, c(2:6, 38)]
+doses <- output[, c(grep(pattern = "n_pev", x = names(output), value = T), "month")]
 doses <- aggregate(cbind(doses[1:5]),
                     by = list(doses$month), 
                     FUN = sum)
diff --git a/vignettes/Variation.Rmd b/vignettes/Variation.Rmd
index 7146d47f..fb563ca0 100644
--- a/vignettes/Variation.Rmd
+++ b/vignettes/Variation.Rmd
@@ -29,7 +29,7 @@ plot_prev <- function(output, ylab = TRUE, ylim = c(0,1)){
   if (ylab == TRUE) {
     ylab = "Prevalence in children aged 2-10 years"
     } else {ylab = ""}
-  plot(x = output$timestep, y = output$n_detect_730_3650 / output$n_730_3650, 
+  plot(x = output$timestep, y = output$n_detect_lm_730_3650 / output$n_730_3650, 
        type = "l", col = cols[3], ylim = ylim,
        xlab = "Time (days)", ylab = ylab, lwd = 1,
        xaxs = "i", yaxs = "i")
@@ -49,7 +49,7 @@ plot_inci <- function(output, ylab = TRUE, ylim){
 
 aggregate_function <- function(df){
   tmp <- aggregate(
-    df$n_detect_730_3650,
+    df$n_detect_lm_730_3650,
     by=list(df$timestep),
     FUN = function(x) {
       c(median = median(x),
@@ -108,7 +108,7 @@ simparams_big <- set_equilibrium(parameters = simparams_big, init_EIR = 50)
 
 ### Simulations
 
-The `n_detect_730_3650` output below shows the total number of individuals in the age group rendered (here, 730-3650 timesteps or 2-10 years) who have microscopy-detectable malaria. Notice that the output is smoother at a higher population.  
+The `n_detect_lm_730_3650` output below shows the total number of individuals in the age group rendered (here, 730-3650 timesteps or 2-10 years) who have microscopy-detectable malaria. Notice that the output is smoother at a higher population.  
 
 Some outcomes will be more sensitive than others to stochastic variation even with the same population size. In the plots below, prevalence is smoother than incidence even at the same population. This is because prevalence is a measure of existing infection, while incidence is recording new cases per timestep.
 
diff --git a/vignettes/VectorControl_Bednets.Rmd b/vignettes/VectorControl_Bednets.Rmd
index ae6a62c4..c1c801a8 100644
--- a/vignettes/VectorControl_Bednets.Rmd
+++ b/vignettes/VectorControl_Bednets.Rmd
@@ -30,11 +30,11 @@ We can create a few plotting functions to visualise the output.
 ```{r}
 # Plotting functions
 plot_prev <- function() {
-  plot(x = output$timestep, y = output$n_detect_730_3650 / output$n_730_3650, 
+  plot(x = output$timestep, y = output$n_detect_lm_730_3650 / output$n_730_3650, 
        type = "l", col = cols[3], lwd = 1,
        xlab = "Time (days)", ylab = expression(paste(italic(Pf),"PR"[2-10])),
        xaxs = "i", yaxs = "i", ylim = c(0, 1))
-  lines(x = output_control$timestep, y = output_control$n_detect_730_3650 / output_control$n_730_3650,
+  lines(x = output_control$timestep, y = output_control$n_detect_lm_730_3650 / output_control$n_730_3650,
         col = cols[5], lwd = 1)
   abline(v = bednetstimesteps, col = "black", lty = 2, lwd = 1)
   text(x = bednetstimesteps + 10, y = 0.95, labels = "Bed net int.", adj = 0, cex = 0.8)
diff --git a/vignettes/VectorControl_IRS.Rmd b/vignettes/VectorControl_IRS.Rmd
index 3cf8ab7e..4dc924a7 100644
--- a/vignettes/VectorControl_IRS.Rmd
+++ b/vignettes/VectorControl_IRS.Rmd
@@ -30,11 +30,11 @@ We will create a few plotting functions to visualise the output.
 ```{r}
 # Plotting functions
 plot_prev <- function() {
-  plot(x = output$timestep, y = output$n_detect_730_3650 / output$n_730_3650, 
+  plot(x = output$timestep, y = output$n_detect_lm_730_3650 / output$n_730_3650, 
        type = "l", col = cols[3], lwd = 1,
        xlab = "Time (days)", ylab = expression(paste(italic(Pf),"PR"[2-10])),
        xaxs = "i", yaxs = "i", ylim = c(0,1))
-  lines(x = output_control$timestep, y = output_control$n_detect_730_3650 / output_control$n_730_3650,
+  lines(x = output_control$timestep, y = output_control$n_detect_lm_730_3650 / output_control$n_730_3650,
         col = cols[5], lwd = 1)
   abline(v = sprayingtimesteps, lty = 2, lwd = 1, col = "black")
   text(x = sprayingtimesteps + 10, y = 0.9, labels = "Spraying\nint.", adj = 0, cex = 0.8)