Added in functionality to start at equilibrium, based on Michael Whit…

…e's equilibrium solution. This is accessed during the initial equilibrium calculation in compatability and during the variable creation. The equilibrium solution is found in malariaEquilibriumVivax.

Note that we draw individuals from a joint probability distribution over age and het group to generate human state, immunities and number of hypnozoite batches.

tests included to check correct equilibrium is generated and maintained.

I have also added the An. koliensis parameters as found in Michael White's repo in order to compare equilibrium solutions.

A few errors were identified. Total infection rates were not being calculated correctly in human_infection.R. Lagged infectivity needs to be saved before the current infectivity is extracted during the biting_process.R. Also, when immunity is boosted through infection, this should incorporate the decrease in immunity from decay (there may be better ways of achieving this).

DESCRIPTION

@@ -68,10 +68,12 @@ Encoding: UTF-8
 LazyData: true
 Remotes:
   mrc-ide/malariaEquilibrium,
+  mrc-ide/malariaEquilibriumVivax,
   mrc-ide/individual
 Imports:
     individual (>= 0.1.17),
     malariaEquilibrium (>= 1.0.1),
+    malariaEquilibriumVivax,
     Rcpp,
     statmod,
     MASS,

NAMESPACE

@@ -10,6 +10,7 @@ export(gamb_params)
 export(get_correlation_parameters)
 export(get_init_carrying_capacity)
 export(get_parameters)
+export(kol_params)
 export(parameterise_mosquito_equilibrium)
 export(parameterise_total_M)
 export(peak_season_offset)

R/biting_process.R

@@ -169,14 +169,14 @@ simulate_bites <- function(
       }
     }
 
+    lagged_infectivity$save(sum(human_infectivity * .pi), timestep)
+
     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)

R/compatibility.R

@@ -4,6 +4,10 @@ inverse_param <- function(name, new_name) {
   function(params) { list(new_name, 1 / params[[name]]) }
 }
 
