Antiparasite immunity function shapes the duration of subpatent infec…

…tion/rate of recovery from subpatent infection in the Pv model. These values are now assigned during the variable creation and calculation of disease progression rates. Antiparasite immunity is labelled IAA and IAM for acquired and maternal immunities, where their decay, boosting, maternal immunity and reset at death are now included. The decay of the immunity occurs daily and has an impact of disease progression rate, which is now calculated and updated daily.

We must also adjust the prevalence rendering, which for p.f is based on detectability immunity, and for p.v is explicit in the human states.

`du` and ID variables and parameters are removed for vivax.

Tests are added for du rates and for prevalence rendering.

The vignette is extended to include the removal of ID and the sub-patent duration calculation.

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
     )
   }