diff --git a/NAMESPACE b/NAMESPACE
index 46953fd3..267d2003 100644
--- a/NAMESPACE
+++ b/NAMESPACE
@@ -8,23 +8,21 @@ export(fun_params)
 export(gamb_params)
 export(get_correlation_parameters)
 export(get_parameters)
-export(parameterise_human_equilibrium)
 export(parameterise_mosquito_equilibrium)
 export(parameterise_total_M)
 export(peak_season_offset)
-export(remove_unused_jamie)
 export(run_simulation)
 export(run_simulation_with_repetitions)
 export(set_bednets)
 export(set_clinical_treatment)
 export(set_drugs)
+export(set_equilibrium)
 export(set_mda)
 export(set_rtss)
 export(set_smc)
 export(set_species)
 export(set_spraying)
 export(set_tbv)
-export(translate_jamie)
 importFrom(MASS,mvrnorm)
 importFrom(Rcpp,sourceCpp)
 importFrom(stats,dweibull)
diff --git a/R/RcppExports.R b/R/RcppExports.R
index 658e9057..0bd73cad 100644
--- a/R/RcppExports.R
+++ b/R/RcppExports.R
@@ -13,6 +13,18 @@ create_adult_solver <- function(model, init) {
     .Call(`_malariasimulation_create_adult_solver`, model, init)
 }
 
+create_history <- function(size, default_value) {
+    .Call(`_malariasimulation_create_history`, size, default_value)
+}
+
+history_at <- function(history, timestep) {
+    .Call(`_malariasimulation_history_at`, history, timestep)
+}
+
+history_push <- function(history, value, timestep) {
+    invisible(.Call(`_malariasimulation_history_push`, history, value, timestep))
+}
+
 carrying_capacity <- function(timestep, model_seasonality, g0, g, h, K0, R_bar) {
     .Call(`_malariasimulation_carrying_capacity`, timestep, model_seasonality, g0, g, h, K0, R_bar)
 }
@@ -25,19 +37,6 @@ rainfall <- function(t, g0, g, h) {
     .Call(`_malariasimulation_rainfall`, t, g0, g, h)
 }
 
-#' @title Mosquito emergence process
-#' @description Move mosquitos from NonExistent to Sm in line with the number of
-#' pupals in the ODE models
-#'
-#' @param solvers a list of solver objects for each species of mosquito
-#' @param state the variable for the mosquito state
-#' @param species the variable for the mosquito species
-#' @param species_names a vector of category names for the species variable
-#' @param dpl the delay for pupal growth (in timesteps)
-create_mosquito_emergence_process_cpp <- function(solvers, state, species, species_names, dpl) {
-    .Call(`_malariasimulation_create_mosquito_emergence_process_cpp`, solvers, state, species, species_names, dpl)
-}
-
 create_mosquito_model <- function(beta, de, mue, K0, gamma, dl, mul, dp, mup, total_M, model_seasonality, g0, g, h, R_bar) {
     .Call(`_malariasimulation_create_mosquito_model`, beta, de, mue, K0, gamma, dl, mul, dp, mup, total_M, model_seasonality, g0, g, h, R_bar)
 }
