diff --git a/.github/workflows/R-CMD-check.yaml b/.github/workflows/R-CMD-check.yaml
index 4b6a405e..b0d19a8d 100644
--- a/.github/workflows/R-CMD-check.yaml
+++ b/.github/workflows/R-CMD-check.yaml
@@ -35,7 +35,7 @@ jobs:
         with:
           r-version: ${{ matrix.config.r }}
 
-      - uses: r-lib/actions/setup-pandoc@v1
+      - uses: r-lib/actions/setup-pandoc@v2
 
       - name: Query dependencies
         run: |
diff --git a/DESCRIPTION b/DESCRIPTION
index 49cc359f..e60e26c9 100644
--- a/DESCRIPTION
+++ b/DESCRIPTION
@@ -1,6 +1,6 @@
 Package: malariasimulation
 Title: An individual based model for malaria
-Version: 1.3.0
+Version: 1.4.0
 Authors@R: c(
   person(
     given = "Giovanni",
@@ -10,27 +10,27 @@ Authors@R: c(
   ),
   person(
     given = "Peter",
-    family = "Windskill",
+    family = "Winskill",
     role = c('aut'),
-    email = 'p.windskill@ic.ac.uk'
+    email = 'p.winskill@ic.ac.uk'
   ),
   person(
     given = "Hillary",
     family = "Topazian",
     role = c('aut'),
-    email = 'h.topazian@ic.ac.uk'
+    email = 'h.topazian@imperial.ac.uk'
   ),
   person(
     given = "Joseph",
     family = "Challenger",
     role = c('aut'),
-    email = 'j.challenger@ic.ac.uk'
+    email = 'j.challenger@imperial.ac.uk'
   ),
   person(
     given = "Richard",
     family = "Fitzjohn",
     role = c('aut'),
-    email = 'r.fitzjohn@ic.ac.uk'
+    email = 'r.fitzjohn@imperial.ac.uk'
   ),
   person(
     given = "Imperial College of Science, Technology and Medicine",
@@ -66,7 +66,7 @@ Suggests:
   ggplot2,
   covr,
   mgcv 
-RoxygenNote: 7.1.2
+RoxygenNote: 7.2.1
 Roxygen: list(markdown = TRUE)
 LinkingTo:
   Rcpp,
diff --git a/NAMESPACE b/NAMESPACE
index e8f26fc0..c0d7331a 100644
--- a/NAMESPACE
+++ b/NAMESPACE
@@ -12,6 +12,7 @@ export(get_parameters)
 export(parameterise_mosquito_equilibrium)
 export(parameterise_total_M)
 export(peak_season_offset)
+export(run_metapop_simulation)
 export(run_simulation)
 export(run_simulation_with_repetitions)
 export(set_bednets)
diff --git a/R/biting_process.R b/R/biting_process.R
index 97d801f0..9dbb2d41 100644
--- a/R/biting_process.R
+++ b/R/biting_process.R
@@ -8,8 +8,13 @@
 #' @param variables a list of all of the model variables
 #' @param events a list of all of the model events
 #' @param parameters model pararmeters
-#' @param lagged_infectivity a LaggedValue class with historical sums of infectivity
+#' @param lagged_infectivity a list of LaggedValue objects with historical sums
+#' of infectivity, one for every metapopulation
 #' @param lagged_eir a LaggedValue class with historical EIRs
+#' @param mixing a vector of mixing coefficients for the lagged_infectivity
+#' values (default: 1)
+#' @param mixing_index an index for this population's position in the
+#' lagged_infectivity list (default: 1)
 #' @noRd
 create_biting_process <- function(
   renderer,
@@ -19,7 +24,9 @@ create_biting_process <- function(
   events,
   parameters,
   lagged_infectivity,
-  lagged_eir
+  lagged_eir,
+  mixing = 1,
+  mixing_index = 1
   ) {
   function(timestep) {
     # Calculate combined EIR
@@ -35,7 +42,9 @@ create_biting_process <- function(
       parameters,
       timestep,
       lagged_infectivity,
-      lagged_eir
+      lagged_eir,
+      mixing,
+      mixing_index
     )
 
     simulate_infection(
@@ -61,7 +70,9 @@ simulate_bites <- function(
   parameters,
   timestep,
   lagged_infectivity,
-  lagged_eir
+  lagged_eir,
+  mixing = 1,
+  mixing_index = 1
   ) {
   bitten_humans <- individual::Bitset$new(parameters$human_population)
 
@@ -117,8 +128,17 @@ simulate_bites <- function(
       n_infectious <- calculate_infectious_compartmental(solver_states)
     }
 
-    lagged_eir[[s_i]]$save(n_infectious * a, timestep)
-    species_eir <- lagged_eir[[s_i]]$get(timestep - parameters$de)
+    # store the current population's EIR for later
+    lagged_eir[[mixing_index]][[s_i]]$save(n_infectious * a, timestep)
+
+    # calculated the EIR for this timestep after mixing
+    species_eir <- sum(
+      vnapply(
+        lagged_eir,
+        function(l) l[[s_i]]$get(timestep - parameters$de)
+      ) * mixing
+    )
+
     renderer$render(paste0('EIR_', species_name), species_eir, timestep)
     EIR <- EIR + species_eir
     expected_bites <- species_eir * mean(psi)
@@ -131,9 +151,15 @@ simulate_bites <- function(
       }
     }
 
-    infectivity <- lagged_infectivity$get(timestep - parameters$delay_gam)
-    lagged_infectivity$save(sum(human_infectivity * .pi), timestep)
-    foim <- calculate_foim(a, infectivity)
+    infectivity <- vnapply(
+      lagged_infectivity,
+      function(l) l$get(timestep - parameters$delay_gam)
+    )
+    lagged_infectivity[[mixing_index]]$save(
+      sum(human_infectivity * .pi),
+      timestep
+    )
+    foim <- calculate_foim(a, infectivity, mixing)
     renderer$render(paste0('FOIM_', species_name), foim, timestep)
     mu <- death_rate(f, W, Z, s_i, parameters)
     renderer$render(paste0('mu_', species_name), mu, timestep)
@@ -276,8 +302,10 @@ unique_biting_rate <- function(age, parameters) {
 #' @title Calculate the force of infection towards mosquitoes
 #'
 #' @param a human blood meal rate
-#' @param infectivity_sum a sum of infectivity weighted by relative biting rate
+#' @param infectivity_sum a vector of sums of infectivity weighted by relative
+#' biting rate for each population
+#' @param mixing a vector of mixing coefficients for each population
 #' @noRd
-calculate_foim <- function(a, infectivity_sum) {
-  a * infectivity_sum
+calculate_foim <- function(a, infectivity_sum, mixing) {
+  a * sum(infectivity_sum * mixing)
 }
diff --git a/R/disease_progression.R b/R/disease_progression.R
index e7d8a654..05367252 100644
--- a/R/disease_progression.R
+++ b/R/disease_progression.R
@@ -7,31 +7,50 @@
 #' @param infectivity the handle for the infectivity variable
 #' @param new_infectivity the new infectivity of the progressed individuals
 #' @noRd
-create_infection_update_listener <- function(
- state,
- to_state,
- infectivity,
- new_infectivity
- ) {
-  function(timestep, to_move) {
-    state$queue_update(to_state, to_move)
-    infectivity$queue_update(new_infectivity, to_move)
-  }
+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(event, state, from_state, rate) {
+create_progression_process <- function(
+  state,
+  from_state,
+  to_state,
+  rate,
+  infectivity,
+  new_infectivity
+  ) {
   function(timestep) {
-    event$schedule(state$get_index_of(from_state)$sample(1/rate), 0)
+    to_move <- state$get_index_of(from_state)$sample(1/rate)
+    update_infection(
+      state,
+      to_state,
+      infectivity,
+      new_infectivity,
+      to_move
+    )
   }
 }
 
-create_rate_listener <- function(from_state, to_state, renderer) {
-  renderer$set_default(paste0('rate_', from_state, '_', to_state), 0)
-  function(timestep, target) {
-    renderer$render(
-      paste0('rate_', from_state, '_', to_state),
-      target$size(),
-      timestep
+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
     )
   }
 }
@@ -43,22 +62,25 @@ create_rate_listener <- function(from_state, to_state, renderer) {
 #' @param variables the available human variables
 #' @param parameters model parameters
 #' @noRd
-create_asymptomatic_update_listener <- function(variables, parameters) {
-  function(timestep, to_move) {
-    if (to_move$size() > 0) {
-      variables$state$queue_update('A', to_move)
-      new_infectivity <- asymptomatic_infectivity(
-        get_age(
-          variables$birth$get_values(to_move),
-          timestep
-        ),
-        variables$id$get_values(to_move),
-        parameters
-      )
-      variables$infectivity$queue_update(
-        new_infectivity,
-        to_move
-      )
-    }
+update_to_asymptomatic_infection <- function(
+  variables,
+  parameters,
+  timestep,
+  to_move
+  ) {
+  if (to_move$size() > 0) {
+    variables$state$queue_update('A', to_move)
+    new_infectivity <- asymptomatic_infectivity(
+      get_age(
+        variables$birth$get_values(to_move),
+        timestep
+      ),
+      variables$id$get_values(to_move),
+      parameters
+    )
+    variables$infectivity$queue_update(
+      new_infectivity,
+      to_move
+    )
   }
 }
diff --git a/R/events.R b/R/events.R
index e8e5704c..2c17e30d 100644
--- a/R/events.R
+++ b/R/events.R
@@ -1,14 +1,5 @@
 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),
-
     # MDA events
     mda_administer = individual::Event$new(),
     smc_administer = individual::Event$new(),
@@ -100,65 +91,6 @@ attach_event_listeners <- function(
   correlations,
   renderer
   ) {
-
-  # =============
-  # State updates
-  # =============
-  # When infection events fire, update the corresponding states and infectivity
-  # variables
-
-  # Infection events
-  events$clinical_infection$add_listener(
-    create_infection_update_listener(
-      variables$state,
-      'D',
-      variables$infectivity,
-      parameters$cd
-    )
-  )
-
-  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,
-      parameters
-    )
-  )
-
-  # Recovery events
-  events$subpatent_progression$add_listener(
-    create_infection_update_listener(
-      variables$state,
-      'U',
-      variables$infectivity,
-      parameters$cu
-    )
-  )
-
-  events$subpatent_progression$add_listener(
-    create_rate_listener('A', 'U', renderer)
-  )
-
-  events$recovery$add_listener(
-    create_infection_update_listener(
-      variables$state,
-      'S',
-      variables$infectivity,
-      0
-    )
-  )
-  events$recovery$add_listener(
-    create_rate_listener('U', 'S', renderer)
-  )
-
   # ===========
   # Progression
   # ===========
@@ -226,8 +158,8 @@ attach_event_listeners <- function(
       parameters$mda_drug,
       parameters$mda_timesteps,
       parameters$mda_coverages,
-      parameters$mda_min_age,
-      parameters$mda_max_age,
+      parameters$mda_min_ages,
+      parameters$mda_max_ages,
       correlations,
       'mda',
       parameters,
@@ -242,8 +174,8 @@ attach_event_listeners <- function(
       parameters$smc_drug,
       parameters$smc_timesteps,
       parameters$smc_coverages,
-      parameters$smc_min_age,
-      parameters$smc_max_age,
+      parameters$smc_min_ages,
+      parameters$smc_max_ages,
       correlations,
       'smc',
       parameters,
diff --git a/R/human_infection.R b/R/human_infection.R
index bd76ac97..a7474382 100644
--- a/R/human_infection.R
+++ b/R/human_infection.R
@@ -73,7 +73,6 @@ simulate_infection <- function(
   treated <- calculate_treated(
     variables,
     clinical_infections,
-    events$recovery,
     parameters,
     timestep,
     renderer
@@ -83,11 +82,12 @@ simulate_infection <- function(
   renderer$render('n_infections', infected_humans$size(), timestep)
 
   schedule_infections(
-    events,
+    variables,
     clinical_infections,
     treated,
     infected_humans,
-    parameters
+    parameters,
+    timestep
   )
 }
 
@@ -255,7 +255,6 @@ update_severe_disease <- function(
 #' Sample treated humans from the clinically infected
 #' @param variables a list of all of the model variables
 #' @param clinical_infections a bitset of clinically infected humans
-#' @param recovery the recovery event
 #' @param parameters model parameters
 #' @param timestep the current timestep
 #' @param renderer simulation renderer
@@ -263,7 +262,6 @@ update_severe_disease <- function(
 calculate_treated <- function(
   variables,
   clinical_infections,
-  recovery,
   parameters,
   timestep,
   renderer
@@ -319,12 +317,12 @@ calculate_treated <- function(
 #' @param parameters model parameters
 #' @noRd
 schedule_infections <- function(
-  events,
+  variables,
   clinical_infections,
   treated,
   infections,
   parameters,
-  asymptomatics
+  timestep
   ) {
   included <- treated$not(TRUE)
 
@@ -334,13 +332,22 @@ schedule_infections <- function(
   )
 
   if(to_infect$size() > 0) {
-    infection_times <- log_uniform(to_infect$size(), parameters$de)
-    events$clinical_infection$schedule(to_infect, 0)
+    update_infection(
+      variables$state,
+      'D',
+      variables$infectivity,
+      parameters$cd,
+      to_infect
+    )
   }
 
   if(to_infect_asym$size() > 0) {
-    infection_times <- log_uniform(to_infect_asym$size(), parameters$de)
-    events$asymptomatic_infection$schedule(to_infect_asym, 0)
+    update_to_asymptomatic_infection(
+      variables,
+      parameters,
+      timestep,
+      to_infect_asym
+    )
   }
 }
 
diff --git a/R/mda_parameters.R b/R/mda_parameters.R
index 166f8a1c..287ab61b 100644
--- a/R/mda_parameters.R
+++ b/R/mda_parameters.R
@@ -2,10 +2,10 @@
 #' @param parameters a list of parameters to modify
 #' @param drug the index of the drug to administer
 #' @param timesteps a vector of timesteps for each round of mda
-#' @param coverages the proportion of the target population who recieve each
+#' @param coverages a vector of the proportion of the target population who receive each
 #' round
-#' @param min_age the minimum age of the target population exclusive (in timesteps)
-#' @param max_age the maximum age of the target population exclusive (in timesteps)
+#' @param min_ages a vector of minimum ages of the target population for each round exclusive (in timesteps)
+#' @param max_ages a vector of maximum ages of the target population for each round exclusive (in timesteps)
 #' drug
 #' @export
 set_mda <- function(
@@ -13,15 +13,26 @@ set_mda <- function(
   drug,
   timesteps,
   coverages,
-  min_age,
-  max_age
+  min_ages,
+  max_ages
   ) {
+  
+  if(length(coverages) != length(timesteps)){
+    stop("coverages and timesteps do no align")
+  }
+  if(length(min_ages) != length(timesteps)){
+    stop("minimum ages and timesteps do no align")
+  }
+  if(length(max_ages) != length(timesteps)){
+    stop("maximum ages and timesteps do no align")
+  }
+  
   parameters$mda <- TRUE
   parameters$mda_drug <- drug
   parameters$mda_timesteps <- timesteps
   parameters$mda_coverages <- coverages
-  parameters$mda_min_age <- min_age
-  parameters$mda_max_age <- max_age
+  parameters$mda_min_ages <- min_ages
+  parameters$mda_max_ages <- max_ages
   parameters
 }
 
@@ -29,10 +40,10 @@ set_mda <- function(
 #' @param parameters a list of parameters to modify
 #' @param drug the index of the drug to administer
 #' @param timesteps a vector of timesteps for each round of smc
-#' @param coverages the proportion of the target population who recieve each
+#' @param coverages a vector of the proportion of the target population who receive each
 #' round
-#' @param min_age the minimum age of the target population exclusive (in timesteps)
-#' @param max_age the maximum age of the target population exclusive (in timesteps)
+#' @param min_ages a vector of minimum ages of the target population for each round exclusive (in timesteps)
+#' @param max_ages a vector of maximum ages of the target population for each round exclusive (in timesteps)
 #' drug
 #' @export
 set_smc <- function(
@@ -40,14 +51,25 @@ set_smc <- function(
   drug,
   timesteps,
   coverages,
-  min_age,
-  max_age
+  min_ages,
+  max_ages
   ) {
+  
+  if(length(coverages) != length(timesteps)){
+    stop("coverages and timesteps do no align")
+  }
+  if(length(min_ages) != length(timesteps)){
+    stop("minimum ages and timesteps do no align")
+  }
+  if(length(max_ages) != length(timesteps)){
+    stop("maximum ages and timesteps do no align")
+  }
+  
   parameters$smc <- TRUE
   parameters$smc_drug <- drug
   parameters$smc_timesteps <- timesteps
   parameters$smc_coverages <- coverages
-  parameters$smc_min_age <- min_age
-  parameters$smc_max_age <- max_age
+  parameters$smc_min_ages <- min_ages
+  parameters$smc_max_ages <- max_ages
   parameters
 }
diff --git a/R/mda_processes.R b/R/mda_processes.R
index 3833a152..348e224b 100644
--- a/R/mda_processes.R
+++ b/R/mda_processes.R
@@ -4,8 +4,8 @@
 #' @param drug the drug to administer
 #' @param timesteps timesteps for each round
 #' @param coverages the coverage for each round
-#' @param min_age minimum age for the target population
-#' @param max_age maximum age for the target population
+#' @param min_ages minimum age for the target population for each round
+#' @param max_ages maximum age for the target population for each round
 #' @param correlations correlation parameters
 #' @param int_name the name of this intervention (either 'smc' or 'mda')
 #' @param parameters the model parameters
@@ -18,8 +18,8 @@ create_mda_listeners <- function(
   drug,
   timesteps,
   coverages,
-  min_age,
-  max_age,
+  min_ages,
+  max_ages,
   correlations,
   int_name,
   parameters,
@@ -30,7 +30,7 @@ create_mda_listeners <- function(
     coverage <- coverages[[time_index]]
     age <- get_age(variables$birth$get_values(), timestep)
 
-    in_age <- which((age > min_age) & (age < max_age))
+    in_age <- which((age > min_ages[[time_index]]) & (age < max_ages[[time_index]]))
     target <- in_age[sample_intervention(in_age, int_name, coverage, correlations)]
 
     successful_treatments <- bernoulli(
diff --git a/R/model.R b/R/model.R
index 6c88e551..86262e49 100644
--- a/R/model.R
+++ b/R/model.R
@@ -71,6 +71,7 @@
 #' subpatent
 #'  * rate_U_S: rate that humans transition from subpatent to
 #' susceptible
+#'  * net_usage: the number people protected by a bed net
 #'  * mosquito_deaths: number of adult female mosquitoes who die this timestep
 #'
 #' @param timesteps the number of timesteps to run the simulation for (in days)
@@ -111,7 +112,9 @@ run_simulation <- function(
       parameters,
       vector_models,
       solvers,
-      correlations
+      correlations,
+      list(create_lagged_eir(variables, solvers, parameters)),
+      list(create_lagged_infectivity(variables, parameters))
     ),
     variables = variables,
     events = unlist(events),
@@ -120,6 +123,102 @@ run_simulation <- function(
   renderer$to_dataframe()
 }
 
+#' @title Run a metapopulation model
+#'
+#' @param timesteps the number of timesteps to run the simulation for (in days)
+#' @param parameters a list of model parameter lists for each population
+#' @param correlations a list of correlation parameters for each population
+#' (default: NULL)
+#' @param mixing matrix of mixing coefficients for infectivity towards
+#' mosquitoes. Each element must be between 0 and 1 and all rows and columns must sum to 1.
+#' @return a list of dataframe of results
+#' @export
+run_metapop_simulation <- function(
+  timesteps,
+  parameters,
+  correlations = NULL,
+  mixing
+  ) {
+  random_seed(ceiling(runif(1) * .Machine$integer.max))
+  if (nrow(mixing) != ncol(mixing)) {
+    stop('mixing matrix must be square')
+  }
+  if (nrow(mixing) != length(parameters)) {
+    stop('mixing matrix rows must match length of parameters')
+  }
+  if (!all(round(rowSums(mixing), 1) == 1)) {
+    stop('all mixing matrix rows must sum to 1')
+  }
+  if (!all(round(colSums(mixing), 1) == 1)) {
+    stop('all mixing matrix columns must sum to 1')
+  }
+  if (is.null(correlations)) {
+    correlations <- lapply(parameters, get_correlation_parameters)
+  }
+  variables <- lapply(parameters, create_variables)
+  events <- lapply(parameters, create_events)
+  renderer <- lapply(parameters, function(.) individual::Render$new(timesteps))
+  for (i in seq_along(parameters)) {
+    # NOTE: forceAndCall is necessary here to make sure i refers to the current
+    # iteration
+    forceAndCall(
+      3,
+      initialise_events,
+      events[[i]],
+      variables[[i]],
+      parameters[[i]]
+    )
+    forceAndCall(
+      5,
+      attach_event_listeners,
+      events[[i]],
+      variables[[i]],
+      parameters[[i]],
+      correlations[[i]],
+      renderer[[i]]
+    )
+  }
+  vector_models <- lapply(parameters, parameterise_mosquito_models)
+  solvers <- lapply(
+    seq_along(parameters),
+    function(i) parameterise_solvers(vector_models[[i]], parameters[[i]])
+  )
+  lagged_eir <- lapply(
+    seq_along(parameters),
+    function(i) create_lagged_eir(variables[[i]], solvers[[i]], parameters[[i]])
+  )
+  lagged_infectivity <- lapply(
+    seq_along(parameters),
+    function(i) create_lagged_infectivity(variables[[i]], parameters[[i]])
+  )
+  processes <- lapply(
+    seq_along(parameters),
+    function(i) {
+      create_processes(
+        renderer[[i]],
+        variables[[i]],
+        events[[i]],
+        parameters[[i]],
+        vector_models[[i]],
+        solvers[[i]],
+        correlations[[i]],
+        lagged_eir,
+        lagged_infectivity,
+        mixing[i,],
+        i
+      )
+    }
+  )
+  individual::simulation_loop(
+    processes = unlist(processes),
+    variables = unlist(variables),
+    events = unlist(events),
+    timesteps = timesteps
+  )
+  
+  lapply(renderer, function(r) r$to_dataframe())
+}
+
 #' @title Run the simulation with repetitions
 #'
 #' @param timesteps the number of timesteps to run the simulation for
diff --git a/R/mortality_processes.R b/R/mortality_processes.R
index 4abc2053..5884c928 100644
--- a/R/mortality_processes.R
+++ b/R/mortality_processes.R
@@ -78,12 +78,6 @@ reset_target <- function(variables, events, target, state, timestep) {
   if (target$size() > 0) {
     # clear events
     to_clear <- c(
-      'asymptomatic_progression',
-      'subpatent_progression',
-      'recovery',
-      'clinical_infection',
-      'asymptomatic_infection',
-      'detection',
       'throw_away_net',
       'rtss_mass_doses',
       'rtss_mass_booster',
diff --git a/R/parameters.R b/R/parameters.R
index 44b4a512..309b0265 100644
--- a/R/parameters.R
+++ b/R/parameters.R
@@ -364,14 +364,14 @@ get_parameters <- function(overrides = list()) {
     mda_drug = 0,
     mda_timesteps = NULL,
     mda_coverages = NULL,
-    mda_min_age = -1,
-    mda_max_age = -1,
+    mda_min_ages = -1,
+    mda_max_ages = -1,
     smc = FALSE,
     smc_drug = 0,
     smc_timesteps = NULL,
     smc_coverages = NULL,
-    smc_min_age = -1,
-    smc_max_age = -1,
+    smc_min_ages = -1,
+    smc_max_ages = -1,
     # tbv
     tbv = FALSE,
     tbv_mt = 35,
diff --git a/R/processes.R b/R/processes.R
index fc7749ee..beedff6c 100644
--- a/R/processes.R
+++ b/R/processes.R
@@ -10,6 +10,14 @@
 #' @param models a list of vector models, one for each species
 #' @param solvers a list of ode solvers, one for each species
 #' @param correlations the intervention correlations object
+#' @param lagged_eir a list of list of LaggedValue objects for EIR for each
+#' population and species in the simulation
+#' @param lagged_infectivity a list of LaggedValue objects for FOIM for each population
+#' in the simulation
+#' @param mixing a vector of mixing coefficients for the lagged transmission
+#' values (default: 1)
+#' @param mixing_index an index for this population's position in the
+#' lagged transmission lists (default: 1)
 #' @noRd
 create_processes <- function(
   renderer,
@@ -18,7 +26,11 @@ create_processes <- function(
   parameters,
   models,
   solvers,
-  correlations
+  correlations,
+  lagged_eir,
+  lagged_infectivity,
+  mixing = 1,
+  mixing_index = 1
   ) {
   # ========
   # Immunity
@@ -54,31 +66,6 @@ create_processes <- function(
   # schedule infections for humans and set last_boosted_*
   # move mosquitoes into incubating state
   # kill mosquitoes caught in vector control
-  lagged_eir <- lapply(
-    seq_along(parameters$species),
-    function(species) {
-      LaggedValue$new(
-        max_lag = parameters$de + 1,
-        default = calculate_eir(
-          species,
-          solvers,
-          variables,
-          parameters,
-          0
-        )
-      )
-    }
-  )
-
-  age <- get_age(variables$birth$get_values(), 0)
-  psi <- unique_biting_rate(age, parameters)
-  .pi <- human_pi(psi, variables$zeta$get_values())
-  init_infectivity <- sum(.pi * variables$infectivity$get_values())
-  lagged_infectivity <- LaggedValue$new(
-    max_lag = parameters$delay_gam + 2,
-    default = init_infectivity
-  )
-
   processes <- c(
     processes,
     create_biting_process(
@@ -89,32 +76,40 @@ create_processes <- function(
       events,
       parameters,
       lagged_infectivity,
-      lagged_eir
+      lagged_eir,
+      mixing,
+      mixing_index
     ),
     create_mortality_process(variables, events, renderer, parameters),
-    create_progression_process(
-      events$asymptomatic_progression,
+    create_asymptomatic_progression_process(
       variables$state,
-      'D',
-      parameters$dd
+      parameters$dd,
+      variables,
+      parameters
     ),
     create_progression_process(
-      events$subpatent_progression,
       variables$state,
       'A',
-      parameters$da
+      'U',
+      parameters$da,
+      variables$infectivity,
+      parameters$cu
     ),
     create_progression_process(
-      events$recovery,
       variables$state,
       'U',
-      parameters$du
+      'S',
+      parameters$du,
+      variables$infectivity,
+      0
     ),
     create_progression_process(
-      events$recovery,
       variables$state,
       'Tr',
-      parameters$dt
+      'S',
+      parameters$dt,
+      variables$infectivity,
+      0
     )
   )
 
@@ -231,3 +226,43 @@ create_exponential_decay_process <- function(variable, rate) {
   decay_rate <- exp(-1/rate)
   function(timestep) variable$queue_update(variable$get_values() * decay_rate)
 }
+
+#' @title Create and initialise lagged_infectivity object
+#'
+#' @param variables model variables for initialisation
+#' @param parameters model parameters 
+#' @noRd
+create_lagged_infectivity <- function(variables, parameters) {
+  age <- get_age(variables$birth$get_values(), 0)
+  psi <- unique_biting_rate(age, parameters)
+  .pi <- human_pi(psi, variables$zeta$get_values())
+  init_infectivity <- sum(.pi * variables$infectivity$get_values())
+  LaggedValue$new(
+    max_lag = parameters$delay_gam + 2,
+    default = init_infectivity
+  )
+}
+
+#' @title Create and initialise a list of lagged_eir objects per species
+#'
+#' @param variables model variables for initialisation
+#' @param solvers model differential equation solvers
+#' @param parameters model parameters 
+#' @noRd
+create_lagged_eir <- function(variables, solvers, parameters) {
+  lapply(
+    seq_along(parameters$species),
+    function(species) {
+      LaggedValue$new(
+        max_lag = parameters$de + 1,
+        default = calculate_eir(
+          species,
+          solvers,
+          variables,
+          parameters,
+          0
+        )
+      )
+    }
+  )
+}
diff --git a/R/variables.R b/R/variables.R
index 09d4f56e..95acb018 100644
--- a/R/variables.R
+++ b/R/variables.R
@@ -307,6 +307,11 @@ create_variables <- function(parameters) {
   variables
 }
 
+
+create_export_variable <- function(metapop_params) {
+  individual::DoubleVariable$new(rep(0, length(metapop_params$x)))
+}
+
 # =========
 # Utilities
 # =========
diff --git a/man/CorrelationParameters.Rd b/man/CorrelationParameters.Rd
index b1d22578..c2e6ada7 100644
--- a/man/CorrelationParameters.Rd
+++ b/man/CorrelationParameters.Rd
@@ -14,17 +14,17 @@ Describes an event in the simulation
 \section{Methods}{
 \subsection{Public methods}{
 \itemize{
-\item \href{#method-new}{\code{CorrelationParameters$new()}}
-\item \href{#method-inter_round_rho}{\code{CorrelationParameters$inter_round_rho()}}
-\item \href{#method-inter_intervention_rho}{\code{CorrelationParameters$inter_intervention_rho()}}
-\item \href{#method-sigma}{\code{CorrelationParameters$sigma()}}
-\item \href{#method-mvnorm}{\code{CorrelationParameters$mvnorm()}}
-\item \href{#method-clone}{\code{CorrelationParameters$clone()}}
+\item \href{#method-CorrelationParameters-new}{\code{CorrelationParameters$new()}}
+\item \href{#method-CorrelationParameters-inter_round_rho}{\code{CorrelationParameters$inter_round_rho()}}
+\item \href{#method-CorrelationParameters-inter_intervention_rho}{\code{CorrelationParameters$inter_intervention_rho()}}
+\item \href{#method-CorrelationParameters-sigma}{\code{CorrelationParameters$sigma()}}
+\item \href{#method-CorrelationParameters-mvnorm}{\code{CorrelationParameters$mvnorm()}}
+\item \href{#method-CorrelationParameters-clone}{\code{CorrelationParameters$clone()}}
 }
 }
 \if{html}{\out{<hr>}}
-\if{html}{\out{<a id="method-new"></a>}}
-\if{latex}{\out{\hypertarget{method-new}{}}}
+\if{html}{\out{<a id="method-CorrelationParameters-new"></a>}}
+\if{latex}{\out{\hypertarget{method-CorrelationParameters-new}{}}}
 \subsection{Method \code{new()}}{
 initialise correlation parameters
 \subsection{Usage}{
@@ -40,8 +40,8 @@ initialise correlation parameters
 }
 }
 \if{html}{\out{<hr>}}
-\if{html}{\out{<a id="method-inter_round_rho"></a>}}
-\if{latex}{\out{\hypertarget{method-inter_round_rho}{}}}
+\if{html}{\out{<a id="method-CorrelationParameters-inter_round_rho"></a>}}
+\if{latex}{\out{\hypertarget{method-CorrelationParameters-inter_round_rho}{}}}
 \subsection{Method \code{inter_round_rho()}}{
 Add rho between rounds
 \subsection{Usage}{
@@ -60,8 +60,8 @@ the intervention}
 }
 }
 \if{html}{\out{<hr>}}
-\if{html}{\out{<a id="method-inter_intervention_rho"></a>}}
-\if{latex}{\out{\hypertarget{method-inter_intervention_rho}{}}}
+\if{html}{\out{<a id="method-CorrelationParameters-inter_intervention_rho"></a>}}
+\if{latex}{\out{\hypertarget{method-CorrelationParameters-inter_intervention_rho}{}}}
 \subsection{Method \code{inter_intervention_rho()}}{
 Add rho between interventions
 \subsection{Usage}{
@@ -83,8 +83,8 @@ the intervention}
 }
 }
 \if{html}{\out{<hr>}}
-\if{html}{\out{<a id="method-sigma"></a>}}
-\if{latex}{\out{\hypertarget{method-sigma}{}}}
+\if{html}{\out{<a id="method-CorrelationParameters-sigma"></a>}}
+\if{latex}{\out{\hypertarget{method-CorrelationParameters-sigma}{}}}
 \subsection{Method \code{sigma()}}{
 Standard deviation of each intervention between rounds
 \subsection{Usage}{
@@ -93,8 +93,8 @@ Standard deviation of each intervention between rounds
 
 }
 \if{html}{\out{<hr>}}
-\if{html}{\out{<a id="method-mvnorm"></a>}}
-\if{latex}{\out{\hypertarget{method-mvnorm}{}}}
+\if{html}{\out{<a id="method-CorrelationParameters-mvnorm"></a>}}
+\if{latex}{\out{\hypertarget{method-CorrelationParameters-mvnorm}{}}}
 \subsection{Method \code{mvnorm()}}{
 multivariate norm draws for these parameters
 \subsection{Usage}{
@@ -103,8 +103,8 @@ multivariate norm draws for these parameters
 
 }
 \if{html}{\out{<hr>}}
-\if{html}{\out{<a id="method-clone"></a>}}
-\if{latex}{\out{\hypertarget{method-clone}{}}}
+\if{html}{\out{<a id="method-CorrelationParameters-clone"></a>}}
+\if{latex}{\out{\hypertarget{method-CorrelationParameters-clone}{}}}
 \subsection{Method \code{clone()}}{
 The objects of this class are cloneable with this method.
 \subsection{Usage}{
diff --git a/man/run_metapop_simulation.Rd b/man/run_metapop_simulation.Rd
new file mode 100644
index 00000000..83b83968
--- /dev/null
+++ b/man/run_metapop_simulation.Rd
@@ -0,0 +1,25 @@
+% Generated by roxygen2: do not edit by hand
+% Please edit documentation in R/model.R
+\name{run_metapop_simulation}
+\alias{run_metapop_simulation}
+\title{Run a metapopulation model}
+\usage{
+run_metapop_simulation(timesteps, parameters, correlations = NULL, mixing)
+}
+\arguments{
+\item{timesteps}{the number of timesteps to run the simulation for (in days)}
+
+\item{parameters}{a list of model parameter lists for each population}
+
+\item{correlations}{a list of correlation parameters for each population
+(default: NULL)}
+
+\item{mixing}{matrix of mixing coefficients for infectivity towards
+mosquitoes}
+}
+\value{
+a list of dataframe of results
+}
+\description{
+Run a metapopulation model
+}
diff --git a/man/run_simulation.Rd b/man/run_simulation.Rd
index 8c371f8e..7e4f795c 100644
--- a/man/run_simulation.Rd
+++ b/man/run_simulation.Rd
@@ -88,6 +88,7 @@ asymptomatic
 subpatent
 \item rate_U_S: rate that humans transition from subpatent to
 susceptible
+\item net_usage: the proportion of the population protected by a bed net
 \item mosquito_deaths: number of adult female mosquitoes who die this timestep
 }
 }
diff --git a/man/set_equilibrium.Rd b/man/set_equilibrium.Rd
index f41905bd..f1c83c54 100644
--- a/man/set_equilibrium.Rd
+++ b/man/set_equilibrium.Rd
@@ -9,7 +9,8 @@ set_equilibrium(parameters, init_EIR, eq_params = NULL)
 \arguments{
 \item{parameters}{model parameters to update}
 
-\item{init_EIR}{the desired initial EIR}
+\item{init_EIR}{the desired initial EIR (infectious bites per day over the entire human
+population)}
 
 \item{eq_params}{parameters from the malariaEquilibrium package, if null.
 The default malariaEquilibrium parameters will be used}
diff --git a/man/set_mda.Rd b/man/set_mda.Rd
index 68db2958..0b6f81c0 100644
--- a/man/set_mda.Rd
+++ b/man/set_mda.Rd
@@ -4,7 +4,7 @@
 \alias{set_mda}
 \title{Parameterise a Mass Drug Administration}
 \usage{
-set_mda(parameters, drug, timesteps, coverages, min_age, max_age)
+set_mda(parameters, drug, timesteps, coverages, min_ages, max_ages)
 }
 \arguments{
 \item{parameters}{a list of parameters to modify}
@@ -13,12 +13,12 @@ set_mda(parameters, drug, timesteps, coverages, min_age, max_age)
 
 \item{timesteps}{a vector of timesteps for each round of mda}
 
-\item{coverages}{the proportion of the target population who recieve each
+\item{coverages}{a vector of the proportion of the target population who receive each
 round}
 
-\item{min_age}{the minimum age of the target population exclusive (in timesteps)}
+\item{min_ages}{a vector of minimum ages of the target population for each round exclusive (in timesteps)}
 
-\item{max_age}{the maximum age of the target population exclusive (in timesteps)
+\item{max_ages}{a vector of maximum ages of the target population for each round exclusive (in timesteps)
 drug}
 }
 \description{
diff --git a/man/set_smc.Rd b/man/set_smc.Rd
index adf66701..95004a51 100644
--- a/man/set_smc.Rd
+++ b/man/set_smc.Rd
@@ -4,7 +4,7 @@
 \alias{set_smc}
 \title{Parameterise a Seasonal Malaria Chemoprevention}
 \usage{
-set_smc(parameters, drug, timesteps, coverages, min_age, max_age)
+set_smc(parameters, drug, timesteps, coverages, min_ages, max_ages)
 }
 \arguments{
 \item{parameters}{a list of parameters to modify}
@@ -13,12 +13,12 @@ set_smc(parameters, drug, timesteps, coverages, min_age, max_age)
 
 \item{timesteps}{a vector of timesteps for each round of smc}
 
-\item{coverages}{the proportion of the target population who recieve each
+\item{coverages}{a vector of the proportion of the target population who receive each
 round}
 
-\item{min_age}{the minimum age of the target population exclusive (in timesteps)}
+\item{min_ages}{a vector of minimum ages of the target population for each round exclusive (in timesteps)}
 
-\item{max_age}{the maximum age of the target population exclusive (in timesteps)
+\item{max_ages}{a vector of maximum ages of the target population for each round exclusive (in timesteps)
 drug}
 }
 \description{
diff --git a/tests/testthat/test-biting-integration.R b/tests/testthat/test-biting-integration.R
index 2b82a1ba..26e2da6f 100644
--- a/tests/testthat/test-biting-integration.R
+++ b/tests/testthat/test-biting-integration.R
@@ -45,7 +45,9 @@ test_that('biting_process integrates mosquito effects and human infection', {
     parameters,
     timestep,
     lagged_foim,
-    lagged_eir
+    lagged_eir,
+    1,
+    1
   )
 
   mockery::expect_args(
@@ -98,8 +100,8 @@ test_that('simulate_bites integrates eir calculation and mosquito side effects',
   mockery::stub(simulate_bites, 'aquatic_mosquito_model_update', eqs_update)
   models <- parameterise_mosquito_models(parameters)
   solvers <- parameterise_solvers(models, parameters)
-  lagged_foim <- LaggedValue$new(12.5, .001)
-  lagged_eir <- list(LaggedValue$new(12, 10))
+  lagged_foim <- list(LaggedValue$new(12.5, .001))
+  lagged_eir <- list(list(LaggedValue$new(12, 10)))
   bitten <- simulate_bites(
     renderer,
     solvers,
@@ -142,3 +144,42 @@ test_that('simulate_bites integrates eir calculation and mosquito side effects',
     c(.5, .5, .5, .5)
   )
 })
+
+
+test_that('simulate_bites works with mixed populations', {
+  population <- 4
+  timestep <- 5
+  renderer <- individual::Render$new(5)
+  parameters <- get_parameters(
+    list(human_population = population)
+  )
+  events <- create_events(parameters)
+  variables <- create_variables(parameters)
+
+  mock_foim <- mockery::mock(1)
+  mock_a <- mockery::mock(.3)
+
+  mockery::stub(simulate_bites, 'calculate_foim', mock_foim)
+  mockery::stub(simulate_bites, '.human_blood_meal_rate', mock_a)
+  models <- parameterise_mosquito_models(parameters)
+  solvers <- parameterise_solvers(models, parameters)
+  lagged_foim <- list(LaggedValue$new(12.5, .001), LaggedValue$new(12.5, .01))
+  lagged_eir <- list(list(LaggedValue$new(12, 10)), list(LaggedValue$new(12, 10)))
+  age <- c(20, 24, 5, 39) * 365
+  bitten <- simulate_bites(
+    renderer,
+    solvers,
+    models,
+    variables,
+    events,
+    age,
+    parameters,
+    timestep,
+    lagged_foim,
+    lagged_eir,
+    c(0.2, 0.8),
+    2
+  )
+
+  mockery::expect_args(mock_foim, 1, .3, c(.001, .01), c(.2, .8))
+})
diff --git a/tests/testthat/test-biting.R b/tests/testthat/test-biting.R
index e6c19fb9..91a60be4 100644
--- a/tests/testthat/test-biting.R
+++ b/tests/testthat/test-biting.R
@@ -10,3 +10,14 @@ test_that('gonotrophic_cycle cannot be negative', {
   vparams$blood_meal_rates <- 5
   expect_error(set_species(params, list(vparams), 1))
 })
+
+test_that('calculate_foim mixes correctly', {
+  expect_equal( # check 1 case
+    calculate_foim(a=0.5, infectivity=6, mixing=1),
+    3
+  )
+  expect_equal( # check 2 cases
+    calculate_foim(a=0.5, infectivity=c(6, 4), mixing=c(.5, .5)),
+    2.5
+  )
+})
diff --git a/tests/testthat/test-infection-integration.R b/tests/testthat/test-infection-integration.R
index 819c5a69..9ff16b3a 100644
--- a/tests/testthat/test-infection-integration.R
+++ b/tests/testthat/test-infection-integration.R
@@ -16,9 +16,6 @@ test_that('simulate_infection integrates different types of infection and schedu
     state = list(get_index_of = mockery::mock(asymptomatics))
   )
 
-  total_eir <- 5
-  eir <- rep(total_eir / population, population)
-
   bitten <- individual::Bitset$new(population)$insert(c(1, 3, 5, 7))
   boost_immunity_mock <- mockery::mock()
   infected <- individual::Bitset$new(population)$insert(c(1, 3, 5))
@@ -94,7 +91,6 @@ test_that('simulate_infection integrates different types of infection and schedu
     1,
     variables,
     clinical,
-    events$recovery,
     parameters,
     timestep,
     renderer
@@ -103,11 +99,12 @@ test_that('simulate_infection integrates different types of infection and schedu
   mockery::expect_args(
     schedule_mock,
     1,
-    events,
+    variables,
     clinical,
     treated,
     infected,
-    parameters
+    parameters,
+    timestep
   )
 })
 
@@ -268,7 +265,6 @@ test_that('calculate_treated correctly samples treated and updates the drug stat
   calculate_treated(
     variables,
     clinical_infections,
-    events$recovery,
     parameters,
     timestep,
     mock_render(timestep)
@@ -300,44 +296,39 @@ test_that('calculate_treated correctly samples treated and updates the drug stat
 })
 
 test_that('schedule_infections correctly schedules new infections', {
-  parameters <- get_parameters()
-
-  clinical_mock <- mockery::mock()
-  asym_mock <- mockery::mock()
-
-  events <- list(
-    clinical_infection = list(
-      schedule = clinical_mock
-    ),
-    asymptomatic_infection = list(
-      schedule = asym_mock
-    ),
-    subpatent_infection = list(clear_schedule = mockery::mock()),
-    recovery = list(clear_schedule = mockery::mock())
-  )
+  parameters <- get_parameters(list(human_population = 20))
+  variables <- create_variables(parameters)
 
   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)
 
+  infection_mock <- mockery::mock()
+  asymp_mock <- mockery::mock()
+
+  mockery::stub(schedule_infections, 'update_infection', infection_mock)
+  mockery::stub(schedule_infections, 'update_to_asymptomatic_infection', asymp_mock)
+
   schedule_infections(
-    events,
+    variables,
     clinical_infections,
     treated,
     infections,
-    parameters
+    parameters,
+    42 
   )
 
-  expect_bitset_schedule(
-    clinical_mock,
-    c(5, 6, 13, 14, 15),
-    0
+  actual_infected <- mockery::mock_args(infection_mock)[[1]][[5]]$to_vector()
+  actual_asymp_infected <- mockery::mock_args(asymp_mock)[[1]][[4]]$to_vector()
+
+  expect_equal(
+    actual_infected,
+    c(5, 6, 13, 14, 15)
   )
 
-  expect_bitset_schedule(
-    asym_mock,
-    c(1, 2, 3, 4, 16, 17, 18, 19, 20),
-    0,
+  expect_equal(
+    actual_asymp_infected,
+    c(1, 2, 3, 4, 16, 17, 18, 19, 20)
   )
 })
 
diff --git a/tests/testthat/test-mda.R b/tests/testthat/test-mda.R
index d9d65ada..034501ec 100644
--- a/tests/testthat/test-mda.R
+++ b/tests/testthat/test-mda.R
@@ -1,3 +1,36 @@
+test_that('MDA inputs align', {
+  parameters <- get_parameters()
+  expect_error(
+    parameters <- set_mda(
+      parameters,
+      drug = 1,
+      timesteps = c(50, 150),
+      coverages= 1,
+      min_ages = c(5* 365, 5 * 365),
+      max_ages = c(10 * 365, 10 * 365)
+    ), "coverages and timesteps do no align")
+  
+  expect_error(
+    parameters <- set_mda(
+      parameters,
+      drug = 1,
+      timesteps = c(50, 150),
+      coverages = c(1, 1),
+      min_ages = 5 * 365,
+      max_ages = c(10 * 365, 10 * 365)
+    ), "minimum ages and timesteps do no align")
+  
+  expect_error(
+    parameters <- set_mda(
+      parameters,
+      drug = 1,
+      timesteps = c(50, 150),
+      coverages= c(1, 1),
+      min_ages = c(5* 365, 5 * 365),
+      max_ages = 10 * 365
+    ), "maximum ages and timesteps do no align")
+ 
+})
 
 test_that('MDA moves the diseased and non-diseased population correctly', {
   timestep <- 50
@@ -9,11 +42,11 @@ test_that('MDA moves the diseased and non-diseased population correctly', {
     drug = 1,
     timesteps = c(50, 150),
     coverages= c(1, 1),
-    min_age = 5 * 365,
-    max_age = 10 * 365
+    min_ages = c(5 * 365, 0),
+    max_ages = c(10 * 365, 100 * 365)
   )
   events <- create_events(parameters)
-
+  
   variables <- list(
     state = mock_category(
       c('D', 'S', 'A', 'U', 'Tr'),
@@ -24,28 +57,28 @@ test_that('MDA moves the diseased and non-diseased population correctly', {
     drug_time = mock_double(c(1, 2, 3, 4)),
     drug = mock_double(c(1, 2, 1, 2))
   )
-
+  
   events$mda_administer <- mock_event(events$mda_administer)
-
+  
   listener <- create_mda_listeners(
     variables,
     events$mda_administer,
     parameters$mda_drug,
     parameters$mda_timesteps,
     parameters$mda_coverages,
-    parameters$mda_min_age,
-    parameters$mda_max_age,
+    parameters$mda_min_ages,
+    parameters$mda_max_ages,
     get_correlation_parameters(parameters),
     'mda',
     parameters,
     renderer
   )
-
+  
   mockery::stub(listener, 'bernoulli', mockery::mock(c(TRUE, TRUE)))
   mock_correlation <- mockery::mock(c(TRUE, TRUE))
   mockery::stub(listener, 'sample_intervention', mock_correlation)
   listener(timestep)
-
+  
   expect_equal(
     mockery::mock_args(mock_correlation)[[1]][[1]],
     c(3, 4)
@@ -55,35 +88,118 @@ test_that('MDA moves the diseased and non-diseased population correctly', {
     'Tr',
     3
   )
-
+  
   expect_bitset_update(
     variables$state$queue_update_mock(),
     'S',
     4,
     call = 2
   )
-
+  
   expect_bitset_update(
     variables$infectivity$queue_update_mock(),
     c(.3, .4) * SP_AQ_params[[2]],
     c(3, 4)
   )
-
+  
   expect_bitset_update(
     variables$drug$queue_update_mock(),
     1,
     c(3, 4)
   )
-
+  
   expect_bitset_update(
     variables$drug_time$queue_update_mock(),
     50,
     c(3, 4)
   )
-
+  
   mockery::expect_args(
     events$mda_administer$schedule,
     1,
     100
   )
 })
+
+test_that('MDA moves the diseased and non-diseased population correctly - second round, varying age range', {
+  timestep <- 150
+  renderer <- individual::Render$new(timestep)
+  parameters <- get_parameters(list(human_population = 4))
+  parameters <- set_drugs(parameters, list(SP_AQ_params))
+  parameters <- set_mda(
+    parameters,
+    drug = 1,
+    timesteps = c(50, 150),
+    coverages= c(1, 1),
+    min_ages = c(5 * 365, 0),
+    max_ages = c(10 * 365, 100 * 365)
+  )
+  events <- create_events(parameters)
+  
+  variables <- list(
+    state = mock_category(
+      c('D', 'S', 'A', 'U', 'Tr'),
+      c('D', 'S', 'A', 'U')
+    ),
+    birth = mock_double(-365 * c(2, 20, 5, 7)),
+    infectivity = mock_double(c(.1, .2, .3, .4)),
+    drug_time = mock_double(c(1, 2, 3, 4)),
+    drug = mock_double(c(1, 2, 1, 2))
+  )
+  
+  events$mda_administer <- mock_event(events$mda_administer)
+  
+  listener <- create_mda_listeners(
+    variables,
+    events$mda_administer,
+    parameters$mda_drug,
+    parameters$mda_timesteps,
+    parameters$mda_coverages,
+    parameters$mda_min_ages,
+    parameters$mda_max_ages,
+    get_correlation_parameters(parameters),
+    'mda',
+    parameters,
+    renderer
+  )
+  
+  mockery::stub(listener, 'bernoulli', mockery::mock(c(TRUE, TRUE, TRUE, TRUE)))
+  mock_correlation <- mockery::mock(c(TRUE, TRUE, TRUE, TRUE))
+  mockery::stub(listener, 'sample_intervention', mock_correlation)
+  listener(timestep)
+  
+  expect_equal(
+    mockery::mock_args(mock_correlation)[[1]][[1]],
+    c(1, 2, 3, 4)
+  )
+  expect_bitset_update(
+    variables$state$queue_update_mock(),
+    'Tr',
+    c(1, 3)
+  )
+  
+  expect_bitset_update(
+    variables$state$queue_update_mock(),
+    'S',
+    c(2, 4),
+    call = 2
+  )
+  
+  expect_bitset_update(
+    variables$infectivity$queue_update_mock(),
+    c(.1, .2, .3, .4) * SP_AQ_params[[2]],
+    c(1, 2, 3, 4)
+  )
+  
+  expect_bitset_update(
+    variables$drug$queue_update_mock(),
+    1,
+    c(1, 2, 3, 4)
+  )
+  
+  expect_bitset_update(
+    variables$drug_time$queue_update_mock(),
+    150,
+    c(1, 2, 3, 4)
+  )
+})
diff --git a/tests/testthat/test-mortality-integration.R b/tests/testthat/test-mortality-integration.R
index b4ce6e04..138f3d19 100644
--- a/tests/testthat/test-mortality-integration.R
+++ b/tests/testthat/test-mortality-integration.R
@@ -7,8 +7,8 @@ test_that('mortality_process resets humans correctly', {
     drug = 1,
     timesteps = c(50, 150),
     coverages= c(1, 1),
-    min_age = 5 * 365,
-    max_age = 10 * 365
+    min_ages = c(5 * 365, 5 * 365),
+    max_ages = c(10 * 365, 10 * 365)
   )
   events <- create_events(parameters)
   variables <- create_variables(parameters)
diff --git a/tests/testthat/test-progression.R b/tests/testthat/test-progression.R
index d4ee3459..a8b7415c 100644
--- a/tests/testthat/test-progression.R
+++ b/tests/testthat/test-progression.R
@@ -2,16 +2,17 @@
 test_that('Asymptomatic infection listener updates infectivity correctly', {
   parameters <- get_parameters(list(human_population = 4))
   variables <- create_variables(parameters)
-  listener <- create_asymptomatic_update_listener(variables, parameters)
-
-  mockery::stub(listener, 'asymptomatic_infectivity', mockery::mock(c(.2, .5)))
 
   update_mock <- mockery::mock()
   variables$infectivity <- list(queue_update = update_mock)
 
-  to_move <- individual::Bitset$new(4)
-  to_move$insert(c(2, 4))
-  listener(1, to_move)
+  to_move <- individual::Bitset$new(4)$insert(c(2, 4))
+  mockery::stub(
+    update_to_asymptomatic_infection,
+    'asymptomatic_infectivity',
+    mockery::mock(c(.2, .5))
+  )
+  update_to_asymptomatic_infection(variables, parameters, 1, to_move)
 
   expect_bitset_update(
     update_mock,
diff --git a/tests/testthat/test-simulation-e2e.R b/tests/testthat/test-simulation-e2e.R
index f5f1b8ab..0e7f7a6c 100644
--- a/tests/testthat/test-simulation-e2e.R
+++ b/tests/testthat/test-simulation-e2e.R
@@ -2,3 +2,28 @@ test_that('Simulation runs for a few timesteps', {
   sim <- run_simulation(100)
   expect_equal(nrow(sim), 100)
 })
+
+test_that('run_metapop_simulation fails with incorrect mixing matrix', {
+  population <- 4
+  timestep <- 5
+  parameters <- get_parameters(list(human_population = population))
+  paramlist <- list(parameters, parameters)
+  # incorrect params
+  mixing <- matrix(c(1, 1), nrow = 1, ncol = 2)
+  expect_error(
+    run_metapop_simulation(timesteps, parameters, NULL, mixing)
+  )
+})
+
+test_that('run_metapop_simulation integrates two models correctly', {
+  population <- 4
+  timesteps <- 5
+  parameters <- get_parameters(list(human_population = population))
+  parametersets <- list(parameters, parameters)
+  mixing <- diag(1, nrow = 2, ncol = 2)
+
+  outputs <- run_metapop_simulation(timesteps, parametersets, NULL, mixing)
+  expect_equal(length(outputs), 2)
+  expect_equal(nrow(outputs[[1]]), 5)
+  expect_equal(nrow(outputs[[2]]), 5)
+})
diff --git a/vignettes/MDA.Rmd b/vignettes/MDA.Rmd
index 1b1bba09..7ac3f20c 100644
--- a/vignettes/MDA.Rmd
+++ b/vignettes/MDA.Rmd
@@ -109,8 +109,8 @@ mdaparams <- set_mda(
   drug = 1,
   timesteps = mda_events$timestep,
   coverages = rep(.8, 2),
-  min_age = 0,
-  max_age = 200 * 365
+  min_ages = rep(0, length(mda_events$timestep)),
+  max_ages = rep(200 * 365, length(mda_events$timestep))
 )
 
 output <- run_simulation(sim_length, mdaparams)
@@ -126,7 +126,7 @@ add_mda_lines(plot_prevalence(output), mda_events)
 
 ## SMC
 
-This is a dose of SP-AQ to 90% of 2 - 11 year olds once a year a month before the peak season for mosquitos.
+This is a dose of SP-AQ to 90% of 2 - 11 year olds once a year a month before the peak season for mosquitoes.
 
 ```{r}
 smcparams <- simparams
@@ -142,8 +142,8 @@ smcparams <- set_smc(
   drug = 1,
   timesteps = smc_events$timestep,
   coverages = rep(.9, 2),
-  min_age = 2 * 365 - 1,
-  max_age = 11 * 365
+  min_ages = rep(2 * 365 - 1, length(smc_events$timestep)),
+  max_ages = rep(11 * 365, length(smc_events$timestep))
 )
 
 output <- run_simulation(sim_length, smcparams)
diff --git a/vignettes/Metapopulation.Rmd b/vignettes/Metapopulation.Rmd
new file mode 100644
index 00000000..c3674a90
--- /dev/null
+++ b/vignettes/Metapopulation.Rmd
@@ -0,0 +1,124 @@
+---
+title: "Metapopulation"
+output: rmarkdown::html_vignette
+vignette: >
+  %\VignetteIndexEntry{Metapopulation}
+  %\VignetteEngine{knitr::rmarkdown}
+  %\VignetteEncoding{UTF-8}
+---
+
+```{r, include = FALSE}
+knitr::opts_chunk$set(
+  collapse = TRUE,
+  comment = "#>"
+)
+```
+
+```{r setup}
+suppressPackageStartupMessages(library(ggplot2))
+library(malariasimulation)
+library(malariaEquilibrium)
+```
+
+# Parameterisation
+
+The metapopulation model runs the individual-based model for multiple parameterized units (countries, regions, admins, etc.) simultaneously. The inputted mixing matrix allows for the transmission of one unit to affect the transmission of another unit. 'Mixing' in the model occurs through the variables `foim_m` and `EIR`. 
+
+Here we will set up a case study of three distinct units and compare output with various transmission mixing patterns. 
+
+```{r}
+# set variables
+year <- 365
+human_population <- 1000
+sim_length <- 5 * year
+EIR_vector <- c(1, 2, 5)
+
+# get parameters
+ms_parameterize <- function(x){ # index of EIR
+  
+  params <- get_parameters(list(human_population = human_population,
+                                model_seasonality = FALSE,
+                                individual_mosquitoes = FALSE))
+  
+  # setting treatment
+  params <- set_drugs(params, list(AL_params))
+  params <- set_clinical_treatment(params, drug = 1, timesteps = 1, coverages = 0.40)
+
+  params <- set_equilibrium(params, init_EIR = EIR_vector[x])
+  
+  return(params)
+  
+}
+
+# creating a list of three parameter lists
+paramslist <- lapply(seq(1, length(EIR_vector), 1), ms_parameterize)
+
+```
+
+# Modelling 
+
+Our parameters for the three distinct units are stored in the object `paramslist`. Next we will run the metapopulation model with these parameters. We will plug in three mixing matrices - A) isolated, B) semi-mixed, C) perfectly mixed. 
+
+```{r, warning = F}
+# isolated
+mix_1 <- diag(length(EIR_vector))
+
+# semi-mixed
+mix_2 <- matrix(c(0.8, 0.1, 0.1, 
+                  0.1, 0.8, 0.1,
+                  0.1, 0.1, 0.8),
+               nrow = 3, ncol = 3)
+
+# perfectly-mixed
+mix_3 <- matrix(rep(.33, 9), nrow = 3, ncol = 3)
+
+# run model
+set.seed(123)
+
+metapop_loop <- function(mixing, mixnam){ # mixing matrix 
+
+  output <- run_metapop_simulation(timesteps = sim_length,
+                                   parameters = paramslist,
+                                   correlations = NULL,
+                                   mixing = mixing)
+  
+  # convert to dataframe and label EIR and mixing matrix type
+  output <-  do.call('rbind', output)
+  output$EIR <- c(sort(rep(EIR_vector, sim_length)))
+  
+  return(output)
+  
+}
+
+output1 <- metapop_loop(mix_1)
+output1$mix <- 'isolated'
+  
+output2 <- metapop_loop(mix_2)
+output2$mix <- 'semi-mixed'
+  
+output3 <- metapop_loop(mix_3)
+output3$mix <- 'perfectly-mixed'
+  
+output <- rbind(output1, output2, output3)
+
+# get mean PfPR 2-10 by year
+output$prev2to10 = output$p_detect_730_3650 / output$n_730_3650
+output$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)
+
+```
+
+Now let's visualize the results of mixing on PfPR2-10:
+
+```{r}
+# plot
+ggplot(data = output) + 
+  geom_line(aes(x = year, y = prev2to10, color = factor(EIR))) +
+  facet_wrap(~ mix) +
+  scale_y_continuous(limits = c(0, 0.35)) + 
+  labs(x = 'time (years)',
+       y = 'PfPR 2-10 (month)',
+       color = 'EIR') +
+  theme_classic()
+```