From 5f5260a38d0fb5460c86fc25cdc6c5709c04952c Mon Sep 17 00:00:00 2001
From: RJSheppard <rjs11@ic.ac.uk>
Date: Tue, 17 Sep 2024 19:02:42 +0100
Subject: [PATCH] Hynozoite batches are now modelled, including formation via
 bite infections, decay of hypnozoite batches, new infections via relapse and
 the reseting of batch number to 0 for new births.

A hypnozoite batch variable is created.

The vivax model considers the number of new bites an individual receives (in contrast to whether an individual has been bitten regardless of number of bites, as currently in the P. falciparum model) to calculate the infection rates. This is now captured in bitten_humans, which is now a list that also includes n_bites_each.

The vivax model does not factor the eir by the population level biting heterogeneity (although it does use the heterogeneity to distribute biting, but retains the overall total biting according to the EIR). biting_rate.R

In vivax, we track bites in all individuals, not just SAU as in falciparum. This is because new hypnozoite batches can occur in any human state, regardless of whether this results in a new infection.

The number of bites is considered in calculating probability of a new infection.

To calculate the infection rates for vivax, we take the biting infection rates and add them to the relapse rate. We also save the relative rates in the competing hazards object to resolve these competing hazards later.

A new function to resolve between bites and relapse infections is created. For each bite we increase the number of batches - with capacity for prophylaxis to be added at a future step. We also cap the total number of batches.