diff --git a/R/biting_process.R b/R/biting_process.R
index 49200fb4..14647139 100644
--- a/R/biting_process.R
+++ b/R/biting_process.R
@@ -15,7 +15,9 @@ create_biting_process <- function(
   models,
   variables,
   events,
-  parameters
+  parameters,
+  lagged_foim,
+  lagged_eir
   ) {
   function(timestep) {
     # Calculate combined EIR
@@ -29,7 +31,9 @@ create_biting_process <- function(
       events,
       age,
       parameters,
-      timestep
+      timestep,
+      lagged_foim,
+      lagged_eir
     )
 
     simulate_infection(
@@ -53,7 +57,9 @@ simulate_bites <- function(
   events,
   age,
   parameters,
-  timestep
+  timestep,
+  lagged_foim,
+  lagged_eir
   ) {
   bitten_humans <- individual::Bitset$new(parameters$human_population)
 
@@ -70,7 +76,8 @@ simulate_bites <- function(
 
   # Calculate pi (the relative biting rate for each human)
   psi <- unique_biting_rate(age, parameters)
-  .pi <- human_pi(variables$zeta$get_values(), psi)
+  zeta <- variables$zeta$get_values()
+  .pi <- human_pi(zeta, psi)
 
   # Get some indices for later
   if (parameters$individual_mosquitoes) {
@@ -82,43 +89,22 @@ simulate_bites <- function(
   EIR <- 0
 
   for (s_i in seq_along(parameters$species)) {
-    if (!parameters$individual_mosquitoes) {
-      solver_states <- solver_get_states(solvers[[s_i]])
-    }
-
-    if (parameters$individual_mosquitoes) {
-      species_index <- variables$species$get_index_of(
-        parameters$species[[s_i]]
-      )$and(adult_index)
-      n_infectious <- infectious_index$copy()$and(species_index)$size()
-    } else {
-      n_infectious <- solver_states[[ADULT_ODE_INDICES['Im']]]
-    }
-
+    solver_states <- solver_get_states(solvers[[s_i]])
     p_bitten <- prob_bitten(timestep, variables, s_i, parameters)
-
     Q0 <- parameters$Q0[[s_i]]
-    Z <- average_p_repelled(p_bitten$prob_repelled, .pi, Q0)
     W <- average_p_successful(p_bitten$prob_bitten_survives, .pi, Q0)
-    renderer$render(
-      paste0('p_repelled_', parameters$species[[s_i]]),
-      Z,
-      timestep
-    )
-    renderer$render(
-      paste0('p_feed_survives_', parameters$species[[s_i]]),
-      W,
-      timestep
-    )
+    Z <- average_p_repelled(p_bitten$prob_repelled, .pi, Q0)
     f <- blood_meal_rate(s_i, Z, parameters)
 
-    lambda <- effective_biting_rate(
-      .pi,
-      age,
+    lambda <- .effective_biting_rates(
       s_i,
+      psi,
+      .pi,
+      zeta,
       p_bitten,
-      f,
       W,
+      Z,
+      f,
       parameters
     )
 
@@ -128,8 +114,26 @@ simulate_bites <- function(
       timestep
     )
 
-    n_bites <- rpois(1, n_infectious * sum(lambda))
-    EIR <- EIR + n_bites
+    if (parameters$individual_mosquitoes) {
+      species_index <- variables$species$get_index_of(
+        parameters$species[[s_i]]
+      )$and(adult_index)
+      n_infectious <- calculate_infectious_individual(
+        s_i,
+        variables,
+        infectious_index,
+        adult_index,
+        species_index,
+        parameters
+      )
+    } else {
+      n_infectious <- calculate_infectious_compartmental(solver_states)
+    }
+
+    lagged_eir[[s_i]]$save(n_infectious * sum(lambda), timestep)
+    species_eir <- lagged_eir[[s_i]]$get(timestep - parameters$de)
+    EIR <- EIR + species_eir / sum(psi)
+    n_bites <- rpois(1, species_eir)
     if (n_bites > 0) {
       bitten_humans$insert(
         sample.int(
@@ -142,6 +146,7 @@ simulate_bites <- function(
     }
 
     foim <- calculate_foim(human_infectivity, lambda)
+    lagged_foim$save(foim, timestep)
     renderer$render(paste0('FOIM_', s_i), foim, timestep)
     mu <- death_rate(f, W, Z, s_i, parameters)
     renderer$render(paste0('mu_', s_i), mu, timestep)
@@ -155,7 +160,7 @@ simulate_bites <- function(
 
       biting_effects_individual(
         variables,
-        foim,
+        lagged_foim$get(timestep - parameters$delay_gam),
         events,
         s_i,
         susceptible_species_index,
@@ -168,7 +173,7 @@ simulate_bites <- function(
       adult_mosquito_model_update(
         models[[s_i]],
         mu,
-        foim,
+        lagged_foim$get(timestep - parameters$delay_gam),
         solver_states[[ADULT_ODE_INDICES['Sm']]],
         f
       )
@@ -180,25 +185,82 @@ simulate_bites <- function(
   bitten_humans
 }
 
+
 # =================
 # Utility functions
 # =================
 
-#' @title Calculate the effective biting rate for a species on each human given
-#' vector control interventions
-#' @description
-#' Implemented from Griffin et al 2010 S2 page 6
-#' @param .pi relative biting rate for each human
-#' @param age of each human (timesteps)
-#' @param species to model
-#' @param p_bitten the probabilities of feeding given vector controls
-#' @param f blood meal rate
-#' @param W average probability of a successful bite
-#' @param parameters of the model
-#' @noRd
-effective_biting_rate <- function(.pi, age, species, p_bitten, f, W, parameters) {
+calculate_eir <- function(species, solvers, variables, parameters, timestep) {
+  lambda <- effective_biting_rates(species, variables, parameters, timestep)
+  infectious <- calculate_infectious(species, solvers, variables, parameters)
+  infectious * sum(lambda)
+}
+
+effective_biting_rates <- function(species, variables, parameters, timestep) {
+  age <- get_age(variables$birth$get_values(), timestep)
+  psi <- unique_biting_rate(age, parameters)
+  zeta <- variables$zeta$get_values()
+  p_bitten <- prob_bitten(timestep, variables, species, parameters)
+  .pi <- human_pi(zeta, psi)
+  Q0 <- parameters$Q0[[species]]
+  W <- average_p_successful(p_bitten$prob_bitten_survives, .pi, Q0)
+  Z <- average_p_repelled(p_bitten$prob_repelled, .pi, Q0)
+  f <- blood_meal_rate(species, Z, parameters)
+  .effective_biting_rates(species, psi, .pi, zeta, p_bitten, W, Z, f, parameters)
+}
+
+.effective_biting_rates <- function(
+  species,
+  psi,
+  .pi,
+  zeta,
+  p_bitten,
+  W,
+  Z,
+  f,
+  parameters
+  ) {
   a <- human_blood_meal_rate(f, species, W, parameters)
-  a * .pi * p_bitten$prob_bitten / sum(.pi * p_bitten$prob_bitten_survives)
+  a * intervention_coefficient(p_bitten) * zeta * psi
+}
+
+calculate_infectious <- function(species, solvers, variables, parameters) {
+  if (parameters$individual_mosquitoes) {
+    adult_index <- variables$mosquito_state$get_index_of('NonExistent')$not()
+    return(
+      calculate_infectious_individual(
+        species,
+        variables,
+        variables$mosquito_state$get_index_of('Im'),
+        adult_index,
+        variables$species$get_index_of(
+          parameters$species[[species]]
+        )$and(adult_index),
+        parameters
+      )
+    )
+  }
+  calculate_infectious_compartmental(solver_get_states(solvers[[species]]))
+}
+
+calculate_infectious_individual <- function(
+  species,
+  variables,
+  infectious_index,
+  adult_index,
+  species_index,
+  parameters
+  ) {
+  
+  infectious_index$copy()$and(species_index)$size()
+}
+
+calculate_infectious_compartmental <- function(solver_states) {
+  solver_states[[ADULT_ODE_INDICES['Im']]]
+}
+
+intervention_coefficient <- function(p_bitten) {
+   p_bitten$prob_bitten / sum(p_bitten$prob_bitten_survives)
 }
 
 human_pi <- function(zeta, psi) {
diff --git a/R/compatibility.R b/R/compatibility.R
index b6920d13..e14a64c8 100644
--- a/R/compatibility.R
+++ b/R/compatibility.R
@@ -1,8 +1,15 @@
+EQUILIBRIUM_AGES <- 0:999 / 10
 
 inverse_param <- function(name, new_name) {
   function(params) { list(new_name, 1 / params[[name]]) }
 }
 
+mean_param <- function(new_name, name, weights) {
+  function(params) {
+    list(new_name, weighted.mean(params[[name]], params[[weights]]))
+  }
+}
+
 translations = list(
   eta = inverse_param('eta', 'average_age'),
   rho = 'rho',
@@ -43,14 +50,59 @@ translations = list(
   PM = 'pcm',
   tau = 'dem',
   Q0 = 'Q0',
-  f = 'blood_meal_rates'
+  f = 'blood_meal_rates',
+  tl = 'delay_gam'
 )
 
-#' @title translate parameter keys from jamie's format to ones compatible
-#' with this IBM 
-#' @param params with keys in the jamie's format
-#' @export
-translate_jamie <- function(params) {
+back_translations = list(
+  average_age = inverse_param('average_age', 'eta'),
+  rho = 'rho',
+  a0  = 'a0',
+  da = inverse_param('da', 'rA'),
+  dd = inverse_param('dd', 'rD'),
+  du  = inverse_param('du', 'rU'),
+  dt  = inverse_param('dt', 'rT'),
+  de  = 'dE',
+  cd  = 'cD',
+  cu  = 'cU',
+  ct  = 'cT',
+  d1  = 'd1',
+  rid = 'dd',
+  id0 = 'ID0',
+  kd  = 'kd',
+  ud  = 'ud',
+  ad  = 'ad0',
+  gammad  = 'gd',
+  b0  = 'b0',
+  b1  = 'b1',
+  d1  = 'd1',
+  rb  = 'db',
+  ib0 = 'IB0',
+  kb  = 'kb',
+  ub  = 'ub',
+  phi0= 'phi0',
+  phi1= 'phi1',
+  rc  = 'dc',
+  ic0 = 'IC0',
+  kc  = 'kc',
+  uc  = 'uc',
+  rm  = 'dm',
+  mum  = mean_param('mu', 'mum', 'species_proportions'),
+  sigma_squared = 's2',
+  fd0 = 'fd0',
+  gamma1 = 'g_inf',
+  pcm = 'PM',
+  dem = 'tau',
+  Q0 = mean_param('Q0', 'Q0', 'species_proportions'),
+  blood_meal_rates = mean_param('f', 'blood_meal_rates', 'species_proportions'),
+  delay_gam = 'tl'
+)
+
+#' @description translate parameter keys from the malariaEquilibrium format
+#' to ones compatible with this IBM 
+#' @param params with keys in the malariaEquilibrium format
+#' @noRd
+translate_equilibrium <- function(params) {
   translated <- list()
   for (name in names(params)) {
     if(!name %in% names(translations)) {
@@ -68,20 +120,75 @@ translate_jamie <- function(params) {
   translated
 }
 
-#' @title remove parameter keys from jamie's format that are not used
+#' @description translate model parameters to the malariaEquilibrium format
+#' @param params model params
+#' @noRd
+translate_parameters <- function(params) {
+  translated <- malariaEquilibrium::load_parameter_set()
+  for (name in names(params)) {
+    if(name %in% names(back_translations)) {
+      translation <- back_translations[[name]]
+      if (is.character(translation)) {
+        translated[[translation]] <- params[[name]]
+      }
+      if (is.function(translation)) {
+        t <- translation(params)
+        translated[[t[[1]]]] <- t[[2]]
+      }
+    }
+  }
+  translated
+}
+
+#' @title remove parameter keys from the malariaEquilibrium format that are not used
 #' in this IBM 
-#' @param params with keys in the jamie's format
-#' @export
-remove_unused_jamie <- function(params) {
+#' @param params with keys in the malariaEquilibrium format
+#' @noRd
+remove_unused_equilibrium <- function(params) {
   remove_keys(
     params,
     c(
       'rP', # Prophylaxis state is no longer used, see `drug_parameters.R`
-      'tl', # unused!
-      'aA', # used for microscopy and pcr calculations
-      'aU', # used for microscopy and pcr calculations
+      'aA', # used for pcr calculations
+      'aU', # used for pcr calculations
       'cd_w', # unused!
       'cd_p'  # unused!
     )
   )
 }
+
+#' @title Set equilibrium
+#' @description This will update the IBM parameters to match the
+#' equilibrium parameters and set up the initial human and mosquito population
+#' to acheive init_EIR
+#' @param parameters model parameters to update
+#' @param init_EIR the desired initial EIR
+#' @param eq_params parameters from the malariaEquilibrium package, if null.
+#' The default malariaEquilibrium parameters will be used
+#' @export
+set_equilibrium <- function(parameters, init_EIR, eq_params = NULL) {
+  if (is.null(eq_params)) {
+    eq_params <- translate_parameters(parameters)
+  } else {
+    parameters <- c(
+      translate_equilibrium(remove_unused_equilibrium(eq_params)),
+      parameters
+    )
+  }
+  eq <- malariaEquilibrium::human_equilibrium(
+    EIR = init_EIR,
+    ft = sum(get_treatment_coverages(parameters, 1)),
+    p = eq_params,
+    age = EQUILIBRIUM_AGES,
+    h = malariaEquilibrium::gq_normal(parameters$n_heterogeneity_groups)
+  )
+  parameters <- c(
+    list(
+      init_foim = eq$FOIM,
+      init_EIR = init_EIR,
+      eq_params = eq_params
+    ),
+    parameters
+  )
+  parameterise_mosquito_equilibrium(parameters, init_EIR)
+}
diff --git a/R/disease_progression.R b/R/disease_progression.R
index 95be6c5a..c030637c 100644
--- a/R/disease_progression.R
+++ b/R/disease_progression.R
@@ -19,30 +19,19 @@ create_infection_update_listener <- function(
   }
 }
 
-#' @title Schedule progression of human disease at the start of the simulation
-#'
-#' @param event the event to schedule
-#' @param state the human state variable
-#' @param from_state the state this event applies to
-#' @param rate the average time spent in this state
-#' @noRd
-initialise_progression <- function(event, state, from_state, rate) {
-  target <- state$get_index_of(from_state)
-  event$schedule(
-    target,
-    log_uniform(target$size(), rate)
-  )
+create_progression_process <- function(event, state, from_state, rate) {
+  function(timestep) {
+    event$schedule(state$get_index_of(from_state)$sample(1/rate), 0)
+  }
 }
 
-#' @title Modelling the progression of the human disease
-#' @description schedules follow up infection events with the log uniform dist.
-#'
-#' @param event the event to schedule
-#' @param rate the average time spent in this state
-#' @noRd
-create_progression_listener <- function(event, rate) {
+create_rate_listener <- function(from_state, to_state, renderer) {
   function(timestep, target) {
-    event$schedule(target, log_uniform(target$size(), rate))
+    renderer$render(
+      paste0('rate_', from_state, '_', to_state),
+      target$size(),
+      timestep
+    )
   }
 }
 
diff --git a/R/events.R b/R/events.R
index f1d9790a..d1a37982 100644
--- a/R/events.R
+++ b/R/events.R
@@ -1,12 +1,16 @@
 create_events <- function(parameters) {
   events <- list(
+    # Disease progression events
+    asymptomatic_progression = individual::TargetedEvent$new(parameters$human_population),
+    subpatent_progression = individual::TargetedEvent$new(parameters$human_population),
+    recovery = individual::TargetedEvent$new(parameters$human_population),
+
     # Human infection events
     clinical_infection = individual::TargetedEvent$new(parameters$human_population),
     asymptomatic_infection = individual::TargetedEvent$new(parameters$human_population),
+
     # whether the infection is detected
     detection = individual::TargetedEvent$new(parameters$human_population), 
-    subpatent_infection = individual::TargetedEvent$new(parameters$human_population),
-    recovery = individual::TargetedEvent$new(parameters$human_population),
 
     # Vaccination events
     rtss_vaccination = individual::Event$new(),
@@ -33,31 +37,6 @@ create_events <- function(parameters) {
 }
 
 initialise_events <- function(events, variables, parameters) {
-  initialise_progression(
-    events$asymptomatic_infection,
-    variables$state,
-    'D',
-    parameters$de
-  )
-  initialise_progression(
-    events$subpatent_infection,
-    variables$state,
-    'A',
-    parameters$da
-  )
-  initialise_progression(
-    events$recovery,
-    variables$state,
-    'U',
-    parameters$du
-  )
-  initialise_progression(
-    events$recovery,
-    variables$state,
-    'Tr',
-    parameters$dt
-  )
-
   if (parameters$individual_mosquitoes) {
     incubating <- variables$mosquito_state$get_index_of('Pm')
     events$mosquito_infection$schedule(
@@ -113,6 +92,15 @@ attach_event_listeners <- function(
     )
   )
 
+  events$asymptomatic_progression$add_listener(
+    create_asymptomatic_update_listener(
+      variables,
+      parameters
+    )
+  )
+  events$asymptomatic_progression$add_listener(
+    create_rate_listener('D', 'A', renderer)
+  )
   events$asymptomatic_infection$add_listener(
     create_asymptomatic_update_listener(
       variables,
@@ -121,7 +109,7 @@ attach_event_listeners <- function(
   )
 
   # Recovery events
-  events$subpatent_infection$add_listener(
+  events$subpatent_progression$add_listener(
     create_infection_update_listener(
       variables$state,
       'U',
@@ -129,6 +117,10 @@ attach_event_listeners <- function(
       parameters$cu
     )
   )
+  events$subpatent_progression$add_listener(
+    create_rate_listener('A', 'U', renderer)
+  )
+
   events$recovery$add_listener(
     create_infection_update_listener(
       variables$state,
@@ -137,18 +129,15 @@ attach_event_listeners <- function(
       0
     )
   )
+  events$recovery$add_listener(
+    create_rate_listener('U', 'S', renderer)
+  )
 
   # ===========
   # Progression
   # ===========
   # When infection events fire, schedule the next stages of infection
 
-  events$clinical_infection$add_listener(
-    create_progression_listener(
-      events$asymptomatic_infection,
-      parameters$dd
-    )
-  )
   events$clinical_infection$add_listener(
     create_clinical_incidence_renderer(
       variables$birth,
@@ -157,19 +146,6 @@ attach_event_listeners <- function(
     )
   )
 
-  events$asymptomatic_infection$add_listener(
-    create_progression_listener(
-      events$subpatent_infection,
-      parameters$da
-    )
-  )
-  events$subpatent_infection$add_listener(
-    create_progression_listener(
-      events$recovery,
-      parameters$du
-    )
-  )
-
   events$detection$add_listener(
     create_incidence_renderer(
       variables$birth,
diff --git a/R/human_infection.R b/R/human_infection.R
index db50e299..a02f71f6 100644
--- a/R/human_infection.R
+++ b/R/human_infection.R
@@ -81,6 +81,7 @@ simulate_infection <- function(
   )
 
   renderer$render('n_treated', treated$size(), timestep)
+  renderer$render('n_infections', infected_humans$size(), timestep)
 
   schedule_infections(
     events,
@@ -263,10 +264,6 @@ calculate_treated <- function(
       timestep,
       treated_index
     )
-    recovery$schedule(
-      treated_index,
-      log_uniform(treated_index$size(), parameters$dt)
-    )
     detection$schedule(treated_index, 0)
   }
   treated_index
@@ -286,36 +283,27 @@ schedule_infections <- function(
   clinical_infections,
   treated,
   infections,
-  parameters
+  parameters,
+  asymptomatics
   ) {
-  included <- events$clinical_infection$get_scheduled()$or(
-    events$asymptomatic_infection$get_scheduled()
-  )$or(treated)$not()
+  included <- treated$not()
 
   to_infect <- clinical_infections$and(included)
-  to_infect_asym <- clinical_infections$not()$and(infections)$and(included)
+  to_infect_asym <- clinical_infections$not()$and(infections)$and(
+    included
+  )
 
   # change to symptomatic
   if(to_infect$size() > 0) {
     infection_times <- log_uniform(to_infect$size(), parameters$de)
-    events$clinical_infection$schedule(
-      to_infect,
-      infection_times
-    )
-    events$detection$schedule(to_infect, infection_times)
-    events$subpatent_infection$clear_schedule(to_infect)
-    events$recovery$clear_schedule(to_infect)
+    events$clinical_infection$schedule(to_infect, 0)
+    events$detection$schedule(to_infect, 0)
   }
 
   if(to_infect_asym$size() > 0) {
     infection_times <- log_uniform(to_infect_asym$size(), parameters$de)
-    events$asymptomatic_infection$schedule(
-      to_infect_asym,
-      infection_times
-    )
-    events$detection$schedule(to_infect_asym, infection_times)
-    events$subpatent_infection$clear_schedule(to_infect_asym)
-    events$recovery$clear_schedule(to_infect_asym)
+    events$asymptomatic_infection$schedule(to_infect_asym, 0)
+    events$detection$schedule(to_infect_asym, 0)
   }
 }
 
diff --git a/R/lag.R b/R/lag.R
new file mode 100644
index 00000000..2892e042
--- /dev/null
+++ b/R/lag.R
@@ -0,0 +1,19 @@
+LaggedValue <- R6::R6Class(
+  'LaggedValue',
+  private = list(
+    history = NULL
+  ),
+  public = list(
+    initialize = function(max_lag, default) {
+      private$history <- create_history(max_lag, default)
+    },
+
+    save = function(value, timestep) {
+      history_push(private$history, value, timestep)
+    },
+    
+    get = function(timestep) {
+      history_at(private$history, timestep)
+    }
+  )
+)
diff --git a/R/mortality_processes.R b/R/mortality_processes.R
index 5ec16c01..5fe18474 100644
--- a/R/mortality_processes.R
+++ b/R/mortality_processes.R
@@ -10,10 +10,10 @@
 #' @noRd
 create_mortality_process <- function(variables, events, renderer, parameters) {
   function(timestep) {
-    age <- get_age(
+    age <- trunc(get_age(
       variables$birth$get_values(),
       timestep
-    ) / 365
+    ) / 365)
 
     died <- individual::Bitset$new(parameters$human_population)
     died <- sample_bitset(died$not(), 1 / parameters$average_age)
@@ -36,7 +36,7 @@ create_mortality_process <- function(variables, events, renderer, parameters) {
     if (died$size() > 0) {
       # inherit immunity from parent in group
       sampleable <- individual::Bitset$new(parameters$human_population)
-      sampleable$insert(which(age >= 15 | age <= 35))
+      sampleable$insert(which(age == 20))
       for (group in seq(parameters$n_heterogeneity_groups)) {
 
         # get the individuals who died in this group
@@ -46,6 +46,9 @@ create_mortality_process <- function(variables, events, renderer, parameters) {
         if (died_group$size() > 0) {
           # find their mothers
           potential_mothers <- group_index$and(sampleable)$to_vector()
+          if (length(potential_mothers) == 0) {
+            potential_mothers = seq(parameters$human_population)
+          }
           mothers <- potential_mothers[
             sample.int(
               length(potential_mothers),
@@ -65,7 +68,8 @@ create_mortality_process <- function(variables, events, renderer, parameters) {
       events$detection$clear_schedule(died)
       events$clinical_infection$clear_schedule(died)
       events$asymptomatic_infection$clear_schedule(died)
-      events$subpatent_infection$clear_schedule(died)
+      events$asymptomatic_progression$clear_schedule(died)
+      events$subpatent_progression$clear_schedule(died)
       events$recovery$clear_schedule(died)
       events$throw_away_net$clear_schedule(died)
 
@@ -92,7 +96,6 @@ create_mortality_process <- function(variables, events, renderer, parameters) {
       variables$rtss_boosted$queue_update(-1, died)
       variables$tbv_vaccinated$queue_update(-1, died)
 
-
       # onwards infectiousness
       variables$infectivity$queue_update(0, died)
 
diff --git a/R/mosquito_biology.R b/R/mosquito_biology.R
index 8235e388..62c88105 100644
--- a/R/mosquito_biology.R
+++ b/R/mosquito_biology.R
@@ -197,3 +197,50 @@ biting_effects_individual <- function(
 
   events$mosquito_death$schedule(died, 0)
 }
+#' @title Mosquito emergence process
+#' @description Move mosquitos from NonExistent to Sm in line with the number of
+#' pupals in the ODE models
+#'
+#' @param solvers a list of solver objects for each species of mosquito
+#' @param state the variable for the mosquito state
+#' @param species the variable for the mosquito species
+#' @param species_names a character vector of species names for each solver
+#' @param dpl the delay for pupal growth (in timesteps)
+#' @noRd
+create_mosquito_emergence_process <- function(
+  solvers,
+  state,
+  species,
+  species_names,
+  dpl
+  ) {
+  rate <- .5 * 1 / dpl
+  function(timestep) {
+    p_counts <- vnapply(
+      solvers,
+      function(solver) {
+        solver_get_states(solver)[[ODE_INDICES[['P']]]]
+      }
+    )
+    n <- sum(p_counts) * rate
+    available <- state$get_size_of('NonExistent')
+    if (n > available) {
+      stop(paste0(
+        'Not enough mosquitoes (short by ',
+        n - available,
+        '). Please raise parameters$mosquito_limit. ',
+        'If you have used parameterise_mosquito_equilibrium,',
+        'your seasonality parameters lead to more mosquitoes than expected.'
+      ))
+    }
+    non_existent <- state$get_index_of('NonExistent')
+    latest <- 1
+    for (i in seq_along(species_names)) {
+      to_hatch <- p_counts[[i]] * rate
+      hatched <- bitset_at(non_existent, seq(latest, latest + to_hatch))
+      state$queue_update('Sm', hatched)
+      species$queue_update(species_names[[i]], hatched)
+      latest <- latest + to_hatch + 1
+    }
+  }
+}
diff --git a/R/parameters.R b/R/parameters.R
index 191cc0dd..5a85d3e3 100644
--- a/R/parameters.R
+++ b/R/parameters.R
@@ -126,8 +126,9 @@
 #'
 #' incubation periods:
 #'
-#' * de - delay for infection
-#' * dem - delay for infection in mosquitoes
+#' * de - Duration of the human latent period of infection
+#' * delay_gam - Lag from parasites to infectious gametocytes
+#' * dem - Extrinsic incubation period in mosquito population model
 #'
 #' vector biology:
 #' species specific values are vectors
@@ -232,13 +233,15 @@
 #' simulation
 #' * individual_mosquitoes - boolean whether adult mosquitoes are modelled
 #' individually or compartmentaly
+#' * enable_heterogeneity - boolean whether to include heterogeneity in biting
+#' rates
 #'
 #' @export
 get_parameters <- function(overrides = list()) {
   parameters <- list(
     dd    = 5,
     dt    = 5,
-    da    = 195,
+    da    = 200,
     du    = 110,
     del   = 6.64,
     dl    = 3.72,
@@ -255,10 +258,10 @@ get_parameters <- function(overrides = list()) {
     rva   = 30 * 365,
     rid   = 10 * 365,
     # blood immunity parameters
-    b0    = 0.590076,
+    b0    = 0.59,
     b1    = 0.5,
-    ib0   = 43.8787,
-    kb    = 2.15506,
+    ib0   = 43.9,
+    kb    = 2.16,
     # immunity boost grace periods
     ub    = 7.2,
     uc    = 6.06,
@@ -288,8 +291,9 @@ get_parameters <- function(overrides = list()) {
     gammav  = 2.91282,
     iv0     = 1.09629,
     # delay for infection
-    de      = 12,
-    dem     = 10,
+    de        = 12,
+    delay_gam = 12.5,
+    dem       = 10,
     # asymptomatic immunity parameters
     fd0   = 0.007055,
     ad    = 21.9 * 365,
@@ -414,7 +418,8 @@ get_parameters <- function(overrides = list()) {
     # misc
     human_population = 100,
     mosquito_limit   = 100 * 1000,
-    individual_mosquitoes = TRUE
+    individual_mosquitoes = TRUE,
+    enable_heterogeneity = TRUE
   )
 
   # Override parameters with any client specified ones
@@ -444,26 +449,6 @@ get_parameters <- function(overrides = list()) {
   parameters
 }
 
-#' @title Parameterise equilibrium proportions
-#' @description parameterise equilibrium proportions from a list
-#'
-#' @param parameters the model parameters
-#' @param eq the equilibrium solution output
-#' @export
-parameterise_human_equilibrium <- function(parameters, eq) {
-  state_props <- colSums(eq$states[,c('S', 'D', 'A', 'U')])
-  parameters$s_proportion <- state_props[['S']]
-  parameters$d_proportion <- state_props[['D']]
-  parameters$a_proportion <- state_props[['A']]
-  parameters$u_proportion <- state_props[['U']]
-  parameters$init_ica <- eq$states[,'ICA']
-  parameters$init_icm <- eq$states[,'ICM']
-  parameters$init_ib <- eq$states[,'IB']
-  parameters$init_id <- eq$states[,'ID']
-  parameters$init_foim <- eq$FOIM
-  parameters
-}
-
 #' @title Parameterise total_M and carrying capacity for mosquitos from EIR
 #'
 #' @description NOTE: the inital EIR is likely to change unless the rest of the
diff --git a/R/processes.R b/R/processes.R
index 2e94da6c..249d3eb4 100644
--- a/R/processes.R
+++ b/R/processes.R
@@ -38,10 +38,10 @@ create_processes <- function(
   if (parameters$individual_mosquitoes) {
     processes <- c(
       processes,
-      create_mosquito_emergence_process_cpp(
+      create_mosquito_emergence_process(
         solvers,
-        variables$mosquito_state$.variable,
-        variables$species$.variable,
+        variables$mosquito_state,
+        variables$species,
         parameters$species,
         parameters$dpl
       )
@@ -54,6 +54,22 @@ create_processes <- function(
   # schedule infections for humans and set last_boosted_*
   # move mosquitoes into incubating state
   # kill mosquitoes caught in vector control
+  init_eir <- lapply(
+    seq_along(parameters$species),
+    function(species) {
+      LaggedValue$new(
+        parameters$de + 1,
+        calculate_eir(
+          species,
+          solvers,
+          variables,
+          parameters,
+          0
+        )
+      )
+    }
+  )
+
   processes <- c(
     processes,
     create_biting_process(
@@ -62,10 +78,38 @@ create_processes <- function(
       models,
       variables,
       events,
-      parameters
+      parameters,
+      LaggedValue$new(
+        parameters$delay_gam + 1,
+        parameters$init_foim
+      ), # lagged FOIM
+      init_eir # lagged EIR
     ),
-
-    create_mortality_process(variables, events, renderer, parameters)
+    create_mortality_process(variables, events, renderer, parameters),
+    create_progression_process(
+      events$asymptomatic_progression,
+      variables$state,
+      'D',
+      parameters$dd
+    ),
+    create_progression_process(
+      events$subpatent_progression,
+      variables$state,
+      'A',
+      parameters$da
+    ),
+    create_progression_process(
+      events$recovery,
+      variables$state,
+      'U',
+      parameters$du
+    ),
+    create_progression_process(
+      events$recovery,
+      variables$state,
+      'Tr',
+      parameters$dt
+    )
   )
 
   # ===============
diff --git a/R/utils.R b/R/utils.R
index 15a6c865..958426ad 100644
--- a/R/utils.R
+++ b/R/utils.R
@@ -1,5 +1,6 @@
 vnapply <- function(X, FUN, ...) vapply(X, FUN, ..., numeric(1))
 vlapply <- function(X, FUN, ...) vapply(X, FUN, ..., logical(1))
+vcapply <- function(X, FUN, ...) vapply(X, FUN, ..., character(1))
 
 #' @importFrom stats rbinom 
 bernoulli <- function(size, p) sample.int(size, rbinom(1, size, min(p, 1)))
@@ -19,19 +20,6 @@ log_uniform <- function(size, rate) -rate * log(runif(size))
 
 approx_sum <- function(X, n) abs(sum(X) - n) < sqrt(.Machine$double.eps)
 
-discretise_normal <- function(values, n_groups) {
-  quads <- statmod::gauss.quad.prob(n_groups, dist='normal')
-  breaks <- quads$nodes
-  breaks[[1]] <- min(min(values), min(breaks))
-  breaks <- c(breaks, max(max(values), max(breaks) + 1))
-  as.numeric(cut(
-    values,
-    breaks,
-    include.lowest = TRUE,
-    right = FALSE
-  ))
-}
-
 get_age <- function(birth_timesteps, current_timestep) {
   current_timestep - birth_timesteps
 }
diff --git a/R/variables.R b/R/variables.R
index e7e39fb5..9378a509 100644
--- a/R/variables.R
+++ b/R/variables.R
@@ -1,14 +1,3 @@
-initial_immunity <- function(parameter, age) {
-  if (length(parameter) == 1) {
-    return(rep(parameter, length(age)))
-  } else if (length(parameter) == 100) {
-    age <- floor(age / 365)
-    age[age > 99] <- 99
-    return(parameter[age + 1])
-  }
-  stop('unsupported immunity parameter')
-}
-
 #' @title Define model variables
 #' @description
 #' create_variables creates the human and mosquito variables for
@@ -56,27 +45,81 @@ initial_immunity <- function(parameter, age) {
 create_variables <- function(parameters) {
   size <- parameters$human_population
 
-  initial_age <- floor(rexp(size, rate=1/parameters$average_age))
+  initial_age <- round(rexp(size, rate=1/parameters$average_age))
 
-  initial_counts <- calculate_initial_counts(parameters)
+  if (parameters$enable_heterogeneity) {
+    #zeta_norm <- rnorm(size)
+    #het_nodes <- gaussian_quadrature(parameters$n_heterogeneity_groups)
+    quads <- statmod::gauss.quad.prob(
+      parameters$n_heterogeneity_groups,
+      dist='normal'
+    )
+    groups <- sample.int(
+      parameters$n_heterogeneity_groups,
+      size,
+      replace = TRUE,
+      prob = quads$weights
+    )
+    zeta_norm <- quads$nodes[groups]
+    zeta <- individual::DoubleVariable$new(
+      calculate_zeta(zeta_norm, parameters)
+    )
+    zeta_group <- individual::CategoricalVariable$new(
+      to_char_vector(seq(parameters$n_heterogeneity_groups)),
+      to_char_vector(groups)
+    )
+    if (!is.null(parameters$init_EIR)) {
+      eq <- calculate_eq(quads$nodes, parameters)
+    } else {
+      eq <- NULL
+    }
+  } else {
+    zeta <- individual::DoubleVariable$new(rep(1, size))
+    groups <- rep(1, size)
+    zeta_group <- individual::CategoricalVariable$new(
+      to_char_vector(seq(parameters$n_heterogeneity_groups)),
+      to_char_vector(groups)
+    )
+    if (!is.null(parameters$init_EIR)) {
+      eq <-	list(
+        malariaEquilibrium::human_equilibrium_no_het(
+          parameters$init_EIR,
+          sum(get_treatment_coverages(parameters, 0)),
+          parameters$eq_params,
+          EQUILIBRIUM_AGES
+        )
+      )
+    } else {
+      eq <- NULL
+    }
+  }
 
-  # Define variables
   states <- c('S', 'D', 'A', 'U', 'Tr')
   state <- individual::CategoricalVariable$new(
     states,
-    rep(states, times = initial_counts)
+    initial_state(parameters, initial_age, groups, eq)
   )
   birth <- individual::DoubleVariable$new(-initial_age)
-  last_boosted_ib <- individual::DoubleVariable$new(rep(-1, size) )
+  last_boosted_ib <- individual::DoubleVariable$new(rep(-1, size))
   last_boosted_ica <- individual::DoubleVariable$new(rep(-1, size))
   last_boosted_iva <- individual::DoubleVariable$new(rep(-1, size))
   last_boosted_id <- individual::DoubleVariable$new(rep(-1, size))
 
-  is_severe <- individual::CategoricalVariable$new(c('yes', 'no'), rep('no', size))
+  is_severe <- individual::CategoricalVariable$new(
+    c('yes', 'no'),
+    rep('no', size)
+  )
 
   # Maternal immunity
   icm <- individual::DoubleVariable$new(
-    initial_immunity(parameters$init_icm, initial_age)
+    initial_immunity(
+      parameters$init_icm,
+      initial_age,
+      groups,
+      eq,
+      parameters,
+      'ICM'
+    )
   )
 
   ivm <- individual::DoubleVariable$new(
@@ -85,11 +128,25 @@ create_variables <- function(parameters) {
 
   # Pre-erythoctic immunity
   ib  <- individual::DoubleVariable$new(
-    initial_immunity(parameters$init_ib, initial_age)
+    initial_immunity(
+      parameters$init_ib,
+      initial_age,
+      groups,
+      eq,
+      parameters,
+      'IB'
+    )
   )
   # Acquired immunity to clinical disease
   ica <- individual::DoubleVariable$new(
-    initial_immunity(parameters$init_ica, initial_age)
+    initial_immunity(
+      parameters$init_ica,
+      initial_age,
+      groups,
+      eq,
+      parameters,
+      'ICA'
+    )
   )
   # Acquired immunity to severe disease
   iva <- individual::DoubleVariable$new(
@@ -97,21 +154,16 @@ create_variables <- function(parameters) {
   )
   # Acquired immunity to detectability
   id <- individual::DoubleVariable$new(
-    initial_immunity(parameters$init_id, initial_age)
-  )
-
-  zeta_norm <- rnorm(size)
-  zeta <- individual::DoubleVariable$new(
-    exp(
-      zeta_norm * sqrt(parameters$sigma_squared) - parameters$sigma_squared/2
+    initial_immunity(
+      parameters$init_id,
+      initial_age,
+      groups,
+      eq,
+      parameters,
+      'ID'
     )
   )
 
-  zeta_group <- individual::CategoricalVariable$new(
-    to_char_vector(seq(parameters$n_heterogeneity_groups)),
-    to_char_vector(discretise_normal(zeta_norm, parameters$n_heterogeneity_groups))
-  )
-
   # Initialise infectiousness of humans -> mosquitoes
   # NOTE: not yet supporting initialisation of infectiousness of Treated individuals
   infectivity_values <- rep(0, parameters$human_population)
@@ -244,6 +296,53 @@ create_variables <- function(parameters) {
   variables
 }
 
+# =========
+# Utilities
+# =========
+
+initial_immunity <- function(
+  parameter,
+  age,
+  groups = NULL,
+  eq = NULL,
+  parameters = NULL,
+  eq_name = NULL
+  ) {
+  if (!is.null(eq)) {
+    age <- age / 365
+    return(vnapply(
+      seq_along(age),
+      function(i) {
+        g <- groups[[i]]
+        a <- age[[i]]
+        eq[[g]][which.max(a < eq[[g]][, 'age']), eq_name]
+      }
+    ))
+  }
+  rep(parameter, length(age))
+}
+
+initial_state <- function(parameters, age, groups, eq) {
+  ibm_states <- c('S', 'A', 'D', 'U', 'Tr')
+  if (!is.null(eq)) {
+    eq_states <- c('S', 'A', 'D', 'U', 'T')
+    age <- age / 365
+    return(vcapply(
+      seq_along(age),
+      function(i) {
+        g <- groups[[i]]
+        a <- age[[i]]
+        sample(
+          ibm_states,
+          size = 1,
+          prob = eq[[g]][which.max(a < eq[[g]][, 'age']), eq_states]
+        )
+      }
+    ))
+  }
+  rep(ibm_states, times = calculate_initial_counts(parameters))
+}
+
 calculate_initial_counts <- function(parameters) {
   initial_counts <- round(
     c(
@@ -258,3 +357,24 @@ calculate_initial_counts <- function(parameters) {
   initial_counts[[1]] <- initial_counts[[1]] + left_over
   initial_counts
 }
+
+calculate_eq <- function(het_nodes, parameters) {
+  ft <- sum(get_treatment_coverages(parameters, 0))
+	lapply(
+		het_nodes,
+		function(n) {
+			malariaEquilibrium::human_equilibrium_no_het(
+				parameters$init_EIR * calculate_zeta(n, parameters),
+				ft,
+				parameters$eq_params,
+        EQUILIBRIUM_AGES
+			)
+		}
+	)
+}
+
+calculate_zeta <- function(zeta_norm, parameters) {
+  exp(
+    zeta_norm * sqrt(parameters$sigma_squared) - parameters$sigma_squared/2
+  )
+}
diff --git a/man/create_mosquito_emergence_process_cpp.Rd b/man/create_mosquito_emergence_process_cpp.Rd
deleted file mode 100644
index 71d7a606..00000000
--- a/man/create_mosquito_emergence_process_cpp.Rd
+++ /dev/null
@@ -1,29 +0,0 @@
-% Generated by roxygen2: do not edit by hand
-% Please edit documentation in R/RcppExports.R
-\name{create_mosquito_emergence_process_cpp}
-\alias{create_mosquito_emergence_process_cpp}
-\title{Mosquito emergence process}
-\usage{
-create_mosquito_emergence_process_cpp(
-  solvers,
-  state,
-  species,
-  species_names,
-  dpl
-)
-}
-\arguments{
-\item{solvers}{a list of solver objects for each species of mosquito}
-
-\item{state}{the variable for the mosquito state}
-
-\item{species}{the variable for the mosquito species}
-
-\item{species_names}{a vector of category names for the species variable}
-
-\item{dpl}{the delay for pupal growth (in timesteps)}
-}
-\description{
-Move mosquitos from NonExistent to Sm in line with the number of
-pupals in the ODE models
-}
diff --git a/man/get_parameters.Rd b/man/get_parameters.Rd
index 6c040849..772c9bca 100644
--- a/man/get_parameters.Rd
+++ b/man/get_parameters.Rd
@@ -141,8 +141,9 @@ initial immunity values:
 
 incubation periods:
 \itemize{
-\item de - delay for infection
-\item dem - delay for infection in mosquitoes
+\item de - Duration of the human latent period of infection
+\item delay_gam - Lag from parasites to infectious gametocytes
+\item dem - Extrinsic incubation period in mosquito population model
 }
 
 vector biology:
@@ -254,6 +255,8 @@ miscellaneous:
 simulation
 \item individual_mosquitoes - boolean whether adult mosquitoes are modelled
 individually or compartmentaly
+\item enable_heterogeneity - boolean whether to include heterogeneity in biting
+rates
 }}
 }
 \description{
diff --git a/man/parameterise_human_equilibrium.Rd b/man/parameterise_human_equilibrium.Rd
deleted file mode 100644
index 15d25765..00000000
--- a/man/parameterise_human_equilibrium.Rd
+++ /dev/null
@@ -1,16 +0,0 @@
-% Generated by roxygen2: do not edit by hand
-% Please edit documentation in R/parameters.R
-\name{parameterise_human_equilibrium}
-\alias{parameterise_human_equilibrium}
-\title{Parameterise equilibrium proportions}
-\usage{
-parameterise_human_equilibrium(parameters, eq)
-}
-\arguments{
-\item{parameters}{the model parameters}
-
-\item{eq}{the equilibrium solution output}
-}
-\description{
-parameterise equilibrium proportions from a list
-}
diff --git a/man/remove_unused_jamie.Rd b/man/remove_unused_jamie.Rd
deleted file mode 100644
index cd1e9215..00000000
--- a/man/remove_unused_jamie.Rd
+++ /dev/null
@@ -1,16 +0,0 @@
-% Generated by roxygen2: do not edit by hand
-% Please edit documentation in R/compatibility.R
-\name{remove_unused_jamie}
-\alias{remove_unused_jamie}
-\title{remove parameter keys from jamie's format that are not used
-in this IBM}
-\usage{
-remove_unused_jamie(params)
-}
-\arguments{
-\item{params}{with keys in the jamie's format}
-}
-\description{
-remove parameter keys from jamie's format that are not used
-in this IBM
-}
diff --git a/man/set_equilibrium.Rd b/man/set_equilibrium.Rd
new file mode 100644
index 00000000..f41905bd
--- /dev/null
+++ b/man/set_equilibrium.Rd
@@ -0,0 +1,21 @@
+% Generated by roxygen2: do not edit by hand
+% Please edit documentation in R/compatibility.R
+\name{set_equilibrium}
+\alias{set_equilibrium}
+\title{Set equilibrium}
+\usage{
+set_equilibrium(parameters, init_EIR, eq_params = NULL)
+}
+\arguments{
+\item{parameters}{model parameters to update}
+
+\item{init_EIR}{the desired initial EIR}
+
+\item{eq_params}{parameters from the malariaEquilibrium package, if null.
+The default malariaEquilibrium parameters will be used}
+}
+\description{
+This will update the IBM parameters to match the
+equilibrium parameters and set up the initial human and mosquito population
+to acheive init_EIR
+}
diff --git a/man/translate_jamie.Rd b/man/translate_jamie.Rd
deleted file mode 100644
index 7481b61b..00000000
--- a/man/translate_jamie.Rd
+++ /dev/null
@@ -1,16 +0,0 @@
-% Generated by roxygen2: do not edit by hand
-% Please edit documentation in R/compatibility.R
-\name{translate_jamie}
-\alias{translate_jamie}
-\title{translate parameter keys from jamie's format to ones compatible
-with this IBM}
-\usage{
-translate_jamie(params)
-}
-\arguments{
-\item{params}{with keys in the jamie's format}
-}
-\description{
-translate parameter keys from jamie's format to ones compatible
-with this IBM
-}
diff --git a/src/RcppExports.cpp b/src/RcppExports.cpp
index 0cdfd5d8..ae0cebc9 100644
--- a/src/RcppExports.cpp
+++ b/src/RcppExports.cpp
@@ -46,6 +46,42 @@ BEGIN_RCPP
     return rcpp_result_gen;
 END_RCPP
 }
+// create_history
+Rcpp::XPtr<History> create_history(size_t size, double default_value);
+RcppExport SEXP _malariasimulation_create_history(SEXP sizeSEXP, SEXP default_valueSEXP) {
+BEGIN_RCPP
+    Rcpp::RObject rcpp_result_gen;
+    Rcpp::RNGScope rcpp_rngScope_gen;
+    Rcpp::traits::input_parameter< size_t >::type size(sizeSEXP);
+    Rcpp::traits::input_parameter< double >::type default_value(default_valueSEXP);
+    rcpp_result_gen = Rcpp::wrap(create_history(size, default_value));
+    return rcpp_result_gen;
+END_RCPP
+}
+// history_at
+double history_at(Rcpp::XPtr<History> history, double timestep);
+RcppExport SEXP _malariasimulation_history_at(SEXP historySEXP, SEXP timestepSEXP) {
+BEGIN_RCPP
+    Rcpp::RObject rcpp_result_gen;
+    Rcpp::RNGScope rcpp_rngScope_gen;
+    Rcpp::traits::input_parameter< Rcpp::XPtr<History> >::type history(historySEXP);
+    Rcpp::traits::input_parameter< double >::type timestep(timestepSEXP);
+    rcpp_result_gen = Rcpp::wrap(history_at(history, timestep));
+    return rcpp_result_gen;
+END_RCPP
+}
+// history_push
+void history_push(Rcpp::XPtr<History> history, double value, double timestep);
+RcppExport SEXP _malariasimulation_history_push(SEXP historySEXP, SEXP valueSEXP, SEXP timestepSEXP) {
+BEGIN_RCPP
+    Rcpp::RNGScope rcpp_rngScope_gen;
+    Rcpp::traits::input_parameter< Rcpp::XPtr<History> >::type history(historySEXP);
+    Rcpp::traits::input_parameter< double >::type value(valueSEXP);
+    Rcpp::traits::input_parameter< double >::type timestep(timestepSEXP);
+    history_push(history, value, timestep);
+    return R_NilValue;
+END_RCPP
+}
 // carrying_capacity
 double carrying_capacity(const size_t timestep, const bool model_seasonality, const double g0, const std::vector<double>& g, const std::vector<double>& h, const double K0, const double R_bar);
 RcppExport SEXP _malariasimulation_carrying_capacity(SEXP timestepSEXP, SEXP model_seasonalitySEXP, SEXP g0SEXP, SEXP gSEXP, SEXP hSEXP, SEXP K0SEXP, SEXP R_barSEXP) {
@@ -90,21 +126,6 @@ BEGIN_RCPP
     return rcpp_result_gen;
 END_RCPP
 }
-// create_mosquito_emergence_process_cpp
-Rcpp::XPtr<process_t> create_mosquito_emergence_process_cpp(Rcpp::List solvers, Rcpp::XPtr<CategoricalVariable> state, Rcpp::XPtr<CategoricalVariable> species, std::vector<std::string> species_names, double dpl);
-RcppExport SEXP _malariasimulation_create_mosquito_emergence_process_cpp(SEXP solversSEXP, SEXP stateSEXP, SEXP speciesSEXP, SEXP species_namesSEXP, SEXP dplSEXP) {
-BEGIN_RCPP
-    Rcpp::RObject rcpp_result_gen;
-    Rcpp::RNGScope rcpp_rngScope_gen;
-    Rcpp::traits::input_parameter< Rcpp::List >::type solvers(solversSEXP);
-    Rcpp::traits::input_parameter< Rcpp::XPtr<CategoricalVariable> >::type state(stateSEXP);
-    Rcpp::traits::input_parameter< Rcpp::XPtr<CategoricalVariable> >::type species(speciesSEXP);
-    Rcpp::traits::input_parameter< std::vector<std::string> >::type species_names(species_namesSEXP);
-    Rcpp::traits::input_parameter< double >::type dpl(dplSEXP);
-    rcpp_result_gen = Rcpp::wrap(create_mosquito_emergence_process_cpp(solvers, state, species, species_names, dpl));
-    return rcpp_result_gen;
-END_RCPP
-}
 // create_mosquito_model
 Rcpp::XPtr<MosquitoModel> create_mosquito_model(double beta, double de, double mue, double K0, double gamma, double dl, double mul, double dp, double mup, size_t total_M, bool model_seasonality, double g0, std::vector<double> g, std::vector<double> h, double R_bar);
 RcppExport SEXP _malariasimulation_create_mosquito_model(SEXP betaSEXP, SEXP deSEXP, SEXP mueSEXP, SEXP K0SEXP, SEXP gammaSEXP, SEXP dlSEXP, SEXP mulSEXP, SEXP dpSEXP, SEXP mupSEXP, SEXP total_MSEXP, SEXP model_seasonalitySEXP, SEXP g0SEXP, SEXP gSEXP, SEXP hSEXP, SEXP R_barSEXP) {
@@ -194,10 +215,12 @@ static const R_CallMethodDef CallEntries[] = {
     {"_malariasimulation_create_adult_mosquito_model", (DL_FUNC) &_malariasimulation_create_adult_mosquito_model, 4},
     {"_malariasimulation_adult_mosquito_model_update", (DL_FUNC) &_malariasimulation_adult_mosquito_model_update, 5},
     {"_malariasimulation_create_adult_solver", (DL_FUNC) &_malariasimulation_create_adult_solver, 2},
+    {"_malariasimulation_create_history", (DL_FUNC) &_malariasimulation_create_history, 2},
+    {"_malariasimulation_history_at", (DL_FUNC) &_malariasimulation_history_at, 2},
+    {"_malariasimulation_history_push", (DL_FUNC) &_malariasimulation_history_push, 3},
     {"_malariasimulation_carrying_capacity", (DL_FUNC) &_malariasimulation_carrying_capacity, 7},
     {"_malariasimulation_eggs_laid", (DL_FUNC) &_malariasimulation_eggs_laid, 3},
     {"_malariasimulation_rainfall", (DL_FUNC) &_malariasimulation_rainfall, 4},
-    {"_malariasimulation_create_mosquito_emergence_process_cpp", (DL_FUNC) &_malariasimulation_create_mosquito_emergence_process_cpp, 5},
     {"_malariasimulation_create_mosquito_model", (DL_FUNC) &_malariasimulation_create_mosquito_model, 15},
     {"_malariasimulation_mosquito_model_update", (DL_FUNC) &_malariasimulation_mosquito_model_update, 4},
     {"_malariasimulation_create_solver", (DL_FUNC) &_malariasimulation_create_solver, 2},
diff --git a/src/history.cpp b/src/history.cpp
new file mode 100644
index 00000000..43fa7c47
--- /dev/null
+++ b/src/history.cpp
@@ -0,0 +1,72 @@
+/*
+ * history.cpp
+ *
+ *  Created on: 08 Apr 2021
+ *      Author: gc1610
+ */
+
+#include <Rcpp.h>
+#include <algorithm>
+#include "history.h"
+
+History::History() : max_size(-1), has_default(false) {}
+History::History(size_t max_size) : max_size(max_size), has_default(false) {}
+History::History(size_t max_size, double default_value)
+    : max_size(max_size), has_default(true), default_value(default_value) {}
+
+void History::push(double value, double time) {
+    values.insert({time, value});
+    if (max_size != -1) {
+        while(values.size() > max_size) {
+            values.erase(values.begin());
+        }
+    }
+}
+
+double History::at(double time) const {
+    auto it = values.lower_bound(time);
+    if (it == values.end()) {
+        if (has_default) {
+            return default_value;
+        }
+        Rcpp::stop("`time` is later than the stored history");
+    }
+
+    // Check if we've landed on the exact time
+    if (it->first == time) {
+        return it->second;
+    }
+
+    // Find the element before
+    auto after_element = *it;
+    while(it->first > time) {
+        // Check if we're at the start of the history
+        if (it == values.begin()) {
+            if (has_default) {
+                return default_value;
+            }
+            Rcpp::stop("`time` is before our stored history");
+        }
+        it--;
+    }
+
+    // Interpolate
+    return it->second + (
+        (time - it->first) / (after_element.first - it->first)
+    ) * (after_element.second - it->second);
+}
+
+//[[Rcpp::export]]
+Rcpp::XPtr<History> create_history(size_t size, double default_value) {
+    return Rcpp::XPtr<History>(new History(size, default_value), true);
+}
+
+//[[Rcpp::export]]
+double history_at(Rcpp::XPtr<History> history, double timestep) {
+    return history->at(timestep);
+}
+
+//[[Rcpp::export]]
+void history_push(Rcpp::XPtr<History> history, double value, double timestep) {
+    return history->push(value, timestep);
+}
diff --git a/src/history.h b/src/history.h
new file mode 100644
index 00000000..5a5fd618
--- /dev/null
+++ b/src/history.h
@@ -0,0 +1,28 @@
+/*
+ * history.h
+ *
+ *  Created on: 08 Apr 2021
+ *      Author: gc1610
+ */
+
+#ifndef SRC_HISTORY_H_
+#define SRC_HISTORY_H_
+
+#include <map>
+
+class History {
+private:
+    std::map<double,double> values;
+    size_t max_size;
+    void clean();
+    bool has_default;
+    double default_value;
+public:
+    History();
+    History(size_t);
+    History(size_t, double);
+    void push(double, double);
+    double at(double) const;
+};
+
+#endif /* SRC_HISTORY_H_ */
diff --git a/src/malariasimulation_types.h b/src/malariasimulation_types.h
index ca5c8f94..5c00ddbd 100644
--- a/src/malariasimulation_types.h
+++ b/src/malariasimulation_types.h
@@ -12,5 +12,6 @@
 #include "mosquito_ode.h"
 #include "adult_mosquito_ode.h"
 #include "solver.h"
+#include "history.h"
 
 #endif /* MALARIASIMULATION_TYPES_H_ */
diff --git a/src/mosquito_emergence.cpp b/src/mosquito_emergence.cpp
deleted file mode 100644
index 937fad47..00000000
--- a/src/mosquito_emergence.cpp
+++ /dev/null
@@ -1,66 +0,0 @@
-/*
- * mosquito_emergence.cpp
- *
- *  Created on: 4 Aug 2020
- *      Author: gc1610
- */
-
-#include "mosquito_emergence.h"
-#include "mosquito_ode.h"
-#include <sstream>
-
-//' @title Mosquito emergence process
-//' @description Move mosquitos from NonExistent to Sm in line with the number of
-//' pupals in the ODE models
-//'
-//' @param solvers a list of solver objects for each species of mosquito
-//' @param state the variable for the mosquito state
-//' @param species the variable for the mosquito species
-//' @param species_names a vector of category names for the species variable
-//' @param dpl the delay for pupal growth (in timesteps)
-//[[Rcpp::export]]
-Rcpp::XPtr<process_t> create_mosquito_emergence_process_cpp(
-    Rcpp::List solvers,
-    Rcpp::XPtr<CategoricalVariable> state,
-    Rcpp::XPtr<CategoricalVariable> species,
-    std::vector<std::string> species_names,
-    double dpl
-    ) {
-    auto rate = .5 * 1./dpl;
-    return Rcpp::XPtr<process_t>(
-        new process_t([=] (size_t t) {
-            auto n = 0u;
-            for (Rcpp::XPtr<Solver> solver: solvers) {
-                n += solver->state[get_idx(ODEState::P)] * rate;
-            }
-            auto source = state->get_index_of(std::vector<std::string>{"NonExistent"});
-            if (source.size() < n) {
-                std::stringstream m;
-                m << "Not enough mosquitos (short by ";
-                m << n - source.size();
-                m << "). Please raise the value of parameters$mosquito_limit. ";
-                m << "If you have used `parameterise_mosquito_equilibrium`, ";
-                m << "your seasonality parameters may have lead to more mosquitoes than expected.";
-                Rcpp::stop(m.str());
-            }
-            if (n > 0) {
-                auto species_i = 0u;
-                auto sourceit = source.begin();
-                for (Rcpp::XPtr<Solver> solver : solvers) {
-                    auto to_set = static_cast<size_t>(
-                        solver->state[get_idx(ODEState::P)] * rate
-                    );
-                    auto target = individual_index_t(state->size);
-                    for (auto i = 0u; i < to_set; ++i) {
-                        target.insert(*sourceit);
-                        ++sourceit;
-                    }
-                    state->queue_update("Sm", target);
-                    species->queue_update(species_names[species_i], target);
-                    ++species_i;
-                }
-            }
-        }),
-        true
-    );
-}
diff --git a/src/mosquito_emergence.h b/src/mosquito_emergence.h
deleted file mode 100644
index 42d67fef..00000000
--- a/src/mosquito_emergence.h
+++ /dev/null
@@ -1,21 +0,0 @@
-/*
- * mosquito_emergence.h
- *
- *  Created on: 4 Aug 2020
- *      Author: gc1610
- */
-
-#ifndef SRC_MOSQUITO_EMERGENCE_H_
-#define SRC_MOSQUITO_EMERGENCE_H_
-
-#include <individual.h>
-
-Rcpp::XPtr<process_t> create_mosquito_emergence_process_cpp(
-    Rcpp::List,
-    Rcpp::XPtr<CategoricalVariable>,
-    Rcpp::XPtr<CategoricalVariable>,
-    std::vector<std::string>,
-    double
-);
-
-#endif /* SRC_MOSQUITO_EMERGENCE_H_ */
diff --git a/src/test-emergence.cpp b/src/test-emergence.cpp
deleted file mode 100644
index 8d5b1366..00000000
--- a/src/test-emergence.cpp
+++ /dev/null
@@ -1,148 +0,0 @@
-// All test files should include the <testthat.h>
-// header file.
-
-#include "solver.h"
-#include <testthat.h>
-#include "test-mock.h"
-#include "mosquito_emergence.h"
-
-context("Emergence works") {
-
-    test_that("emergence process fails when there are not enough individuals") {
-        //mock a solver
-        Solver solver(state_t{1000, 500, 100}, mock_integration);
-
-        //mock my variables
-        auto state_values = std::vector<std::string>(2000);
-        for (auto n = 0; n < 2000; ++n) {
-            state_values[n] = "Sm";
-        }
-        MockCategory state(
-            {"Sm", "NonExistent"},
-            state_values
-        );
-        MockCategory species(
-            {"gamb"},
-            std::vector<std::string>(2000, "gamb")
-        );
-
-        auto process = create_mosquito_emergence_process_cpp(
-            Rcpp::List::create(Rcpp::XPtr<Solver>(&solver, false)),
-            Rcpp::XPtr<CategoricalVariable>(&state, false),
-            Rcpp::XPtr<CategoricalVariable>(&species, false),
-            {"gamb"},
-            2
-        );
-        CATCH_REQUIRE_THROWS((*process)(1));
-    }
-
-    test_that("emergence process fails at the boundary") {
-        //mock a solver
-        Solver solver(state_t{1000, 500, 100}, mock_integration);
-
-        //mock my variables
-        auto state_values = std::vector<std::string>(2000);
-        for (auto n = 0; n < 2000; ++n) {
-            state_values[n] = "Sm";
-        }
-        for (auto n = 1000; n < 1024; ++n) {
-            state_values[n] = "NonExistent";
-        }
-        MockCategory state(
-            {"Sm", "NonExistent"},
-            state_values
-        );
-        MockCategory species(
-            {"gamb"},
-            std::vector<std::string>(2000, "gamb")
-        );
-
-        auto process = create_mosquito_emergence_process_cpp(
-            Rcpp::List::create(Rcpp::XPtr<Solver>(&solver, false)),
-            Rcpp::XPtr<CategoricalVariable>(&state, false),
-            Rcpp::XPtr<CategoricalVariable>(&species, false),
-            {"gamb"},
-            2
-        );
-        CATCH_REQUIRE_THROWS((*process)(1));
-    }
-
-    test_that("emergence process creates the correct number of susceptibles") {
-        //mock a solver
-        Solver solver(state_t{1000, 500, 100}, mock_integration);
-
-        //mock my variables
-        auto state_values = std::vector<std::string>(2000);
-        for (auto n = 0; n < 1000; ++n) {
-            state_values[n] = "Sm";
-        }
-        for (auto n = 1000; n < 2000; ++n) {
-            state_values[n] = "NonExistent";
-        }
-        MockCategory state(
-            {"Sm", "NonExistent"},
-            state_values
-        );
-        MockCategory species(
-            {"gamb"},
-            std::vector<std::string>(2000, "gamb")
-        );
-
-        auto process = create_mosquito_emergence_process_cpp(
-            Rcpp::List::create(Rcpp::XPtr<Solver>(&solver, false)),
-            Rcpp::XPtr<CategoricalVariable>(&state, false),
-            Rcpp::XPtr<CategoricalVariable>(&species, false),
-            {"gamb"},
-            2
-        );
-        
-        auto expected_to_emerge = 25u;
-        individual_index_t target(2000);
-        for (auto i = 0ul; i < expected_to_emerge; ++i) {
-            target.insert(i + 1000);
-        }
-
-        REQUIRE_CALL(state, queue_update("Sm", target));
-        REQUIRE_CALL(species, queue_update("gamb", target));
-        (*process)(1);
-    }
-
-    test_that("emergence process works at the boundary") {
-        //mock a solver
-        Solver solver(state_t{1000, 500, 100}, mock_integration);
-
-        //mock my variables
-        auto state_values = std::vector<std::string>(2000);
-        for (auto n = 0; n < 2000; ++n) {
-            state_values[n] = "Sm";
-        }
-        auto expected_to_emerge = 25u;
-        for (auto n = 1000; n < 1000 + expected_to_emerge; ++n) {
-            state_values[n] = "NonExistent";
-        }
-        MockCategory state(
-            {"Sm", "NonExistent"},
-            state_values
-        );
-        MockCategory species(
-            {"gamb"},
-            std::vector<std::string>(2000, "gamb")
-        );
-
-        auto process = create_mosquito_emergence_process_cpp(
-            Rcpp::List::create(Rcpp::XPtr<Solver>(&solver, false)),
-            Rcpp::XPtr<CategoricalVariable>(&state, false),
-            Rcpp::XPtr<CategoricalVariable>(&species, false),
-            {"gamb"},
-            2
-        );
-        individual_index_t target(2000);
-        for (auto i = 0ul; i < expected_to_emerge; ++i) {
-            target.insert(i + 1000);
-        }
-
-        REQUIRE_CALL(state, queue_update("Sm", target));
-        REQUIRE_CALL(species, queue_update("gamb", target));
-        (*process)(1);
-    }
-}
diff --git a/src/test_history.cpp b/src/test_history.cpp
new file mode 100644
index 00000000..644680de
--- /dev/null
+++ b/src/test_history.cpp
@@ -0,0 +1,47 @@
+/*
+ * test-history.cpp
+ *
+ *  Created on: 08 Apr 2021
+ *      Author: gc1610
+ */
+
+#include <testthat.h>
+#include "history.h"
+
+context("History") {
+    test_that("Can get history at exact time point") {
+        History h;
+        h.push(1., 1.);
+        h.push(1.5, 4.);
+        expect_true(h.at(1.) == 1.);
+    }
+
+    test_that("Can get history at interpolated time point") {
+        History h;
+        h.push(0., 1.);
+        h.push(4., 2.);
+        expect_true(h.at(1.5) == 2.);
+    }
+
+    test_that("Errors gracefully when I search before the history") {
+        History h;
+        h.push(0., 1.);
+        h.push(4., 2.);
+        expect_error(h.at(0.));
+    }
+
+    test_that("Errors gracefully when I search after the history") {
+        History h;
+        h.push(0., 1.);
+        h.push(4., 2.);
+        expect_error(h.at(4.));
+    }
+
+    test_that("Can restrict size of history") {
+        History h(2);
+        h.push(3, 2.5);
+        h.push(0., 1.);
+        h.push(4., 2.);
+        expect_error(h.at(1));
+    }
+}
diff --git a/tests/testthat/helper-integration.R b/tests/testthat/helper-integration.R
index 9be9ba25..e60a800c 100644
--- a/tests/testthat/helper-integration.R
+++ b/tests/testthat/helper-integration.R
@@ -22,7 +22,9 @@ mock_category <- function(...) {
   v <- individual::CategoricalVariable$new(...)
   list(
     get_index_of = v$get_index_of,
-    queue_update = mockery::mock()
+    get_size_of = v$get_size_of,
+    queue_update = mockery::mock(),
+    get_categories = v$get_categories
   )
 }
 
diff --git a/tests/testthat/test-biology.R b/tests/testthat/test-biology.R
index d57aa1ff..9ffe4415 100644
--- a/tests/testthat/test-biology.R
+++ b/tests/testthat/test-biology.R
@@ -1,50 +1,24 @@
 test_that('total_M and EIR functions are consistent with equilibrium EIR', {
   EIR <- 5
-  jamie_parameters <- malariaEquilibrium::load_parameter_set()
-  h_eq <- malariaEquilibrium::human_equilibrium_no_het(
-    EIR,
-    0,
-    jamie_parameters,
-    0:99
-  )
-  foim <- sum(h_eq[,'inf'] * h_eq[,'psi'])
   population <- 100000
-  parameters <- get_parameters(c(
-    translate_jamie(remove_unused_jamie(jamie_parameters)),
+  parameters <- get_parameters(
     list(
-      init_foim = foim,
-      species = 'all',
-      species_proportions = 1,
-      human_population = population
+      human_population = population,
+      enable_heterogeneity = FALSE
     )
-  ))
-  parameters <- parameterise_mosquito_equilibrium(parameters, EIR)
-  m_eq <- initial_mosquito_counts(parameters, 1, foim, parameters$total_M)
+  )
+  parameters <- set_equilibrium(parameters, EIR)
 
   #set up arguments for EIR calculation
   variables <- create_variables(parameters)
-  age <- get_age(variables$birth$get_values(), 0)
-  zeta <- rep(1, length(age))
-  variables$zeta <- individual::DoubleVariable$new(zeta)
-  p_bitten <- prob_bitten(0, variables, 1, parameters)
-  psi <- unique_biting_rate(age, parameters)
-  .pi <- human_pi(zeta, psi)
-  Z <- average_p_repelled(p_bitten$prob_repelled, .pi, parameters$Q0)
-  W <- average_p_successful(p_bitten$prob_bitten_survives, .pi, parameters$Q0)
-  f <- blood_meal_rate(1, Z, parameters)
-  estimated_eir <- m_eq[[6]] * effective_biting_rate(
-    .pi,
-    age,
-    1,
-    p_bitten,
-    f,
-    W,
-    parameters
-  )
+  vector_models <- parameterise_mosquito_models(parameters)
+  solvers <- parameterise_solvers(vector_models, parameters)
+  omega <- sum(unique_biting_rate(-variables$birth$get_values(), parameters))
+  estimated_eir <- calculate_eir(1, solvers, variables, parameters, 0) / omega
   expect_equal(
     mean(estimated_eir) * 365,
     EIR,
-    tolerance = 1
+    tolerance = 1e-2
   )
 })
 
@@ -53,53 +27,147 @@ test_that('total_M and EIR functions are consistent with equilibrium EIR (with h
 
   EIR <- 50
 
-  jamie_parameters <- malariaEquilibrium::load_parameter_set()
+  parameters <- get_parameters(list(human_population = population))
+  parameters <- set_equilibrium(parameters, EIR)
 
-  h_eq <- malariaEquilibrium::human_equilibrium(
+  variables <- create_variables(parameters)
+  vector_models <- parameterise_mosquito_models(parameters)
+  solvers <- parameterise_solvers(vector_models, parameters)
+  omega <- sum(unique_biting_rate(-variables$birth$get_values(), parameters))
+  estimated_eir <- calculate_eir(1, solvers, variables, parameters, 0) / omega
+  expect_equal(
+    mean(estimated_eir) * 365,
+    EIR,
+    tolerance = 1
+  )
+})
+
+test_that('FOIM is consistent with equilibrium', {
+  population <- 1e5
+
+  EIR <- 5
+
+  eq_params <- malariaEquilibrium::load_parameter_set()
+
+  ages <- 0:999 / 10
+  foim <- malariaEquilibrium::human_equilibrium(
     EIR,
     0,
-    jamie_parameters,
-    0:99
+    eq_params,
+    ages
+  )$FOIM
+
+  parameters <- get_parameters(c(list(human_population = population)))
+  parameters <- set_equilibrium(parameters, EIR)
+
+  variables <- create_variables(parameters)
+  vector_models <- parameterise_mosquito_models(parameters)
+  solvers <- parameterise_solvers(vector_models, parameters)
+  lambda <- effective_biting_rates(1, variables, parameters, 0)
+  expect_equal(
+    calculate_foim(variables$infectivity$get_values(), lambda),
+    foim,
+    tolerance=1e-2
   )
-  foim <- h_eq$FOIM
+})
 
-  parameters <- get_parameters(c(
-    translate_jamie(remove_unused_jamie(jamie_parameters)),
-    list(
-      init_foim = foim,
-      species = 'all',
-      species_proportions = 1,
-      human_population = population
-    )
-  ))
-  parameters <- parameterise_mosquito_equilibrium(parameters, EIR)
-  m_eq <- initial_mosquito_counts(parameters, 1, foim, parameters$total_M)
+test_that('FOI is consistent with equilibrium', {
+  population <- 1e5
+
+  EIR <- 5
+
+  eq_params <- malariaEquilibrium::load_parameter_set()
+
+  ages <- 0:999 / 10
+  eq <- malariaEquilibrium::human_equilibrium(
+    EIR,
+    0,
+    eq_params,
+    ages
+  )
+  foi <- weighted.mean(eq$states[,'FOI'], eq$states[,'prop'])
+
+  parameters <- get_parameters(list(human_population = population))
+  parameters <- set_equilibrium(parameters, EIR)
 
-  #set up arguments for EIR calculation
   variables <- create_variables(parameters)
-  age <- get_age(variables$birth$get_values(), 0)
-  zeta <- variables$zeta$get_values()
-  p_bitten <- prob_bitten(0, variables, 1, parameters)
-  psi <- unique_biting_rate(age, parameters)
-  .pi <- human_pi(zeta, psi)
-  Z <- average_p_repelled(p_bitten$prob_repelled, .pi, parameters$Q0)
-  W <- average_p_successful(p_bitten$prob_bitten_survives, .pi, parameters$Q0)
-  f <- blood_meal_rate(1, Z, parameters)
-  estimated_eir <- m_eq[[6]] * effective_biting_rate(
-    .pi,
-    age,
-    1,
-    p_bitten,
-    f,
-    W,
-    parameters
+  lambda <- effective_biting_rates(1, variables, parameters, 0)
+  ib <- blood_immunity(variables$ib$get_values(), parameters)
+  expect_equal(
+    mean(lambda * ib),
+    foi,
+    tolerance=1e-2
   )
+})
+
+test_that('phi is consistent with equilibrium at high EIR (no het)', {
+  population <- 1e5
+
+  EIR <- 100
+  ages <- 0:999 / 10
+
+  parameters <- get_parameters(list(
+    human_population = population,
+    enable_heterogeneity = FALSE
+  ))
+  parameters <- set_equilibrium(parameters, EIR)
+  eq_params <- malariaEquilibrium::load_parameter_set()
 
+  eq <-  malariaEquilibrium::human_equilibrium_no_het(
+    EIR,
+    0,
+    eq_params,
+    ages
+  )
+
+  variables <- create_variables(parameters)
   expect_equal(
-    mean(estimated_eir) * 365,
+    mean(
+      clinical_immunity(
+        variables$ica$get_values(),
+        variables$icm$get_values(),
+        parameters
+      )
+    ),
+    weighted.mean(eq[,'phi'], eq[,'prop']),
+    tolerance=1e-2
+  )
+})
+
+test_that('phi is consistent with equilibrium at high EIR', {
+  population <- 1e5
+
+  EIR <- 100
+  ages <- 0:999 / 10
+
+  parameters <- get_parameters(list(human_population = population))
+  parameters <- set_equilibrium(parameters, EIR)
+  eq_params <- malariaEquilibrium::load_parameter_set()
+
+  eq <-  malariaEquilibrium::human_equilibrium(
     EIR,
-    tolerance = 1
+    0,
+    eq_params,
+    ages
   )
+
+  variables <- create_variables(parameters)
+  expect_equal(
+    mean(
+      clinical_immunity(
+        variables$ica$get_values(),
+        variables$icm$get_values(),
+        parameters
+      )
+    ),
+    weighted.mean(eq$states[,'phi'], eq$states[,'prop']),
+    tolerance=1e-2
+  )
+})
+
+test_that('mosquito_limit is set to 1 for 0 EIR', {
+  parameters <- parameterise_mosquito_equilibrium(get_parameters(), 0)
+  expect_equal(parameters$mosquito_limit, 1)
 })
 
 test_that('mosquito_limit is set to 1 for 0 EIR', {
diff --git a/tests/testthat/test-biting-integration.R b/tests/testthat/test-biting-integration.R
index 585056be..3b06d5d5 100644
--- a/tests/testthat/test-biting-integration.R
+++ b/tests/testthat/test-biting-integration.R
@@ -9,14 +9,20 @@ test_that('biting_process integrates mosquito effects and human infection', {
   events <- mockery::mock()
   age <- c(20, 24, 5, 39) * 365
   variables <- list(birth = individual::DoubleVariable$new((-age + timestep)))
+  lagged_foim <- LaggedValue$new(1, 1)
+  lagged_eir <- LaggedValue$new(1, 1)
+  models <- parameterise_mosquito_models(parameters)
+  solvers <- parameterise_solvers(models, parameters)
 
   biting_process <- create_biting_process(
     renderer,
-    list(), #solvers
-    list(), #models
+    solvers,
+    models,
     variables,
     events,
-    parameters
+    parameters,
+    lagged_foim,
+    lagged_eir
   )
 
   bitten <- individual::Bitset$new(parameters$human_population)
@@ -31,13 +37,15 @@ test_that('biting_process integrates mosquito effects and human infection', {
     bites_mock,
     1,
     renderer,
-    list(),
-    list(),
+    solvers,
+    models,
     variables,
     events,
     age,
     parameters,
-    timestep
+    timestep,
+    lagged_foim,
+    lagged_eir
   )
 
   mockery::expect_args(
@@ -77,53 +85,42 @@ test_that('simulate_bites integrates eir calculation and mosquito side effects',
     rep('All', 100)
   )
 
-  lambda_mock <- mockery::mock(.5, .5, .5)
+  lambda_mock <- mockery::mock(c(.5, .5, .5, .5))
   mosquito_effects_mock <- mockery::mock()
   ode_update <- mockery::mock()
+  sample_mock <- mockery::mock(c(2, 3))
+  pois_mock <- mockery::mock(2)
 
-  mockery::stub(simulate_bites, 'effective_biting_rate', lambda_mock)
   mockery::stub(simulate_bites, 'biting_effects_individual', mosquito_effects_mock)
-  mockery::stub(simulate_bites, 'rpois', mockery::mock(2))
-  mockery::stub(simulate_bites, 'sample.int', mockery::mock(c(2, 3)))
-  .pi <- rep(1 / population, population)
-  mockery::stub(simulate_bites, 'human_pi', mockery::mock(.pi))
+  mockery::stub(simulate_bites, 'rpois', pois_mock)
+  mockery::stub(simulate_bites, 'sample.int', sample_mock)
+  mockery::stub(simulate_bites, '.effective_biting_rates', lambda_mock)
   mockery::stub(simulate_bites, 'mosquito_model_update', ode_update)
-  ode_model <- mockery::mock()
+  models <- parameterise_mosquito_models(parameters)
+  solvers <- parameterise_solvers(models, parameters)
+  lagged_foim <- LaggedValue$new(12.5, .001)
+  lagged_eir <- list(LaggedValue$new(12, 10))
   bitten <- simulate_bites(
     renderer,
-    list(), #solvers
-    list(ode_model), #models
+    solvers,
+    models,
     variables,
     events,
     age,
     parameters,
-    timestep
+    timestep,
+    lagged_foim,
+    lagged_eir
   )
 
   expect_equal(bitten$to_vector(), c(2, 3))
 
   f <- parameters$blood_meal_rates[[1]]
 
-  mockery::expect_args(
-    lambda_mock,
-    1,
-    .pi,
-    age,
-    1,
-    list(
-      prob_bitten_survives = rep(1, population),
-      prob_bitten = rep(1, population),
-      prob_repelled = rep(0, population)
-    ),
-    f,
-    1,
-    parameters
-  )
-
   effects_args <- mockery::mock_args(mosquito_effects_mock)
 
   expect_equal(effects_args[[1]][[1]], variables)
-  expect_equal(effects_args[[1]][[2]], .55)
+  expect_equal(effects_args[[1]][[2]], .001)
   expect_equal(effects_args[[1]][[3]], events)
   expect_equal(effects_args[[1]][[4]], 1)
   expect_equal(effects_args[[1]][[5]]$to_vector(), 11:25)
@@ -132,5 +129,14 @@ test_that('simulate_bites integrates eir calculation and mosquito side effects',
   expect_equal(effects_args[[1]][[8]], parameters)
   expect_equal(effects_args[[1]][[9]], timestep)
 
-  mockery::expect_args(ode_update, 1, ode_model, 25, f, parameters$mum)
+  mockery::expect_args(ode_update, 1, models[[1]], 25, f, parameters$mum)
+  mockery::expect_args(pois_mock, 1, 1, 10)
+  mockery::expect_args(
+    sample_mock,
+    1,
+    parameters$human_population,
+    2,
+    replace=TRUE,
+    prob=c(.5, .5, .5, .5)
+  )
 })
diff --git a/tests/testthat/test-compatibility.R b/tests/testthat/test-compatibility.R
new file mode 100644
index 00000000..5aa44e90
--- /dev/null
+++ b/tests/testthat/test-compatibility.R
@@ -0,0 +1,10 @@
+test_that('default params can be translated to equilibrium', {
+  malariaEquilibrium::human_equilibrium(
+    5,
+    0,
+    translate_parameters(get_parameters()),
+    0:100
+  )
+
+  expect_true(TRUE)
+})
diff --git a/tests/testthat/test-emergence-integration.R b/tests/testthat/test-emergence-integration.R
new file mode 100644
index 00000000..400e0e18
--- /dev/null
+++ b/tests/testthat/test-emergence-integration.R
@@ -0,0 +1,80 @@
+test_that('emergence process fails when there are not enough individuals', {
+  parameters <- get_parameters()
+  state <- individual::CategoricalVariable$new(
+    c('Sm', 'Pm', 'Im', 'NonExistent'),
+    c(rep('Im', 1000), rep('Sm', 1000))
+  )
+  species <- individual::CategoricalVariable$new(
+    c('All'),
+    rep('All', 2000)
+  )
+  emergence_process <- create_mosquito_emergence_process(
+    list(),
+    state,
+    species,
+    c('All'),
+    parameters$dpl
+  )
+  mockery::stub(
+    emergence_process,
+    'solver_get_states', 
+    mockery::mock(
+      c(1000, 500, 100),
+      c(1000, 500, 100)
+    )
+  )
+  expect_error(emergence_process(0), '*')
+})
+
+test_that('emergence_process creates the correct number of susceptables', {
+  parameters <- get_parameters()
+  state <- mock_category(
+    c('Sm', 'Pm', 'Im', 'NonExistent'),
+    c(rep('Im', 1000), rep('Sm', 1000), rep('NonExistent', 100000))
+  )
+  species <- mock_category(
+    c('a', 'b'),
+    c('a', 'b')
+  )
+  emergence_process <- create_mosquito_emergence_process(
+    list(mockery::mock(), mockery::mock()),
+    state,
+    species,
+    c('a', 'b'),
+    parameters$dpl
+  )
+
+  mockery::stub(
+    emergence_process,
+    'solver_get_states', 
+    mockery::mock(
+      c(100000, 50000, 10000),
+      c(10000, 5000, 1000)
+    )
+  )
+
+  emergence_process(0)
+
+  expect_bitset_update(
+    state$queue_update,
+    'Sm',
+    seq(7777) + 2000
+  )
+  expect_bitset_update(
+    species$queue_update,
+    'a',
+    seq(7777) + 2000
+  )
+  expect_bitset_update(
+    state$queue_update,
+    'Sm',
+    seq(778) + 2000 + 7777,
+    call = 2
+  )
+  expect_bitset_update(
+    species$queue_update,
+    'b',
+    seq(778) + 2000 + 7777,
+    call = 2
+  )
+})
diff --git a/tests/testthat/test-infection-integration.R b/tests/testthat/test-infection-integration.R
index 4485f6ef..a2ee2d09 100644
--- a/tests/testthat/test-infection-integration.R
+++ b/tests/testthat/test-infection-integration.R
@@ -297,11 +297,6 @@ test_that('calculate_treated correctly samples treated and updates the drug stat
   expect_bitset_update(variables$drug$queue_update, c(2, 1), c(1, 4))
   expect_bitset_update(variables$drug_time$queue_update, 5, c(1, 4))
 
-  expect_bitset_schedule(
-    recovery_mock,
-    c(1, 4),
-    c(3, 4)
-  )
   expect_bitset_schedule(
     detection_mock,
     c(1, 4),
@@ -312,7 +307,6 @@ test_that('calculate_treated correctly samples treated and updates the drug stat
 test_that('schedule_infections correctly schedules new infections', {
   parameters <- get_parameters()
 
-  scheduled <- individual::Bitset$new(20)$insert(c(1, 3, 7, 15))
   clinical_mock <- mockery::mock()
   asym_mock <- mockery::mock()
   detection_mock <- mockery::mock()
@@ -320,32 +314,15 @@ test_that('schedule_infections correctly schedules new infections', {
   events <- list(
     detection = list(schedule = detection_mock),
     clinical_infection = list(
-      schedule = clinical_mock,
-      get_scheduled = mockery::mock(individual::Bitset$new(20)$insert(c(1, 3)))
+      schedule = clinical_mock
     ),
     asymptomatic_infection = list(
-      schedule = asym_mock,
-      get_scheduled = mockery::mock(individual::Bitset$new(20)$insert(c(7, 15)))
+      schedule = asym_mock
     ),
     subpatent_infection = list(clear_schedule = mockery::mock()),
     recovery = list(clear_schedule = mockery::mock())
   )
 
-  mockery::stub(
-    schedule_infections,
-    'log_uniform',
-    mockery::mock(
-      c(5, 6, 13, 14),
-      c(2, 4, 16, 17, 18, 19, 20)
-    )
-  )
-
-  mockery::stub(
-    schedule_infections,
-    'bernoulli_multi_p',
-    mockery::mock(c(1, 3, 6))
-  )
-
   infections <- individual::Bitset$new(20)$insert(1:20)
   clinical_infections <- individual::Bitset$new(20)$insert(5:15)
   treated <- individual::Bitset$new(20)$insert(7:12)
@@ -360,26 +337,26 @@ test_that('schedule_infections correctly schedules new infections', {
 
   expect_bitset_schedule(
     clinical_mock,
-    c(5, 6, 13, 14),
-    c(5, 6, 13, 14)
+    c(5, 6, 13, 14, 15),
+    0
   )
 
   expect_bitset_schedule(
     detection_mock,
-    c(5, 6, 13, 14),
-    c(5, 6, 13, 14)
+    c(5, 6, 13, 14, 15),
+    0
   )
 
   expect_bitset_schedule(
     asym_mock,
-    c(2, 4, 16, 17, 18, 19, 20),
-    c(2, 4, 16, 17, 18, 19, 20)
+    c(1, 2, 3, 4, 16, 17, 18, 19, 20),
+    0,
   )
 
   expect_bitset_schedule(
     detection_mock,
-    c(2, 4, 16, 17, 18, 19, 20),
-    c(2, 4, 16, 17, 18, 19, 20),
+    c(1, 2, 3, 4, 16, 17, 18, 19, 20),
+    0,
     call = 2
   )
 })
diff --git a/tests/testthat/test-lag.R b/tests/testthat/test-lag.R
new file mode 100644
index 00000000..cfd9cad9
--- /dev/null
+++ b/tests/testthat/test-lag.R
@@ -0,0 +1,60 @@
+test_that('lag gives default when empty', {
+  expect_equal(
+    LaggedValue$new(12.5, 10)$get(1),
+    10
+  )
+})
+
+test_that('lag gives default before', {
+  expect_equal(
+    LaggedValue$new(12.5, 10)$get(1),
+    10
+  )
+})
+
+test_that('lag gives default before', {
+  lagged <- LaggedValue$new(12.5, 10)
+  lagged$save(2, 2)
+  expect_equal(
+    lagged$get(1),
+    10
+  )
+})
+
+test_that('lag gives default after', {
+  lagged <- LaggedValue$new(12.5, 10)
+  lagged$save(2, 2)
+  expect_equal(
+    lagged$get(3),
+    10
+  )
+})
+
+test_that('lag interpolates inbetween', {
+  lagged <- LaggedValue$new(12.5, 10)
+  lagged$save(2, 2)
+  lagged$save(4, 4)
+  expect_equal(
+    lagged$get(3),
+    3
+  )
+})
+
+test_that('returns on the boundary', {
+  lagged <- LaggedValue$new(12.5, 10)
+  lagged$save(2, 2)
+  expect_equal(
+    lagged$get(2),
+    2
+  )
+})
+
+test_that('returns on the boundary (multiple points)', {
+  lagged <- LaggedValue$new(12.5, 10)
+  lagged$save(2, 2)
+  lagged$save(4, 4)
+  expect_equal(
+    lagged$get(4),
+    4
+  )
+})
diff --git a/tests/testthat/test-mortality-integration.R b/tests/testthat/test-mortality-integration.R
index f7068e21..60be84b9 100644
--- a/tests/testthat/test-mortality-integration.R
+++ b/tests/testthat/test-mortality-integration.R
@@ -66,8 +66,8 @@ test_that('mortality_process resets humans correctly', {
   expect_bitset_update(variables$is_severe$queue_update, 'no', c(2, 4))
   expect_bitset_update(variables$icm$queue_update, parameters$pcm, 2)
   expect_bitset_update(variables$ivm$queue_update, parameters$pvm, 2)
-  expect_bitset_update(variables$icm$queue_update, 3 * parameters$pcm, 4, call = 2)
-  expect_bitset_update(variables$ivm$queue_update, 3 * parameters$pvm, 4, call = 2)
+  expect_bitset_update(variables$icm$queue_update, parameters$pcm, 4, call = 2)
+  expect_bitset_update(variables$ivm$queue_update, parameters$pvm, 4, call = 2)
 })
 
 test_that('died from severe samples correctly', {
diff --git a/tests/testthat/test-mortality.R b/tests/testthat/test-mortality.R
deleted file mode 100644
index 56534c4b..00000000
--- a/tests/testthat/test-mortality.R
+++ /dev/null
@@ -1,6 +0,0 @@
-test_that('discretise normal works outside the boundaries', {
-  expect_equal(
-    discretise_normal(c(-100, 2, 100), 5),
-    c(1, 4, 5)
-  )
-})
diff --git a/vignettes/MDA.Rmd b/vignettes/MDA.Rmd
index 9e0c3c8b..b0317ae4 100644
--- a/vignettes/MDA.Rmd
+++ b/vignettes/MDA.Rmd
@@ -31,11 +31,7 @@ sim_length <- 3 * year
 human_population <- 1000
 starting_EIR <- 50
 
-jamie_params <- load_parameter_set()
-eq <- human_equilibrium(EIR = starting_EIR, ft = 0, p = jamie_params, age = 0:99)
-
-simparams <- get_parameters(c(
-  translate_jamie(remove_unused_jamie(jamie_params)),
+simparams <- get_parameters(
   list(
     human_population = human_population,
     model_seasonality = TRUE, # Let's try a bi-modal model
@@ -43,11 +39,9 @@ simparams <- get_parameters(c(
     g = c(0.20636, -0.0740318, -0.0009293),
     h = c(0.173743, -0.0730962, -0.116019)
   )
-))
-
+)
+simparams <- set_equilibrium(simparams, starting_EIR)
 simparams <- set_drugs(simparams, list(SP_AQ_params))
-simparams <- parameterise_human_equilibrium(simparams, eq)
-simparams <- parameterise_mosquito_equilibrium(simparams, starting_EIR)
 
 # Plotting functions
 plot_prevalence <- function(output) {
diff --git a/vignettes/Treatment.Rmd b/vignettes/Treatment.Rmd
index a7ca9ff7..bd37375b 100644
--- a/vignettes/Treatment.Rmd
+++ b/vignettes/Treatment.Rmd
@@ -27,14 +27,7 @@ year <- 365
 sim_length <- 1 * year
 human_population <- 1000
 
-jamie_params <- load_parameter_set()
-
-simparams <- get_parameters(c(
-  translate_jamie(remove_unused_jamie(jamie_params)),
-  list(
-    human_population = human_population
-  )
-))
+simparams <- get_parameters(list(human_population = human_population))
 simparams <- set_drugs(simparams, list(AL_params, DHA_PQP_params))
 
 # NOTE: this equilibrium is calculated using a different prophylaxis model
@@ -52,11 +45,10 @@ starting_EIR <- 10
 fts <- c(0, .25, .5, .75, 1.)
 
 for (ft in fts) {
-  eq <- get_equilibrium(starting_EIR, ft)
-  simparams <- parameterise_human_equilibrium(simparams, eq)
-  simparams <- parameterise_mosquito_equilibrium(simparams, starting_EIR)
+  eq_params <- load_parameter_set()
   simparams <- set_clinical_treatment(simparams, 1, 1, ft / 2)
   simparams <- set_clinical_treatment(simparams, 2, 1, ft / 2)
+  simparams <- set_equilibrium(simparams, starting_EIR, eq_params)
   output <- run_simulation(sim_length, simparams)
   outputs[[length(outputs) + 1]] <- output
 }
diff --git a/vignettes/Vaccines.Rmd b/vignettes/Vaccines.Rmd
index fb05cfd0..2809a27b 100644
--- a/vignettes/Vaccines.Rmd
+++ b/vignettes/Vaccines.Rmd
@@ -31,12 +31,7 @@ sim_length <- 3 * year
 human_population <- 1000
 starting_EIR <- 50
 
-jamie_params <- load_parameter_set()
-eq <- human_equilibrium(EIR = starting_EIR, ft = 0, p = jamie_params, age = 0:99)
-
-simparams <- get_parameters(c(
-  translate_jamie(remove_unused_jamie(jamie_params)),
-  list(
+simparams <- get_parameters(list(
     human_population = human_population,
     model_seasonality = TRUE, # Let's try a bi-modal model
     g0 = 0.28605,
@@ -48,10 +43,9 @@ simparams <- get_parameters(c(
     severe_incidence_rendering_min_ages = 0,
     severe_incidence_rendering_max_ages = 18 * year
   )
-))
+)
 
-simparams <- parameterise_human_equilibrium(simparams, eq)
-simparams <- parameterise_mosquito_equilibrium(simparams, starting_EIR)
+simparams <- set_equilibrium(simparams, starting_EIR)
 
 # Plotting functions
 plot_prevalence <- function(output) {
diff --git a/vignettes/VectorControl.Rmd b/vignettes/VectorControl.Rmd
index 58ea390b..a2567aa9 100644
--- a/vignettes/VectorControl.Rmd
+++ b/vignettes/VectorControl.Rmd
@@ -32,11 +32,7 @@ sim_length <- 3 * year
 human_population <- 1000
 starting_EIR <- 50
 
-jamie_params <- load_parameter_set()
-eq <- human_equilibrium(EIR = starting_EIR, ft = 0, p = jamie_params, age = 0:99)
-
-simparams <- get_parameters(c(
-  translate_jamie(remove_unused_jamie(jamie_params)),
+simparams <- get_parameters(
   list(
     human_population = human_population,
     model_seasonality = TRUE, # Let's try a bi-modal model
@@ -54,12 +50,9 @@ simparams <- get_parameters(c(
     gammas = .5 * 365,
     endophily = .813
   )
-))
-
-simparams <- set_drugs(simparams, list(SP_AQ_params))
-simparams <- parameterise_human_equilibrium(simparams, eq)
-simparams <- parameterise_mosquito_equilibrium(simparams, starting_EIR)
+)
 
+simparams <- set_equilibrium(simparams, starting_EIR)
 
 # Plotting functions
 plot_prevalence <- function(output) {