Coding in capacity to adapt housing to increase barrier to mosquitoes…

… and induce mortality with waning effect.

NAMESPACE

@@ -24,6 +24,7 @@ export(set_clinical_treatment)
 export(set_demography)
 export(set_drugs)
 export(set_equilibrium)
+export(set_housing)
 export(set_mass_pev)
 export(set_mda)
 export(set_parameter_draw)

R/processes.R

@@ -238,7 +238,12 @@ create_processes <- function(
   if (parameters$housing) {
     processes <- c(
       processes,
-      housing_improvement(variables$house_time, parameters, correlations)
+      housing_improvement(
+        variables,
+        parameters,
+        correlations
+      ),
+      house_usage_renderer(variables$house_time, renderer)
     )
   }
 

R/variables.R

@@ -26,7 +26,7 @@
 #' haven't been any
 #' * net_time - The timestep when a net was last put up (-1 if never)
 #' * spray_time - The timestep when the house was last sprayed (-1 if never)
-#' * house_time - The timestep when the house was last improved (-1 if never)
+#' * house_time - The timestep for adaptation to the house to reduce entry/kill vectors (-1 if never)
 #' * infectivity - The onward infectiousness to mosquitos
 #' * drug - The last prescribed drug
 #' * drug_time - The timestep of the last drug

R/vector_control.R

@@ -65,53 +65,56 @@ prob_bitten <- function(
       js_prime,
       parameters$k0
     )
+    spray_on = 1
     rs_comp <- 1 - rs
     ss <- rep(1, n)
     ss[protected_index] <- prob_survives_spraying(
       ks_prime,
       parameters$k0
     )
   } else {
-    # phi_indoors <- 0 ## if housing on too this has to be phi_indoors
+    spray_on = 0
     rs <- 0
     rs_comp <- 1
     ss <- 1
   }
 
   if (parameters$housing) {
+    phi_housing <- parameters$phi_housing[[species]]  ## if a change in behaviour is caused by the housing changes, we can increase outdoor biting proportion
+    phi_indoors <- parameters$phi_indoors[[species]]
     house_time <- variables$house_time$get_values()
-    since_house <- timestep - house_time
-    matches <- match(house_time, parameters$house_timesteps)
-    phi_housing <- parameters$phi_housing[[matches]]
-    rn_house <- parameters$rn_house[[matches]]
-    # phi_housing <- parameters$phi_housing[[matches]]
-    # rn_house <- prob_repelled_house(matches, since_house, species, parameters)
-    sn <- 1 - rn
+    since_housing <- timestep - house_time
+    matches <- match(house_time, parameters$housing_timesteps)
+    rh <- prob_repelled_house(matches, since_housing, species, parameters) ## if housing prevents entry to house, we increase proportion needing to repeat foraging
+    sh <- prob_survives_house(rh, matches, since_housing, species, parameters)
     unused <- house_time == -1
+    sh[unused] <- 1
+    rh[unused] <- 0
   } else {
     phi_housing <- 1
-    rn_house <- 0
+    rh <- 0
+    sh <- 1
   }