+product_param <- function(new_name, name_1, name_2) {
+  function(params) { list(new_name, params[[name_1]] * params[[name_2]]) }
+}
+
 mean_param <- function(new_name, name, weights) {
   function(params) {
     list(new_name, weighted.mean(params[[name]], params[[weights]]))
@@ -120,6 +124,56 @@ back_translations = list(
   kv = 'kv'
 )
 
+vivax_translations = list(
+
+  mean_age  = 'average_age',
+  rho_age   = 'rho',
+  age_0     = 'a0',
+  N_het     = 'n_heterogeneity_groups',
+
+  bb  = 'b',      ## mosquito -> human transmission probability
+
+  c_PCR = 'cu',   ## human -> mosquito transmission probability (PCR)
+  c_LM  = 'ca',   ## human -> mosquito transmission probability (LM-detectable)
+  c_D   = 'cd',   ## human -> mosquito transmission probability (disease state)
+  c_T   = 'ct',   ## human -> mosquito transmission probability (treatment)
+
+  d_E   = 'de',   ## duration of liver-stage latency
+  r_D   = inverse_param('dd', 'r_D'),   ## duraton of disease = 1/rate
+  r_T   = inverse_param('dt', 'r_T'),   ## duraton of prophylaxis = 1/rate
+
+  r_par = inverse_param('ra', 'r_par'),    ## rate of decay of anti-parasite immunity
+  r_clin= inverse_param('rc', 'r_clin'),    ## rate of decay of clinical immunity
+
+  mu_M  = mean_param('mum', 'mum', 'species_proportions'),  ## mosquito death rate = 1/(mosquito life expectancy)
+  Q0    = mean_param('Q0', 'Q0', 'species_proportions'),
+  blood_meal_rates = mean_param('blood_meal_rates', 'blood_meal_rates', 'species_proportions'),
+  tau_M = 'dem', 	## duration of sporogony
+
+  ff      = 'f',      ## relapse rate
+  gamma_L = 'gammal', ## duration of liver-stage carriage
+  K_max   = 'kmax', ## maximum hypnozoite batches
+
+  u_par        = 'ua',        ## refractory period for anti-parasite immune boosting
+  phi_LM_max   = 'philm_max', ## probability of LM_detectable infection with no immunity
+  phi_LM_min   = 'philm_min', ## probability of LM_detectable infection with maximum immunity
+  A_LM_50pc    = 'alm50',     ## blood-stage immunity scale parameter
+  K_LM         = 'klm',       ## blood-stage immunity shape parameter
+  u_clin       = 'uc',        ## refractory period for clinical immune boosting
+  phi_D_max    = 'phi0',      ## probability of clinical episode with no immunity
+  phi_D_min    = product_param("phi_D_min", "phi0", "phi1"),    ## probability of clinical episode with maximum immunity
+  A_D_50pc     = 'ic0',       ## clinical immunity scale parameter
+  K_D          = 'kc',        ## clinical immunity shape parameter
+  A_d_PCR_50pc = 'apcr50',    ## scale parameter for effect of anti-parasite immunity on PCR-detectable infection
+  K_d_PCR      = 'kpcr',      ## shape parameter for effect of anti-parasite immunity on PCR-detectable infection
+  d_PCR_max    = 'dpcr_max',  ## maximum duration on PCR-detectable infection
+  d_PCR_min    = 'dpcr_min',  ## maximum duration of PCR-detectable infection
+  d_LM         = 'da',        ## duration of LM-detectable infection
+  P_MI         = 'pcm',       ## Proportion of immunity acquired maternally
+  d_MI         = 'rm'         ## Rate of waning of maternal immunity
+
+)
+
 #' @description translate parameter keys from the malariaEquilibrium format
 #' to ones compatible with this IBM 
 #' @param params with keys in the malariaEquilibrium format
@@ -162,6 +216,23 @@ translate_parameters <- function(params) {
   translated
 }
 
+#' @description translate parameter keys from the malariaVivaxEquilibrium format
+#' to ones compatible with this IBM
+#' @param params with keys in the malariaVivaxEquilibrium format
+#' @noRd
+translate_vivax_parameters <- function(params) {
+  translated <- params
+  for (i in 1:length(vivax_translations)) {
+    if (is.character(vivax_translations[[i]])) {
+      translated[[names(vivax_translations)[i]]] <- params[[vivax_translations[[i]]]]
+    }
+    if (is.function(vivax_translations[[i]])) {
+      translated[[names(vivax_translations)[i]]] <- vivax_translations[[i]](params)[[2]]
+    }
+  }
+  translated
+}
+
 #' @title remove parameter keys from the malariaEquilibrium format that are not used
 #' in this IBM 
 #' @param params with keys in the malariaEquilibrium format
@@ -182,36 +253,58 @@ remove_unused_equilibrium <- function(params) {
 #' @title Set equilibrium
 #' @description This will update the IBM parameters to match the
 #' equilibrium parameters and set up the initial human and mosquito population
-#' to acheive init_EIR
+#' to achieve init_EIR
 #' @param parameters model parameters to update
 #' @param init_EIR the desired initial EIR (infectious bites per person per day over the entire human
 #' population)
 #' @param eq_params parameters from the malariaEquilibrium package, if null.
 #' The default malariaEquilibrium parameters will be used
 #' @export
 set_equilibrium <- function(parameters, init_EIR, eq_params = NULL) {
-  if (is.null(eq_params)) {
-    eq_params <- translate_parameters(parameters)
-  } else {
+  if(parameters$parasite == "falciparum"){
+    if (is.null(eq_params)) {
+      eq_params <- translate_parameters(parameters)
+    } else {
+      parameters <- c(
+        translate_equilibrium(remove_unused_equilibrium(eq_params)),
+        parameters
+      )
+    }
+    eq <- malariaEquilibrium::human_equilibrium(
+      EIR = init_EIR,
+      ft = sum(get_treatment_coverages(parameters, 1)),
+      p = eq_params,
+      age = EQUILIBRIUM_AGES,
+      h = malariaEquilibrium::gq_normal(parameters$n_heterogeneity_groups)
+    )
+    parameters <- c(
+      list(
+        init_foim = eq$FOIM,
+        init_EIR = init_EIR,
+        eq_params = eq_params
+      ),
+      parameters
+    )
+  } else if (parameters$parasite == "vivax"){
+
+    if (!(is.null(eq_params))) {
+      stop("Importing MalariaEquilibriumVivax parameters is not supported")
+    }
+
+    eq <- malariaEquilibriumVivax::vivax_equilibrium(
+      EIR = init_EIR,
+      ft = sum(get_treatment_coverages(parameters, 1)),
+      p = translate_vivax_parameters(parameters),
+      age = EQUILIBRIUM_AGES
+    )
+
     parameters <- c(
-      translate_equilibrium(remove_unused_equilibrium(eq_params)),
+      list(
+        init_foim = eq$FOIM,
+        init_EIR = init_EIR
+      ),
       parameters
     )
   }
-  eq <- malariaEquilibrium::human_equilibrium(
-    EIR = init_EIR,
-    ft = sum(get_treatment_coverages(parameters, 1)),
-    p = eq_params,
-    age = EQUILIBRIUM_AGES,
-    h = malariaEquilibrium::gq_normal(parameters$n_heterogeneity_groups)
-  )
-  parameters <- c(
-    list(
-      init_foim = eq$FOIM,
-      init_EIR = init_EIR,
-      eq_params = eq_params
-    ),
-    parameters
-  )
   parameterise_mosquito_equilibrium(parameters, init_EIR)
 }

R/human_infection.R

@@ -119,7 +119,6 @@ calculate_infections <- function(
       )
     }
 
-    # infection_rates <- rep(0, length = parameters$human_population)
     if(parameters$parasite == "falciparum"){
 
       ## p.f models blood immunity
@@ -205,15 +204,16 @@ infection_outcome_process <- function(
       timestep
     )
 
-    boost_immunity(
-      variables$ica,
-      infected_humans,
-      variables$last_boosted_ica,
-      timestep,
-      parameters$uc
-    )
-
     if(parameters$parasite == "falciparum"){
+
+      boost_immunity(
+        variables$ica,
+        infected_humans,
+        variables$last_boosted_ica,
+        timestep,
+        parameters$uc
+      )
+
       boost_immunity(
         variables$id,
         infected_humans,
@@ -261,13 +261,23 @@ infection_outcome_process <- function(
         renderer,
         timestep
       )
-
+
+      boost_immunity(
+        variables$ica,
+        infected_humans,
+        variables$last_boosted_ica,
+        timestep,
+        parameters$uc,
+        1/parameters$rc
+      )
+
       boost_immunity(
         variables$iaa,
         infected_humans,
         variables$last_boosted_iaa,
         timestep,
-        parameters$ua
+        parameters$ua,
+        1/parameters$ra
       )
 
       ## Only S and U infections are considered in generating lm-det infections
@@ -744,7 +754,8 @@ boost_immunity <- function(
     exposed_index,
     last_boosted_variable,
     timestep,
-    delay
+    delay,
+    decay_rate = 0
 ) {
   # record who can be boosted
   exposed_index_vector <- exposed_index$to_vector()
@@ -754,7 +765,7 @@ boost_immunity <- function(
   if (sum(to_boost) > 0) {
     # boost the variable
     immunity_variable$queue_update(
-      immunity_variable$get_values(exposed_to_boost) + 1,
+      immunity_variable$get_values(exposed_to_boost) * (1-decay_rate) + 1,
       exposed_to_boost
     )
     # record last boosted

R/processes.R

@@ -47,6 +47,27 @@ create_processes <- function(
     immunity_process = create_exponential_decay_process(variables$ica,
                                                         parameters$rc)
   )
+
+  if(parameters$parasite == "falciparum"){
+    processes <- c(
+      processes,
+      # Blood immunity
+      create_exponential_decay_process(variables$ib, parameters$rb),
+      # Immunity to severe disease
+      create_exponential_decay_process(variables$ivm, parameters$rvm),
+      create_exponential_decay_process(variables$iva, parameters$rva),
+      # Immunity to detectability
+      create_exponential_decay_process(variables$id, parameters$rid)
+    )
+  } else if (parameters$parasite == "vivax"){
+    processes <- c(
+      processes,
+      # Anti-parasite immunity
+      create_exponential_decay_process(variables$iam, parameters$rm),
+      create_exponential_decay_process(variables$iaa, parameters$ra),
+      create_hypnozoite_batch_decay_process(variables$hypnozoites, parameters$gammal)
+    )
+  }
 
   if(parameters$parasite == "falciparum"){
     processes <- c(