diff --git a/DESCRIPTION b/DESCRIPTION
index 5b1ded5a..c3f62550 100644
--- a/DESCRIPTION
+++ b/DESCRIPTION
@@ -1,6 +1,6 @@
 Package: malariasimulation
 Title: An individual based model for malaria
-Version: 1.0.1
+Version: 1.0.2
 Authors@R: c(
   person(
     given = "Giovanni",
diff --git a/R/biting_process.R b/R/biting_process.R
index 20b18fd7..b67e5feb 100644
--- a/R/biting_process.R
+++ b/R/biting_process.R
@@ -8,6 +8,8 @@
 #' @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_eir a LaggedValue class with historical EIRs
 #' @noRd
 create_biting_process <- function(
   renderer,
@@ -16,7 +18,7 @@ create_biting_process <- function(
   variables,
   events,
   parameters,
-  lagged_foim,
+  lagged_infectivity,
   lagged_eir
   ) {
   function(timestep) {
@@ -32,7 +34,7 @@ create_biting_process <- function(
       age,
       parameters,
       timestep,
-      lagged_foim,
+      lagged_infectivity,
       lagged_eir
     )
 
@@ -58,7 +60,7 @@ simulate_bites <- function(
   age,
   parameters,
   timestep,
-  lagged_foim,
+  lagged_infectivity,
   lagged_eir
   ) {
   bitten_humans <- individual::Bitset$new(parameters$human_population)
@@ -95,30 +97,8 @@ simulate_bites <- function(
     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(s_i, Z, parameters)
-
-    lambda <- .effective_biting_rates(
-      s_i,
-      psi,
-      .pi,
-      zeta,
-      p_bitten,
-      W,
-      Z,
-      f,
-      parameters
-    )
-
-    renderer$render(
-      paste0('lambda_', parameters$species[[s_i]]),
-      sum(lambda),
-      timestep
-    )
-
-    renderer$render(
-      paste0('normal_lambda_', parameters$species[[s_i]]),
-      sum(lambda) / mean(psi),
-      timestep
-    )
+    a <- .human_blood_meal_rate(f, s_i, W, parameters)
+    lambda <- effective_biting_rates(a, .pi, p_bitten)
 
     if (parameters$individual_mosquitoes) {
       species_index <- variables$species$get_index_of(
@@ -136,18 +116,19 @@ simulate_bites <- function(
       n_infectious <- calculate_infectious_compartmental(solver_states)
     }
 
-    lagged_eir[[s_i]]$save(n_infectious * sum(lambda), timestep)
+    lagged_eir[[s_i]]$save(n_infectious * a, timestep)
     species_eir <- lagged_eir[[s_i]]$get(timestep - parameters$de)
-    EIR <- EIR + species_eir / sum(psi)
-    n_bites <- rpois(1, species_eir)
+    EIR <- EIR + species_eir
+    n_bites <- rpois(1, species_eir * mean(psi))
     if (n_bites > 0) {
       bitten_humans$insert(
         fast_weighted_sample(n_bites, lambda)
       )
     }
 
-    foim <- calculate_foim(human_infectivity, lambda, psi)
-    lagged_foim$save(foim, timestep)
+    infectivity <- lagged_infectivity$get(timestep - parameters$delay_gam)
+    lagged_infectivity$save(sum(human_infectivity * .pi), timestep)
+    foim <- calculate_foim(a, infectivity)
     renderer$render(paste0('FOIM_', s_i), foim, timestep)
     mu <- death_rate(f, W, Z, s_i, parameters)
     renderer$render(paste0('mu_', s_i), mu, timestep)
@@ -161,7 +142,7 @@ simulate_bites <- function(
 
       biting_effects_individual(
         variables,
-        lagged_foim$get(timestep - parameters$delay_gam),
+        foim,
         events,
         s_i,
         susceptible_species_index,
@@ -174,7 +155,7 @@ simulate_bites <- function(
       adult_mosquito_model_update(
         models[[s_i]],
         mu,
-        lagged_foim$get(timestep - parameters$delay_gam),
+        foim,
         solver_states[[ADULT_ODE_INDICES['Sm']]],
         f
       )
@@ -192,37 +173,15 @@ simulate_bites <- function(
 # =================
 
 calculate_eir <- function(species, solvers, variables, parameters, timestep) {
-  lambda <- effective_biting_rates(species, variables, parameters, timestep)
+  a <- human_blood_meal_rate(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)
+  infectious * a * mean(psi)
 }
 
-.effective_biting_rates <- function(
-  species,
-  psi,
-  .pi,
-  zeta,
-  p_bitten,
-  W,
-  Z,
-  f,
-  parameters
-  ) {
-  a <- human_blood_meal_rate(f, species, W, parameters)
-  a * intervention_coefficient(p_bitten) * zeta * psi
+effective_biting_rates <- function(a, .pi, p_bitten) {
+  a * .pi * p_bitten$prob_bitten / sum(.pi * p_bitten$prob_bitten_survives)
 }
 
 calculate_infectious <- function(species, solvers, variables, parameters) {
@@ -274,7 +233,20 @@ blood_meal_rate <- function(v, z, parameters) {
   1 / (interrupted_foraging_time + gonotrophic_cycle)
 }
 
-human_blood_meal_rate <- function(f, v, W, parameters) {
+human_blood_meal_rate <- 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)
+  .human_blood_meal_rate(f, species, W, parameters)
+}
+
+.human_blood_meal_rate <- function(f, v, W, parameters) {
   Q <- 1 - (1 - parameters$Q0[[v]]) / W
   Q * f
 }
@@ -294,10 +266,9 @@ unique_biting_rate <- function(age, parameters) {
 
 #' @title Calculate the force of infection towards mosquitoes
 #'
-#' @param human_infectivity a vector of infectivities for each human
-#' @param lambda a vector of biting rates for each human
-#' @param psi a vector of age-based relative biting rates for each human
+#' @param a human blood meal rate
+#' @param infectivity_sum a sum of infectivity weighted by relative biting rate
 #' @noRd
-calculate_foim <- function(human_infectivity, lambda, psi) {
-  sum(human_infectivity * lambda) / mean(psi)
+calculate_foim <- function(a, infectivity_sum) {
+  a * infectivity_sum
 }
diff --git a/R/parameters.R b/R/parameters.R
index 1fdb8a8e..7664f01c 100644
--- a/R/parameters.R
+++ b/R/parameters.R
@@ -146,17 +146,10 @@
 #' please set vector control strategies using `set_betnets` and `set_spraying`
 #'
 #' * bednets - boolean for if bednets are enabled; default = FALSE 
-#' * rn - probability mosquito is repelled by the bednet; default = 0.56
-#' * rnm - minimum probability mosquito is repelled by the bednet ; default = 0.24
-#' * dn0 - probability killed by the bednet; default = 0.533
+#' * phi_bednets - proportion of bites taken in bed; default = 0.85
+#' * k0 - proportion of females bloodfed with no net; default = 0.699
 #' * spraying - boolean for if indoor spraying is enabled; default = FALSE 
-#' * rs - probability repelled by indoor spraying; default = 0.2
-#' * phi_indoors - proportion of bites taken indoors; default = 0.97
-#' * phi_bednets - proportion of bites taken in bed; default = 0.89
-#' * endophily - proportion of mosquitoes resting indoors after feeding with no
-#' intervention; default = 0.813
-#' * gammas - the half-life of spraying efficacy (timesteps); default = 91.25
-#' * gamman - the half-life of bednet efficacy (timesteps); default = 963.6
+#' * phi_indoors - proportion of bites taken indoors; default = 0.90
 #'
 #' treatment parameters:
 #' please set treatment parameters with the convenience functions in
@@ -337,20 +330,14 @@ get_parameters <- function(overrides = list()) {
     species_proportions = 1,
     blood_meal_rates    = 1/3,
     Q0                  = .92,
-    endophily           = .813,
     foraging_time       = .69,
     # bed nets
     bednets = FALSE,
-    rn = .56,
-    rnm = .24,
-    dn0 = .533,
-    phi_bednets = .89,
-    gamman = 2.64 * 365,
+    phi_bednets = .85,
+    k0 = .699,
     # indoor spraying
     spraying = FALSE,
-    rs = .2,
-    phi_indoors = .97,
-    gammas = .25 * 365,
+    phi_indoors = .90,
     # treatment
     drug_efficacy          = numeric(0),
     drug_rel_c             = numeric(0),
diff --git a/R/processes.R b/R/processes.R
index cf2871b1..6236236f 100644
--- a/R/processes.R
+++ b/R/processes.R
@@ -54,12 +54,12 @@ 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(
+  lagged_eir <- lapply(
     seq_along(parameters$species),
     function(species) {
       LaggedValue$new(
-        parameters$de + 1,
-        calculate_eir(
+        max_lag = parameters$de + 1,
+        default = calculate_eir(
           species,
           solvers,
           variables,
@@ -70,6 +70,15 @@ create_processes <- function(
     }
   )
 
+  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(
@@ -79,11 +88,8 @@ create_processes <- function(
       variables,
       events,
       parameters,
-      LaggedValue$new(
-        parameters$delay_gam + 2,
-        parameters$init_foim
-      ), # lagged FOIM
-      init_eir # lagged EIR
+      lagged_infectivity,
+      lagged_eir
     ),
     create_mortality_process(variables, events, renderer, parameters),
     create_progression_process(
diff --git a/R/variables.R b/R/variables.R
index a6885734..7635da83 100644
--- a/R/variables.R
+++ b/R/variables.R
@@ -48,8 +48,6 @@ create_variables <- function(parameters) {
   initial_age <- round(rexp(size, rate=1/parameters$average_age))
 
   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'
@@ -209,8 +207,8 @@ create_variables <- function(parameters) {
   tbv_vaccinated <- individual::DoubleVariable$new(rep(-1, size))
 
   # Init vector controls
-  net_time <- individual::DoubleVariable$new(rep(-1, size))
-  spray_time <- individual::DoubleVariable$new(rep(-1, size))
+  net_time <- individual::IntegerVariable$new(rep(-1, size))
+  spray_time <- individual::IntegerVariable$new(rep(-1, size))
 
   variables <- list(
     state = state,
diff --git a/R/vector_control.R b/R/vector_control.R
index 2629e472..dc39dd1c 100644
--- a/R/vector_control.R
+++ b/R/vector_control.R
@@ -1,13 +1,17 @@
-
-#' @title return the probability of being bitten given vector controls
+#' @title The probability of being bitten given vector controls
 #' @param timestep current timestep
 #' @param variables a list of available variables
 #' @param species the species to calculate for
 #' @param parameters model parameters
 #' @noRd
-prob_bitten <- function(timestep, variables, species, parameters) {
+prob_bitten <- function(
+  timestep,
+  variables,
+  species,
+  parameters
+  ) {
+  n <- parameters$human_population
   if (!(parameters$bednets || parameters$spraying)) {
-    n <- parameters$human_population
     return(
       list(
         prob_bitten_survives = rep(1, n),
@@ -21,8 +25,9 @@ prob_bitten <- function(timestep, variables, species, parameters) {
     phi_bednets <- parameters$phi_bednets[[species]]
     net_time <- variables$net_time$get_values()
     since_net <- timestep - net_time
-    rn <- prob_repelled_bednets(since_net, species, parameters)
-    sn <- prob_survives_bednets(rn, since_net, species, parameters)
+    matches <- match(net_time, parameters$bednet_timesteps)
+    rn <- prob_repelled_bednets(matches, since_net, species, parameters)
+    sn <- prob_survives_bednets(rn, matches, since_net, species, parameters)
     unused <- net_time == -1
     sn[unused] <- 1
     rn[unused] <- 0
@@ -34,22 +39,39 @@ prob_bitten <- function(timestep, variables, species, parameters) {
 
   if (parameters$spraying) {
     phi_indoors <- parameters$phi_indoors[[species]]
-    spray_time <- variables$spray_time$get_values()
+    unprotected <- variables$spray_time$get_index_of(set=-1)
+    protected <- unprotected$not()
+    spray_time <- variables$spray_time$get_values(protected)
+    matches <- match(spray_time, parameters$spraying_timesteps)
+    ls_theta <- parameters$spraying_ls_theta[matches, species]
+    ls_gamma <- parameters$spraying_ls_gamma[matches, species]
+    ks_theta <- parameters$spraying_ks_theta[matches, species]
+    ks_gamma <- parameters$spraying_ks_gamma[matches, species]
+    ms_theta <- parameters$spraying_ms_theta[matches, species]
+    ms_gamma <- parameters$spraying_ms_gamma[matches, species]
     since_spray <- timestep - spray_time
-    unused <- spray_time == -1
-    rs <- prob_spraying_repels(
-      since_spray,
-      parameters$rs[[species]],
-      parameters$gammas
+    ls <- spraying_decay(since_spray, ls_theta, ls_gamma)
+    ks <- parameters$k0 * spraying_decay(since_spray, ks_theta, ks_gamma)
+    ms <- spraying_decay(since_spray, ms_theta, ms_gamma)
+    js <- 1 - ls - ks
+    ms_comp <- (1 - ms)
+    ls_prime <- ls * ms_comp
+    ks_prime <- ks * ms_comp
+    js_prime <- js * ms_comp + ms
+    protected_index <- protected$to_vector()
+    rs <- rep(0, n)
+    rs[protected_index] <- prob_spraying_repels(
+      ls_prime,
+      ks_prime,
+      js_prime,
+      parameters$k0
     )
-    rs[unused] <- 0
     rs_comp <- 1 - rs
-    ss <- prob_survives_spraying(
-      since_spray,
-      parameters$endophily[[species]],
-      parameters$gammas
+    ss <- rep(1, n)
+    ss[protected_index] <- prob_survives_spraying(
+      ks_prime,
+      parameters$k0
     )
-    ss[unused] <- 1
   } else {
     phi_indoors <- 0
     rs <- 0
@@ -57,7 +79,6 @@ prob_bitten <- function(timestep, variables, species, parameters) {
     ss <- 1
   }
 
-  
   list(
     prob_bitten_survives = (
       1 - phi_indoors +
@@ -121,7 +142,10 @@ distribute_nets <- function(variables, throw_away_net, parameters, correlations)
         correlations
       ))
       variables$net_time$queue_update(timestep, target)
-      throw_away_net$schedule(target, log_uniform(length(target), parameters$bednet_retention))
+      throw_away_net$schedule(
+        target,
+        log_uniform(length(target), parameters$bednet_retention)
+      )
     }
   }
 }
@@ -135,26 +159,30 @@ throw_away_nets <- function(variables) {
 # =================
 # Utility functions
 # =================
-prob_spraying_repels <- function(t, rs0, gammas) {
-  rs0 * vector_control_decay(t, gammas)
+prob_spraying_repels <- function(ls_prime, ks_prime, js_prime, k0) {
+  (1 - ks_prime / k0) * (js_prime / (ls_prime + js_prime))
 }
 
-prob_survives_spraying <- function(t, endophily, gammas) {
-  ds <- endophily * vector_control_decay(t, gammas)
-  1 - ds
+prob_survives_spraying <- function(ks_prime, k0) {
+  ks_prime / k0
 }
 
-prob_repelled_bednets <- function(t, species, parameters) {
-  rnm <- parameters$rnm[[species]]
-  gamman <- parameters$gamman
-  (parameters$rn[[species]] - rnm) * vector_control_decay(t, gamman) + rnm
+prob_repelled_bednets <- function(matches, dt, species, parameters) {
+  rnm <- parameters$bednet_rnm[matches, species]
+  gamman <- parameters$bednet_gamman[matches]
+  (parameters$bednet_rn[matches, species] - rnm) * bednet_decay(dt, gamman) + rnm
 }
 
-prob_survives_bednets <- function(rn, t, species, parameters) {
-  dn <- parameters$dn0[[species]] * vector_control_decay(t, parameters$gamman)
+prob_survives_bednets <- function(rn, matches, dt, species, parameters) {
+  dn0 <- parameters$bednet_dn0[matches, species]
+  dn <- dn0 * bednet_decay(dt, parameters$bednet_gamman[matches])
   1 - rn - dn
 }
 
-vector_control_decay <- function(t, gamma) {
+bednet_decay <- function(t, gamma) {
   exp(-t / gamma)
 }
+
+spraying_decay <- function(t, theta, gamma) {
+  1 / (1 + exp(-(theta + gamma * t)))
+}
diff --git a/R/vector_control_parameters.R b/R/vector_control_parameters.R
index 3509d1f1..412b64e1 100644
--- a/R/vector_control_parameters.R
+++ b/R/vector_control_parameters.R
@@ -1,36 +1,63 @@
 #' @title Parameterise a bed net strategy
 #'
-#' @description The model will distribute bednets at `timesteps` to a random
+#' @description The model will distribute bed nets at `timesteps` to a random
 #' sample of the entire human population. The sample size will be a proportion
 #' of the human population taken from the corresponding `coverages`.
-#' The sample _can_ contain humans who already have bednets.
+#' The sample _can_ contain humans who already have bed nets.
 #'
-#' All of the sample "use" their bednets on the timestep after they are
+#' All of the sample "use" their bed nets on the timestep after they are
 #' distributed. Incomplete usage is not part of this model.
 #'
-#' If a human in the sample already has a bednet, their bednet will be replaced
+#' If a human in the sample already has a bed net, their bed net will be replaced
 #' by a new one.
 #'
-#' Bednets will be randomly removed each timestep with a rate of `1 -
+#' Bed nets will be randomly removed each timestep with a rate of `1 -
 #' exp(-1/retention)`
 #'
+#' The structure for the bed net model is documented in the
+#' S.I. of 10.1038/s41467-018-07357-w
+#'
 #' @param parameters a list of parameters to modify
-#' @param timesteps the timesteps at which to distribute bednets
-#' @param coverages the proportion of the human population who receive bednets
+#' @param timesteps the timesteps at which to distribute bed nets
+#' @param coverages the proportion of the human population who receive bed nets
 #' @param retention the average number of timesteps a net is kept for
+#' @param dn0 a matrix of death probabilities for each species over time.
+#' With dimensions length(species) x length(timesteps)
+#' @param rn a matrix of repelling probabilities for each species over time
+#' With dimensions length(species) x length(timesteps)
+#' @param rnm a matrix of minimum repelling probabilities for each species over time
+#' With dimensions length(species) x length(timesteps)
+#' @param gamman a vector of bednet half-lives for each distribution timestep
 #' @export
 set_bednets <- function(
   parameters,
   timesteps,
   coverages,
-  retention
+  retention,
+  dn0,
+  rn,
+  rnm,
+  gamman
   ) {
-  if (length(timesteps) != length(coverages)) {
-    stop('timesteps and coverages must align')
+  lengths <- vnapply(list(coverages, gamman), length)
+  if (!all(lengths == length(timesteps))) {
+    stop('timesteps and time-varying parameters must align')
+  }
+  for (x in list(dn0, rn, rnm)) {
+    if (ncol(x) != length(parameters$species)) {
+      stop('death and repelling probabilities rows need to align with species')
+    }
+    if (nrow(x) != length(timesteps)) {
+      stop('death and repelling probabilities columns need to align with timesteps')
+    }
   }
   parameters$bednets <- TRUE
   parameters$bednet_timesteps <- timesteps
   parameters$bednet_coverages <- coverages
+  parameters$bednet_dn0 <- dn0
+  parameters$bednet_rn <- rn
+  parameters$bednet_rnm <- rnm
+  parameters$bednet_gamman <- gamman
   parameters$bednet_retention <- retention
   parameters
 }
@@ -40,26 +67,69 @@ set_bednets <- function(
 #' @description The model will apply indoor spraying at `timesteps` to a random
 #' sample of the entire human population. The sample size will be a proportion
 #' of the human population taken from the corresponding `coverages`.
-#' The sample _can_ contain humans who have already benefitted from spraying.
+#' The sample _can_ contain humans who have already benefited from spraying.
 #'
 #' If a human in the sample lives in a sprayed house, the efficacy of the
 #' spraying will be returned to the maximum.
 #'
+#' The structure for the indoor residual spraying model is documented in the
+#' S.I. of 10.1038/s41467-018-07357-w
+#'
 #' @param parameters a list of parameters to modify
 #' @param timesteps the timesteps at which to spray
 #' @param coverages the proportion of the population who get indoor
 #' spraying
+#' @param ls_theta matrix of mortality parameters
+#' With dimensions length(species) x length(timesteps)
+#' @param ls_gamma matrix of mortality parameters per timestep
+#' With dimensions length(species) x length(timesteps)
+#' @param ks_theta matrix of feeding success parameters per timestep 
+#' With dimensions length(species) x length(timesteps)
+#' @param ks_gamma matrix of feeding success parameters per timestep 
+#' With dimensions length(species) x length(timesteps)
+#' @param ms_theta matrix of deterrence parameters per timestep 
+#' With dimensions length(species) x length(timesteps)
+#' @param ms_gamma matrix of deterrence parameters per timestep 
+#' With dimensions length(species) x length(timesteps)
 #' @export
 set_spraying <- function(
   parameters,
   timesteps,
-  coverages
+  coverages,
+  ls_theta,
+  ls_gamma,
+  ks_theta,
+  ks_gamma,
+  ms_theta,
+  ms_gamma
   ) {
-  if (length(timesteps) != length(coverages)) {
-    stop('timesteps and coverages must align')
+  if (length(coverages) != length(timesteps)) {
+    stop('coverages and timesteps must must align')
+  }
+  decays <- list(
+    ls_theta,
+    ls_gamma,
+    ks_theta,
+    ks_gamma,
+    ms_theta,
+    ms_gamma
+  )
+  for (x in decays) {
+    if (ncol(x) != length(parameters$species)) {
+      stop('theta and gamma rows need to align with species')
+    }
+    if (nrow(x) != length(timesteps)) {
+      stop('theta and gamma cols need to align with timesteps')
+    }
   }
   parameters$spraying <- TRUE
   parameters$spraying_timesteps <- timesteps
   parameters$spraying_coverages <- coverages
+  parameters$spraying_ls_theta <- ls_theta
+  parameters$spraying_ls_gamma <- ls_gamma
+  parameters$spraying_ks_theta <- ks_theta
+  parameters$spraying_ks_gamma <- ks_gamma
+  parameters$spraying_ms_theta <- ms_theta
+  parameters$spraying_ms_gamma <- ms_gamma
   parameters
 }
diff --git a/R/vector_parameters.R b/R/vector_parameters.R
index 760178d2..888477e1 100644
--- a/R/vector_parameters.R
+++ b/R/vector_parameters.R
@@ -3,26 +3,19 @@
 #' species: "gamb"
 #' blood_meal_rates: 0.3333333
 #' Q0: 0.92
-#' endophily: 0.813
-#' rn: 0.56
-#' rnm: 0.24
-#' dn0: 0.533
-#' phi_bednets: 0.89
-#' rs: 0.2
-#' phi_indoors: 0.97
-#' mum: 0.132"
+#' phi_bednets: 0.85
+#' phi_indoors: 0.90
+#' mum: 0.132
+#'
+#' parameters from:
+#' https://www.pnas.org/content/pnas/early/2019/07/02/1820646116.full.pdf
 #' @export
 gamb_params <- list(
   species = 'gamb',
   blood_meal_rates = 1/3,
   Q0 = .92,
-  endophily = .813,
-  rn = .56,
-  rnm = .24,
-  dn0 = .533,
-  phi_bednets = .89,
-  rs = .2,
-  phi_indoors = .97,
+  phi_bednets = .85,
+  phi_indoors = .90,
   mum = .132
 )
 
@@ -31,26 +24,19 @@ gamb_params <- list(
 #' species: "arab"
 #' blood_meal_rates: 0.3333333
 #' Q0: 0.71
-#' endophily: 0.422
-#' rn: 0.46
-#' rnm: 0.1
-#' dn0: 0.533
-#' phi_bednets: 0.9
-#' rs: 0.2
-#' phi_indoors: 0.96
+#' phi_bednets: 0.8
+#' phi_indoors: 0.86
 #' mum: 0.132
+#'
+#' parameters from:
+#' https://www.pnas.org/content/pnas/early/2019/07/02/1820646116.full.pdf
 #' @export
 arab_params <- list(
   species = 'arab',
   blood_meal_rates = 1/3,
   Q0 = .71,
-  endophily = .422,
-  rn = .46,
-  rnm = .1,
-  dn0 = .533,
-  phi_bednets = .9,
-  rs = .2,
-  phi_indoors = .96,
+  phi_bednets = .8,
+  phi_indoors = .86,
   mum = .132
 )
 
@@ -59,26 +45,19 @@ arab_params <- list(
 #' species: "fun"
 #' blood_meal_rates: 0.3333333
 #' Q0: 0.94
-#' endophily: 0.813
-#' rn: 0.56
-#' rnm: 0.1
-#' dn0: 0.533
-#' phi_bednets: 0.9
-#' rs: 0.2
-#' phi_indoors: 0.98
+#' phi_bednets: 0.78
+#' phi_indoors: 0.87
 #' mum: 0.112
+#'
+#' parameters from:
+#' https://www.pnas.org/content/pnas/early/2019/07/02/1820646116.full.pdf
 #' @export
 fun_params <- list(
   species = 'fun',
   blood_meal_rates = 1/3,
   Q0 = .94,
-  endophily = .813,
-  rn = .56,
-  rnm = .1,
-  dn0 = .533,
-  phi_bednets = .9,
-  rs = .2,
-  phi_indoors = .98,
+  phi_bednets = .78,
+  phi_indoors = .87,
   mum = .112
 )
 
@@ -99,12 +78,7 @@ set_species <- function(parameters, species, proportions) {
     'species',
     'blood_meal_rates',
     'Q0',
-    'endophily',
-    'rn',
-    'rnm',
-    'dn0',
     'phi_bednets',
-    'rs',
     'phi_indoors',
     'mum'
   )
diff --git a/man/arab_params.Rd b/man/arab_params.Rd
index 625c992c..b5655b25 100644
--- a/man/arab_params.Rd
+++ b/man/arab_params.Rd
@@ -5,7 +5,7 @@
 \alias{arab_params}
 \title{Preset parameters for the An. arabiensis vector}
 \format{
-An object of class \code{list} of length 11.
+An object of class \code{list} of length 6.
 }
 \usage{
 arab_params
@@ -18,13 +18,11 @@ Default parameters:
 species: "arab"
 blood_meal_rates: 0.3333333
 Q0: 0.71
-endophily: 0.422
-rn: 0.46
-rnm: 0.1
-dn0: 0.533
-phi_bednets: 0.9
-rs: 0.2
-phi_indoors: 0.96
+phi_bednets: 0.8
+phi_indoors: 0.86
 mum: 0.132
+
+parameters from:
+https://www.pnas.org/content/pnas/early/2019/07/02/1820646116.full.pdf
 }
 \keyword{datasets}
diff --git a/man/fun_params.Rd b/man/fun_params.Rd
index 222abcd5..b5d0eaff 100644
--- a/man/fun_params.Rd
+++ b/man/fun_params.Rd
@@ -5,7 +5,7 @@
 \alias{fun_params}
 \title{Preset parameters for the An. funestus vector}
 \format{
-An object of class \code{list} of length 11.
+An object of class \code{list} of length 6.
 }
 \usage{
 fun_params
@@ -18,13 +18,11 @@ Default parameters:
 species: "fun"
 blood_meal_rates: 0.3333333
 Q0: 0.94
-endophily: 0.813
-rn: 0.56
-rnm: 0.1
-dn0: 0.533
-phi_bednets: 0.9
-rs: 0.2
-phi_indoors: 0.98
+phi_bednets: 0.78
+phi_indoors: 0.87
 mum: 0.112
+
+parameters from:
+https://www.pnas.org/content/pnas/early/2019/07/02/1820646116.full.pdf
 }
 \keyword{datasets}
diff --git a/man/gamb_params.Rd b/man/gamb_params.Rd
index e22451c4..434c46fa 100644
--- a/man/gamb_params.Rd
+++ b/man/gamb_params.Rd
@@ -5,7 +5,7 @@
 \alias{gamb_params}
 \title{Preset parameters for the An. gambiae s.s vector}
 \format{
-An object of class \code{list} of length 11.
+An object of class \code{list} of length 6.
 }
 \usage{
 gamb_params
@@ -18,13 +18,11 @@ Default parameters:
 species: "gamb"
 blood_meal_rates: 0.3333333
 Q0: 0.92
-endophily: 0.813
-rn: 0.56
-rnm: 0.24
-dn0: 0.533
-phi_bednets: 0.89
-rs: 0.2
-phi_indoors: 0.97
-mum: 0.132"
+phi_bednets: 0.85
+phi_indoors: 0.90
+mum: 0.132
+
+parameters from:
+https://www.pnas.org/content/pnas/early/2019/07/02/1820646116.full.pdf
 }
 \keyword{datasets}
diff --git a/man/get_parameters.Rd b/man/get_parameters.Rd
index 85eee836..584271e9 100644
--- a/man/get_parameters.Rd
+++ b/man/get_parameters.Rd
@@ -163,17 +163,10 @@ feeding cycle:
 please set vector control strategies using \code{set_betnets} and \code{set_spraying}
 \itemize{
 \item bednets - boolean for if bednets are enabled; default = FALSE
-\item rn - probability mosquito is repelled by the bednet; default = 0.56
-\item rnm - minimum probability mosquito is repelled by the bednet ; default = 0.24
-\item dn0 - probability killed by the bednet; default = 0.533
+\item phi_bednets - proportion of bites taken in bed; default = 0.85
+\item k0 - proportion of females bloodfed with no net; default = 0.699
 \item spraying - boolean for if indoor spraying is enabled; default = FALSE
-\item rs - probability repelled by indoor spraying; default = 0.2
-\item phi_indoors - proportion of bites taken indoors; default = 0.97
-\item phi_bednets - proportion of bites taken in bed; default = 0.89
-\item endophily - proportion of mosquitoes resting indoors after feeding with no
-intervention; default = 0.813
-\item gammas - the half-life of spraying efficacy (timesteps); default = 91.25
-\item gamman - the half-life of bednet efficacy (timesteps); default = 963.6
+\item phi_indoors - proportion of bites taken indoors; default = 0.90
 }
 
 treatment parameters:
diff --git a/man/set_bednets.Rd b/man/set_bednets.Rd
index db77a8ce..7dfbe1a5 100644
--- a/man/set_bednets.Rd
+++ b/man/set_bednets.Rd
@@ -4,28 +4,42 @@
 \alias{set_bednets}
 \title{Parameterise a bed net strategy}
 \usage{
-set_bednets(parameters, timesteps, coverages, retention)
+set_bednets(parameters, timesteps, coverages, retention, dn0, rn, rnm, gamman)
 }
 \arguments{
 \item{parameters}{a list of parameters to modify}
 
-\item{timesteps}{the timesteps at which to distribute bednets}
+\item{timesteps}{the timesteps at which to distribute bed nets}
 
-\item{coverages}{the proportion of the human population who receive bednets}
+\item{coverages}{the proportion of the human population who receive bed nets}
 
 \item{retention}{the average number of timesteps a net is kept for}
+
+\item{dn0}{a matrix of death probabilities for each species over time.
+With dimensions length(species) x length(timesteps)}
+
+\item{rn}{a matrix of repelling probabilities for each species over time
+With dimensions length(species) x length(timesteps)}
+
+\item{rnm}{a matrix of minimum repelling probabilities for each species over time
+With dimensions length(species) x length(timesteps)}
+
+\item{gamman}{a vector of bednet half-lives for each distribution timestep}
 }
 \description{
-The model will distribute bednets at \code{timesteps} to a random
+The model will distribute bed nets at \code{timesteps} to a random
 sample of the entire human population. The sample size will be a proportion
 of the human population taken from the corresponding \code{coverages}.
-The sample \emph{can} contain humans who already have bednets.
+The sample \emph{can} contain humans who already have bed nets.
 
-All of the sample "use" their bednets on the timestep after they are
+All of the sample "use" their bed nets on the timestep after they are
 distributed. Incomplete usage is not part of this model.
 
-If a human in the sample already has a bednet, their bednet will be replaced
+If a human in the sample already has a bed net, their bed net will be replaced
 by a new one.
 
-Bednets will be randomly removed each timestep with a rate of \code{1 - exp(-1/retention)}
+Bed nets will be randomly removed each timestep with a rate of \code{1 - exp(-1/retention)}
+
+The structure for the bed net model is documented in the
+S.I. of 10.1038/s41467-018-07357-w
 }
diff --git a/man/set_spraying.Rd b/man/set_spraying.Rd
index c750d9d2..7a6bdda0 100644
--- a/man/set_spraying.Rd
+++ b/man/set_spraying.Rd
@@ -4,7 +4,17 @@
 \alias{set_spraying}
 \title{Parameterise an indoor spraying strategy}
 \usage{
-set_spraying(parameters, timesteps, coverages)
+set_spraying(
+  parameters,
+  timesteps,
+  coverages,
+  ls_theta,
+  ls_gamma,
+  ks_theta,
+  ks_gamma,
+  ms_theta,
+  ms_gamma
+)
 }
 \arguments{
 \item{parameters}{a list of parameters to modify}
@@ -13,13 +23,34 @@ set_spraying(parameters, timesteps, coverages)
 
 \item{coverages}{the proportion of the population who get indoor
 spraying}
+
+\item{ls_theta}{matrix of mortality parameters
+With dimensions length(species) x length(timesteps)}
+
+\item{ls_gamma}{matrix of mortality parameters per timestep
+With dimensions length(species) x length(timesteps)}
+
+\item{ks_theta}{matrix of feeding success parameters per timestep
+With dimensions length(species) x length(timesteps)}
+
+\item{ks_gamma}{matrix of feeding success parameters per timestep
+With dimensions length(species) x length(timesteps)}
+
+\item{ms_theta}{matrix of deterrence parameters per timestep
+With dimensions length(species) x length(timesteps)}
+
+\item{ms_gamma}{matrix of deterrence parameters per timestep
+With dimensions length(species) x length(timesteps)}
 }
 \description{
 The model will apply indoor spraying at \code{timesteps} to a random
 sample of the entire human population. The sample size will be a proportion
 of the human population taken from the corresponding \code{coverages}.
-The sample \emph{can} contain humans who have already benefitted from spraying.
+The sample \emph{can} contain humans who have already benefited from spraying.
 
 If a human in the sample lives in a sprayed house, the efficacy of the
 spraying will be returned to the maximum.
+
+The structure for the indoor residual spraying model is documented in the
+S.I. of 10.1038/s41467-018-07357-w
 }
diff --git a/tests/testthat/test-biology.R b/tests/testthat/test-biology.R
index 8f8948c9..1d450353 100644
--- a/tests/testthat/test-biology.R
+++ b/tests/testthat/test-biology.R
@@ -1,5 +1,6 @@
 test_that('total_M and EIR functions are consistent with equilibrium EIR', {
-  EIR <- 5
+  skip_on_ci()
+  expected_EIR <- c(1, 5, 10, 100)
   population <- 100000
   parameters <- get_parameters(
     list(
@@ -7,101 +8,102 @@ test_that('total_M and EIR functions are consistent with equilibrium EIR', {
       enable_heterogeneity = FALSE
     )
   )
-  parameters <- set_equilibrium(parameters, EIR)
+  actual_EIR <- vnapply(expected_EIR, function(EIR) {
+    parameters <- set_equilibrium(parameters, EIR)
+
+    #set up arguments for EIR calculation
+    variables <- create_variables(parameters)
+    vector_models <- parameterise_mosquito_models(parameters)
+    solvers <- parameterise_solvers(vector_models, parameters)
+    estimated_eir <- calculate_eir(1, solvers, variables, parameters, 0)
+    age <- get_age(variables$birth$get_values(), 0)
+    psi <- unique_biting_rate(age, parameters)
+    omega <- mean(psi)
+    mean(estimated_eir) / omega / population * 365
+  })
 
-  #set up arguments for EIR calculation
-  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,
+    actual_EIR,
+    expected_EIR,
     tolerance = 1e-2
   )
 })
 
 test_that('total_M and EIR functions are consistent with equilibrium EIR (with het)', {
+  skip_on_ci()
   population <- 100000
 
-  EIR <- 50
+  expected_EIR <- c(1, 5, 10, 100)
 
   parameters <- get_parameters(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)
-  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,
-    eq_params,
-    ages
-  )$FOIM
+  actual_EIR <- vnapply(expected_EIR, function(EIR) {
+    parameters <- set_equilibrium(parameters, EIR)
 
-  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)
+    estimated_eir <- calculate_eir(1, solvers, variables, parameters, 0)
+    age <- get_age(variables$birth$get_values(), 0)
+    psi <- unique_biting_rate(age, parameters)
+    omega <- mean(psi)
+    mean(estimated_eir) / omega / population * 365
+  })
 
-  variables <- create_variables(parameters)
-  vector_models <- parameterise_mosquito_models(parameters)
-  solvers <- parameterise_solvers(vector_models, parameters)
-  lambda <- effective_biting_rates(1, variables, parameters, 0)
-  psi <- unique_biting_rate(-variables$birth$get_values(), parameters)
   expect_equal(
-    calculate_foim(variables$infectivity$get_values(), lambda, psi),
-    foim,
-    tolerance=1e-2
+    actual_EIR,
+    expected_EIR,
+    tolerance = 1e-2
   )
 })
 
-test_that('FOI is consistent with equilibrium', {
-  population <- 1e5
+test_that('FOIM is consistent with equilibrium', {
+  skip_on_ci()
+  population <- 100000
 
-  EIR <- 5
+  EIRs <- c(1, 5, 10, 100)
 
   eq_params <- malariaEquilibrium::load_parameter_set()
 
   ages <- 0:999 / 10
-  eq <- malariaEquilibrium::human_equilibrium(
-    EIR,
-    0,
-    eq_params,
-    ages
+  expected_foim <- vnapply(
+    EIRs,
+    function(EIR) {
+      malariaEquilibrium::human_equilibrium(
+        EIR,
+        0,
+        eq_params,
+        ages
+      )$FOIM
+    }
   )
-  foi <- weighted.mean(eq$states[,'FOI'], eq$states[,'prop'])
-
-  parameters <- get_parameters(list(human_population = population))
-  parameters <- set_equilibrium(parameters, EIR)
 
-  variables <- create_variables(parameters)
-  lambda <- effective_biting_rates(1, variables, parameters, 0)
-  ib <- blood_immunity(variables$ib$get_values(), parameters)
+  actual_foim <- vnapply(
+    EIRs,
+    function(EIR) {
+      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)
+      a <- human_blood_meal_rate(1, variables, parameters, 0)
+      age <- get_age(variables$birth$get_values(), 0)
+      psi <- unique_biting_rate(age, parameters)
+      zeta <- variables$zeta$get_values()
+      .pi <- human_pi(psi, zeta)
+      calculate_foim(a, sum(.pi * variables$infectivity$get_values()))
+    }
+  )
   expect_equal(
-    mean(lambda * ib),
-    foi,
-    tolerance=1e-2
+    expected_foim,
+    actual_foim,
+    tolerance = 1e-4
   )
 })
 
 test_that('phi is consistent with equilibrium at high EIR (no het)', {
+  skip_on_ci()
   population <- 1e5
 
   EIR <- 100
@@ -136,6 +138,7 @@ test_that('phi is consistent with equilibrium at high EIR (no het)', {
 })
 
 test_that('phi is consistent with equilibrium at high EIR', {
+  skip_on_ci()
   population <- 1e5
 
   EIR <- 100
@@ -177,8 +180,6 @@ test_that('mosquito_limit is set to 1 for 0 EIR', {
 })
 
 test_that('mosquito_limit is set to a sensible level', {
-  #EIRs <- c(5, 50, 1000)
-  #EIRs <- 1000
   EIRs <- 5
 
   seasonalities <- list(
diff --git a/tests/testthat/test-biting-integration.R b/tests/testthat/test-biting-integration.R
index bd76959e..37a6a30c 100644
--- a/tests/testthat/test-biting-integration.R
+++ b/tests/testthat/test-biting-integration.R
@@ -94,7 +94,7 @@ test_that('simulate_bites integrates eir calculation and mosquito side effects',
   mockery::stub(simulate_bites, 'biting_effects_individual', mosquito_effects_mock)
   mockery::stub(simulate_bites, 'rpois', pois_mock)
   mockery::stub(simulate_bites, 'fast_weighted_sample', sample_mock)
-  mockery::stub(simulate_bites, '.effective_biting_rates', lambda_mock)
+  mockery::stub(simulate_bites, 'effective_biting_rates', lambda_mock)
   mockery::stub(simulate_bites, 'mosquito_model_update', ode_update)
   models <- parameterise_mosquito_models(parameters)
   solvers <- parameterise_solvers(models, parameters)
@@ -120,7 +120,6 @@ test_that('simulate_bites integrates eir calculation and mosquito side effects',
   effects_args <- mockery::mock_args(mosquito_effects_mock)
 
   expect_equal(effects_args[[1]][[1]], variables)
-  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)
@@ -130,7 +129,12 @@ test_that('simulate_bites integrates eir calculation and mosquito side effects',
   expect_equal(effects_args[[1]][[9]], timestep)
 
   mockery::expect_args(ode_update, 1, models[[1]], 25, f, parameters$mum)
-  mockery::expect_args(pois_mock, 1, 1, 10)
+  mockery::expect_args(
+    pois_mock,
+    1,
+    1,
+    10 * mean(unique_biting_rate(age, parameters))
+  )
   mockery::expect_args(
     sample_mock,
     1,
diff --git a/tests/testthat/test-vector-control.R b/tests/testthat/test-vector-control.R
index 53c53e05..c7112355 100644
--- a/tests/testthat/test-vector-control.R
+++ b/tests/testthat/test-vector-control.R
@@ -1,14 +1,48 @@
-test_that('set_bednets validates parameters', {
+test_that('set_bednets validates coverages', {
   parameters <- get_parameters()
   expect_error(
-    set_bednets(parameters, c(5, 50), .8, 40),
-    '*'
+    set_bednets(
+      parameters,
+      timesteps = c(5, 50),
+      coverages = c(.5),
+      retention = 40,
+      dn0 = matrix(c(.533, .533), nrow=2, ncol=1),
+      rn = matrix(c(.56, .56), nrow=2, ncol=1),
+      rnm = matrix(c(.24, .24), nrow=2, ncol=1),
+      gamman = c(963.6, 963.6)
+    )
+  )
+})
+
+test_that('set_bednets validates matrices', {
+  parameters <- get_parameters()
+  parameters <- set_species(parameters, list(gamb_params, fun_params), c(.1, .9))
+  expect_error(
+    set_bednets(
+      parameters,
+      timesteps = c(5, 50),
+      coverages = c(.5, .9),
+      retention = 40,
+      dn0 = matrix(c(.533, .533), nrow=2, ncol=1),
+      rn = matrix(c(.56, .56), nrow=2, ncol=1),
+      rnm = matrix(c(.24, .24), nrow=2, ncol=1),
+      gamman = c(963.6, 963.6)
+    )
   )
 })
 
 test_that('set_bednets sets parameters', {
   parameters <- get_parameters()
-  parameters <- set_bednets(parameters, c(5, 50), c(.5, .9), 40)
+  parameters <- set_bednets(
+    parameters,
+    timesteps = c(5, 50),
+    coverages = c(.5, .9),
+    retention = 40,
+    dn0 = matrix(c(.533, .533), nrow=2, ncol=1),
+    rn = matrix(c(.56, .56), nrow=2, ncol=1),
+    rnm = matrix(c(.24, .24), nrow=2, ncol=1),
+    gamman = c(963.6, 963.6)
+  )
   expect_true(parameters$bednets)
   expect_equal(parameters$bednet_timesteps, c(5, 50))
   expect_equal(parameters$bednet_coverages, c(.5, .9))
@@ -18,14 +52,33 @@ test_that('set_bednets sets parameters', {
 test_that('set_spraying validates parameters', {
   parameters <- get_parameters()
   expect_error(
-    set_spraying(parameters, c(5, 50), .8),
-    '*'
+    set_spraying(
+      parameters,
+      timesteps = c(5, 50),
+      coverages = c(.5, .9),
+      ls_theta = matrix(c(2.025, 2.025), nrow=2, ncol=1),
+      ls_gamma = matrix(c(-0.009, -0.009), nrow=2, ncol=1),
+      ks_theta = matrix(c(-2.222, -2.222), nrow=2, ncol=1),
+      ks_gamma = matrix(c(0.008, 0.008), nrow=2, ncol=1),
+      ms_theta = matrix(c(-1.232, -1.232), nrow=2, ncol=1),
+      ms_gamma = matrix(c(-0.009), nrow=1, ncol=1)
+    )
   )
 })
 
 test_that('set_spraying sets parameters', {
   parameters <- get_parameters()
-  parameters <- set_spraying(parameters, c(5, 50), c(.5, .9))
+  parameters <- set_spraying(
+    parameters,
+    timesteps = c(5, 50),
+    coverages = c(.5, .9),
+    ls_theta = matrix(c(2.025, 2.025), nrow=2, ncol=1),
+    ls_gamma = matrix(c(-0.009, -0.009), nrow=2, ncol=1),
+    ks_theta = matrix(c(-2.222, -2.222), nrow=2, ncol=1),
+    ks_gamma = matrix(c(0.008, 0.008), nrow=2, ncol=1),
+    ms_theta = matrix(c(-1.232, -1.232), nrow=2, ncol=1),
+    ms_gamma = matrix(c(-0.009, -0.009), nrow=2, ncol=1)
+  )
   expect_true(parameters$spraying)
   expect_equal(parameters$spraying_timesteps, c(5, 50))
   expect_equal(parameters$spraying_coverages, c(.5, .9))
@@ -34,7 +87,16 @@ test_that('set_spraying sets parameters', {
 test_that('distribute_bednets process sets net_time correctly', {
   timestep <- 50
   parameters <- get_parameters(list(human_population = 4))
-  parameters <- set_bednets(parameters, c(5, 50), c(.5, .9), 40)
+  parameters <- set_bednets(
+    parameters,
+    timesteps = c(5, 50),
+    coverages = c(.5, .9),
+    retention = 40,
+    dn0 = matrix(c(.533, .533), nrow=2, ncol=1),
+    rn = matrix(c(.56, .56), nrow=2, ncol=1),
+    rnm = matrix(c(.24, .24), nrow=2, ncol=1),
+    gamman = c(963.6, 963.6)
+  )
   events <- create_events(parameters)
   variables <- create_variables(parameters)
   variables$net_time <- mock_double(rep(0, 4))
@@ -71,7 +133,16 @@ test_that('distribute_bednets process sets net_time correctly', {
 test_that('throw_away_bednets process resets net_time correctly', {
   timestep <- 1
   parameters <- get_parameters(list(human_population = 4))
-  parameters <- set_bednets(parameters, c(5, 50), c(.5, .9), 40)
+  parameters <- set_bednets(
+    parameters,
+    timesteps = c(5, 50),
+    coverages = c(.5, .9),
+    retention = 40,
+    dn0 = matrix(c(.533, .533), nrow=2, ncol=1),
+    rn = matrix(c(.56, .56), nrow=2, ncol=1),
+    rnm = matrix(c(.24, .24), nrow=2, ncol=1),
+    gamman = c(963.6, 963.6)
+  )
   events <- create_events(parameters)
   variables <- create_variables(parameters)
   variables$net_time <- mock_double(rep(0, 4))
@@ -89,7 +160,17 @@ test_that('throw_away_bednets process resets net_time correctly', {
 test_that('indoor_spraying process sets spray_time correctly', {
   timestep <- 50
   parameters <- get_parameters(list(human_population = 4))
-  parameters <- set_spraying(parameters, c(5, 50), c(.5, .9))
+  parameters <- set_spraying(
+    parameters,
+    timesteps = c(5, 50),
+    coverages = c(.5, .9),
+    ls_theta = matrix(c(2.025, 2.025), nrow=2, ncol=1),
+    ls_gamma = matrix(c(-0.009, -0.009), nrow=2, ncol=1),
+    ks_theta = matrix(c(-2.222, -2.222), nrow=2, ncol=1),
+    ks_gamma = matrix(c(0.008, 0.008), nrow=2, ncol=1),
+    ms_theta = matrix(c(-1.232, -1.232), nrow=2, ncol=1),
+    ms_gamma = matrix(c(-0.009, -0.009), nrow=2, ncol=1)
+  )
   spray_time <- mock_double(rep(0, 4))
   correlations <- get_correlation_parameters(parameters)
   process <- indoor_spraying(
@@ -131,10 +212,17 @@ test_that('prob_bitten defaults to 1 with no protection', {
 
 test_that('prob_bitten correctly calculates net only probabilities', {
   timestep <- 100
-  parameters <- get_parameters(list(
-    bednets = TRUE,
-    gamman = 25
-  ))
+  parameters <- get_parameters()
+  parameters <- set_bednets(
+    parameters,
+    timesteps = c(5, 50, 100),
+    coverages = c(.5, .9, .2),
+    retention = 40,
+    dn0 = matrix(rep(.533, 3), nrow=3, ncol=1),
+    rn = matrix(rep(.56, 3), nrow=3, ncol=1),
+    rnm = matrix(rep(.24, 3), nrow=3, ncol=1),
+    gamman = rep(25, 3)
+  )
   variables <- create_variables(parameters)
   variables$net_time <- individual::DoubleVariable$new(
     c(-1, 5, 50, 100)
@@ -144,9 +232,9 @@ test_that('prob_bitten correctly calculates net only probabilities', {
   expect_equal(
     prob_bitten(timestep, variables, 1, parameters),
     list(
-      prob_bitten_survives = c(1, 0.7694168, 0.6836575, 0.0272300),
-      prob_bitten = c(1, 0.7694168, 0.6836575, 0.0272300),
-      prob_repelled = c(0, 0.2199712, 0.2521435, 0.4984000)
+      prob_bitten_survives = c(1, 0.7797801, 0.6978752, 0.0709500),
+      prob_bitten = c(1, 0.7797801, 0.6978752, 0.0709500),
+      prob_repelled = c(0, 0.2100848, 0.2408112, 0.4760000)
     ),
     tolerance = 1e-5
   )
@@ -154,20 +242,31 @@ test_that('prob_bitten correctly calculates net only probabilities', {
 
 test_that('prob_bitten correctly calculates spraying only probabilities', {
   timestep <- 100
-  parameters <- get_parameters(list(spraying = TRUE, gammas = 25))
+  parameters <- get_parameters(list(human_population = 4))
+  parameters <- set_spraying(
+    parameters,
+    timesteps = c(5, 50, 100),
+    coverages = c(.5, .9, .2),
+    ls_theta = matrix(rep(2.025, 3), nrow=3, ncol=1),
+    ls_gamma = matrix(rep(-0.009, 3), nrow=3, ncol=1),
+    ks_theta = matrix(rep(-2.222, 3), nrow=3, ncol=1),
+    ks_gamma = matrix(rep(0.008, 3), nrow=3, ncol=1),
+    ms_theta = matrix(rep(-1.232, 3), nrow=3, ncol=1),
+    ms_gamma = matrix(rep(-0.009, 3), nrow=3, ncol=1)
+  )
   variables <- create_variables(parameters)
 
-  variables$net_time <- individual::DoubleVariable$new(rep(-1, 4))
-  variables$spray_time <- individual::DoubleVariable$new(
+  variables$net_time <- individual::IntegerVariable$new(rep(-1, 4))
+  variables$spray_time <- individual::IntegerVariable$new(
     c(-1, 5, 50, 100)
   )
 
   expect_equal(
     prob_bitten(timestep, variables, 1, parameters),
     list(
-      prob_bitten_survives = c(1, 0.9780972, 0.8699070, 0.1751120),
-      prob_bitten = c(1, 0.9956601, 0.9737450, 0.8060000),
-      prob_repelled = c(0, 0.00433993, 0.02625504, 0.19400000)
+      prob_bitten_survives = c(1, 0.2216652, 0.1833448, 0.1506359),
+      prob_bitten = c(1, 0.8268157, 0.8101383, 0.7688352),
+      prob_repelled = c(0, 0.1731843, 0.1898617, 0.2311648)
     ),
     tolerance = 1e-5
   )
@@ -175,21 +274,42 @@ test_that('prob_bitten correctly calculates spraying only probabilities', {
 
 test_that('prob_bitten correctly combines spraying and net probabilities', {
   timestep <- 100
-  parameters <- get_parameters(list(bednets = TRUE, spraying = TRUE, gamman = 25))
+  parameters <- get_parameters(list(human_population = 4))
+  parameters <- set_bednets(
+    parameters,
+    timesteps = c(5, 50, 100),
+    coverages = c(.5, .9, .2),
+    retention = 40,
+    dn0 = matrix(rep(.533, 3), nrow=3, ncol=1),
+    rn = matrix(rep(.56, 3), nrow=3, ncol=1),
+    rnm = matrix(rep(.24, 3), nrow=3, ncol=1),
+    gamman = rep(25, 3)
+  )
+  parameters <- set_spraying(
+    parameters,
+    timesteps = c(5, 50, 100),
+    coverages = c(.5, .9, .2),
+    ls_theta = matrix(rep(2.025, 3), nrow=3, ncol=1),
+    ls_gamma = matrix(rep(-0.009, 3), nrow=3, ncol=1),
+    ks_theta = matrix(rep(-2.222, 3), nrow=3, ncol=1),
+    ks_gamma = matrix(rep(0.008, 3), nrow=3, ncol=1),
+    ms_theta = matrix(rep(-1.232, 3), nrow=3, ncol=1),
+    ms_gamma = matrix(rep(-0.009, 3), nrow=3, ncol=1)
+  )
   variables <- create_variables(parameters)
-  variables$net_time <- individual::DoubleVariable$new(
+  variables$net_time <- individual::IntegerVariable$new(
     c(100, 50, 5, -1)
   )
-  variables$spray_time <- individual::DoubleVariable$new(
+  variables$spray_time <- individual::IntegerVariable$new(
     c(-1, 5, 50, 100)
   )
 
   expect_equal(
     prob_bitten(timestep, variables, 1, parameters),
     list(
-      prob_bitten_survives = c(0.0272300, 0.4631212, 0.3765613, 0.1751120),
-      prob_bitten = c(0.0272300, 0.6375005, 0.6839200, 0.8060000),
-      prob_repelled = c(0.4984000, 0.3028339, 0.3066950, 0.1940000)
+      prob_bitten_survives = c(0.0709500, 0.1808229, 0.1629512, 0.1506359),
+      prob_bitten = c(0.0709500, 0.5828278, 0.6363754, 0.7688352),
+      prob_repelled = c(0.4760000, 0.3676569, 0.3556276, 0.2311648)
     ),
     tolerance=1e-4
   )
diff --git a/vignettes/VectorControl.Rmd b/vignettes/VectorControl.Rmd
index a2567aa9..e4ff54c8 100644
--- a/vignettes/VectorControl.Rmd
+++ b/vignettes/VectorControl.Rmd
@@ -38,17 +38,7 @@ simparams <- get_parameters(
     model_seasonality = TRUE, # Let's try a bi-modal model
     g0 = 0.28605,
     g = c(0.20636, -0.0740318, -0.0009293),
-    h = c(0.173743, -0.0730962, -0.116019),
-    # vector control parameters
-    rn = .31,
-    rnm = .24,
-    dn0 = .51,
-    phi_bednets = .89,
-    gamman = 2.64 * 365,
-    rs = .2,
-    phi_indoors = .97,
-    gammas = .5 * 365,
-    endophily = .813
+    h = c(0.173743, -0.0730962, -0.116019)
   )
 )
 
@@ -81,21 +71,26 @@ Then we can run the simulation for a variety of vector control strategies:
 
 ## Bed nets
 
-We can distribute bed nets once a year for 4 years, 3 months before peak season:
+We can distribute bed nets once a year for two years, changing the
+characteristics of the bed nets on each distribution...
 
 ```{r}
 bednetparams <- simparams
 
 bednet_events = data.frame(
   timestep = c(1, 2) * year,
-  name=c("Bednets start", "Bednets end")
+  name=c("Bednets 1", "Bednets 2")
 )
 
 bednetparams <- set_bednets(
   bednetparams,
-  timesteps = c(1, 2) * year,
-  coverages = rep(.8, 2),
-  retention = 5 * year
+  timesteps = bednet_events$timestep,
+  coverages = c(.8, .8),
+  retention = 5 * year,
+  dn0 = matrix(c(.533, .45), nrow=2, ncol=1),
+  rn = matrix(c(.56, .5), nrow=2, ncol=1),
+  rnm = matrix(c(.24, .24), nrow=2, ncol=1),
+  gamman = rep(2.64 * 365, 2)
 )
 
 output <- run_simulation(sim_length, bednetparams)
@@ -105,7 +100,7 @@ add_intervention_lines(plot_prevalence(output), bednet_events)
 
 ## Indoor spraying
 
-We can spray indoors once a year for 4 years, 3 months before peak season:
+We can do the same for IRS...
 
 ```{r}
 sprayingparams <- simparams
@@ -113,14 +108,19 @@ sprayingparams <- simparams
 peak <- peak_season_offset(sprayingparams)
 spraying_events = data.frame(
   timestep = c(1, 2) * year + peak - 3 * month,
-  name=c("Spraying starts", "Spraying ends")
+  name=c("Spraying 1", "Spraying 2")
 )
 
-sprayingparams <- set_bednets(
+sprayingparams <- set_spraying(
   sprayingparams,
-  timesteps = c(1, 2) * year,
+  timesteps = spraying_events$timestep,
   coverages = rep(.8, 2),
-  retention = 5 * year
+  ls_theta = matrix(c(2.025, 2.025), nrow=2, ncol=1),
+  ls_gamma = matrix(c(-0.009, -0.009), nrow=2, ncol=1),
+  ks_theta = matrix(c(-2.222, -2.222), nrow=2, ncol=1),
+  ks_gamma = matrix(c(0.008, 0.008), nrow=2, ncol=1),
+  ms_theta = matrix(c(-1.232, -1.232), nrow=2, ncol=1),
+  ms_gamma = matrix(c(-0.009, -0.009), nrow=2, ncol=1)
 )
 
 output <- run_simulation(sim_length, sprayingparams)