-  
-  if (!(parameters$housing & parameters$spraying)) {
+
+  if ((!parameters$housing & !parameters$spraying)) {
     phi_indoors <- 0 ## we want phi_indoors to be applied if housing is on
   }
-
+  
   list(
     prob_bitten_survives = (
-      1 - phi_indoors * phi_housing +
-        (1 - rn_house) * (phi_bednets * phi_housing * rs_comp * sn * ss) +
-        (1 - rn_house) * ((phi_indoors * phi_housing - phi_bednets * phi_housing) * rs_comp * ss)
+      1 - phi_indoors * phi_housing + ##
+        (1 - rh) * (phi_bednets * phi_housing * rs_comp * sn * ss * sh^2) +  ## * sh if some mortality from housing
+        (1 - rh) * ((phi_indoors * phi_housing - phi_bednets * phi_housing) * rs_comp * ss * sh^2) ## * sh if some mortality from housing
     ),
     prob_bitten = (
       1 - phi_indoors * phi_housing +
-        (1 - rn_house) * (phi_bednets * phi_housing * rs_comp * sn) +
-        (1 - rn_house) * ((phi_indoors * phi_housing - phi_bednets * phi_housing) * rs_comp)
+        (1 - rh) * (phi_bednets * phi_housing * rs_comp * sn * sh) + ## * sh if some mortality from housing
+        (1 - rh) * ((phi_indoors * phi_housing - phi_bednets * phi_housing) * rs_comp * sh) ## * sh if some mortality from housing
     ),
     prob_repelled = (
-      rn_house * phi_indoors * phi_housing + 
-        (1 - rn_house) * phi_bednets * phi_housing * rs_comp * rn + 
-        (1 - rn_house) * phi_indoors * phi_housing * rs
+      phi_bednets * phi_housing * rs_comp * sh * rn + 
+        phi_indoors * phi_housing * rs * sh * spray_on + 
+        phi_indoors * phi_housing * rh
     )
   )
 }