We add a hypnozoite batch decay process. Note that this process differs from the immunity decay processes. We work in discrete outputs (you can't have a fractional batch). We also work in the probability that an individual experiences a single loss, rather than loss at the hypnozoite level.

When an individual dies, we reset the hypnozoite batches back to 0.

We now render number of relapses, average hypnozoite batch number and the number of individuals with hypnozoite batches. We may also output these by specified age groups.

Added a test to check relapses occur as expected and expanded the p.v vignette.
---
 R/biting_process.R                          |  30 ++-
 R/competing_hazards.R                       |  11 +-
 R/human_infection.R                         | 267 ++++++++++++++------
 R/mortality_processes.R                     |   1 +
 R/parameters.R                              |   4 +-
 R/processes.R                               |  68 +++--
 R/render.R                                  |  40 +++
 R/variables.R                               |   9 +-
 man/get_parameters.Rd                       |   4 +-
 tests/testthat/test-biting-integration.R    |  10 +-
 tests/testthat/test-infection-integration.R |  11 +-
 tests/testthat/test-pev.R                   |   1 +
 tests/testthat/test-vivax.R                 |  60 +++++
 vignettes/Plasmodium_vivax.Rmd              |   6 +
 14 files changed, 394 insertions(+), 128 deletions(-)

diff --git a/R/biting_process.R b/R/biting_process.R
index df494604..17b75366 100644
--- a/R/biting_process.R
+++ b/R/biting_process.R
@@ -34,7 +34,7 @@ create_biting_process <- function(
   function(timestep) {
     # Calculate combined EIR
     age <- get_age(variables$birth$get_values(), timestep)
-    bitten_humans <- simulate_bites(
+    bitten <- simulate_bites(
       renderer,
       solvers,
       models,
@@ -52,7 +52,8 @@ create_biting_process <- function(
     simulate_infection(
       variables,
       events,
-      bitten_humans,
+      bitten$bitten_humans,
+      bitten$n_bites_per_person,
       age,
       parameters,
       timestep,
@@ -78,6 +79,7 @@ simulate_bites <- function(
   mixing_index = 1
   ) {
   bitten_humans <- individual::Bitset$new(parameters$human_population)
+  n_bites_per_person <- numeric(0)
   
   human_infectivity <- variables$infectivity$get_values()
   if (parameters$tbv) {
@@ -146,13 +148,24 @@ simulate_bites <- function(
 
     renderer$render(paste0('EIR_', species_name), species_eir, timestep)
     EIR <- EIR + species_eir
-    expected_bites <- species_eir * mean(psi)
+    if(parameters$parasite == "falciparum"){
+      # p.f model factors eir by psi
+      expected_bites <- species_eir * mean(psi)
+    } else if (parameters$parasite == "vivax"){
+      # p.v model standardises biting rate het to eir
+      expected_bites <- species_eir
+    }
+    
     if (expected_bites > 0) {
       n_bites <- rpois(1, expected_bites)
       if (n_bites > 0) {
-        bitten_humans$insert(
-          fast_weighted_sample(n_bites, lambda)
-        )
+        bitten <- fast_weighted_sample(n_bites, lambda)
+        bitten_humans$insert(bitten)
+        renderer$render('n_bitten', bitten_humans$size(), timestep)
+        if(parameters$parasite == "vivax"){
+          # p.v must pass through the number of bites per person
+          n_bites_per_person <- tabulate(bitten, nbins = length(lambda))
+        }
       }
     }
 
@@ -202,9 +215,8 @@ simulate_bites <- function(
       )
     }
   }
-  
-  renderer$render('n_bitten', bitten_humans$size(), timestep)
-  bitten_humans
+
+  list(bitten_humans = bitten_humans, n_bites_per_person = n_bites_per_person)
 }
 
 
diff --git a/R/competing_hazards.R b/R/competing_hazards.R
index 13c522ad..5a44e3e6 100644
--- a/R/competing_hazards.R
+++ b/R/competing_hazards.R
@@ -16,6 +16,7 @@ CompetingOutcome <- R6::R6Class(
 
       self$target <- individual::Bitset$new(size)
       self$rates <- NULL
+      self$relative_rates <- NULL
     },
     set_rates = function(target, rates){
       stopifnot(target$size() == length(rates))
@@ -23,15 +24,23 @@ CompetingOutcome <- R6::R6Class(
       self$target$copy_from(target)
       self$rates <- rates
     },
+    set_relative_rates = function(target, relative_rates){
+      stopifnot(target$size() == length(relative_rates))
+      
+      self$target$copy_from(target)
+      self$relative_rates <- relative_rates
+    },
     execute = function(t, target){
       private$targeted_process(t, target)
     },
     reset = function() {
       self$target$clear()
       self$rates <- NULL
+      self$relative_rates <- NULL
     },
     target = NULL,
-    rates = NULL
+    rates = NULL,
+    relative_rates = NULL
   )
 )
 
diff --git a/R/human_infection.R b/R/human_infection.R
index 33e8614d..f1af20f5 100644
--- a/R/human_infection.R
+++ b/R/human_infection.R
@@ -5,6 +5,7 @@
 #' @param variables a list of all of the model variables
 #' @param events a list of all of the model events
 #' @param bitten_humans a bitset of bitten humans
+#' @param n_bites_per_person vector of number of bites each person receives (p.v only)
 #' @param age of each human (timesteps)
 #' @param parameters of the model
 #' @param timestep current timestep
@@ -15,6 +16,7 @@ simulate_infection <- function(
     variables,
     events,
     bitten_humans,
+    n_bites_per_person,
     age,
     parameters,
     timestep,
@@ -37,6 +39,7 @@ simulate_infection <- function(
   calculate_infections(
     variables,
     bitten_humans,
+    n_bites_per_person,
     parameters,
     renderer,
     timestep,
@@ -48,6 +51,7 @@ simulate_infection <- function(
 #' @description Infection rates are stored in the infection outcome competing hazards object
 #' @param variables a list of all of the model variables
 #' @param bitten_humans bitset of bitten humans
+#' @param n_bites_per_person vector of number of bites each person receives (p.v only)
 #' @param parameters model parameters
 #' @param renderer model render object
 #' @param timestep current timestep
@@ -55,84 +59,116 @@ simulate_infection <- function(
 calculate_infections <- function(
     variables,
     bitten_humans,
+    n_bites_per_person,
     parameters,
     renderer,
     timestep,
     infection_outcome
 ) {
-  source_humans <- variables$state$get_index_of(
-    c('S', 'A', 'U'))$and(bitten_humans)
 
-  if(parameters$parasite == "falciparum"){
-    ## p.f models blood immunity
-    b <- blood_immunity(variables$ib$get_values(source_humans), parameters)
+    if(parameters$parasite == "falciparum"){
+      source_humans <- variables$state$get_index_of(c('S','A','U'))$and(bitten_humans)
+    } else if (parameters$parasite == "vivax"){
+      ## source_humans must include individuals with hypnozoites which may be impacted by prophylaxis/vaccination
+      hypnozoites_humans <- variables$hypnozoites$get_index_of(set = 0)$not(T)
+      source_humans <- bitten_humans$or(hypnozoites_humans)
+    }
+  
+  if(source_humans$size() > 0){
+   
+    source_vector <- source_humans$to_vector()
     
-  } else if (parameters$parasite == "vivax"){
-    ## p.v does not model blood immunity
-    b <- parameters$b
-  }
-
-  source_vector <- source_humans$to_vector()
-
-  # calculate prophylaxis
-  prophylaxis <- rep(0, length(source_vector))
-  drug <- variables$drug$get_values(source_vector)
-  medicated <- (drug > 0)
-  if (any(medicated)) {
-    drug <- drug[medicated]
-    drug_time <- variables$drug_time$get_values(source_vector[medicated])
-    prophylaxis[medicated] <- weibull_survival(
-      timestep - drug_time,
-      parameters$drug_prophylaxis_shape[drug],
-      parameters$drug_prophylaxis_scale[drug]
-    )
-  }
-
-  # calculate vaccine efficacy
-  vaccine_efficacy <- rep(0, length(source_vector))
-  vaccine_times <- variables$last_eff_pev_timestep$get_values(source_vector)
-  pev_profile <- variables$pev_profile$get_values(source_vector)
-  # get vector of individuals who have received their 3rd dose
-  vaccinated <- vaccine_times > -1
-  pev_profile <- pev_profile[vaccinated]
-  if (length(vaccinated) > 0) {
-    antibodies <- calculate_pev_antibodies(
-      timestep - vaccine_times[vaccinated],
-      exp(sample_pev_param(pev_profile, parameters$pev_profiles, 'cs')),
-      invlogit(sample_pev_param(pev_profile, parameters$pev_profiles, 'rho')),
-      exp(sample_pev_param(pev_profile, parameters$pev_profiles, 'ds')),
-      exp(sample_pev_param(pev_profile, parameters$pev_profiles, 'dl')),
-      parameters
-    )
-    vmax <- vnapply(parameters$pev_profiles, function(p) p$vmax)
-    beta <- vnapply(parameters$pev_profiles, function(p) p$beta)
-    alpha <- vnapply(parameters$pev_profiles, function(p) p$alpha)
-    vaccine_efficacy[vaccinated] <- calculate_pev_efficacy(
-      antibodies,
-      vmax[pev_profile],
-      beta[pev_profile],
-      alpha[pev_profile]
+    # calculate prophylaxis
+    prophylaxis <- rep(0, length(source_vector))
+    drug <- variables$drug$get_values(source_vector)
+    medicated <- (drug > 0)
+    if (any(medicated)) {
+      drug <- drug[medicated]
+      drug_time <- variables$drug_time$get_values(source_vector[medicated])
+      prophylaxis[medicated] <- weibull_survival(
+        timestep - drug_time,
+        parameters$drug_prophylaxis_shape[drug],
+        parameters$drug_prophylaxis_scale[drug]
+      )
+    }
+    
+    # calculate vaccine efficacy
+    vaccine_efficacy <- rep(0, length(source_vector))
+    vaccine_times <- variables$last_eff_pev_timestep$get_values(source_vector)
+    pev_profile <- variables$pev_profile$get_values(source_vector)
+    # get vector of individuals who have received their 3rd dose
+    vaccinated <- vaccine_times > -1
+    pev_profile <- pev_profile[vaccinated]
+    if (length(vaccinated) > 0) {
+      antibodies <- calculate_pev_antibodies(
+        timestep - vaccine_times[vaccinated],
+        exp(sample_pev_param(pev_profile, parameters$pev_profiles, 'cs')),
+        invlogit(sample_pev_param(pev_profile, parameters$pev_profiles, 'rho')),
+        exp(sample_pev_param(pev_profile, parameters$pev_profiles, 'ds')),
+        exp(sample_pev_param(pev_profile, parameters$pev_profiles, 'dl')),
+        parameters
+      )
+      vmax <- vnapply(parameters$pev_profiles, function(p) p$vmax)
+      beta <- vnapply(parameters$pev_profiles, function(p) p$beta)
+      alpha <- vnapply(parameters$pev_profiles, function(p) p$alpha)
+      vaccine_efficacy[vaccinated] <- calculate_pev_efficacy(
+        antibodies,
+        vmax[pev_profile],
+        beta[pev_profile],
+        alpha[pev_profile]
+      )
+    }
+    
+    # infection_rates <- rep(0, length = parameters$human_population)
+    if(parameters$parasite == "falciparum"){
+      
+      ## p.f models blood immunity
+      b <- blood_immunity(variables$ib$get_values(source_humans), parameters)
+      
+    } else if (parameters$parasite == "vivax"){
+      
+      ## p.v does not model blood immunity but must take into account multiple bites per person
+      b <- 1 - (1 - parameters$b)^n_bites_per_person[bitten_humans$to_vector()]
+      
+      ## get infection rates for all bitten or with hypnozoites
+      infection_rates <- rep(0, length = source_humans$size())
+      infection_rates[bitset_index(source_humans, bitten_humans)] <- prob_to_rate(b)
+      
+      # Add relapse rates for individuals with hypnozoites
+      if(hypnozoites_humans$size()>0){
+        relapse_rates <- variables$hypnozoites$get_values(hypnozoites_humans) * parameters$f
+        hyp_source_index <- bitset_index(source_humans, hypnozoites_humans)
+        infection_rates[hyp_source_index] <- infection_rates[hyp_source_index] + relapse_rates
+        # Get and store relative rates for bite/relapse competing hazards resolution
+        relative_rates <- relapse_rates/infection_rates[hyp_source_index]
+        infection_outcome$set_relative_rates(
+          hypnozoites_humans, 
+          relative_rates)
+      }
+      
+      prob <- rate_to_prob(infection_rates)
+    }
+    
+    prob <- b * (1 - prophylaxis) * (1 - vaccine_efficacy)
+    infection_rates <- prob_to_rate(prob)
+    
+    ## probability of incidence must be rendered at each timestep
+    incidence_probability_renderer(
+      variables$birth,
+      renderer,
+      source_humans,
+      prob,
+      "inc_",
+      parameters$incidence_min_ages,
+      parameters$incidence_max_ages,
+      timestep
     )
+    
+    ## capture infection rates to resolve in competing hazards
+    infection_outcome$set_rates(
+      source_humans, 
+      infection_rates)
   }
-
-  prob <- b * (1 - prophylaxis) * (1 - vaccine_efficacy)
-
-  ## probability of incidence must be rendered at each timestep
-  incidence_probability_renderer(
-    variables$birth,
-    renderer,
-    source_humans,
-    prob,
-    "inc_",
-    parameters$incidence_min_ages,
-    parameters$incidence_max_ages,
-    timestep
-  )
-  
-  ## capture infection rates to resolve in competing hazards
-  infection_outcome$set_rates(
-    source_humans,
-    prob_to_rate(prob))
 }
 
 #' @title Assigns infections to appropriate human states
@@ -145,6 +181,7 @@ calculate_infections <- function(
 #' @param renderer model render object
 #' @param parameters model parameters
 #' @param prob vector of population probabilities of infection
+#' @param relative_rates relative rates of hypnozoite relapse relative to total infection rate
 #' @noRd
 infection_outcome_process <- function(
     timestep,
@@ -152,7 +189,8 @@ infection_outcome_process <- function(
     variables,
     renderer,
     parameters,
-    prob){
+    prob,
+    relative_rates = NULL){
   
   if (infected_humans$size() > 0) {
     
@@ -215,6 +253,15 @@ infection_outcome_process <- function(
       
     } else if (parameters$parasite == "vivax"){
       
+      relapse_bite_infection_hazard_resolution(
+        infected_humans,
+        relative_rates,
+        variables,
+        parameters,
+        renderer,
+        timestep
+      )
+
       boost_immunity(
         variables$iaa,
         infected_humans,
@@ -227,9 +274,7 @@ infection_outcome_process <- function(
       lm_detectable <- calculate_lm_det_infections(
         variables,
         variables$state$get_index_of(c("S","U"))$and(infected_humans),
-        parameters,
-        renderer,
-        timestep
+        parameters
       )
       
       # Lm-detectable level infected S and U, and all A infections may receive clinical infections
@@ -270,31 +315,91 @@ infection_outcome_process <- function(
   }
 }
 
+#' @title Relapse/bite infection competing hazard resolution (p.v only)
+#' @description
+#' Resolves competing hazards of bite and hypnozoite relapse infections. 
+#' For bite infections we increase the batch number and factor in drug prophylaxis.
+#' 
+#' @param variables a list of all of the model variables
+#' @param infected_humans bitset of infected humans
+#' @param relative_rates relative rate of relapse infection
+#' @param variables model variables
+#' @param parameters model parameters
+#' @param renderer model renderer
+#' @param timestep current timestep
+#' @noRd
+relapse_bite_infection_hazard_resolution <- function(
+    infected_humans,
+    relative_rates,
+    variables,
+    parameters,
+    renderer,
+    timestep
+){
+  
+  if(variables$hypnozoites$get_index_of(0)$not(T)$and(infected_humans)$size()>0){
+    
+    hypnozoite_humans <- variables$hypnozoites$get_index_of(0)$not(T)
+    potential_relapse_index <- bitset_index(hypnozoite_humans, infected_humans)
+    relapse_infections <- bitset_at(
+      hypnozoite_humans$and(infected_humans),
+      bernoulli_multi_p(relative_rates[potential_relapse_index])
+    )
+    
+    renderer$render('n_relapses', relapse_infections$size(), timestep)
+    # render relapse infections by age
+    incidence_renderer(
+      variables$birth,
+      renderer,
+      relapse_infections,
+      'inc_relapse_',
+      parameters$incidence_relapse_rendering_min_ages,
+      parameters$incidence_relapse_rendering_max_ages,
+      timestep
+    )
+    
+    # get bite infections
+    bite_infections <- infected_humans$copy()$and(relapse_infections$not(inplace = F))
+    
+  } else {
+    bite_infections <- infected_humans
+  }
+  
+  ## all bitten humans with an infectious bite (incorporating prophylaxis) get a new batch of hypnozoites
+  if(bite_infections$size()>0){
+    
+    # make sure batches are capped
+    current_batches <- variables$hypnozoites$get_values(bite_infections)
+    new_batch_number <- pmin(current_batches + 1, parameters$kmax)
+    
+    variables$hypnozoites$queue_update(
+      new_batch_number,
+      bite_infections
+    )
+  }
+}
+
 #' @title Calculate light microscopy detectable infections (p.v only)
 #' @description
 #' Sample light microscopy detectable infections from all infections
 #' @param variables a list of all of the model variables
 #' @param infections bitset of infected humans
 #' @param parameters model parameters
-#' @param renderer model render
-#' @param timestep current timestep
 #' @noRd
 calculate_lm_det_infections <- function(
     variables,
     infections,
-    parameters,
-    renderer,
-    timestep
+    parameters
 ) {
   
   iaa <- variables$iaa$get_values(infections)
   iam <- variables$iam$get_values(infections)
-
+  
   philm <- anti_parasite_immunity(
     min = parameters$philm_min, max = parameters$philm_max, a50 = parameters$alm50,
     k = parameters$klm, iaa = iaa, iam = iam)
 
-  lm_det_infections <- bitset_at(infections, bernoulli_multi_p(philm))
+  bitset_at(infections, bernoulli_multi_p(philm))
 }
 
 #' @title Calculate clinical infections
@@ -603,7 +708,7 @@ schedule_infections <- function(
         'A',
         variables$infectivity,
         parameters$ca,
-        variables$recovery_rates,
+        variables$progression_rates,
         1/parameters$da,
         to_A
       )
diff --git a/R/mortality_processes.R b/R/mortality_processes.R
index b41aca32..9b530c56 100644
--- a/R/mortality_processes.R
+++ b/R/mortality_processes.R
@@ -112,6 +112,7 @@ reset_target <- function(variables, events, target, state, parameters, timestep)
     } else if (parameters$parasite == "vivax"){
       variables$last_boosted_iaa$queue_update(-1, target)
       variables$iaa$queue_update(0, target)
+      variables$hypnozoites$queue_update(0, target)
     }
 
     # treatment
diff --git a/R/parameters.R b/R/parameters.R
index 79ef7629..f113d908 100644
--- a/R/parameters.R
+++ b/R/parameters.R
@@ -189,9 +189,9 @@
 #' * delay_gam - lag from parasites to infectious gametocytes; default = 12.5
 #' * dem - extrinsic incubation period in mosquito population model; default = 10
 #' 
-#' hypnozoite parameters (p.v only):
+#' hypnozoite batch parameters (p.v only):
 #'
-#' * f - hypnozoite relapse rate; default = 0.02439024
+#' * f - hypnozoite batch relapse rate; default = 0.02439024
 #' * gammal - hypnozoite batch clearance rate; default = 0.002610966
 #' * init_hyp - initial hypnozoite batch number; default = 0
 #' * kmax - maximum number of hypnozoite batches for use in the equilibrium solution; default = 10
diff --git a/R/processes.R b/R/processes.R
index 3330f2b4..3acee329 100644
--- a/R/processes.R
+++ b/R/processes.R
@@ -48,26 +48,6 @@ create_processes <- function(
                                                         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)
-    )
-  }
-
   if(parameters$parasite == "falciparum"){
     processes <- c(
       processes,
@@ -81,14 +61,17 @@ create_processes <- function(
       immunity_process = create_exponential_decay_process(variables$iva,
                                                           parameters$rva),
       # Immunity to detectability
-      immunity_process = create_exponential_decay_process(variables$id, parameters$rid)
+      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)
+      immunity_process = create_exponential_decay_process(variables$iam, 
+                                                          parameters$rm),
+      immunity_process = create_exponential_decay_process(variables$iaa, 
+                                                          parameters$ra)
     )
   }
 
@@ -113,7 +96,8 @@ create_processes <- function(
     targeted_process = function(timestep, target){
       infection_outcome_process(timestep, target, 
                                 variables, renderer, parameters, 
-                                prob = rate_to_prob(infection_outcome$rates))
+                                prob = rate_to_prob(infection_outcome$rates),
+                                relative_rates = infection_outcome$relative_rates)
     },
     size = parameters$human_population
   )
@@ -219,7 +203,23 @@ create_processes <- function(
   if(parameters$parasite == "falciparum"){
     imm_var_names <- c(imm_var_names, 'ib', 'iva', 'ivm', 'id')
   } else if (parameters$parasite == "vivax"){
-    imm_var_names <- c(imm_var_names, 'iaa', 'iam')
+    imm_var_names <- c(imm_var_names, 'iaa', 'iam', 'hypnozoites')
+
+    ## hypnozoite infection prevalence rendering
+    processes <- c(
+      processes,
+      hypnozoite_prevalence_renderer = create_n_with_hypnozoites_renderer_process(
+        renderer,
+        variables$hypnozoites,
+        parameters
+      ),
+      hypnozoite_prevalence_age_group_renderer = create_n_with_hypnozoites_age_renderer_process(
+        variables$hypnozoites,
+        variables$birth,
+        parameters,
+        renderer
+      )
+    )
   }
 
   processes <- c(
@@ -396,3 +396,21 @@ create_lagged_eir <- function(variables, solvers, parameters) {
     }
   )
 }
+#' @title Hypnozoite decay function
+#' @description
+#' calulates the number of individuals in whom a batch decay occurs
+#'
+#' @param variable the hypnozoite variable to update
+#' @param rate the hypnozoite decay rate
+#' @noRd
+create_hypnozoite_batch_decay_process <- function(hypnozoites, gammal){
+  function(timestep){
+    to_decay <- bernoulli_multi_p(p = rate_to_prob(hypnozoites$get_values() * gammal))
+    if(length(to_decay) > 0){
+      hypnozoites$queue_update(
+        hypnozoites$get_values(to_decay) - 1,
+        to_decay
+      )
+    }
+  }
+}
\ No newline at end of file
diff --git a/R/render.R b/R/render.R
index 34f4bd99..1357decf 100644
--- a/R/render.R
+++ b/R/render.R
@@ -266,6 +266,11 @@ populate_incidence_rendering_columns <- function(renderer, parameters){
     renderer$set_default('n_early_treatment_failure', 0)
     renderer$set_default('n_slow_parasite_clearance', 0)
   }
+
+  # relapses only render for the vivax model
+  if(parameters$parasite == "vivax"){
+    renderer$set_default('n_relapses', 0)
+  }
   
   if(length(parameters$incidence_rendering_min_ages)>0){
     render_initial_incidence(renderer,
@@ -303,3 +308,38 @@ populate_metapopulation_incidence_rendering_columns <- function(renderer, parame
     populate_incidence_rendering_columns(renderer[[i]], parameters[[i]])
   }
 }
+
+create_n_with_hypnozoites_renderer_process <- function(
+    renderer,
+    hypnozoites,
+    parameters
+) {
+  function(timestep) {
+    renderer$render(
+      paste0("n_with_hypnozoites"),
+      hypnozoites$size() - hypnozoites$get_size_of(0),
+      timestep
+    )
+  }
+}
+
+create_n_with_hypnozoites_age_renderer_process <- function(
+    hypnozoites,
+    birth,
+    parameters,
+    renderer
+) {
+  function(timestep) {
+    for (i in seq_along(parameters$n_with_hypnozoites_rendering_min_ages)) {
+      lower <- parameters$n_with_hypnozoites_rendering_min_ages[[i]]
+      upper <- parameters$n_with_hypnozoites_rendering_max_ages[[i]]
+      in_age <- in_age_range(birth, timestep, lower, upper)
+      renderer$render(paste0('n_', lower, '_', upper), in_age$size(), timestep)
+      renderer$render(
+        paste0('n_with_hypnozoites_', lower, '_', upper),
+        sum(hypnozoites$get_values(index = in_age)!=0),
+        timestep
+      )
+    }
+  }
+}
diff --git a/R/variables.R b/R/variables.R
index 024b695d..b6bd1f80 100644
--- a/R/variables.R
+++ b/R/variables.R
@@ -18,6 +18,7 @@
 #' * ICA  - Acquired immunity to clinical disease
 #' * IVA  - Acquired immunity to severe disease (p.f only)
 #' * ID - Acquired immunity to detectability (p.f only)
+#' * hypnozoites - Hypnozoite batch number (p.v 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
@@ -98,7 +99,10 @@ 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)
-  
+  if(parameters$parasite == "vivax"){
+    hypnozoites <- individual::IntegerVariable$new(rep(parameters$init_hyp, parameters$human_population))
+  }
+
   # Maternal immunity to clinical disease
   icm <- individual::DoubleVariable$new(
     initial_immunity(
@@ -299,7 +303,8 @@ create_variables <- function(parameters) {
     variables <- c(variables,
                    last_boosted_iaa = last_boosted_iaa,
                    iaa = iaa,
-                   iam = iam
+                   iam = iam,
+                   hypnozoites = hypnozoites
     )
   }
   
diff --git a/man/get_parameters.Rd b/man/get_parameters.Rd
index 76e2db6a..76ce99a8 100644
--- a/man/get_parameters.Rd
+++ b/man/get_parameters.Rd
@@ -209,9 +209,9 @@ parasite incubation periods:
 \item dem - extrinsic incubation period in mosquito population model; default = 10
 }
 
-hypnozoite parameters (p.v only):
+hypnozoite batch parameters (p.v only):
 \itemize{
-\item f - hypnozoite relapse rate; default = 0.02439024
+\item f - hypnozoite batch relapse rate; default = 0.02439024
 \item gammal - hypnozoite batch clearance rate; default = 0.002610966
 \item init_hyp - initial hypnozoite batch number; default = 0
 \item kmax - maximum number of hypnozoite batches for use in the equilibrium solution; default = 10
diff --git a/tests/testthat/test-biting-integration.R b/tests/testthat/test-biting-integration.R
index fbfd1fd4..9ab38a76 100644
--- a/tests/testthat/test-biting-integration.R
+++ b/tests/testthat/test-biting-integration.R
@@ -33,8 +33,9 @@ test_that('biting_process integrates mosquito effects and human infection', {
     infection_outcome=infection_outcome
   )
   
-  bitten <- individual::Bitset$new(parameters$human_population)
-  bites_mock <- mockery::mock(bitten)
+  bitten <- list(bitten_humans = individual::Bitset$new(parameters$human_population),
+                 n_bites_per_person = numeric(0))
+  bites_mock <- mockery::mock(bitten, cycle = T)
   infection_mock <- mockery::mock()
 
   mockery::stub(biting_process, 'simulate_bites', bites_mock)
@@ -63,7 +64,8 @@ test_that('biting_process integrates mosquito effects and human infection', {
     1,
     variables,
     events,
-    bitten,
+    bitten$bitten,
+    bitten$n_bites_per_person,
     age,
     parameters,
     timestep,
@@ -125,7 +127,7 @@ test_that('simulate_bites integrates eir calculation and mosquito side effects',
     lagged_eir
   )
 
-  expect_equal(bitten$to_vector(), c(2, 3))
+  expect_equal(bitten$bitten_humans$to_vector(), c(2, 3))
 
   f <- parameters$blood_meal_rates[[1]]
 
diff --git a/tests/testthat/test-infection-integration.R b/tests/testthat/test-infection-integration.R
index ad83c35f..c29420c9 100644
--- a/tests/testthat/test-infection-integration.R
+++ b/tests/testthat/test-infection-integration.R
@@ -17,6 +17,7 @@ test_that('simulate_infection integrates different types of infection and schedu
   )
   
   bitten <- individual::Bitset$new(population)$insert(c(1, 3, 5, 7))
+  n_bites_per_person <- numeric(0)
   boost_immunity_mock <- mockery::mock()
   infected <- individual::Bitset$new(population)$insert(c(1, 3, 5))
   infection_mock <- mockery::mock(infected)
@@ -35,6 +36,7 @@ test_that('simulate_infection integrates different types of infection and schedu
     variables,
     events,
     bitten,
+    n_bites_per_person,
     age,
     parameters,
     timestep,
@@ -57,6 +59,7 @@ test_that('simulate_infection integrates different types of infection and schedu
     1,
     variables,
     bitten,
+    n_bites_per_person,
     parameters,
     renderer,
     timestep,
@@ -64,7 +67,6 @@ test_that('simulate_infection integrates different types of infection and schedu
   )
 })
 
-
 test_that('simulate_infection integrates different types of infection and scheduling', {
   population <- 8
   timestep <- 5
@@ -156,6 +158,8 @@ test_that('simulate_infection integrates different types of infection and schedu
     renderer
   )
   
+  mockery::mock_args(schedule_mock)
+  
   mockery::expect_args(
     schedule_mock,
     1,
@@ -220,6 +224,7 @@ test_that('calculate_infections works various combinations of drug and vaccinati
   )
   
   bitten_humans <- individual::Bitset$new(4)$insert(c(1, 2, 3, 4))
+  n_bites_per_person <- numeric(0)
   
   infection_outcome <- CompetingOutcome$new(
     targeted_process = function(timestep, target){
@@ -231,6 +236,7 @@ test_that('calculate_infections works various combinations of drug and vaccinati
   infections <- calculate_infections(
     variables,
     bitten_humans, 
+    n_bites_per_person,
     parameters,
     mock_render(timestep),
     timestep,
@@ -699,10 +705,10 @@ test_that('prophylaxis is considered for medicated humans', {
   )
   
   bitten <- individual::Bitset$new(4)$insert(seq(4))
+  n_bites_per_person <- numeric(0)
   m <- mockery::mock(seq(3))
   mockery::stub(calculate_infections, 'bernoulli_multi_p', m)
   
-  
   infection_outcome <- CompetingOutcome$new(
     targeted_process = function(timestep, target){
       infection_outcome_process(timestep, target, variables, renderer, parameters)
@@ -713,6 +719,7 @@ test_that('prophylaxis is considered for medicated humans', {
   infection_rates <- calculate_infections(
     variables,
     bitten,
+    n_bites_per_person,
     parameters,
     mock_render(timestep),
     timestep,
diff --git a/tests/testthat/test-pev.R b/tests/testthat/test-pev.R
index f3c3d793..ce41ac20 100644
--- a/tests/testthat/test-pev.R
+++ b/tests/testthat/test-pev.R
@@ -174,6 +174,7 @@ test_that('Infection considers pev efficacy', {
   infection_rates <- calculate_infections(
     variables = variables,
     bitten_humans = individual::Bitset$new(4)$insert(seq(4)),
+    n_bites_per_person = numeric(0),
     parameters = parameters,
     renderer = mock_render(timestep),
     timestep = timestep,
diff --git a/tests/testthat/test-vivax.R b/tests/testthat/test-vivax.R
index 8fc77956..4f762728 100644
--- a/tests/testthat/test-vivax.R
+++ b/tests/testthat/test-vivax.R
@@ -255,3 +255,63 @@ test_that('vivax schedule_infections correctly schedules new infections', {
     c(1)
   )
 })
+
+test_that('relapses are recognised with division between bite infections and relapses', {
+  timestep <- 50
+  parameters <- get_parameters(parasite = "vivax", overrides = list(human_population = 4))
+  
+  variables <- list(
+    state = individual::CategoricalVariable$new(
+      c('D', 'S', 'A', 'U', 'Tr'),
+      c('D', 'S', 'A', 'U')
+    ),
+    infectivity = individual::DoubleVariable$new(rep(0, 4)),
+    progression_rates = individual::DoubleVariable$new(rep(0, 4)),
+    drug = individual::DoubleVariable$new(rep(0, 4)),
+    drug_time = individual::DoubleVariable$new(rep(-1, 4)),
+    iaa = individual::DoubleVariable$new(rep(0, 4)),
+    iam = individual::DoubleVariable$new(rep(0, 4)),
+    ica = individual::DoubleVariable$new(rep(0, 4)),
+    icm = individual::DoubleVariable$new(rep(0, 4)),
+    last_boosted_iaa = individual::DoubleVariable$new(rep(-1, 4)),
+    last_boosted_ica = individual::DoubleVariable$new(rep(-1, 4)),
+    last_eff_pev_timestep = individual::DoubleVariable$new(rep(-1, 4)),
+    pev_profile = individual::IntegerVariable$new(rep(-1, 4)),
+    hypnozoites = individual::IntegerVariable$new(c(0, 1, 2, 3))
+  )
+  
+  bernoulli_mock <- mockery::mock(c(1, 3), 2, cycle = TRUE)
+  calc_mock <- mockery::mock(individual::Bitset$new(4)$insert(2))
+  mockery::stub(infection_outcome_process, 'bernoulli_multi_p', bernoulli_mock, depth = 2)
+  mockery::stub(infection_outcome_process, 'calculate_clinical_infections', calc_mock)
+  
+  renderer <- mock_render(1)
+  infected_humans <- individual::Bitset$new(4)$insert(c(1, 2, 3, 4))
+  
+  infection_outcome_process(
+    timestep = timestep,
+    infected_humans,
+    variables,
+    renderer,
+    parameters,
+    prob = c(rep(0.5, 4)),
+    relative_rate = c(rep(0.5, 3))
+  )
+
+  mockery::expect_args(
+    renderer$render_mock(),
+    1,
+    'n_infections',
+    4,
+    timestep
+  )
+
+  mockery::expect_args(
+    renderer$render_mock(),
+    2,
+    'n_relapses',
+    2,
+    timestep
+  )
+  
+})
diff --git a/vignettes/Plasmodium_vivax.Rmd b/vignettes/Plasmodium_vivax.Rmd
index 87c56671..67c0227e 100644
--- a/vignettes/Plasmodium_vivax.Rmd
+++ b/vignettes/Plasmodium_vivax.Rmd
@@ -55,6 +55,8 @@ The *P. vivax* model follows a similar structure to the *P. falciparum* model, a
 
 ### New Infections
 
+The rate of infection for an individual who has been bitten increases with the number of bites they have received in the *P. vivax* model. In contrast, the rate of infection in the *P. falciparum* model is independent of the number of bites an individual has received.
+
 Newly infected individuals in the  *P. falciparum* model can move into either the clinically diseased or asymptomatic infection compartment. In addition to these compartments, the *P. vivax* model allows infection to the PCR-detectable compartment (U), where the assignment of light-miscroscopy detectable infections and PCR-detectable infections is mediated by anti-parasite immunity.
 
 ### Immunity
@@ -69,6 +71,10 @@ Anti-parasite immunity has effects in two places. The first is in the separation
 
 While the *P. falciparum* model calculates the onward infectivity of asymptomatic infections (`ca`) using the age and detectability immunity of each individual, the *P. vivax* model uses a constant onwards infectivity for LM-detectable infections (`ca = 0.1`).
 
+### Hypnozoites
+
+The *P. vivax* model tracks the number of hypnozoite batches of each individual which contribute to the rate of new infections. Acquisition of new batches comes from bite infections and are capped (where `kmax = 10` by default). All individuals can acquire new hypnozoite batches via bites, even when these do not result in new blood-stage infection (as in the clinically diseased and treated). Hypnozoite batches are lost probabilistically where an individual's number of batches determines the loss rate of a single batch.
+
 ### Key Model References
 
 White, M. T., Walker, P., Karl, S., Hetzel, M. W., Freeman, T., Waltmann, A., Laman, M., Robinson, L. J., Ghani, A., & Mueller, I. (2018). Mathematical modelling of the impact of expanding levels of malaria control interventions on Plasmodium vivax. Nature Communications, 9(1), 1–10. <https://doi.org/10.1038/s41467-018-05860-8>