Merge pull request #318 from mrc-ide/vivax_competing_hazards_main_sub…

…patent_duration_disease_progression

p.v sub-patent duration is immunity dependent.

R/disease_progression.R

@@ -6,18 +6,34 @@
 #' @param progression_outcome competing hazards object for disease progression rates
 #' @noRd
 create_progression_rates_process <- function(
+  parameters,
   variables,
   progression_outcome
 ) {
   function(timestep){
     target <- variables$state$get_index_of("S")$not()
+    progression_rates <- variables$progression_rates$get_values(target)
+    if (parameters$parasite == "vivax"){
+      # p.v subpatent recovery is immunity-dependent
+      progression_rates <- variables$progression_rates$get_values(target)
+      u_index <- variables$state$get_index_of("U")
+      target_u <- bitset_index(target, u_index)
+      progression_rates[target_u] <-
+        1 / anti_parasite_immunity(
+          min = parameters$dpcr_min, 
+          max = parameters$dpcr_max, 
+          a50 = parameters$apcr50, 
+          k = parameters$kpcr,
+          iaa = variables$iaa$get_values(index = u_index),
+          iam = variables$iam$get_values(index = u_index)
+        )
+    }
     progression_outcome$set_rates(
       target,
-      variables$progression_rates$get_values(target))
+      progression_rates)
   }
 }
 
-
 #' @title Disease progression outcomes
 #' @description Following resolution of competing hazards, update state and
 #' infectivity of sampled individuals

R/human_infection.R

@@ -172,13 +172,23 @@ infection_outcome_process <- function(
       timestep,
       parameters$uc
     )
-    boost_immunity(
-      variables$id,
-      infected_humans,
-      variables$last_boosted_id,
-      timestep,
-      parameters$ud
-    )
+    if(parameters$parasite == "falciparum"){
+      boost_immunity(
+        variables$id,
+        infected_humans,
+        variables$last_boosted_id,
+        timestep,
+        parameters$ud
+      )
+    } else if (parameters$parasite == "vivax"){
+      boost_immunity(
+        variables$iaa,
+        infected_humans,
+        variables$last_boosted_iaa,
+        timestep,
+        parameters$ua
+      )
+    }
   }
 
   clinical_infections <- calculate_clinical_infections(
@@ -624,3 +634,9 @@ blood_immunity <- function(ib, parameters) {
       (1 + (ib / parameters$ib0) ** parameters$kb)
   )
 }
+
+# Implemented from White et al., 2018 - Supplementary Information
+anti_parasite_immunity <- function(min, max, a50, k, iaa, iam){
+  min + (max - min) / (
+    1 + ((iaa + iam) / a50) ** k)
+}

R/mortality_processes.R

@@ -69,6 +69,10 @@ sample_maternal_immunity <- function(variables, target, timestep, parameters) {
         birth_ivm <- variables$iva$get_values(mothers) * parameters$pvm
         variables$icm$queue_update(birth_icm, target_group)
         variables$ivm$queue_update(birth_ivm, target_group)
+        if(parameters$parasite == "vivax"){
+          birth_iam <- variables$iaa$get_values(mothers) * parameters$pcm
+          variables$iam$queue_update(birth_iam, target_group)
+        }
       }
     }
   }
@@ -93,15 +97,19 @@ reset_target <- function(variables, events, target, state, parameters, timestep)
     # non-maternal immunity
     variables$last_boosted_ica$queue_update(-1, target)
     variables$last_boosted_iva$queue_update(-1, target)
-    variables$last_boosted_id$queue_update(-1, target)
     variables$ica$queue_update(0, target)
     variables$iva$queue_update(0, target)
-    variables$id$queue_update(0, target)
     variables$state$queue_update(state, target)
 
     if(parameters$parasite == "falciparum"){
       variables$last_boosted_ib$queue_update(-1, target)
+      variables$last_boosted_id$queue_update(-1, target)
       variables$ib$queue_update(0, target)
+      variables$id$queue_update(0, target)
+
+    } else if (parameters$parasite == "vivax"){
+      variables$last_boosted_iaa$queue_update(-1, target)
+      variables$iaa$queue_update(0, target)
     }
 
     # treatment

R/processes.R

@@ -49,17 +49,24 @@ create_processes <- function(
     immunity_process = create_exponential_decay_process(variables$ica,
                                                         parameters$rc),
     immunity_process = create_exponential_decay_process(variables$iva,
-                                                        parameters$rva),
-    immunity_process = create_exponential_decay_process(variables$id,
-                                                        parameters$rid)
+                                                        parameters$rva)
   )
 
   if(parameters$parasite == "falciparum"){
     processes <- c(
       processes,
       # Blood immunity
       immunity_process = create_exponential_decay_process(variables$ib,
-                                                          parameters$rb)
+                                                          parameters$rb),
+      # Immunity to detectability
+      immunity_process = create_exponential_decay_process(variables$id, parameters$rid)
+    )
+  } else if (parameters$parasite == "vivax"){
+    processes <- c(
+      processes,
+      # Anti-parasite immunity
+      immunity_process = create_exponential_decay_process(variables$iam, parameters$rm),
+      immunity_process = create_exponential_decay_process(variables$iaa, parameters$ra)
     )
   }
 
