diff --git a/DESCRIPTION b/DESCRIPTION
index ab68ccb0..c6d25b2d 100644
--- a/DESCRIPTION
+++ b/DESCRIPTION
@@ -70,7 +70,7 @@ Remotes:
   mrc-ide/malariaEquilibrium,
   mrc-ide/individual
 Imports:
-    individual (>= 0.1.16),
+    individual (>= 0.1.17),
     malariaEquilibrium (>= 1.0.1),
     Rcpp,
     statmod,
diff --git a/R/biting_process.R b/R/biting_process.R
index 5cd594da..df494604 100644
--- a/R/biting_process.R
+++ b/R/biting_process.R
@@ -15,6 +15,8 @@
 #' on the other populations
 #' @param mixing_index an index for this population's position in the
 #' lagged_infectivity list (default: 1)
+#' @param infection_outcome competing hazards object for infection rates
+#' @param timestep the current timestep
 #' @noRd
 create_biting_process <- function(
   renderer,
@@ -26,12 +28,12 @@ create_biting_process <- function(
   lagged_infectivity,
   lagged_eir,
   mixing_fn = NULL,
-  mixing_index = 1
+  mixing_index = 1,
+  infection_outcome
   ) {
   function(timestep) {
     # Calculate combined EIR
     age <- get_age(variables$birth$get_values(), timestep)
-    
     bitten_humans <- simulate_bites(
       renderer,
       solvers,
@@ -54,7 +56,8 @@ create_biting_process <- function(
       age,
       parameters,
       timestep,
-      renderer
+      renderer,
+      infection_outcome
     )
   }
 }
diff --git a/R/competing_hazards.R b/R/competing_hazards.R
new file mode 100644
index 00000000..13c522ad
--- /dev/null
+++ b/R/competing_hazards.R
@@ -0,0 +1,95 @@
+## Define classes to resolve competing hazards
+CompetingOutcome <- R6::R6Class(
+  "CompetingOutcome",
+  private = list(
+    targeted_process = NULL
+  ),
+  public = list(
+    initialize = function(targeted_process, size){
+      if (!is.function(targeted_process)){
+        stop("targeted_process must be a function")
+      }
+      if (!is.numeric(size) || size <= 0){
+        stop("size must be positive integer")
+      }
+      private$targeted_process <- targeted_process
+
+      self$target <- individual::Bitset$new(size)
+      self$rates <- NULL
+    },
+    set_rates = function(target, rates){
+      stopifnot(target$size() == length(rates))
+
+      self$target$copy_from(target)
+      self$rates <- rates
+    },
+    execute = function(t, target){
+      private$targeted_process(t, target)
+    },
+    reset = function() {
+      self$target$clear()
+      self$rates <- NULL
+    },
+    target = NULL,
+    rates = NULL
+  )
+)
+
+CompetingHazard <- R6::R6Class(
+  "CompetingHazard",
+  private = list(
+    size = NULL,
+    outcomes = list(),
+    # RNG is passed in because mockery is not able to stub runif
+    # TODO: change when fixed
+    rng = NULL
+  ),
+  public = list(
+    initialize = function(size, outcomes, rng = runif){
+      if (length(outcomes) == 0){
+        stop("At least one outcome must be provided")
+      }
+      if (!all(sapply(outcomes, function(x) inherits(x, "CompetingOutcome")))){
+        stop("All outcomes must be of class CompetingOutcome")
+      }
+      private$size <- size
+      private$outcomes <- outcomes
+      private$rng <- rng
+    },
+    resolve = function(t){
+      candidates <- individual::Bitset$new(private$size)
+      for (o in private$outcomes) {
+        candidates$or(o$target)
+      }
+
+      rates <- matrix(ncol = length(private$outcomes), nrow = candidates$size(), 0)
+      for (i in seq_along(private$outcomes)) {
+        idx <- bitset_index(
+          candidates,
+          private$outcomes[[i]]$target)
+
+        rates[idx, i] <- private$outcomes[[i]]$rates
+      }
+
+      total_rates <- rowSums(rates)
+      probs <- rate_to_prob(total_rates) * (rates / total_rates)
+      probs[is.na(probs)] <- 0
+
+      rng <- private$rng(candidates$size())
+
+      cumulative <- rep(0, candidates$size())
+
+      for (o in seq_along(private$outcomes)) {
+        next_cumulative <- cumulative + probs[,o]
+        selected <- (rng > cumulative) & (rng <= next_cumulative)
+        cumulative <- next_cumulative
+
+        target <- bitset_at(candidates, selected)
+        if (target$size() > 0) {
+          private$outcomes[[o]]$execute(t, target)
+        }
+        private$outcomes[[o]]$reset()
+      }
+    }
+  )
+)
diff --git a/R/disease_progression.R b/R/disease_progression.R
index d9ac11e4..dde88bb2 100644
--- a/R/disease_progression.R
+++ b/R/disease_progression.R
@@ -1,3 +1,70 @@
+#' @title Calculate recovery rates
+#' @description Calculates recovery rates for each individual in the population
+#' for storage in competing hazards object and subsequent resolution
+#'
+#' @param variables the available human variables
+#' @param recovery_outcome competing hazards object for recovery rates
+#' @noRd
+create_recovery_rates_process <- function(
+  variables,
+  recovery_outcome
+) {
+  function(timestep){
+    target <- variables$state$get_index_of("S")$not()
+    recovery_outcome$set_rates(
+      target,
+      variables$recovery_rates$get_values(target))
+  }
+}
+
+
+#' @title Disease progression outcomes (recovery)
+#' @description Following resolution of competing hazards, update state and
+#' infectivity of sampled individuals
+#'
+#' @param timestep the current timestep
+#' @param target the sampled recovering individuals
+#' @param variables the available human variables
+#' @param parameters model parameters
+#' @param renderer competing hazards object for recovery rates
+#' @noRd
+recovery_outcome_process <- function(
+    timestep,
+    target,
+    variables,
+    parameters,
+    renderer
+){
+  
+  update_to_asymptomatic_infection(
+    variables,
+    parameters,
+    timestep,
+    variables$state$get_index_of("D")$and(target)
+  )
+  
+  update_infection(
+    variables$state,
+    "U",
+    variables$infectivity,
+    parameters$cu,
+    variables$recovery_rates,
+    1/parameters$du,
+    variables$state$get_index_of("A")$and(target)
+  )
+  
+  update_infection(
+    variables$state,
+    "S",
+    variables$infectivity,
+    0,
+    variables$recovery_rates,
+    0,
+    variables$state$get_index_of(c("U","Tr"))$and(target)
+  )
+  
+}
+
 #' @title Update the state of an individual as infection events occur
 #' @description Randomly moves individuals towards the later stages of disease
 #' and updates their infectivity
@@ -8,63 +75,17 @@
 #' @param new_infectivity the new infectivity of the progressed individuals
 #' @noRd
 update_infection <- function(
-  state,
-  to_state,
-  infectivity,
-  new_infectivity,
-  to_move
-  ) {
-  state$queue_update(to_state, to_move)
-  infectivity$queue_update(new_infectivity, to_move)
-}
-
-create_progression_process <- function(
     state,
-    from_state,
     to_state,
-    rate,
     infectivity,
-    new_infectivity
+    new_infectivity,
+    recovery_rates,
+    new_recovery_rate,
+    to_move
 ) {
-  function(timestep) {
-    
-    # Retrieve the indices of all individuals in the to_move state:
-    index <- state$get_index_of(from_state)
-    
-    # If the length of rate is greater than 1 (when it's a variable):
-    if (inherits(rate, "DoubleVariable")) {
-      rate <- rate$get_values(index)
-    }
-    
-    # Sample the individuals to be moved into a new Bitset using the transition rate(s):
-    to_move <- index$sample(1/rate)
-    
-    # Update the infection status of those individuals who are moving:
-    update_infection(
-      state,
-      to_state,
-      infectivity,
-      new_infectivity,
-      to_move
-    )
-  }
-}
-
-create_asymptomatic_progression_process <- function(
-  state,
-  rate,
-  variables,
-  parameters
-  ) {
-  function(timestep) {
-    to_move <- state$get_index_of('D')$sample(1/rate)
-    update_to_asymptomatic_infection(
-      variables,
-      parameters,
-      timestep,
-      to_move
-    )
-  }
+  state$queue_update(to_state, to_move)
+  infectivity$queue_update(new_infectivity, to_move)
+  recovery_rates$queue_update(new_recovery_rate, to_move)
 }
 
 #' @title Modelling the progression to asymptomatic disease
@@ -75,11 +96,11 @@ create_asymptomatic_progression_process <- function(
 #' @param parameters model parameters
 #' @noRd
 update_to_asymptomatic_infection <- function(
-  variables,
-  parameters,
-  timestep,
-  to_move
-  ) {
+    variables,
+    parameters,
+    timestep,
+    to_move
+) {
   if (to_move$size() > 0) {
     variables$state$queue_update('A', to_move)
     new_infectivity <- asymptomatic_infectivity(
@@ -94,5 +115,9 @@ update_to_asymptomatic_infection <- function(
       new_infectivity,
       to_move
     )
+    variables$recovery_rates$queue_update(
+      1/parameters$da,
+      to_move
+    )
   }
 }
diff --git a/R/human_infection.R b/R/human_infection.R
index d62b16f0..6f16a7bb 100644
--- a/R/human_infection.R
+++ b/R/human_infection.R
@@ -1,23 +1,26 @@
 #' @title Simulate malaria infection in humans
 #' @description
-#' Updates human states and variables to represent asymptomatic/clinical/severe
-#' and treated malaria; and resulting boosts in immunity
+#' This function ends with the assignment of rates of infection to the competing
+#' hazard resolution object and boosts immunity given infectious bites.
 #' @param variables a list of all of the model variables
 #' @param events a list of all of the model events
 #' @param bitten_humans a bitset of bitten humans
 #' @param age of each human (timesteps)
 #' @param parameters of the model
 #' @param timestep current timestep
+#' @param renderer the model renderer object
+#' @param infection_outcome competing hazards object for infection rates
 #' @noRd
 simulate_infection <- function(
-  variables,
-  events,
-  bitten_humans,
-  age,
-  parameters,
-  timestep,
-  renderer
-  ) {
+    variables,
+    events,
+    bitten_humans,
+    age,
+    parameters,
+    timestep,
+    renderer,
+    infection_outcome
+) {
   if (bitten_humans$size() > 0) {
     boost_immunity(
       variables$ib,
@@ -29,70 +32,18 @@ simulate_infection <- function(
   }
 
   # Calculate Infected
-  infected_humans <- calculate_infections(
+  calculate_infections(
     variables,
     bitten_humans,
     parameters,
     renderer,
-    timestep
-  )
-
-  if (infected_humans$size() > 0) {
-    boost_immunity(
-      variables$ica,
-      infected_humans,
-      variables$last_boosted_ica,
-      timestep,
-      parameters$uc
-    )
-    boost_immunity(
-      variables$id,
-      infected_humans,
-      variables$last_boosted_id,
-      timestep,
-      parameters$ud
-    )
-  }
-
-  clinical_infections <- calculate_clinical_infections(
-    variables,
-    infected_humans,
-    parameters,
-    renderer,
-    timestep
-  )
-
-  update_severe_disease(
     timestep,
-    infected_humans,
-    variables,
-    parameters,
-    renderer
-  )
-
-  treated <- calculate_treated(
-    variables,
-    clinical_infections,
-    parameters,
-    timestep,
-    renderer
-  )
-
-  renderer$render('n_infections', infected_humans$size(), timestep)
-
-  schedule_infections(
-    variables,
-    clinical_infections,
-    treated,
-    infected_humans,
-    parameters,
-    timestep
+    infection_outcome
   )
 }
 
 #' @title Calculate overall infections for bitten humans
-#' @description
-#' Sample infected humans given prophylaxis and vaccination
+#' @description Infection rates are stored in the infection outcome competing hazards object
 #' @param variables a list of all of the model variables
 #' @param bitten_humans bitset of bitten humans
 #' @param parameters model parameters
@@ -100,12 +51,13 @@ simulate_infection <- function(
 #' @param timestep current timestep
 #' @noRd
 calculate_infections <- function(
-  variables,
-  bitten_humans,
-  parameters,
-  renderer,
-  timestep
-  ) {
+    variables,
+    bitten_humans,
+    parameters,
+    renderer,
+    timestep,
+    infection_outcome
+) {
   source_humans <- variables$state$get_index_of(
     c('S', 'A', 'U'))$and(bitten_humans)
 
@@ -155,21 +107,105 @@ calculate_infections <- function(
   }
 
   prob <- b * (1 - prophylaxis) * (1 - vaccine_efficacy)
-  infected <- bitset_at(source_humans, bernoulli_multi_p(prob))
 
-  incidence_renderer(
+  ## probability of incidence must be rendered at each timestep
+  incidence_probability_renderer(
     variables$birth,
     renderer,
-    infected,
     source_humans,
     prob,
+    "inc_",
+    parameters$incidence_min_ages,
+    parameters$incidence_max_ages,
+    timestep
+  )
+  
+  ## capture infection rates to resolve in competing hazards
+  infection_outcome$set_rates(
+    source_humans,
+    prob_to_rate(prob))
+}
+
+#' @title Assigns infections to appropriate human states
+#' @description
+#' Updates human states and variables to represent asymptomatic/clinical/severe
+#' and treated malaria; and resulting boosts in immunity
+#' @param timestep current timestep
+#' @param infected_humans bitset of infected humans
+#' @param variables a list of all of the model variables
+#' @param renderer model render object
+#' @param parameters model parameters
+#' @param prob vector of population probabilities of infection
+#' @noRd
+infection_outcome_process <- function(
+    timestep,
+    infected_humans,
+    variables,
+    renderer,
+    parameters,
+    prob){
+  
+  incidence_renderer(
+    variables$birth,
+    renderer,
+    infected_humans,
     'inc_',
     parameters$incidence_rendering_min_ages,
     parameters$incidence_rendering_max_ages,
     timestep
   )
-
-  infected
+  
+  if (infected_humans$size() > 0) {
+    boost_immunity(
+      variables$ica,
+      infected_humans,
+      variables$last_boosted_ica,
+      timestep,
+      parameters$uc
+    )
+    boost_immunity(
+      variables$id,
+      infected_humans,
+      variables$last_boosted_id,
+      timestep,
+      parameters$ud
+    )
+  }
+  
+  clinical_infections <- calculate_clinical_infections(
+    variables,
+    infected_humans,
+    parameters,
+    renderer,
+    timestep
+  )
+  
+  update_severe_disease(
+    timestep,
+    infected_humans,
+    variables,
+    parameters,
+    renderer
+  )
+  
+  treated <- calculate_treated(
+    variables,
+    clinical_infections,
+    parameters,
+    timestep,
+    renderer
+  )
+  
+  renderer$render('n_infections', infected_humans$size(), timestep)
+  
+  schedule_infections(
+    variables,
+    clinical_infections,
+    treated,
+    infected_humans,
+    parameters,
+    timestep
+  )
 }
 
 #' @title Calculate clinical infections
@@ -182,12 +218,12 @@ calculate_infections <- function(
 #' @param timestep current timestep
 #' @noRd
 calculate_clinical_infections <- function(
-  variables,
-  infections,
-  parameters,
-  renderer,
-  timestep
-  ) {
+    variables,
+    infections,
+    parameters,
+    renderer,
+    timestep
+) {
   ica <- variables$ica$get_values(infections)
   icm <- variables$icm$get_values(infections)
   phi <- clinical_immunity(ica, icm, parameters)
@@ -196,6 +232,14 @@ calculate_clinical_infections <- function(
     variables$birth,
     renderer,
     clinical_infections,
+    'inc_clinical_',
+    parameters$clinical_incidence_rendering_min_ages,
+    parameters$clinical_incidence_rendering_max_ages,
+    timestep
+  )
+  incidence_probability_renderer(
+    variables$birth,
+    renderer,
     infections,
     phi,
     'inc_clinical_',
@@ -216,12 +260,12 @@ calculate_clinical_infections <- function(
 #' @param renderer model outputs
 #' @noRd
 update_severe_disease <- function(
-  timestep,
-  infections,
-  variables,
-  parameters,
-  renderer
-  ) {
+    timestep,
+    infections,
+    variables,
+    parameters,
+    renderer
+) {
   age <- get_age(variables$birth$get_values(infections), timestep)
   iva <- variables$iva$get_values(infections)
   ivm <- variables$ivm$get_values(infections)
@@ -237,6 +281,14 @@ update_severe_disease <- function(
     variables$birth,
     renderer,
     severe_infections,
+    'inc_severe_',
+    parameters$severe_incidence_rendering_min_ages,
+    parameters$severe_incidence_rendering_max_ages,
+    timestep
+  )
+  incidence_probability_renderer(
+    variables$birth,
+    renderer,
     infections,
     theta,
     'inc_severe_',
@@ -355,14 +407,19 @@ calculate_treated <- function(
       successfully_treated
     )
     if(parameters$antimalarial_resistance) {
-      variables$dt$queue_update(
-        parameters$dt,
+      variables$recovery_rates$queue_update(
+        1/parameters$dt,
         non_slow_parasite_clearance_individuals
       )
-      variables$dt$queue_update(
-        dt_slow_parasite_clearance,
+      variables$recovery_rates$queue_update(
+        1/dt_slow_parasite_clearance,
         slow_parasite_clearance_individuals
       )
+    } else {
+      variables$recovery_rates$queue_update(
+        1/parameters$dt,
+        successfully_treated
+      )
     }
   }
   
@@ -380,13 +437,13 @@ calculate_treated <- function(
 #' @param parameters model parameters
 #' @noRd
 schedule_infections <- function(
-  variables,
-  clinical_infections,
-  treated,
-  infections,
-  parameters,
-  timestep
-  ) {
+    variables,
+    clinical_infections,
+    treated,
+    infections,
+    parameters,
+    timestep
+) {
   included <- treated$not(TRUE)
 
   to_infect <- clinical_infections$and(included)
@@ -400,6 +457,8 @@ schedule_infections <- function(
       'D',
       variables$infectivity,
       parameters$cd,
+      variables$recovery_rates,
+      1/parameters$dd,
       to_infect
     )
   }
@@ -418,12 +477,12 @@ schedule_infections <- function(
 # Utility functions
 # =================
 boost_immunity <- function(
-  immunity_variable,
-  exposed_index,
-  last_boosted_variable,
-  timestep,
-  delay
-  ) {
+    immunity_variable,
+    exposed_index,
+    last_boosted_variable,
+    timestep,
+    delay
+) {
   # record who can be boosted
   exposed_index_vector <- exposed_index$to_vector()
   last_boosted <- last_boosted_variable$get_values(exposed_index)
diff --git a/R/model.R b/R/model.R
index 92e4f7c2..a0d77a20 100644
--- a/R/model.R
+++ b/R/model.R
@@ -120,6 +120,7 @@ run_resumable_simulation <- function(
   events <- create_events(parameters)
   initialise_events(events, variables, parameters)
   renderer <- individual::Render$new(timesteps)
+  populate_incidence_rendering_columns(renderer, parameters)
   attach_event_listeners(
     events,
     variables,
@@ -269,6 +270,7 @@ run_metapop_simulation <- function(
   variables <- lapply(parameters, create_variables)
   events <- lapply(parameters, create_events)
   renderer <- lapply(parameters, function(.) individual::Render$new(timesteps))
+  populate_metapopulation_incidence_rendering_columns(renderer, parameters)
   for (i in seq_along(parameters)) {
     # NOTE: forceAndCall is necessary here to make sure i refers to the current
     # iteration
diff --git a/R/mortality_processes.R b/R/mortality_processes.R
index 37792f6c..e67baf0e 100644
--- a/R/mortality_processes.R
+++ b/R/mortality_processes.R
@@ -113,12 +113,8 @@ reset_target <- function(variables, events, target, state, parameters, timestep)
 
     # onwards infectiousness
     variables$infectivity$queue_update(0, target)
+    variables$recovery_rates$queue_update(0, target)
     
-    # treated compartment residence time:
-    if(!is.null(variables$dt)) {
-      variables$dt$queue_update(parameters$dt, target)
-    }
-
     # zeta and zeta group and vector controls survive rebirth
   }
 }
diff --git a/R/parameters.R b/R/parameters.R
index 3319402c..6da0be95 100644
--- a/R/parameters.R
+++ b/R/parameters.R
@@ -445,7 +445,7 @@ get_parameters <- function(overrides = list()) {
     incidence_rendering_min_ages = numeric(0),
     incidence_rendering_max_ages = numeric(0),
     clinical_incidence_rendering_min_ages = numeric(0),
-    clinical_incidence_rendering_max_ages = 5 * 365,
+    clinical_incidence_rendering_max_ages = numeric(0),
     severe_prevalence_rendering_min_ages = numeric(0),
     severe_prevalence_rendering_max_ages = numeric(0),
     severe_incidence_rendering_min_ages = numeric(0),
diff --git a/R/processes.R b/R/processes.R
index 97aa43ef..9f71948d 100644
--- a/R/processes.R
+++ b/R/processes.R
@@ -22,19 +22,19 @@
 #' lagged transmission lists (default: 1)
 #' @noRd
 create_processes <- function(
-  renderer,
-  variables,
-  events,
-  parameters,
-  models,
-  solvers,
-  correlations,
-  lagged_eir,
-  lagged_infectivity,
-  timesteps,
-  mixing_fn = NULL,
-  mixing_index = 1
-  ) {
+    renderer,
+    variables,
+    events,
+    parameters,
+    models,
+    solvers,
+    correlations,
+    lagged_eir,
+    lagged_infectivity,
+    timesteps,
+    mixing_fn = NULL,
+    mixing_index = 1
+) {
   
   # ========
   # Immunity
@@ -69,11 +69,31 @@ create_processes <- function(
       )
     )
   }
-
+  
+  # =====================================================
+  # Competing Hazard Outcomes (Infections and Recoveries)
+  # =====================================================
+  
+  infection_outcome <- CompetingOutcome$new(
+    targeted_process = function(timestep, target){
+      infection_outcome_process(timestep, target, 
+                                variables, renderer, parameters, 
+                                prob = rate_to_prob(infection_outcome$rates))
+    },
+    size = parameters$human_population
+  )
+  
+  recovery_outcome <- CompetingOutcome$new(
+    targeted_process = function(timestep, target){
+      recovery_outcome_process(timestep, target, variables, parameters, renderer)
+    },
+    size = parameters$human_population
+  )
+  
   # ==============================
   # Biting and mortality processes
   # ==============================
-  # schedule infections for humans and set last_boosted_*
+  # simulate bites, calculates infection rates for bitten humans and set last_boosted_*
   # move mosquitoes into incubating state
   # kill mosquitoes caught in vector control
   processes <- c(
@@ -88,60 +108,29 @@ create_processes <- function(
       lagged_infectivity,
       lagged_eir,
       mixing_fn,
-      mixing_index
-    ),
-    asymptotic_progression_process = create_asymptomatic_progression_process(
-      variables$state,
-      parameters$dd,
-      variables,
-      parameters
-    ),
-    progression_process = create_progression_process(
-      variables$state,
-      'A',
-      'U',
-      parameters$da,
-      variables$infectivity,
-      parameters$cu
-    ),
-    progression_process = create_progression_process(
-      variables$state,
-      'U',
-      'S',
-      parameters$du,
-      variables$infectivity,
-      0
+      mixing_index,
+      infection_outcome
     )
   )
   
-  # =======================
-  # Antimalarial Resistance
-  # =======================
-  # Add an a new process which governs the transition from Tr to S when
-  # antimalarial resistance is simulated. The rate of transition switches
-  # from a parameter to a variable when antimalarial resistance == TRUE.
-  
-  # Assign the dt input to a separate object with the default single parameter value:
-  dt_input <- parameters$dt
-  
-  # If antimalarial resistance is switched on, assign dt variable values to the 
-  if(parameters$antimalarial_resistance) {
-    dt_input <- variables$dt
-  }
+  # ===================
+  # Disease Progression
+  # ===================
   
-  # Create the progression process for Tr --> S specifying dt_input as the rate:
   processes <- c(
     processes,
-    progression_process = create_progression_process(
-      variables$state,
-      'Tr',
-      'S',
-      dt_input,
-      variables$infectivity,
-      0
-    )
+    progression_process = create_recovery_rates_process(
+      variables,
+      recovery_outcome
+    ),
+    
+    # Resolve competing hazards of infection with disease progression
+    hazard_resolution_process = CompetingHazard$new(
+      outcomes = list(infection_outcome, recovery_outcome),
+      size = parameters$human_population
+    )$resolve
   )
-
+  
   # ===============
   # ODE integration
   # ===============
@@ -291,7 +280,11 @@ create_processes <- function(
     )
   }
 
+  # ======================
   # Mortality step
+  # ======================
+  # Mortality is not resolved as a competing hazard
+  
   processes <- c(
     processes,
     mortality_process = create_mortality_process(
diff --git a/R/render.R b/R/render.R
index ac1e5656..5ba9d79d 100644
--- a/R/render.R
+++ b/R/render.R
@@ -18,12 +18,12 @@ in_age_range <- function(birth, timestep, lower, upper) {
 #' 
 #' @noRd
 create_prevalence_renderer <- function(
-  state,
-  birth,
-  immunity,
-  parameters,
-  renderer
-  ) {
+    state,
+    birth,
+    immunity,
+    parameters,
+    renderer
+) {
   function(timestep) {
     asymptomatic <- state$get_index_of('A')
     prob <- probability_of_detection(
@@ -41,11 +41,6 @@ create_prevalence_renderer <- function(
       lower <- parameters$prevalence_rendering_min_ages[[i]]
       upper <- parameters$prevalence_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_detect_lm_', lower, '_', upper),
         in_age$copy()$and(detected)$size(),
@@ -69,13 +64,11 @@ create_prevalence_renderer <- function(
 
 #' @title Render incidence statistics
 #' 
-#' @description renders incidence (new for this timestep) for indivduals
+#' @description renders incidence (new for this timestep) for individuals
 #' 
 #' @param birth variable for birth of the individual
 #' @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
@@ -83,27 +76,54 @@ create_prevalence_renderer <- function(
 #' 
 #' @noRd
 incidence_renderer <- function(
-  birth,
-  renderer,
-  target,
-  source_pop,
-  prob,
-  prefix,
-  lowers,
-  uppers,
-  timestep
-  ) {
+    birth,
+    renderer,
+    target,
+    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 probability of incidence statistics
+#' 
+#' @description renders probability of incidence (new for this timestep) for individuals
+#' 
+#' @param birth variable for birth of the individual
+#' @param renderer object for model outputs
+#' @param source_pop the population which is sampled for infection (bitten SAU for incidence, infecte for clinical/severe)
+#' @param prob probability of infection
+#' @param prefix for model outputs
+#' @param lowers age bounds
+#' @param uppers age bounds
+#' @param timestep current target
+#' 
+#' @noRd
+incidence_probability_renderer <- function(
+    birth,
+    renderer,
+    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('p_', prefix, lower, '_', upper),
       sum(prob[bitset_index(source_pop, in_age)]),
@@ -113,9 +133,9 @@ incidence_renderer <- function(
 }
 
 create_variable_mean_renderer_process <- function(
-  renderer,
-  names,
-  variables
+    renderer,
+    names,
+    variables
 ) {
   function(timestep) {
     for (i in seq_along(variables)) {
@@ -153,12 +173,12 @@ create_age_variable_mean_renderer_process <- function(
 }
 
 create_vector_count_renderer_individual <- function(
-  mosquito_state,
-  species,
-  state,
-  renderer,
-  parameters
-  ) {
+    mosquito_state,
+    species,
+    state,
+    renderer,
+    parameters
+) {
   function(timestep) {
     adult <- mosquito_state$get_index_of('NonExistent')$not(TRUE)
     for (i in seq_along(parameters$species)) {
@@ -192,10 +212,23 @@ create_age_group_renderer <- function(
     parameters,
     renderer
 ) {
+  
+  age_ranges <- rbind(
+    cbind(parameters$prevalence_rendering_min_ages, parameters$prevalence_rendering_max_ages),
+    cbind(parameters$incidence_rendering_min_ages, parameters$incidence_rendering_max_ages),
+    cbind(parameters$clinical_incidence_rendering_min_ages, parameters$clinical_incidence_rendering_max_ages),
+    cbind(parameters$severe_incidence_rendering_min_ages, parameters$severe_incidence_rendering_max_ages),
+    cbind(parameters$age_group_rendering_min_ages, parameters$age_group_rendering_max_ages)
+  )
+  
+  unique_age_combinations <- as.data.frame(unique(age_ranges))
+  ordered_unique_age_combinations <- unique_age_combinations[order(unique_age_combinations$V1, unique_age_combinations$V2), ]
+  
   function(timestep) {
-    for (i in seq_along(parameters$age_group_rendering_min_ages)) {
-      lower <- parameters$age_group_rendering_min_ages[[i]]
-      upper <- parameters$age_group_rendering_max_ages[[i]]
+    
+    for (i in seq_along(ordered_unique_age_combinations$V1)) {
+      lower <- ordered_unique_age_combinations$V1[[i]]
+      upper <- ordered_unique_age_combinations$V2[[i]]
       in_age <- in_age_range(birth, timestep, lower, upper)
       renderer$render(
         paste0('n_age_', lower, '_', upper),
@@ -205,3 +238,60 @@ create_age_group_renderer <- function(
     }
   }
 }
+
+
+populate_incidence_rendering_columns <- function(renderer, parameters){
+  
+  # infections must render in all simulations 
+  renderer$set_default('n_infections', 0)
+  
+  # treatment associated only renders when drugs are used
+  if(sum(unlist(parameters$clinical_treatment_coverages))>0){
+    renderer$set_default('ft', 0)
+    renderer$set_default('n_treated', 0)
+    renderer$set_default('n_drug_efficacy_failures', 0)
+    renderer$set_default('n_successfully_treated', 0)
+  }
+  
+  # ETC, SPC only render when antimalarial resistance is on
+  if(parameters$antimalarial_resistance){
+    renderer$set_default('n_early_treatment_failure', 0)
+    renderer$set_default('n_slow_parasite_clearance', 0)
+  }
+  
+  if(length(parameters$incidence_rendering_min_ages)>0){
+    render_initial_incidence(renderer,
+                             parameters$incidence_rendering_min_ages,
+                             parameters$incidence_rendering_max_ages,
+                             "inc")
+  }
+  
+  if(length(parameters$clinical_incidence_rendering_min_ages)>0){
+    render_initial_incidence(renderer,
+                             parameters$clinical_incidence_rendering_min_ages,
+                             parameters$clinical_incidence_rendering_max_ages,
+                             "inc_clinical")
+  }
+
+  if(length(parameters$severe_incidence_rendering_min_ages)>0){
+    render_initial_incidence(renderer,
+                             parameters$severe_incidence_rendering_min_ages,
+                             parameters$severe_incidence_rendering_max_ages,
+                             "inc_severe")
+  }
+  
+}
+
+render_initial_incidence <- function(renderer, lower_vals, upper_vals, inc_name){
+  for (i in seq_along(lower_vals)){
+    renderer$set_default(paste0('n_', inc_name, "_", lower_vals[i], '_', upper_vals[i]), 0)
+    renderer$set_default(paste0('p_', inc_name, "_", lower_vals[i], '_', upper_vals[i]), 0)
+  }
+}
+  
+
+populate_metapopulation_incidence_rendering_columns <- function(renderer, parameters){
+  for(i in length(parameters)){
+    populate_incidence_rendering_columns(renderer[[i]], parameters[[i]])
+  }
+}
diff --git a/R/utils.R b/R/utils.R
index 03eba932..0c98b099 100644
--- a/R/utils.R
+++ b/R/utils.R
@@ -104,3 +104,17 @@ RandomState <- R6::R6Class(
     }
   )
 )
+
+#'@title Convert probability to a rate
+#'@param prob probability
+#'@noRd
+prob_to_rate <- function(prob){
+  -log(1 - prob)
+}
+
+#'@title Convert rate to a probability
+#'@param rate rate
+#'@noRd
+rate_to_prob <- function(rate){
+  1 - exp(-rate)
+}
diff --git a/R/variables.R b/R/variables.R
index 86a13185..5a6dc749 100644
--- a/R/variables.R
+++ b/R/variables.R
@@ -178,6 +178,7 @@ create_variables <- function(parameters) {
   diseased <- state$get_index_of('D')$to_vector()
   asymptomatic <- state$get_index_of('A')$to_vector()
   subpatent <- state$get_index_of('U')$to_vector()
+  treated <- state$get_index_of('Tr')$to_vector()
 
   # Set the initial infectivity values for each individual
   infectivity_values[diseased] <- parameters$cd
@@ -191,6 +192,16 @@ create_variables <- function(parameters) {
   # Initialise the infectivity variable
   infectivity <- individual::DoubleVariable$new(infectivity_values)
 
+  # Set recovery rate for each individual
+  recovery_values <- rep(0, get_human_population(parameters, 0))
+  recovery_values[diseased] <- 1/parameters$dd
+  recovery_values[asymptomatic] <- 1/parameters$da
+  recovery_values[subpatent] <- 1/parameters$du
+  recovery_values[treated] <- 1/parameters$dt
+
+  # Initialise the recovery rate variable
+  recovery_rates <- individual::DoubleVariable$new(recovery_values)
+  
   drug <- individual::IntegerVariable$new(rep(0, size))
   drug_time <- individual::IntegerVariable$new(rep(-1, size))
 
@@ -220,6 +231,7 @@ create_variables <- function(parameters) {
     zeta = zeta,
     zeta_group = zeta_group,
     infectivity = infectivity,
+    recovery_rates = recovery_rates,
     drug = drug,
     drug_time = drug_time,
     last_pev_timestep = last_pev_timestep,
@@ -230,14 +242,6 @@ create_variables <- function(parameters) {
     spray_time = spray_time
   )
   
-  # Add variables for antimalarial resistance state residency times (dt)
-  if(parameters$antimalarial_resistance) {
-    dt <- individual::DoubleVariable$new(rep(parameters$dt, size))
-    variables <- c(
-      variables,
-      dt = dt
-    )
-  }
 
   # Add variables for individual mosquitoes
   if (parameters$individual_mosquitoes) {
diff --git a/tests/testthat/test-biting-integration.R b/tests/testthat/test-biting-integration.R
index 992ffe42..fbfd1fd4 100644
--- a/tests/testthat/test-biting-integration.R
+++ b/tests/testthat/test-biting-integration.R
@@ -14,6 +14,13 @@ test_that('biting_process integrates mosquito effects and human infection', {
   models <- parameterise_mosquito_models(parameters, timestep)
   solvers <- parameterise_solvers(models, parameters)
 
+  infection_outcome <- CompetingOutcome$new(
+    targeted_process = function(timestep, target){
+      infection_process_resolved_hazard(timestep, target, variables, renderer, parameters)
+    },
+    size = parameters$human_population
+    )
+
   biting_process <- create_biting_process(
     renderer,
     solvers,
@@ -22,9 +29,10 @@ test_that('biting_process integrates mosquito effects and human infection', {
     events,
     parameters,
     lagged_foim,
-    lagged_eir
+    lagged_eir,
+    infection_outcome=infection_outcome
   )
-
+  
   bitten <- individual::Bitset$new(parameters$human_population)
   bites_mock <- mockery::mock(bitten)
   infection_mock <- mockery::mock()
@@ -32,7 +40,7 @@ test_that('biting_process integrates mosquito effects and human infection', {
   mockery::stub(biting_process, 'simulate_bites', bites_mock)
   mockery::stub(biting_process, 'simulate_infection', infection_mock)
   biting_process(timestep)
-
+  
   mockery::expect_args(
     bites_mock,
     1,
@@ -59,7 +67,8 @@ test_that('biting_process integrates mosquito effects and human infection', {
     age,
     parameters,
     timestep,
-    renderer
+    renderer,
+    infection_outcome
   )
 })
 
diff --git a/tests/testthat/test-competing-hazards.R b/tests/testthat/test-competing-hazards.R
new file mode 100644
index 00000000..f7ada3fc
--- /dev/null
+++ b/tests/testthat/test-competing-hazards.R
@@ -0,0 +1,176 @@
+# ==============================
+# Test the CompetingHazard class
+# ==============================
+
+
+test_that("hazard resolves two disjoint outcomes", {
+  size <- 4
+  population <- individual::Bitset$new(size)$not()
+
+  outcome_1_process <- mockery::mock()
+  outcome_1 <- CompetingOutcome$new(
+    targeted_process = outcome_1_process,
+    size = size
+  )
+  outcome_2_process <- mockery::mock()
+  outcome_2 <- CompetingOutcome$new(
+    targeted_process = outcome_2_process,
+    size = size
+  )
+
+  hazard <- CompetingHazard$new(
+    size = size,
+    outcomes = list(outcome_1, outcome_2),
+    rng = mockery::mock(c(.05, .3, .2, .5))
+  )
+  
+  outcome_1$set_rates(population, c(10, 0, 10, 0))
+  outcome_2$set_rates(population, c(0, 10, 0, 10))
+  
+  hazard$resolve(0)
+   
+  mockery::expect_args(outcome_1_process, 1, 0,
+                       individual::Bitset$new(size)$insert(c(1, 3)))
+  mockery::expect_args(outcome_2_process, 1, 0,
+                       individual::Bitset$new(size)$insert(c(2, 4)))
+})
+
+test_that("hazard resolves two competing outcomes", {
+  size <- 4
+  population <- individual::Bitset$new(size)$not()
+
+  outcome_1_process <- mockery::mock()
+  outcome_1 <- CompetingOutcome$new(
+    targeted_process = outcome_1_process,
+    size = size
+  )
+  outcome_2_process <- mockery::mock()
+  outcome_2 <- CompetingOutcome$new(
+    targeted_process = outcome_2_process,
+    size = size
+  )
+
+  hazard <- CompetingHazard$new(
+    size = size,
+    outcomes = list(outcome_1, outcome_2),
+    rng = mockery::mock(c(.7, .3, .2, .6))
+  )
+  
+  outcome_1$set_rates(population, c(5, 5, 5, 5))
+  outcome_2$set_rates(population, c(5, 5, 5, 5))
+  
+  hazard$resolve(0)
+   
+  mockery::expect_args(outcome_1_process, 1, 0,
+                       individual::Bitset$new(size)$insert(c(2, 3)))
+  mockery::expect_args(outcome_2_process, 1, 0,
+                       individual::Bitset$new(size)$insert(c(1, 4)))
+})
+
+test_that("hazard may resolve to neither outcome", {
+  size <- 4
+  population <- individual::Bitset$new(size)$not()
+
+  outcome_1_process <- mockery::mock()
+  outcome_1 <- CompetingOutcome$new(
+    targeted_process = outcome_1_process,
+    size = size
+  )
+  outcome_2_process <- mockery::mock()
+  outcome_2 <- CompetingOutcome$new(
+    targeted_process = outcome_2_process,
+    size = size
+  )
+
+  hazard <- CompetingHazard$new(
+    size = size,
+    outcomes = list(outcome_1, outcome_2),
+    rng = mockery::mock(c(.8, .4, .2, .6))
+  )
+  
+  outcome_1$set_rates(population, prob_to_rate(rep(0.5, size)))
+  outcome_2$set_rates(population, prob_to_rate(rep(0.5, size)))
+  
+  hazard$resolve(0)
+   
+  mockery::expect_args(outcome_1_process, 1, 0,
+                       individual::Bitset$new(size)$insert(c(3)))
+  mockery::expect_args(outcome_2_process, 1, 0,
+                       individual::Bitset$new(size)$insert(c(2, 4)))
+})
+
+test_that("outcomes can define a partial set of rates", {
+  size <- 4
+  population <- individual::Bitset$new(size)$not()
+
+  outcome_1_process <- mockery::mock()
+  outcome_1 <- CompetingOutcome$new(
+    targeted_process = outcome_1_process,
+    size = size
+  )
+  outcome_2_process <- mockery::mock()
+  outcome_2 <- CompetingOutcome$new(
+    targeted_process = outcome_2_process,
+    size = size
+  )
+
+  hazard <- CompetingHazard$new(
+    size = size,
+    outcomes = list(outcome_1, outcome_2),
+
+    # Only individuals 1, 3 and 4 get sampled
+    rng = mockery::mock(c(.2, .3, .6))
+  )
+  
+  outcome_1$set_rates(individual::Bitset$new(size)$insert(c(1,3)), c(5, 5))
+  outcome_2$set_rates(individual::Bitset$new(size)$insert(c(1,4)), c(5, 5))
+  
+  hazard$resolve(0)
+   
+  mockery::expect_args(outcome_1_process, 1, 0,
+                       individual::Bitset$new(size)$insert(c(1, 3)))
+  mockery::expect_args(outcome_2_process, 1, 0,
+                       individual::Bitset$new(size)$insert(c(4))
+  )
+})
+
+test_that("hazard resolves three competing outcomes", {
+  size <- 4
+  population <- individual::Bitset$new(size)$not()
+
+  outcome_1_process <- mockery::mock()
+  outcome_1 <- CompetingOutcome$new(
+    targeted_process = outcome_1_process,
+    size = size
+  )
+  outcome_2_process <- mockery::mock()
+  outcome_2 <- CompetingOutcome$new(
+    targeted_process = outcome_2_process,
+    size = size
+  )
+  outcome_3_process <- mockery::mock()
+  outcome_3 <- CompetingOutcome$new(
+    targeted_process = outcome_3_process,
+    size = size
+  )
+
+  hazard <- CompetingHazard$new(
+    size = size,
+    outcomes = list(outcome_1, outcome_2, outcome_3),
+    rng = mockery::mock(c(.1, .5, .8, .8))
+  )
+  
+  outcome_1$set_rates(population, c(5, 5, 5, 5))
+  outcome_2$set_rates(population, c(5, 5, 5, 5))
+  outcome_3$set_rates(population, c(5, 5, 5, 5))
+  
+  hazard$resolve(0)
+   
+  mockery::expect_args(outcome_1_process, 1, 0,
+                       individual::Bitset$new(size)$insert(c(1)))
+  mockery::expect_args(outcome_2_process, 1, 0,
+                       individual::Bitset$new(size)$insert(c(2)))
+  mockery::expect_args(outcome_3_process, 1, 0,
+                       individual::Bitset$new(size)$insert(c(3, 4)))
+})
+
diff --git a/tests/testthat/test-infection-integration.R b/tests/testthat/test-infection-integration.R
index 83d7def9..17972b6a 100644
--- a/tests/testthat/test-infection-integration.R
+++ b/tests/testthat/test-infection-integration.R
@@ -20,22 +20,17 @@ test_that('simulate_infection integrates different types of infection and schedu
   boost_immunity_mock <- mockery::mock()
   infected <- individual::Bitset$new(population)$insert(c(1, 3, 5))
   infection_mock <- mockery::mock(infected)
-  clinical <- individual::Bitset$new(population)$insert(c(1, 3))
-  clinical_infection_mock <- mockery::mock(clinical)
-  severe <- individual::Bitset$new(population)$insert(1)
-  severe_infection_mock <- mockery::mock(severe)
-  treated <- individual::Bitset$new(population)$insert(3)
-  treated_mock <- mockery::mock(treated)
-  schedule_mock <- mockery::mock()
+  
+  infection_outcome <- CompetingOutcome$new(
+    targeted_process = function(timestep, target){
+      infection_outcome_process(timestep, target, variables, renderer, parameters)
+    },
+    size = parameters$human_population
+  )
   
   mockery::stub(simulate_infection, 'boost_immunity', boost_immunity_mock)
   mockery::stub(simulate_infection, 'calculate_infections', infection_mock)
-  mockery::stub(simulate_infection, 'calculate_clinical_infections', clinical_infection_mock)
-  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,
@@ -43,7 +38,8 @@ test_that('simulate_infection integrates different types of infection and schedu
     age,
     parameters,
     timestep,
-    renderer
+    renderer,
+    infection_outcome
   )
   
   mockery::expect_args(
@@ -63,7 +59,68 @@ test_that('simulate_infection integrates different types of infection and schedu
     bitten,
     parameters,
     renderer,
-    timestep
+    timestep,
+    infection_outcome
+  )
+})
+
+
+test_that('simulate_infection integrates different types of infection and scheduling', {
+  population <- 8
+  timestep <- 5
+  parameters <- get_parameters(list(
+    human_population = population
+  ))
+  events <- create_events(parameters)
+  renderer <- mock_render(timestep)
+  
+  age <- c(20, 24, 5, 39, 20, 24, 5, 39) * 365
+  immunity <- c(.2, .3, .5, .9, .2, .3, .5, .9)
+  # asymptomatics <- mockery::mock()
+  variables <- list(
+    ib = individual::DoubleVariable$new(immunity),
+    id = individual::DoubleVariable$new(immunity),
+    # state = list(get_index_of = mockery::mock(asymptomatics, cycle = T))
+    state = individual::CategoricalVariable$new(categories = c("S","A","U","D","Tr"), initial_values = rep("S", population))#list(get_index_of = mockery::mock(asymptomatics, cycle = T))
+  )
+  prob <- rep(0.5,population)
+  
+  source_humans <- individual::Bitset$new(population)$insert(c(1, 2, 3, 5))
+  infected <- individual::Bitset$new(population)$insert(c(1, 3, 5))
+  clinical <- individual::Bitset$new(population)$insert(c(1, 3))
+  clinical_infection_mock <- mockery::mock(clinical)
+  boost_immunity_mock <- mockery::mock()
+  severe <- individual::Bitset$new(population)$insert(1)
+  severe_infection_mock <- mockery::mock(severe)
+  treated <- individual::Bitset$new(population)$insert(3)
+  treated_mock <- mockery::mock(treated)
+  schedule_mock <- mockery::mock()
+  
+  
+  mockery::stub(infection_outcome_process, 'incidence_renderer', mockery::mock())
+  mockery::stub(infection_outcome_process, 'boost_immunity', boost_immunity_mock)
+  mockery::stub(infection_outcome_process, 'calculate_clinical_infections', clinical_infection_mock)
+  mockery::stub(infection_outcome_process, 'update_severe_disease', severe_infection_mock)
+  mockery::stub(infection_outcome_process, 'calculate_treated', treated_mock)
+  mockery::stub(infection_outcome_process, 'schedule_infections', schedule_mock)
+  
+  infection_outcome_process(
+    timestep,
+    infected,
+    variables,
+    renderer,
+    parameters,
+    prob)
+  
+  
+  mockery::expect_args(
+    boost_immunity_mock,
+    1,
+    variables$ica,
+    infected,
+    variables$last_boosted_ica,
+    5,
+    parameters$uc
   )
   
   mockery::expect_args(
@@ -161,15 +218,23 @@ test_that('calculate_infections works various combinations of drug and vaccinati
   
   bitten_humans <- individual::Bitset$new(4)$insert(c(1, 2, 3, 4))
   
+  infection_outcome <- CompetingOutcome$new(
+    targeted_process = function(timestep, target){
+      infection_outcome_process(timestep, target, variables, renderer, parameters)
+    },
+    size = 4
+  )
+  
   infections <- calculate_infections(
     variables,
     bitten_humans, 
     parameters,
     mock_render(timestep),
-    timestep
+    timestep,
+    infection_outcome
   )
   
-  expect_equal(infections$to_vector(), 3)
+  expect_equal(sum(infections!=0), 3)
   
   mockery::expect_args(immunity_mock, 1, c(.3, .5, .9), parameters)
   mockery::expect_args(
@@ -198,12 +263,6 @@ test_that('calculate_infections works various combinations of drug and vaccinati
     c(rtss_profile$beta, rtss_booster_profile$beta),
     c(rtss_profile$alpha, rtss_booster_profile$alpha)
   )
-  mockery::expect_args(
-    bernoulli_mock,
-    1,
-    c(.2 * .8 * .8, .3 * .7, .4)
-  )
-  
 })
 
 test_that('calculate_clinical_infections correctly samples clinically infected', {
@@ -261,6 +320,7 @@ test_that('calculate_treated correctly samples treated and updates the drug stat
   variables <- list(
     state = list(queue_update = mockery::mock()),
     infectivity = list(queue_update = mockery::mock()),
+    recovery_rates = list(queue_update = mockery::mock()),
     drug = list(queue_update = mockery::mock()),
     drug_time = list(queue_update = mockery::mock())
   )
@@ -350,9 +410,9 @@ test_that('calculate_treated correctly samples treated and updates the drug stat
   variables <- list(
     state = list(queue_update = mockery::mock()),
     infectivity = list(queue_update = mockery::mock()),
+    recovery_rates = list(queue_update = mockery::mock()),
     drug = list(queue_update = mockery::mock()),
-    drug_time = list(queue_update = mockery::mock()),
-    dt = list(queue_update = mockery::mock())
+    drug_time = list(queue_update = mockery::mock())
   )
   renderer <- individual::Render$new(timesteps = timestep)
   
@@ -427,8 +487,8 @@ test_that('calculate_treated correctly samples treated and updates the drug stat
   )
   expect_bitset_update(variables$drug$queue_update, c(2, 1, 1, 1, 2, 2, 2), c(1, 2, 3, 4, 5, 6, 7))
   expect_bitset_update(variables$drug_time$queue_update, 5, c(1, 2, 3, 4, 5, 6, 7))
-  expect_bitset_update(variables$dt$queue_update, 5, c(2, 3, 4, 5, 6, 7), 1)
-  expect_bitset_update(variables$dt$queue_update, 15, c(1), 2)
+  expect_bitset_update(variables$recovery_rates$queue_update, 1/5, c(2, 3, 4, 5, 6, 7), 1)
+  expect_bitset_update(variables$recovery_rates$queue_update, 1/15, c(1), 2)
 })
 
 test_that('calculate_treated correctly samples treated and updates the drug state when resistance not set for all drugs', {
@@ -463,9 +523,9 @@ test_that('calculate_treated correctly samples treated and updates the drug stat
   variables <- list(
     state = list(queue_update = mockery::mock()),
     infectivity = list(queue_update = mockery::mock()),
+    recovery_rates = list(queue_update = mockery::mock()),
     drug = list(queue_update = mockery::mock()),
-    drug_time = list(queue_update = mockery::mock()),
-    dt = list(queue_update = mockery::mock())
+    drug_time = list(queue_update = mockery::mock())
   )
   
   # Create a Bitset of individuals seeking treatment individuals:
@@ -569,15 +629,15 @@ test_that('calculate_treated correctly samples treated and updates the drug stat
   
   # Check that update queued for dt for the non-slow parasite clearance individuals is correct:
   expect_bitset_update(
-    variables$dt$queue_update,
-    parameters$dt,
+    variables$recovery_rates$queue_update,
+    1/parameters$dt,
     c(2, 3, 4, 5, 6, 7), 
     1)
   
   # Check that update queued for dt for the slow parasite clearance individuals is correct:
   expect_bitset_update(
-    variables$dt$queue_update,
-    unlist(parameters$dt_slow_parasite_clearance),
+    variables$recovery_rates$queue_update,
+    1/unlist(parameters$dt_slow_parasite_clearance),
     c(1), 
     2)
   
@@ -605,8 +665,8 @@ test_that('schedule_infections correctly schedules new infections', {
     parameters,
     42 
   )
-  
-  actual_infected <- mockery::mock_args(infection_mock)[[1]][[5]]$to_vector()
+
+  actual_infected <- mockery::mock_args(infection_mock)[[1]][[7]]$to_vector()
   actual_asymp_infected <- mockery::mock_args(asymp_mock)[[1]][[4]]$to_vector()
   
   expect_equal(
@@ -642,16 +702,25 @@ test_that('prophylaxis is considered for medicated humans', {
   m <- mockery::mock(seq(3))
   mockery::stub(calculate_infections, 'bernoulli_multi_p', m)
   
-  calculate_infections(
+  
+  infection_outcome <- CompetingOutcome$new(
+    targeted_process = function(timestep, target){
+      infection_outcome_process(timestep, target, variables, renderer, parameters)
+    },
+    size = 4
+  )
+  
+  infection_rates <- calculate_infections(
     variables,
     bitten,
     parameters,
     mock_render(timestep),
-    timestep
+    timestep,
+    infection_outcome
   )
   
   expect_equal(
-    mockery::mock_args(m)[[1]][[1]],
+    rate_to_prob(infection_rates[infection_rates!=0]),
     c(2.491951e-07, 2.384032e-01, 5.899334e-01),
     tolerance = 1e-3
   )
@@ -774,6 +843,7 @@ test_that('update_severe_disease renders with no infections', {
   
   render_function <- mockery::mock()
   mockery::stub(update_severe_disease, 'incidence_renderer', render_function)
+  mockery::stub(update_severe_disease, 'incidence_probability_renderer', render_function)
   empty <- individual::Bitset$new(population)
   
   update_severe_disease(
@@ -783,13 +853,24 @@ test_that('update_severe_disease renders with no infections', {
     parameters,
     renderer
   )
-  
+
   mockery::expect_args(
     render_function,
     1,
     variables$birth,
     renderer,
     empty,
+    'inc_severe_',
+    0,
+    5 * 365,
+    timestep
+  )
+  
+  mockery::expect_args(
+    render_function,
+    2,
+    variables$birth,
+    renderer,
     empty,
     double(0),
     'inc_severe_',
@@ -972,9 +1053,9 @@ test_that('calculate_treated() successfully adds additional resistance columns t
   variables <- list(
     state = list(queue_update = mockery::mock()),
     infectivity = list(queue_update = mockery::mock()),
+    recovery_rates = list(queue_update = mockery::mock()),
     drug = list(queue_update = mockery::mock()),
-    drug_time = list(queue_update = mockery::mock()),
-    dt = list(queue_update = mockery::mock())
+    drug_time = list(queue_update = mockery::mock())
   )
   renderer <- individual::Render$new(timesteps = 10)
   
diff --git a/tests/testthat/test-output.R b/tests/testthat/test-output.R
index 91d427e4..5cef37b9 100644
--- a/tests/testthat/test-output.R
+++ b/tests/testthat/test-output.R
@@ -17,15 +17,17 @@ test_that('Test age parameter function works', {
   
   sim <- run_simulation(timesteps = 1, parameters)
   
-  expect_true(
-    all(
-      paste0(rep(c("n_age",
-                   "n", "n_inc", "p_inc",
-                   "n","n_inc_clinical","p_inc_clinical",
-                   "n","n_inc_severe","p_inc_severe",
-                   "n","n_detect_lm","p_detect_lm","n_detect_pcr",
-                   "ica_mean", "icm_mean","id_mean","ib_mean","iva_mean","ivm_mean"), each = 3),"_",
-             age_limits[-4]+rep(c(0,rep(c(1:3), each = 3), rep(4, 4), 5:10), each = 3),"_",age_limits[-1]-1+rep(c(0,rep(c(1:3), each = 3), rep(4, 4), 5:10), each = 3)
-      ) %in%
-        names(sim)))
+  prefixes <- c("n_age",
+                "n_inc", "p_inc",
+                "n_inc_clinical","p_inc_clinical",
+                "n_inc_severe","p_inc_severe",
+                "n_detect_lm","p_detect_lm","n_detect_pcr",
+                "ica_mean", "icm_mean","id_mean","ib_mean","iva_mean","ivm_mean")
+  offsets <- c(0, rep(1, 2), rep(2, 2), rep(3, 2), rep(4, 3), 5:10)
+  expect_equal(length(prefixes), length(offsets))
+
+  expected <- paste0(rep(prefixes, each = 3),
+                     "_", age_limits[-4]+rep(offsets, each = 3),
+                     "_", age_limits[-1]-1+rep(offsets, each = 3))
+  expect_in(expected, names(sim))
 })
diff --git a/tests/testthat/test-pev.R b/tests/testthat/test-pev.R
index 09c4cf7d..f3c3d793 100644
--- a/tests/testthat/test-pev.R
+++ b/tests/testthat/test-pev.R
@@ -150,6 +150,13 @@ test_that('Infection considers pev efficacy', {
     rep(.2, 4)
   )
 
+  infection_outcome <- CompetingOutcome$new(
+    targeted_process = function(timestep, target){
+      infection_process_resolved_hazard(timestep, target, variables, renderer, parameters)
+    },
+    size = 4
+  )
+  
   # remove randomness from infection sampling
   bernoulli_mock <- mockery::mock(c(1, 2))
   mockery::stub(calculate_infections, 'bernoulli_multi_p', bernoulli_mock)
@@ -164,16 +171,17 @@ test_that('Infection considers pev efficacy', {
     depth = 4
   )
 
-  calculate_infections(
+  infection_rates <- calculate_infections(
     variables = variables,
     bitten_humans = individual::Bitset$new(4)$insert(seq(4)),
     parameters = parameters,
     renderer = mock_render(timestep),
-    timestep = timestep
+    timestep = timestep,
+    infection_outcome = infection_outcome
   )
 
   expect_equal(
-    mockery::mock_args(bernoulli_mock)[[1]][[1]],
+    rate_to_prob(infection_rates[infection_rates!=0]),
     c(0.590, 0.590, 0.215, 0.244),
     tolerance=1e-3
   )
diff --git a/tests/testthat/test-render.R b/tests/testthat/test-render.R
index f42c15b4..6396d4ca 100644
--- a/tests/testthat/test-render.R
+++ b/tests/testthat/test-render.R
@@ -26,14 +26,6 @@ test_that('that default rendering works', {
   mockery::expect_args(
     renderer$render_mock(),
     1,
-    'n_730_3650',
-    3,
-    timestep
-  )
-  
-  mockery::expect_args(
-    renderer$render_mock(),
-    2,
     'n_detect_lm_730_3650',
     2,
     timestep
@@ -41,7 +33,7 @@ test_that('that default rendering works', {
 
   mockery::expect_args(
     renderer$render_mock(),
-    3,
+    2,
     'p_detect_lm_730_3650',
     1.5,
     timestep
@@ -49,7 +41,7 @@ test_that('that default rendering works', {
 
   mockery::expect_args(
     renderer$render_mock(),
-    4,
+    3,
     'n_detect_pcr_730_3650',
     3,
     timestep
@@ -83,7 +75,7 @@ test_that('that default rendering works when no one is in the age range', {
   mockery::expect_args(
     renderer$render_mock(),
     1,
-    'n_730_3650',
+    'n_detect_lm_730_3650',
     0,
     timestep
   )
@@ -101,6 +93,15 @@ test_that('that clinical incidence rendering works', {
     birth,
     renderer,
     individual::Bitset$new(4)$insert(c(1, 2, 4)),
+    'inc_clinical_',
+    c(0, 2) * year,
+    c(5, 10) * year,
+    timestep
+  )
+
+  incidence_probability_renderer(
+    birth,
+    renderer,
     individual::Bitset$new(4)$insert(seq(4)),
     c(.1, .2, .3, .4),
     'inc_clinical_',
@@ -112,7 +113,7 @@ test_that('that clinical incidence rendering works', {
   mockery::expect_args(
     renderer$render_mock(),
     1,
-    'n_0_1825',
+    'n_inc_clinical_0_1825',
     2,
     timestep
   )
@@ -120,7 +121,7 @@ test_that('that clinical incidence rendering works', {
   mockery::expect_args(
     renderer$render_mock(),
     2,
-    'n_inc_clinical_0_1825',
+    'n_inc_clinical_730_3650',
     2,
     timestep
   )
@@ -132,11 +133,68 @@ test_that('that clinical incidence rendering works', {
     .3,
     timestep
   )
+  
+  mockery::expect_args(
+    renderer$render_mock(),
+    4,
+    'p_inc_clinical_730_3650',
+    .6,
+    timestep
+  )
+})
+
 
+test_that('that age rendering works', {
+  timestep <- 0
+  year <- 365
+  birth <- individual::IntegerVariable$new(
+    -c(0:101) * year
+  )
+  
+  parameters <- get_parameters(overrides = list(prevalence_rendering_min_ages = c(0, 1, 5) * year,
+                                                prevalence_rendering_max_ages = c(2, 5, 8) * year - 1,
+                                                incidence_rendering_min_ages = c(3, 5, 10) * year,
+                                                incidence_rendering_max_ages = c(10, 20, 100) * year - 1,
+                                                clinical_incidence_rendering_min_ages = c(0, 1, 5) * year,
+                                                clinical_incidence_rendering_max_ages = c(2, 5, 8) * year - 1))
+  
+  renderer <- mock_render(1)
+  process <- create_age_group_renderer(
+    birth,
+    parameters, 
+    renderer)
+
+  process(timestep)
+  
+  
+  mockery::expect_args(
+    renderer$render_mock(),
+    1,
+    'n_age_0_729',
+    2,
+    timestep
+  )
+  
   mockery::expect_args(
     renderer$render_mock(),
+    2,
+    'n_age_365_1824',
     4,
-    'n_730_3650',
+    timestep
+  )
+  
+  mockery::expect_args(
+    renderer$render_mock(),
+    3,
+    'n_age_1095_3649',
+    7,
+    timestep
+  )
+  
+  mockery::expect_args(
+    renderer$render_mock(),
+    4,
+    'n_age_1825_2919',
     3,
     timestep
   )
@@ -144,16 +202,16 @@ test_that('that clinical incidence rendering works', {
   mockery::expect_args(
     renderer$render_mock(),
     5,
-    'n_inc_clinical_730_3650',
-    2,
+    'n_age_1825_7299',
+    15,
     timestep
   )
 
   mockery::expect_args(
     renderer$render_mock(),
     6,
-    'p_inc_clinical_730_3650',
-    .6,
+    'n_age_3650_36499',
+    90,
     timestep
   )
 })
diff --git a/vignettes/Antimalarial_Resistance.Rmd b/vignettes/Antimalarial_Resistance.Rmd
index 1c59de93..3bd2fb1a 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_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
+sim_out_baseline$pfpr210 = sim_out_baseline$n_detect_lm_730_3650/sim_out_baseline$n_age_730_3650
+sim_out_clin_treatment$pfpr210 = sim_out_clin_treatment$n_detect_lm_730_3650/sim_out_clin_treatment$n_age_730_3650
+sim_out_resistance$pfpr210 = sim_out_resistance$n_detect_lm_730_3650/sim_out_resistance$n_age_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_lm_730_3650/dynamic_resistance_output$n_730_3650
+dynamic_resistance_output$pfpr210 <- dynamic_resistance_output$n_detect_lm_730_3650/dynamic_resistance_output$n_age_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_lm_730_3650/simulation_outputs$base$n_730_3650,
+     y = simulation_outputs$base$n_detect_lm_730_3650/simulation_outputs$base$n_age_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_lm_730_3650/simulation_outputs$treatment$n_730_3650,
+      y = simulation_outputs$treatment$n_detect_lm_730_3650/simulation_outputs$treatment$n_age_730_3650,
       col = cols[4])
 
 lines(x = simulation_outputs$etf$timestep,
-      y = simulation_outputs$etf$n_detect_lm_730_3650/simulation_outputs$etf$n_730_3650,
+      y = simulation_outputs$etf$n_detect_lm_730_3650/simulation_outputs$etf$n_age_730_3650,
       col = cols[5])
 
 lines(x = simulation_outputs$spc$timestep,
-      y = simulation_outputs$spc$n_detect_lm_730_3650/simulation_outputs$spc$n_730_3650,
+      y = simulation_outputs$spc$n_detect_lm_730_3650/simulation_outputs$spc$n_age_730_3650,
       col = cols[6])
 
 lines(x = simulation_outputs$etf_spc$timestep,
-      y = simulation_outputs$etf_spc$n_detect_lm_730_3650/simulation_outputs$etf_spc$n_730_3650,
+      y = simulation_outputs$etf_spc$n_detect_lm_730_3650/simulation_outputs$etf_spc$n_age_730_3650,
       col = cols[7])
 
 # Add vlines to indicate when SP-AQ were administered:
diff --git a/vignettes/Carrying-capacity.Rmd b/vignettes/Carrying-capacity.Rmd
index db38e54a..e17ac25d 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_lm_730_3650 / s$n_730_3650
+s$pfpr <- s$n_detect_lm_730_3650 / s$n_age_730_3650
 set.seed(123)
 s_seasonal <- run_simulation(timesteps = timesteps, parameters = p_seasonal)
-s_seasonal$pfpr <- s_seasonal$n_detect_lm_730_3650 / s_seasonal$n_730_3650
+s_seasonal$pfpr <- s_seasonal$n_detect_lm_730_3650 / s_seasonal$n_age_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_lm_730_3650 / s_lsm$n_730_3650
+s_lsm$pfpr <- s_lsm$n_detect_lm_730_3650 / s_lsm$n_age_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_lm_730_3650 / s_stephensi$n_730_3650
+s_stephensi$pfpr <- s_stephensi$n_detect_lm_730_3650 / s_stephensi$n_age_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_lm_730_3650 / s_flexible$n_730_3650
+s_flexible$pfpr <- s_flexible$n_detect_lm_730_3650 / s_flexible$n_age_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/Demography.Rmd b/vignettes/Demography.Rmd
index d8c96a93..9de15d9c 100644
--- a/vignettes/Demography.Rmd
+++ b/vignettes/Demography.Rmd
@@ -27,7 +27,7 @@ The dynamics of malaria transmission, and the efficacy of interventions designed
 
 ## Age group rendering
 
-First, we'll establish a base set of parameters using the `get_parameters()` function and accept the default values. The `run_simulation()` function's default behaviour is to output only the number of individuals aged 2-10 years old (`output$n_730_3650`, where 730 and 3650 are the ages in days). However, the user can instruct the model to output the number of individuals in age groups of their choosing using the `age_group_rendering_min_ages` and `age_group_rendering_max_ages` parameters. These arguments take vectors containing the minimum and maximum ages (in daily time steps) of each age group to be rendered. To allow us to see the effect of changing demographic parameters, we'll use this functionality to output the number of individuals in ages groups ranging from 0 to 85 at 5 year intervals. Note that the same is possible for other model outputs using their equivalent min/max age-class rendering arguments (`n_detect` , `p_detect`, `n_severe`, `n_inc`, `p_inc`, `n_inc_clinical`, `p_inc_clinical`, `n_inc_severe`, and `p_inc_severe`, run `?run_simulation()` for more detail).
+First, we'll establish a base set of parameters using the `get_parameters()` function and accept the default values. The `run_simulation()` function's default behaviour is to output only the number of individuals aged 2-10 years old (`output$n_age_730_3650`, where 730 and 3650 are the ages in days). However, the user can instruct the model to output the number of individuals in age groups of their choosing using the `age_group_rendering_min_ages` and `age_group_rendering_max_ages` parameters. These arguments take vectors containing the minimum and maximum ages (in daily time steps) of each age group to be rendered. To allow us to see the effect of changing demographic parameters, we'll use this functionality to output the number of individuals in ages groups ranging from 0 to 85 at 5 year intervals. Note that the same is possible for other model outputs using their equivalent min/max age-class rendering arguments (`n_detect` , `p_detect`, `n_severe`, `n_inc`, `p_inc`, `n_inc_clinical`, `p_inc_clinical`, `n_inc_severe`, and `p_inc_severe`, run `?run_simulation()` for more detail).
 
 We next use the `set_equilibrium()` function to tune the initial parameter set to those required to observe the specified initial entomological inoculation rate (`starting_EIR`) at equilibrium. We now have a set of default parameters ready to use to run simulations.
 
diff --git a/vignettes/EIRprevmatch.Rmd b/vignettes/EIRprevmatch.Rmd
index 7040ed0c..cdf46b21 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_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.
+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_age_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
@@ -161,7 +161,7 @@ malSim_prev <- lapply(
         'n_detect_lm_730_3650'
       ] / output[
         output$timestep %in% seq(4 * 365, 5 * 365),
-        'n_730_3650'
+        'n_age_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_lm_730_3650/x$n_730_3650)
+ prev_2_10 <- mean(x$n_detect_lm_730_3650/x$n_age_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_lm_730_3650 / cali_sim$n_730_3650
-malsim_pfpr2_10 <- malsim_sim$n_detect_lm_730_3650 / malsim_sim$n_730_3650
+cali_pfpr2_10 <- cali_sim$n_detect_lm_730_3650 / cali_sim$n_age_730_3650
+malsim_pfpr2_10 <- malsim_sim$n_detect_lm_730_3650 / malsim_sim$n_age_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 1d71aeb2..2da93a6e 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_lm_730_3650/mda_output$n_730_3650,
+     y = mda_output$n_detect_lm_730_3650/mda_output$n_age_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_lm_730_3650/noint_output$n_730_3650,
+     y = noint_output$n_detect_lm_730_3650/noint_output$n_age_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_lm_730_3650/smc_output$n_730_3650,
+plot(x = smc_output$timestep, y = smc_output$n_detect_lm_730_3650/smc_output$n_age_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_lm_730_3650/smc_output$n_7
 
 # Add the dynamics for no-intervention simulation
 lines(x = noint_output$timestep,
-     y = noint_output$n_detect_lm_730_3650/noint_output$n_730_3650,
+     y = noint_output$n_detect_lm_730_3650/noint_output$n_age_730_3650,
      col = cols[4])
 
 # Add lines to indicate SMC events:
diff --git a/vignettes/Metapopulation.Rmd b/vignettes/Metapopulation.Rmd
index ad624b05..0f3bf776 100644
--- a/vignettes/Metapopulation.Rmd
+++ b/vignettes/Metapopulation.Rmd
@@ -112,7 +112,7 @@ output3$mix <- 'perfectly-mixed'
 output <- rbind(output1, output2, output3)
 
 # get mean PfPR 2-10 by year
-output$prev2to10 = output$p_detect_lm_730_3650 / output$n_730_3650
+output$prev2to10 = output$p_detect_lm_730_3650 / output$n_age_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 2eb55b75..3a9c4a32 100644
--- a/vignettes/Model.Rmd
+++ b/vignettes/Model.Rmd
@@ -153,7 +153,7 @@ The `run_simulation()` function then simulates malaria transmission dynamics and
 -   `id_mean`: mean immunity from detected using microscopy
 -   `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_age_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_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*)
@@ -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_lm_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_age_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_lm_730_3650/test_sim$n_730_3650
+test_sim$PfPR2_10 <- test_sim$n_detect_lm_730_3650/test_sim$n_age_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_lm_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_age_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 9b177b8a..eec04fd7 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_lm_730_3650 / sim$n_730_3650,
+  sim$n_detect_lm_730_3650 / sim$n_age_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_lm_730_3650[1:500] / sim$n_730_3650[1:500],
+  sim$n_detect_lm_730_3650[1:500] / sim$n_age_730_3650[1:500],
   t = "l",
   ylim = c(0, 1),
   ylab = "PfPr",
@@ -78,7 +78,7 @@ for (i in 1:7) {
   
   sim <- run_simulation(timesteps = 500, param_draw)
   
-  lines(sim$timestep, sim$n_detect_lm_730_3650 / sim$n_730_3650, col = cols[i])
+  lines(sim$timestep, sim$n_detect_lm_730_3650 / sim$n_age_730_3650, col = cols[i])
 }
 ```
 
diff --git a/vignettes/SetSpecies.Rmd b/vignettes/SetSpecies.Rmd
index d6a2e61f..762a3ccf 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_lm_730_3650 / output_endophilic$n_730_3650, 
+  plot(x = output_endophilic$timestep, y = output_endophilic$n_detect_lm_730_3650 / output_endophilic$n_age_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_lm_730_3650 / output_exophilic$n_730_3650,
+  lines(x = output_exophilic$timestep, y = output_exophilic$n_detect_lm_730_3650 / output_exophilic$n_age_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/Vaccines.Rmd b/vignettes/Vaccines.Rmd
index 4638b1a8..b36f09ca 100644
--- a/vignettes/Vaccines.Rmd
+++ b/vignettes/Vaccines.Rmd
@@ -37,9 +37,9 @@ plot_inci <- function(type = "not seasonal"){
     comparison <- output_control
     vaccinetime <- 1
   }
-  output$clinical_incidence <- 1000 * output$n_inc_clinical_0_1825 / output$n_0_1825
+  output$clinical_incidence <- 1000 * output$n_inc_clinical_0_1825 / output$n_age_0_1825
   output$time_year <- output$timestep / year
-  comparison$clinical_incidence <- 1000 * comparison$n_inc_clinical_0_1825 / comparison$n_0_1825
+  comparison$clinical_incidence <- 1000 * comparison$n_inc_clinical_0_1825 / comparison$n_age_0_1825
   comparison$time_year <- comparison$timestep / year
 
   plot(x = output$time_year, y = output$clinical_incidence,
@@ -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_lm_730_3650/output$n_730_3650, 
+  plot(x = output$time_year, y = output$n_detect_lm_730_3650/output$n_age_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_lm_730_3650/comparison$n_730_3650,
+  lines(x = comparison$time_year, y = comparison$n_detect_lm_730_3650/comparison$n_age_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)
diff --git a/vignettes/Variation.Rmd b/vignettes/Variation.Rmd
index fb563ca0..10da8a80 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_lm_730_3650 / output$n_730_3650, 
+  plot(x = output$timestep, y = output$n_detect_lm_730_3650 / output$n_age_730_3650, 
        type = "l", col = cols[3], ylim = ylim,
        xlab = "Time (days)", ylab = ylab, lwd = 1,
        xaxs = "i", yaxs = "i")
@@ -40,7 +40,7 @@ plot_inci <- function(output, ylab = TRUE, ylim){
   if (ylab == TRUE) {
     ylab = "Incidence per 1000 children aged 0-5 years"
     } else {ylab = ""}
-  plot(x = output$timestep, y = output$n_inc_clinical_0_1825 / output$n_0_1825 * 1000, 
+  plot(x = output$timestep, y = output$n_inc_clinical_0_1825 / output$n_age_0_1825 * 1000, 
        type = "l", col = cols[5], ylim = ylim,
        xlab = "Time (days)", ylab = ylab, lwd = 1,
        xaxs = "i", yaxs = "i")
diff --git a/vignettes/VectorControl_Bednets.Rmd b/vignettes/VectorControl_Bednets.Rmd
index c1c801a8..f8cf4dd9 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_lm_730_3650 / output$n_730_3650, 
+  plot(x = output$timestep, y = output$n_detect_lm_730_3650 / output$n_age_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_lm_730_3650 / output_control$n_730_3650,
+  lines(x = output_control$timestep, y = output_control$n_detect_lm_730_3650 / output_control$n_age_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 4dc924a7..1db0727d 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_lm_730_3650 / output$n_730_3650, 
+  plot(x = output$timestep, y = output$n_detect_lm_730_3650 / output$n_age_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_lm_730_3650 / output_control$n_730_3650,
+  lines(x = output_control$timestep, y = output_control$n_detect_lm_730_3650 / output_control$n_age_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)