@@ -127,15 +130,15 @@ prob_bitten <- function(
 #' @noRd
 housing_improvement <- function(variables, parameters, correlations) {
   function(timestep) {
-    matches <- timestep == parameters$house_timesteps
+    matches <- timestep == parameters$housing_timesteps
     if (any(matches)) {
       target <- which(sample_intervention(
         seq(parameters$human_population),
         'housing',
-        parameters$house_coverages[matches],
+        parameters$housing_coverages[matches],
         correlations
       ))
-      house_time$queue_update(timestep, target)
+      variables$house_time$queue_update(timestep, target)
     }
   }
 }
@@ -231,10 +234,20 @@ spraying_decay <- function(t, theta, gamma) {
   1 / (1 + exp(-(theta + gamma * t)))
 }
 
-prob_repelled_house <- function(matches, dt, parameters) {
-  rn_house <- parameters$rn_house[matches]
-  rn_house <- dn0 * bednet_decay(dt, 1000)
-  rn_house    ## make this continual through time or we could have a decay as nets but very long lasting?
+house_decay <- function(t, gamma) {
+  exp(-t / gamma)
+}
+
+prob_repelled_house <- function(matches, dt, species, parameters) {
+  rhm <- parameters$house_rhm[matches, species]
+  gammah <- parameters$house_gammah[matches]
+  (parameters$house_rh[matches, species] - rhm) * house_decay(dt, gammah) + rhm
+}
+
+prob_survives_house <- function(rh, matches, dt, species, parameters) {
+  dh0 <- parameters$house_dh0[matches, species]
+  dh <- dh0 * house_decay(dt, parameters$house_gammah[matches])
+  1 - rh - dh
 }
 
 net_usage_renderer <- function(net_time, renderer) {
@@ -246,3 +259,13 @@ net_usage_renderer <- function(net_time, renderer) {
     )
   }
 }
+
+house_usage_renderer <- function(house_time, renderer) {
+  function(t) {
+    renderer$render(
+      'n_use_house_adaptation',
+      house_time$get_index_of(-1)$not(TRUE)$size(),
+      t
+    )
+  }
+}

R/vector_control_parameters.R

@@ -195,39 +195,45 @@ get_init_carrying_capacity <- function(parameters){
 #' Sherrard-Smith et al in prep
 #'
 #' @param parameters a list of parameters to modify
-#' @param timesteps the timesteps at which to improve the house
-#' @param coverages the proportion of the population who get a benefit from housing improvements
-#' @param phi_housing a parameter to potentially increase outdoor biting, default 1
-#' @param rn_house matrix of repellence parameters, will increase repelled, decrease indoor bites
+#' @param timesteps the timesteps at which to distribute housing adaptations
+#' @param coverages the proportion of the human population who reside in protected housing
+#' @param dh0 a matrix of death probabilities for each species over time.
 #' With nrows=length(timesteps), ncols=length(species)
+#' @param rh a matrix of repelling probabilities for each species over time
+#' With nrows=length(timesteps), ncols=length(species)
+#' @param rhm a matrix of minimum repelling probabilities for each species over time
+#' With nrows=length(timesteps), ncols=length(species)
+#' @param gammah a vector of house adaptation half-lives for each distribution timestep
 #' @export
 set_housing <- function(
     parameters,
     timesteps,
     coverages,
     phi_housing,
-    rn_house
+    dh0,
+    rh,
+    rhm,
+    gammah
 ) {
   stopifnot(all(coverages >= 0) && all(coverages <= 1))
-  if (length(coverages) != length(timesteps)) {
-    stop('coverages and timesteps must must align')
+  lengths <- vnapply(list(coverages, gammah), length)
+  if (!all(lengths == length(timesteps))) {
+    stop('timesteps and time-varying parameters must align')
   }
-  houses <- list(
-    phi_housing,
-    rn_house
-  )
-  for (x in houses) {
+  for (x in list(dh0, rh, rhm)) {
     if (ncol(x) != length(parameters$species)) {
-      stop('rn_house rows need to align with species')
+      stop('death and repelling probabilities rows need to align with species')
     }
     if (nrow(x) != length(timesteps)) {
-      stop('rn_house cols need to align with timesteps')
+      stop('death and repelling probabilities columns need to align with timesteps')
     }
   }
   parameters$housing <- TRUE
-  parameters$house_timesteps <- timesteps
-  parameters$house_coverages <- coverages
-  parameters$phi_housing <- phi_housing
-  parameters$rn_house <- rn_house
+  parameters$housing_timesteps <- timesteps
+  parameters$housing_coverages <- coverages
+  parameters$house_dh0 <- dh0
+  parameters$house_rh <- rh
+  parameters$house_rhm <- rhm
+  parameters$house_gammah <- gammah
   parameters
 }

vignettes/VectorControl_Housing.Rmd

@@ -1,5 +1,5 @@
 ---
-title: "Vector Control: Housing improvements"
+title: "Vector Control: Housing adaptations"
 output: 
   rmarkdown::html_vignette:
 vignette: >
@@ -24,7 +24,7 @@ library(malariasimulation)
 cols <- c("#E69F00", "#56B4E9", "#009E73", "#F0E442", "#0072B2", "#D55E00", "#CC79A7")
 ```
 
-Housing improvements are continually occurring in malaria-endemic parts of the world. Improvements could reduce the penetrability of housing to vectors. This will be a longer-lasting protective barrier effect than other 'indoor interventions' like ITNs or spraying. The model provides the flexibility to turn on 'improved' housing using a mechanism that adjusts the mosquito foraging pathway. The female mosquito now has some probability of being repelled by the adjustments to the house prior to any effects of spray or ITNs. The impacts from the three types of interventions are coded to be multiplicative and independent. When housing is 0, there is no impact. The model provides the user with the flexibility to specify parameters that describe the housing campaign (e.g. timing, coverage, and efficacy). It is currently assumed that housing improvements will represent a permanent change. We can illustrate this below:
+Housing adaptations are continually occurring in malaria-endemic parts of the world. Improvements could reduce the penetrability of housing to vectors. This will be a longer-lasting protective barrier effect than other 'indoor interventions' like ITNs or spraying. The model provides the flexibility to turn on 'improved' housing using a mechanism that adjusts the mosquito foraging pathway. The female mosquito now has some probability of being repelled by the adjustments to the house prior to any effects of spray or ITNs. The impacts from the three types of interventions are coded to be multiplicative and independent. When housing is 0, there is no impact. The model provides the user with the flexibility to specify parameters that describe the housing campaign (e.g. timing, coverage, and efficacy). It is currently assumed that housing improvements will represent a permanent change. We can illustrate this below:
 
 We can create a few plotting functions to visualise the output.
 ```{r}
@@ -36,12 +36,19 @@ plot_prev <- function() {
        xaxs = "i", yaxs = "i", ylim = c(0, 1))
   lines(x = output_control$timestep, y = output_control$n_detect_730_3650 / output_control$n_730_3650,
         col = cols[5], lwd = 1)
+  lines(x = output_nets_housing$timestep, y = output_nets_housing$n_detect_730_3650 / output_nets_housing$n_730_3650,
+        col = cols[1], lwd = 1)
   abline(v = bednetstimesteps, col = "black", lty = 2, lwd = 1)
+  abline(v = housingtimesteps, col = "darkred", lty = 2, lwd = 1)
+  
   text(x = bednetstimesteps + 10, y = 0.95, labels = "Bed net int.", adj = 0, cex = 0.8)
+  text(x = housingtimesteps + 10, y = 0.95, labels = "Housing int.", adj = 0, cex = 0.8)
   grid(lty = 2, col = "grey80", lwd = 0.5)
   legend("bottomleft", box.lty = 0, bg = "white",
-         legend = c("Prevalence for housing scenario","Prevalence for control scenario"),
-         col = c(cols[3], cols[5]), lty = c(1,1), lwd = 2, cex = 0.8, y.intersp = 1.3)
+         legend = c("Prevalence for housing scenario",
+                    "Prevalence for nets and housing scenario",
+                    "Prevalence for control scenario"),
+         col = c(cols[3], cols[4],cols[1],cols[5]), lty = c(1,1,1,1), lwd = 2, cex = 0.8, y.intersp = 1.3)
 }
 ```
 
@@ -80,31 +87,33 @@ simparams$species_proportions
 Having established a base set of parameters, we can now create a copy of this parameter list and update it to specify that housing has been altered in some way to reduce access to mosquitoes. In the example below, we simulate that 30% of housing is 'improved' in the 2nd year, and that nets are used by a random 50% of the population every three years. We can change the 'repellence' introduced from households multiple times and differently for each mosquito species potentially.
 
 ```{r, fig.align = 'center', out.width='100%'}
-housingtimesteps <- c(2, 2.5, 3) * year # housing will be improved at year 2 to year 3 in a stepped fashion
-
 housingparams <- set_housing(
   simparams,
   timesteps = housingtimesteps,
-  coverages = c(0.1,0.2,0.4),
-  phi_housing = matrix(c(0.9,0.9,1), nrow = 3, ncol = 1),
-  rn_house = matrix(c(0.2,0.3,0.5), nrow = 3, ncol = 1)
+  coverages = c(0.8,0.8,0.9),
+  phi_housing = c(1,1,1),
+  dh0 = matrix(c(.2, .2), nrow = 3, ncol = 2), 
+  rh = matrix(c(.5, .5), nrow = 3, ncol = 2), 
+  rhm = matrix(c(.3, .3), nrow = 3, ncol = 2), 
+  gammah = rep(0.5 * 365, 3) 
 )
+output <- run_simulation(timesteps = sim_length, parameters = housingparams)
 
 
 bednetstimesteps <- c(1, 4) * year # The bed nets will be distributed at the end of the first and the 4th year. 
 
 bednetparams <- set_bednets(
   housingparams,
   timesteps = bednetstimesteps,
-  coverages = c(.5, .8),  # The first round is distributed to 50% of the population, the second round to 80%.
-  retention = 5 * year, # Nets are kept on average 5 years
-  dn0 = matrix(c(.533, .533), nrow = 2, ncol = 1), # Matrix of death probabilities for each mosquito species over time
-  rn = matrix(c(.56, .56), nrow = 2, ncol = 1), # Matrix of repelling probabilities for each mosquito species over time
-  rnm = matrix(c(.24, .24), nrow = 2, ncol = 1), # Matrix of minimum repelling probabilities for each mosquito species over time
-  gamman = rep(2.64 * 365, 2) # Vector of bed net half-lives for each distribution timestep
+  coverages = c(.5, .8),  
+  retention = 5 * year, 
+  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 = rep(2.64 * 365, 2) 
 )
 
-output <- run_simulation(timesteps = sim_length, parameters = bednetparams)
+output_nets_housing <- run_simulation(timesteps = sim_length, parameters = bednetparams)
 ```
 
 ### Visualisation