Feat/metapopulation model

R/biting_process.R

@@ -11,23 +11,23 @@
 #' @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_fn a function to retrieve the mixed EIR and infectivity based
+#' on the other populations
 #' @param mixing_index an index for this population's position in the
 #' lagged_infectivity list (default: 1)
 #' @noRd
 create_biting_process <- function(
-    renderer,
-    solvers,
-    models,
-    variables,
-    events,
-    parameters,
-    lagged_infectivity,
-    lagged_eir,
-    mixing = 1,
-    mixing_index = 1
-) {
+  renderer,
+  solvers,
+  models,
+  variables,
+  events,
+  parameters,
+  lagged_infectivity,
+  lagged_eir,
+  mixing_fn = NULL,
+  mixing_index = 1
+  ) {
   function(timestep) {
     # Calculate combined EIR
     age <- get_age(variables$birth$get_values(), timestep)
@@ -43,7 +43,7 @@ create_biting_process <- function(
       timestep,
       lagged_infectivity,
       lagged_eir,
-      mixing,
+      mixing_fn,
       mixing_index
     )
 
@@ -61,19 +61,19 @@ create_biting_process <- function(
 
 #' @importFrom stats rpois
 simulate_bites <- function(
-    renderer,
-    solvers,
-    models,
-    variables,
-    events,
-    age,
-    parameters,
-    timestep,
-    lagged_infectivity,
-    lagged_eir,
-    mixing = 1,
-    mixing_index = 1
-) {
+  renderer,
+  solvers,
+  models,
+  variables,
+  events,
+  age,
+  parameters,
+  timestep,
+  lagged_infectivity,
+  lagged_eir,
+  mixing_fn = NULL,
+  mixing_index = 1
+  ) {
   bitten_humans <- individual::Bitset$new(parameters$human_population)
 
   human_infectivity <- variables$infectivity$get_values()
@@ -111,7 +111,7 @@ simulate_bites <- function(
     f <- blood_meal_rate(s_i, Z, parameters)
     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(
         parameters$species[[s_i]]
@@ -129,16 +129,18 @@ simulate_bites <- function(
     }
 
     # 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
+    lagged_eir[[s_i]]$save(
+      n_infectious * a,
+      timestep
     )
-
+
+    # lagged EIR
+    if (is.null(mixing_fn)) {
+      species_eir <- lagged_eir[[s_i]]$get(timestep - parameters$de)
+    } else {
+      species_eir <- mixing_fn(timestep=timestep)$eir[[mixing_index, s_i]]
+    }
+
     renderer$render(paste0('EIR_', species_name), species_eir, timestep)
     EIR <- EIR + species_eir
     expected_bites <- species_eir * mean(psi)
@@ -150,16 +152,16 @@ simulate_bites <- function(
         )
       }
     }
-    
-    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)
+
+    if (is.null(mixing_fn)) {
+      infectivity <- lagged_infectivity$get(timestep - parameters$delay_gam)
+    } else {
+      infectivity <- mixing_fn(timestep=timestep)$inf[[mixing_index]]
+    }
+
+    lagged_infectivity$save(sum(human_infectivity * .pi), timestep)
+
+    foim <- calculate_foim(a, infectivity)
     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)
@@ -237,14 +239,13 @@ calculate_infectious <- function(species, solvers, variables, parameters) {
 }
 
 calculate_infectious_individual <- function(
-    species,
-    variables,
-    infectious_index,
-    adult_index,
-    species_index,
-    parameters
-) {
-
+  species,
+  variables,
+  infectious_index,
+  adult_index,
+  species_index,
+  parameters
+  ) {
   infectious_index$copy()$and(species_index)$size()
 }
 
@@ -300,10 +301,8 @@ unique_biting_rate <- function(age, parameters) {
 #' @title Calculate the force of infection towards mosquitoes
 #'
 #' @param a human blood meal 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
+#' @param infectivity_sum the sum of each individual's infectivity 
 #' @noRd
-calculate_foim <- function(a, infectivity_sum, mixing) {
-  a * sum(infectivity_sum * mixing)
+calculate_foim <- function(a, infectivity_sum) {
+  a * infectivity_sum
 }

R/mixing.R

@@ -0,0 +1,72 @@
+create_transmission_mixer <- function(
+  variables,
+  parameters,
+  lagged_eir,
+  lagged_infectivity,
+  mixing_tt,
+  export_mixing,
+  import_mixing,
+  p_captured_tt,
+  p_captured,
+  p_success
+  ) {
+  function (timestep) {
+    n_pops <- length(variables)
+    rdt_positive <- vnapply(
+      seq_len(n_pops),
+      function(i) {
+        rdt_detectable(variables[[i]], parameters[[i]], timestep)
+      }
+    )
+    mix_t <- match_timestep(mixing_tt, timestep)
+    p_mix_export <- export_mixing[[mix_t]]
+    p_mix_import <- import_mixing[[mix_t]]
+    p_captured_t <- p_captured[[match_timestep(p_captured_tt, timestep)]]
+    n_species <- length(parameters[[1]]$species)
+    species_eir <- t(vapply(
+      seq_along(lagged_eir),
+      function(i) {
+        vnapply(
+          lagged_eir[[i]],
+          function(l) l$get(timestep - parameters[[i]]$de)
+        )
+      },
+      numeric(n_species)
+    ))
+    if (n_species == 1) {
+      species_eir <- t(species_eir)
+    }
+    inf <- vnapply(
+      seq_along(lagged_infectivity),
+      function(i) lagged_infectivity[[i]]$get(timestep - parameters[[i]]$delay_gam)
+    )
+
+    test_and_treat_coeff <- (1 - p_captured_t * rdt_positive * p_success)
+    diag(test_and_treat_coeff) <- 1
+
+    eir <- vapply(
+      seq_len(n_species),
+      function(i) {
+        mixed_eir <- t(species_eir[,i] * p_mix_import)
+        rowSums(mixed_eir * test_and_treat_coeff)
+      },
+      numeric(n_pops)
+    )
+
+    inf <- rowSums(inf * p_mix_export * test_and_treat_coeff)
+
+    list(eir = eir, inf = inf)
+  }
+}
+
+# Estimates RDT prevalence from PCR prevalence,
+# We assume PCR prevalence is close to the true proportion of infectious people
+# in the population
+# values take from Wu et al 2015: https://doi.org/10.1038/nature16039
+rdt_detectable <- function(variables, parameters, timestep) {
+  infectious_prev <- variables$state$get_size_of(
+    c('D', 'A', 'U')) / parameters$human_population
+  logit_prev <- log(infectious_prev / (1 - infectious_prev))
+  logit_rdt <- parameters$rdt_intercept + parameters$rdt_coeff * logit_prev
+  exp(logit_rdt) / (1 + exp(logit_rdt))
+}