@@ -126,6 +133,7 @@ create_processes <- function(
   processes <- c(
     processes,
     progression_process = create_progression_rates_process(
+      parameters,
       variables,
       progression_outcome
     ),
@@ -185,9 +193,11 @@ create_processes <- function(
   # Rendering
   # =========
 
-  imm_var_names <- c('ica', 'icm', 'id', 'iva', 'ivm')
+  imm_var_names <- c('ica', 'icm', 'iva', 'ivm')
   if(parameters$parasite == "falciparum"){
-    imm_var_names <- c(imm_var_names, 'ib')
+    imm_var_names <- c(imm_var_names, 'ib', 'id')
+  } else if (parameters$parasite == "vivax"){
+    imm_var_names <- c(imm_var_names, 'iaa', 'iam')
   }
 
   processes <- c(

R/render.R

@@ -6,9 +6,10 @@ in_age_range <- function(birth, timestep, lower, upper) {
 #' 
 #' @description renders prevalence numerators and denominators for individuals
 #' detected by light microscopy and pcr, where those infected asymptomatically by
-#' P. falciparum have reduced probability of infection due to detectability
+#' p.f has reduced probability of infection due to detectability
 #' immunity (reported as an integer sample: n_detect_lm, and summing over
-#' detection probabilities: p_detect_lm)
+#' detection probabilities: p_detect_lm), whearas p.v human states are defined
+#' explicitly by lm and pcr detectability.
 #' 
 #' @param state human infection state
 #' @param birth variable for birth of the individual
@@ -26,12 +27,17 @@ create_prevalence_renderer <- function(
 ) {
   function(timestep) {
     asymptomatic <- state$get_index_of('A')
-    prob <- probability_of_detection(
-      get_age(birth$get_values(asymptomatic), timestep),
-      immunity$get_values(asymptomatic),
-      parameters
-    )
-    asymptomatic_detected <- bitset_at(asymptomatic, bernoulli_multi_p(prob))
+
+    if(parameters$parasite == "falciparum"){
+      prob <- probability_of_detection(
+        get_age(birth$get_values(asymptomatic), timestep),
+        immunity$get_values(asymptomatic),
+        parameters
+      )
+      asymptomatic_detected <- bitset_at(asymptomatic, bernoulli_multi_p(prob))
+    } else if (parameters$parasite == "vivax") {
+      asymptomatic_detected <- asymptomatic
+    }
 
     clinically_detected <- state$get_index_of(c('Tr', 'D'))
     detected <- clinically_detected$copy()$or(asymptomatic_detected)
@@ -46,13 +52,15 @@ create_prevalence_renderer <- function(
         in_age$copy()$and(detected)$size(),
         timestep
       )
-      renderer$render(
-        paste0('p_detect_lm_', lower, '_', upper),
-        in_age$copy()$and(clinically_detected)$size() + sum(
-          prob[bitset_index(asymptomatic, in_age)]
-        ),
-        timestep
-      )
+      if(parameters$parasite == "falciparum"){
+        renderer$render(
+          paste0('p_detect_lm_', lower, '_', upper),
+          in_age$copy()$and(clinically_detected)$size() + sum(
+            prob[bitset_index(asymptomatic, in_age)]
+          ),
+          timestep
+        )
+      }
       renderer$render(
         paste0('n_detect_pcr_', lower, '_', upper),
         pcr_detected$copy()$and(in_age)$size(),

R/variables.R

@@ -10,12 +10,14 @@
 #' * birth - an integer representing the timestep when this individual was born
 #' * last_boosted_* - the last timestep at which this individual's immunity was
 #' boosted for tracking grace periods in the boost of immunity
+#' * IAM - Maternal anti-parasite immunity (p.v only)
 #' * ICM - Maternal immunity to clinical disease
-#' * IVM - Maternal immunity to severe disease
+#' * IVM - Maternal immunity to severe disease (p.f only)
 #' * IB  - Pre-erythrocytic immunity (p.f only)
+#' * IAA  - Acquired anti-parasite immunity (p.v only)
 #' * ICA  - Acquired immunity to clinical disease
-#' * IVA  - Acquired immunity to severe disease
-#' * ID - Acquired immunity to detectability
+#' * IVA  - Acquired immunity to severe disease (p.f only)
+#' * ID - Acquired immunity to detectability (p.f only)
 #' * zeta - Heterogeneity of human individuals
 #' * zeta_group - Discretised heterogeneity of human individuals
 #' * last_pev_timestep - The timestep of the last pev vaccination (-1 if there
@@ -96,9 +98,9 @@ create_variables <- function(parameters) {
   initial_states <- initial_state(parameters, initial_age, groups, eq, states)
   state <- individual::CategoricalVariable$new(states, initial_states)
   birth <- individual::IntegerVariable$new(-initial_age)
+
   last_boosted_ica <- individual::DoubleVariable$new(rep(-1, size))
   last_boosted_iva <- individual::DoubleVariable$new(rep(-1, size))
-  last_boosted_id <- individual::DoubleVariable$new(rep(-1, size))
 
   # Maternal immunity
   icm <- individual::DoubleVariable$new(
@@ -136,8 +138,47 @@ create_variables <- function(parameters) {
         'IB'
       )
     )
+
+    # Acquired immunity to detectability
+    last_boosted_id <- individual::DoubleVariable$new(rep(-1, size))
+    id <- individual::DoubleVariable$new(
+      initial_immunity(
+        parameters$init_id,
+        initial_age,
+        groups,
+        eq,
+        parameters,
+        'ID'
+      )
+    )
+
+  } else if (parameters$parasite == "vivax"){
+    # Acquired anti-parasite immunity
+    last_boosted_iaa <- individual::DoubleVariable$new(rep(-1, size))
+    iaa <- individual::DoubleVariable$new(
+      initial_immunity(
+        parameters$init_iaa,
+        initial_age,
+        groups,
+        eq,
+        parameters,
+        'IAA'
+      )
+    )
+
+    # Maternal anti-parasite immunity
+    iam <- individual::DoubleVariable$new(
+      initial_immunity(
+        parameters$init_iam,
+        initial_age,
+        groups,
+        eq,
+        parameters,
+        'IAM'
+      )
+    )
   }
-  
+
   # Acquired immunity to clinical disease
   ica <- individual::DoubleVariable$new(
     initial_immunity(
@@ -149,6 +190,7 @@ create_variables <- function(parameters) {
       'ICA'
     )
   )
+
   # Acquired immunity to severe disease
   iva <- individual::DoubleVariable$new(
     initial_immunity(
@@ -160,18 +202,7 @@ create_variables <- function(parameters) {
       'IVA'
     )
   )
-  # Acquired immunity to detectability
-  id <- individual::DoubleVariable$new(
-    initial_immunity(
-      parameters$init_id,
-      initial_age,
-      groups,
-      eq,
-      parameters,
-      'ID'
-    )
-  )
-
+
   # Initialise infectiousness of humans -> mosquitoes
   # NOTE: not yet supporting initialisation of infectiousness of Treated individuals
   infectivity_values <- rep(0, get_human_population(parameters, 0))
@@ -205,7 +236,17 @@ create_variables <- function(parameters) {
   progression_rate_values <- rep(0, get_human_population(parameters, 0))
   progression_rate_values[diseased] <- 1/parameters$dd
   progression_rate_values[asymptomatic] <- 1/parameters$da
-  progression_rate_values[subpatent] <- 1/parameters$du
+  if(parameters$parasite == "falciparum"){
+    # p.f subpatent recovery rate is constant
+    progression_rate_values[subpatent] <- 1/parameters$du
+  } else if (parameters$parasite == "vivax"){
+    # p.v subpatent recovery rate is immunity-dependent
+    progression_rate_values[subpatent] <- 1/anti_parasite_immunity(
+      parameters$dpcr_min, parameters$dpcr_max, parameters$apcr50, parameters$kpcr,
+      iaa$get_values(subpatent),
+      iam$get_values(subpatent)
+    )
+  }
   progression_rate_values[treated] <- 1/parameters$dt
 
   # Initialise the disease progression rate variable
@@ -229,12 +270,10 @@ create_variables <- function(parameters) {
     birth = birth,
     last_boosted_ica = last_boosted_ica,
     last_boosted_iva = last_boosted_iva,
-    last_boosted_id = last_boosted_id,
     icm = icm,
     ivm = ivm,
     ica = ica,
     iva = iva,
-    id = id,
     zeta = zeta,
     zeta_group = zeta_group,
     infectivity = infectivity,
@@ -252,7 +291,15 @@ create_variables <- function(parameters) {
   if(parameters$parasite == "falciparum"){
     variables <- c(variables,
                    last_boosted_ib = last_boosted_ib,
-                   ib = ib
+                   last_boosted_id = last_boosted_id,
+                   ib = ib,
+                   id = id
+    )
+  } else if (parameters$parasite == "vivax"){
+    variables <- c(variables,
+                   last_boosted_iaa = last_boosted_iaa,
+                   iaa = iaa,
+                   iam = iam
     )
   }