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..ea7ffd97 100644
--- a/R/competing_hazards.R
+++ b/R/competing_hazards.R
@@ -16,22 +16,25 @@ CompetingOutcome <- R6::R6Class(
 
       self$target <- individual::Bitset$new(size)
       self$rates <- NULL
+      self$args <- NULL
     },
-    set_rates = function(target, rates){
+    set_rates = function(target, rates, ...){
       stopifnot(target$size() == length(rates))
-
       self$target$copy_from(target)
       self$rates <- rates
+      self$args <- list(...)
     },
     execute = function(t, target){
-      private$targeted_process(t, target)
+      do.call(private$targeted_process, c(t, list(target), self$args))
     },
     reset = function() {
       self$target$clear()
       self$rates <- NULL
+      self$args <- NULL
     },
     target = NULL,
-    rates = NULL
+    rates = NULL,
+    args = NULL
   )
 )
 
diff --git a/R/human_infection.R b/R/human_infection.R
index 36c07897..36aa7b2a 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,14 +16,17 @@ simulate_infection <- function(
     variables,
     events,
     bitten_humans,
+    n_bites_per_person,
     age,
     parameters,
     timestep,
     renderer,
     infection_outcome
 ) {
+  
   if (bitten_humans$size() > 0) {
     if(parameters$parasite == "falciparum"){
+      
       boost_immunity(
         variables$ib,
         bitten_humans,
@@ -30,51 +34,193 @@ simulate_infection <- function(
         timestep,
         parameters$ub
       )
+      
+      # Calculate Infected
+      calculate_falciparum_infections(
+        variables,
+        bitten_humans,
+        parameters,
+        renderer,
+        timestep,
+        infection_outcome
+      )
+      
+    } else if(parameters$parasite == "vivax"){
+      
+      # Calculate Infected
+      calculate_vivax_infections(
+        variables,
+        bitten_humans,
+        n_bites_per_person,
+        parameters,
+        renderer,
+        timestep,
+        infection_outcome
+      )
     }
   }
+  
+}
 
-  # Calculate Infected
-  calculate_infections(
+#' @title Calculate overall falciparum infections for bitten humans
+#' @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 parameters model parameters
+#' @param renderer model render object
+#' @param timestep current timestep
+#' @noRd
+calculate_falciparum_infections <- function(
     variables,
     bitten_humans,
     parameters,
     renderer,
     timestep,
     infection_outcome
-  )
+) {
+  
+  source_humans <- variables$state$get_index_of(c('S','A','U'))$and(bitten_humans)
+  
+  if(source_humans$size() > 0){
+    
+    source_vector <- source_humans$to_vector()
+    
+    # calculate prophylaxis
+    prophylaxis <- treatment_prophylaxis_efficacy(
+      source_vector,
+      variables,
+      parameters,
+      timestep
+    )
+    
+    # calculate vaccine efficacy
+    vaccine_efficacy <- pev_efficacy(
+      source_vector,
+      variables,
+      parameters,
+      timestep)
+    
+    ## p.f models blood immunity
+    prob <- blood_immunity(variables$ib$get_values(source_humans), parameters)
+    prob <- prob * (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)
+  }
 }
 
-#' @title Calculate overall infections for bitten humans
+#' @title Calculate overall vivax infections for bitten humans
 #' @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
 #' @noRd
-calculate_infections <- function(
+calculate_vivax_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)
+  
+  ## 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$copy()$or(hypnozoites_humans)
+  
+  if(source_humans$size() > 0){
+    
+    source_vector <- source_humans$to_vector()
+    
+    # calculate prophylaxis
+    prophylaxis <- treatment_prophylaxis_efficacy(
+      source_vector,
+      variables,
+      parameters,
+      timestep
+    )
+    
+    # calculate vaccine efficacy
+    vaccine_efficacy <- pev_efficacy(
+      source_vector,
+      variables,
+      parameters,
+      timestep)
+    
+    ## 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
+    relative_rates <- NULL
+    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]
+    }
+    
+    prob <- rate_to_prob(infection_rates) * (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
+    )
     
-  } else if (parameters$parasite == "vivax"){
-    ## p.v does not model blood immunity
-    b <- parameters$b
+    ## capture infection rates to resolve in competing hazards
+    infection_outcome$set_rates(
+      source_humans,
+      infection_rates,
+      relative_rates = relative_rates
+    )
   }
+}
 
-  source_vector <- source_humans$to_vector()
-
-  # calculate prophylaxis
+#' @title Calculate protection from infection due to drug prophylaxis
+#' @description This function calculates the probability that drug prophylaxis will 
+#' protect each individual from infection
+#' @param source_vector a vector of individuals with prospective new infections
+#' @param variables a list of all of the model variables
+#' @param parameters model parameters
+#' @param timestep current timestep
+#' @noRd
+treatment_prophylaxis_efficacy <- function(
+    source_vector,
+    variables,
+    parameters,
+    timestep
+){
   prophylaxis <- rep(0, length(source_vector))
   drug <- variables$drug$get_values(source_vector)
   medicated <- (drug > 0)
@@ -87,8 +233,23 @@ calculate_infections <- function(
       parameters$drug_prophylaxis_scale[drug]
     )
   }
+  prophylaxis
+}
 
-  # calculate vaccine efficacy
+#' @title Calculate protection from infection due to vaccination
+#' @description This function calculates the probability that vaccination will 
+#' protect each individual from infection
+#' @param source_vector a vector of individuals with prospective new infections
+#' @param variables a list of all of the model variables
+#' @param parameters model parameters
+#' @param timestep current timestep
+#' @noRd
+pev_efficacy <- function(
+    source_vector,
+    variables,
+    parameters,
+    timestep
+){
   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)
@@ -114,28 +275,97 @@ calculate_infections <- function(
       alpha[pev_profile]
     )
   }
+  vaccine_efficacy
+}
 
-  prob <- b * (1 - prophylaxis) * (1 - vaccine_efficacy)
-
-  ## probability of incidence must be rendered at each timestep
-  incidence_probability_renderer(
-    variables$birth,
+#' @title Assigns falciparum infections to appropriate human states
+#' @description
+#' Updates human states and variables to represent asymptomatic/clinical/severe
+#' and treated malaria; and resulting boosts in immunity
+#' @param timestep current timestep
+#' @param infected_humans bitset of infected humans
+#' @param variables a list of all of the model variables
+#' @param renderer model render object
+#' @param parameters model parameters
+#' @param prob vector of population probabilities of infection
+#' @noRd
+falciparum_infection_outcome_process <- function(
+    timestep,
+    infected_humans,
+    variables,
     renderer,
-    source_humans,
-    prob,
-    "inc_",
-    parameters$incidence_min_ages,
-    parameters$incidence_max_ages,
-    timestep
-  )
+    parameters){
   
-  ## capture infection rates to resolve in competing hazards
-  infection_outcome$set_rates(
-    source_humans,
-    prob_to_rate(prob))
+  if (infected_humans$size() > 0) {
+    
+    renderer$render('n_infections', infected_humans$size(), timestep)
+    incidence_renderer(
+      variables$birth,
+      renderer,
+      infected_humans,
+      'inc_',
+      parameters$incidence_rendering_min_ages,
+      parameters$incidence_rendering_max_ages,
+      timestep
+    )
+    
+    boost_immunity(
+      variables$ica,
+      infected_humans,
+      variables$last_boosted_ica,
+      timestep,
+      parameters$uc
+    )
+    
+    boost_immunity(
+      variables$id,
+      infected_humans,
+      variables$last_boosted_id,
+      timestep,
+      parameters$ud
+    )
+    
+    clinical <- calculate_clinical_infections(
+      variables,
+      infected_humans,
+      parameters,
+      renderer,
+      timestep
+    )
+    
+    treated <- calculate_treated(
+      variables,
+      clinical,
+      parameters,
+      timestep,
+      renderer
+    )
+    
+    update_severe_disease(
+      timestep,
+      infected_humans,
+      variables,
+      parameters,
+      renderer
+    )
+    
+    ## The treated and infected_humans bitsets are re-written so be cautious!
+    to_D <- treated$not(FALSE)$and(clinical)
+    to_A <- infected_humans$and(clinical$not(FALSE))
+    to_U <- NULL
+    
+    schedule_infections(
+      parameters,
+      variables,
+      timestep,
+      to_D,
+      to_A,
+      to_U
+    )
+  }
 }
 
-#' @title Assigns infections to appropriate human states
+#' @title Assigns vivax infections to appropriate human states
 #' @description
 #' Updates human states and variables to represent asymptomatic/clinical/severe
 #' and treated malaria; and resulting boosts in immunity
@@ -144,15 +374,15 @@ calculate_infections <- function(
 #' @param variables a list of all of the model variables
 #' @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(
+vivax_infection_outcome_process <- function(
     timestep,
     infected_humans,
     variables,
     renderer,
     parameters,
-    prob){
+    relative_rates){
   
   if (infected_humans$size() > 0) {
     
@@ -167,6 +397,14 @@ infection_outcome_process <- function(
       timestep
     )
     
+    boost_immunity(
+      variables$iaa,
+      infected_humans,
+      variables$last_boosted_iaa,
+      timestep,
+      parameters$ua
+    )
+    
     boost_immunity(
       variables$ica,
       infected_humans,
@@ -174,90 +412,48 @@ infection_outcome_process <- function(
       timestep,
       parameters$uc
     )
-
-    if(parameters$parasite == "falciparum"){
-      boost_immunity(
-        variables$id,
-        infected_humans,
-        variables$last_boosted_id,
-        timestep,
-        parameters$ud
-      )
-      
-      clinical <- calculate_clinical_infections(
-        variables,
-        infected_humans,
-        parameters,
-        renderer,
-        timestep
-      )
-      
-      treated <- calculate_treated(
-        variables,
-        clinical,
-        parameters,
-        timestep,
-        renderer
-      )
-      
-      update_severe_disease(
-        timestep,
-        infected_humans,
-        variables,
-        parameters,
-        renderer
-      )
-      
-      ## The treated and infected_humans bitsets are re-written so be cautious!
-      to_D <- treated$not(FALSE)$and(clinical)
-      to_A <- infected_humans$and(clinical$not(FALSE))
-      to_U <- NULL
-      
-    } else if (parameters$parasite == "vivax"){
-      
-      boost_immunity(
-        variables$iaa,
-        infected_humans,
-        variables$last_boosted_iaa,
-        timestep,
-        parameters$ua
-      )
-      
-      ## Only S and U infections are considered in generating lm-det infections
-      lm_detectable <- calculate_lm_det_infections(
-        variables,
-        variables$state$get_index_of(c("S","U"))$and(infected_humans),
-        parameters,
-        renderer,
-        timestep
-      )
-      
-      # Lm-detectable level infected S and U, and all A infections may receive clinical infections
-      # There is a different calculation to generate clinical infections, based on current infection level
-      # LM infections must only pass through the clinical calculation, therefore all "A" infections are included
-      # "S" and "U" infections must pass through the lm-detectable calculation prior to and in addition to the clinical
-      # calculation. We therefore consider all "A" infections and only the "S" and "U" infections that are now lm-detectable.
-      clinical <- calculate_clinical_infections(
-        variables,
-        variables$state$get_index_of("A")$and(infected_humans)$or(lm_detectable),
-        parameters,
-        renderer,
-        timestep
-      )
-      
-      treated <- calculate_treated(
-        variables,
-        clinical,
-        parameters,
-        timestep,
-        renderer
-      )
-      
-      ## The infected_humans,lm_detectable and clinical bitsets are re-written so be cautious!
-      to_U <- infected_humans$and(lm_detectable$not(F))$and(variables$state$get_index_of(c("S")))
-      to_A <- lm_detectable$and(clinical$not(F))
-      to_D <- clinical$and(treated$not(F))
-    }
+    
+    relapse_bite_infection_hazard_resolution(
+      infected_humans,
+      relative_rates,
+      variables,
+      parameters,
+      renderer,
+      timestep
+    )
+    
+    ## Only S and U infections are considered in generating lm-det infections
+    lm_detectable <- calculate_lm_det_infections(
+      variables,
+      variables$state$get_index_of(c("S","U"))$and(infected_humans),
+      parameters
+    )
+    
+    # Lm-detectable level infected S and U, and all A infections may receive clinical infections
+    # There is a different calculation to generate clinical infections, based on current infection level
+    # LM infections must only pass through the clinical calculation, therefore all "A" infections are included
+    # "S" and "U" infections must pass through the lm-detectable calculation prior to and in addition to the clinical
+    # calculation. We therefore consider all "A" infections and only the "S" and "U" infections that are now lm-detectable.
+    clinical <- calculate_clinical_infections(
+      variables,
+      variables$state$get_index_of("A")$and(infected_humans)$or(lm_detectable),
+      parameters,
+      renderer,
+      timestep
+    )
+    
+    treated <- calculate_treated(
+      variables,
+      clinical,
+      parameters,
+      timestep,
+      renderer
+    )
+    
+    ## The infected_humans,lm_detectable and clinical bitsets are re-written so be cautious!
+    to_U <- infected_humans$and(lm_detectable$not(F))$and(variables$state$get_index_of(c("S")))
+    to_A <- lm_detectable$and(clinical$not(F))
+    to_D <- clinical$and(treated$not(F))
     
     schedule_infections(
       parameters,
@@ -270,31 +466,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
@@ -574,7 +830,7 @@ schedule_infections <- function(
     to_A,
     to_U
 ) {
-
+  
   if(to_D$size() > 0) {
     update_infection(
       variables$state,
@@ -586,7 +842,7 @@ schedule_infections <- function(
       to_D
     )
   }
-
+  
   if(to_A$size() > 0) {
     if(parameters$parasite == "falciparum"){
       # p.f has immunity-determined asymptomatic infectivity
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 b67dc559..3048ba6e 100644
--- a/R/processes.R
+++ b/R/processes.R
@@ -47,7 +47,7 @@ create_processes <- function(
     immunity_process = create_exponential_decay_process(variables$ica,
                                                         parameters$rc)
   )
-  
+
   if(parameters$parasite == "falciparum"){
     processes <- c(
       processes,
@@ -71,7 +71,9 @@ create_processes <- function(
       immunity_process = create_exponential_decay_process(variables$iam, 
                                                           parameters$rm),
       immunity_process = create_exponential_decay_process(variables$iaa, 
-                                                          parameters$ra)
+                                                          parameters$ra),
+      hypnozoite_process = create_hypnozoite_batch_decay_process(variables$hypnozoites,
+                                                                 parameters$gammal)
     )
   }
 
@@ -92,14 +94,27 @@ create_processes <- function(
   # Competing Hazard Outcomes (infections and disease progression)
   # =====================================================
   
-  infection_outcome <- CompetingOutcome$new(
-    targeted_process = function(timestep, target){
-      infection_outcome_process(timestep, target, 
-                                variables, renderer, parameters, 
-                                prob = rate_to_prob(infection_outcome$rates))
-    },
-    size = parameters$human_population
-  )
+  if(parameters$parasite == "falciparum"){
+    infection_outcome <- CompetingOutcome$new(
+      targeted_process = function(timestep, target){
+        falciparum_infection_outcome_process(timestep, target, 
+                                  variables, renderer, parameters
+        )
+      },
+      size = parameters$human_population
+    )
+    
+  } else if (parameters$parasite == "vivax"){
+    infection_outcome <- CompetingOutcome$new(
+      targeted_process = function(timestep, target, relative_rates){
+        vivax_infection_outcome_process(timestep, target, 
+                                  variables, renderer, parameters, 
+                                  relative_rates
+        )
+      },
+      size = parameters$human_population
+    )
+  }
   
   progression_outcome <- CompetingOutcome$new(
     targeted_process = function(timestep, target){
@@ -202,7 +217,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(
@@ -379,3 +410,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-competing-hazards.R b/tests/testthat/test-competing-hazards.R
index f7ada3fc..21ddefcb 100644
--- a/tests/testthat/test-competing-hazards.R
+++ b/tests/testthat/test-competing-hazards.R
@@ -26,13 +26,15 @@ test_that("hazard resolves two disjoint outcomes", {
   
   outcome_1$set_rates(population, c(10, 0, 10, 0))
   outcome_2$set_rates(population, c(0, 10, 0, 10))
-  
+
   hazard$resolve(0)
-   
+
   mockery::expect_args(outcome_1_process, 1, 0,
                        individual::Bitset$new(size)$insert(c(1, 3)))
+
   mockery::expect_args(outcome_2_process, 1, 0,
                        individual::Bitset$new(size)$insert(c(2, 4)))
+  
 })
 
 test_that("hazard resolves two competing outcomes", {
@@ -130,8 +132,7 @@ test_that("outcomes can define a partial set of rates", {
   mockery::expect_args(outcome_1_process, 1, 0,
                        individual::Bitset$new(size)$insert(c(1, 3)))
   mockery::expect_args(outcome_2_process, 1, 0,
-                       individual::Bitset$new(size)$insert(c(4))
-  )
+                       individual::Bitset$new(size)$insert(c(4)))
 })
 
 test_that("hazard resolves three competing outcomes", {
@@ -173,4 +174,3 @@ test_that("hazard resolves three competing outcomes", {
   mockery::expect_args(outcome_3_process, 1, 0,
                        individual::Bitset$new(size)$insert(c(3, 4)))
 })
-
diff --git a/tests/testthat/test-infection-integration.R b/tests/testthat/test-infection-integration.R
index ad83c35f..a4bcd1ab 100644
--- a/tests/testthat/test-infection-integration.R
+++ b/tests/testthat/test-infection-integration.R
@@ -20,21 +20,23 @@ test_that('simulate_infection integrates different types of infection and schedu
   boost_immunity_mock <- mockery::mock()
   infected <- individual::Bitset$new(population)$insert(c(1, 3, 5))
   infection_mock <- mockery::mock(infected)
+  n_bites_per_person <- NULL
   
   infection_outcome <- CompetingOutcome$new(
-    targeted_process = function(timestep, target){
-      infection_outcome_process(timestep, target, variables, renderer, parameters)
+    targeted_process = function(timestep, target, args){
+      falciparum_infection_outcome_process(timestep, target, variables, renderer, parameters)
     },
     size = parameters$human_population
   )
   
   mockery::stub(simulate_infection, 'boost_immunity', boost_immunity_mock)
-  mockery::stub(simulate_infection, 'calculate_infections', infection_mock)
+  mockery::stub(simulate_infection, 'calculate_falciparum_infections', infection_mock)
   
   simulate_infection(
     variables,
     events,
     bitten,
+    n_bites_per_person,
     age,
     parameters,
     timestep,
@@ -64,7 +66,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
@@ -83,7 +84,6 @@ test_that('simulate_infection integrates different types of infection and schedu
     # state = list(get_index_of = mockery::mock(asymptomatics, cycle = T))
     state = individual::CategoricalVariable$new(categories = c("S","A","U","D","Tr"), initial_values = rep("S", population))#list(get_index_of = mockery::mock(asymptomatics, cycle = T))
   )
-  prob <- rep(0.5,population)
   
   source_humans <- individual::Bitset$new(population)$insert(c(1, 2, 3, 5))
   infected <- individual::Bitset$new(population)$insert(c(1, 3, 5))
@@ -100,20 +100,19 @@ test_that('simulate_infection integrates different types of infection and schedu
   to_A <- infected$and(clinical$not(FALSE))
   to_U <- NULL
 
-  mockery::stub(infection_outcome_process, 'incidence_renderer', mockery::mock())
-  mockery::stub(infection_outcome_process, 'boost_immunity', boost_immunity_mock)
-  mockery::stub(infection_outcome_process, 'calculate_clinical_infections', clinical_infection_mock)
-  mockery::stub(infection_outcome_process, 'update_severe_disease', severe_infection_mock)
-  mockery::stub(infection_outcome_process, 'calculate_treated', treated_mock)
-  mockery::stub(infection_outcome_process, 'schedule_infections', schedule_mock)
-  
-  infection_outcome_process(
+  mockery::stub(falciparum_infection_outcome_process, 'incidence_renderer', mockery::mock())
+  mockery::stub(falciparum_infection_outcome_process, 'boost_immunity', boost_immunity_mock)
+  mockery::stub(falciparum_infection_outcome_process, 'calculate_clinical_infections', clinical_infection_mock)
+  mockery::stub(falciparum_infection_outcome_process, 'update_severe_disease', severe_infection_mock)
+  mockery::stub(falciparum_infection_outcome_process, 'calculate_treated', treated_mock)
+  mockery::stub(falciparum_infection_outcome_process, 'schedule_infections', schedule_mock)
+  
+  falciparum_infection_outcome_process(
     timestep,
     infected,
     variables,
     renderer,
-    parameters,
-    prob)
+    parameters)
 
   mockery::expect_args(
     boost_immunity_mock,
@@ -156,6 +155,8 @@ test_that('simulate_infection integrates different types of infection and schedu
     renderer
   )
   
+  mockery::mock_args(schedule_mock)
+  
   mockery::expect_args(
     schedule_mock,
     1,
@@ -203,32 +204,29 @@ test_that('calculate_infections works various combinations of drug and vaccinati
   vaccine_antibodies_mock <- mockery::mock(c(2, 3))
   vaccine_efficacy_mock <- mockery::mock(c(.2, .3))
   bernoulli_mock <- mockery::mock(2)
-  mockery::stub(calculate_infections, 'blood_immunity', immunity_mock)
-  mockery::stub(calculate_infections, 'weibull_survival', weibull_mock)
-  mockery::stub(calculate_infections, 'calculate_pev_antibodies', vaccine_antibodies_mock)
-  mockery::stub(calculate_infections, 'calculate_pev_efficacy', vaccine_efficacy_mock)
-  mockery::stub(calculate_infections, 'bernoulli_multi_p', bernoulli_mock)
+  mockery::stub(calculate_falciparum_infections, 'blood_immunity', immunity_mock)
+  mockery::stub(calculate_falciparum_infections, 'bernoulli_multi_p', bernoulli_mock)
+  
+  local_mocked_bindings(weibull_survival = weibull_mock)
+  local_mocked_bindings(calculate_pev_antibodies = vaccine_antibodies_mock)
+  local_mocked_bindings(calculate_pev_efficacy = vaccine_efficacy_mock)
   
   # remove randomness from vaccine parameters
-  mockery::stub(
-    calculate_infections,
-    'sample_pev_param',
-    function(index, profiles, name) {
-      vnapply(index, function(i) profiles[[i]][[name]][[1]]) # return mu
-    },
-    depth = 4
-  )
+  local_mocked_bindings(sample_pev_param = function(index, profiles, name) {
+    vnapply(index, function(i) profiles[[i]][[name]][[1]]) # return mu
+  })
   
   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){
-      infection_outcome_process(timestep, target, variables, renderer, parameters)
+      falciparum_infection_outcome_process(timestep, target, variables, renderer, parameters)
     },
     size = 4
   )
   
-  infections <- calculate_infections(
+  infections <- calculate_falciparum_infections(
     variables,
     bitten_humans, 
     parameters,
@@ -236,8 +234,8 @@ test_that('calculate_infections works various combinations of drug and vaccinati
     timestep,
     infection_outcome
   )
-  
-  expect_equal(sum(infections!=0), 3)
+
+  expect_equal(sum(infection_outcome$rates!=0), 3)
   
   mockery::expect_args(immunity_mock, 1, c(.3, .5, .9), parameters)
   mockery::expect_args(
@@ -700,17 +698,16 @@ test_that('prophylaxis is considered for medicated humans', {
   
   bitten <- individual::Bitset$new(4)$insert(seq(4))
   m <- mockery::mock(seq(3))
-  mockery::stub(calculate_infections, 'bernoulli_multi_p', m)
-  
+  mockery::stub(calculate_falciparum_infections, 'bernoulli_multi_p', m)
   
   infection_outcome <- CompetingOutcome$new(
     targeted_process = function(timestep, target){
-      infection_outcome_process(timestep, target, variables, renderer, parameters)
+      falciparum_infection_outcome_process(timestep, target, variables, renderer, parameters)
     },
     size = 4
   )
   
-  infection_rates <- calculate_infections(
+  infection_rates <- calculate_falciparum_infections(
     variables,
     bitten,
     parameters,
@@ -720,7 +717,7 @@ test_that('prophylaxis is considered for medicated humans', {
   )
   
   expect_equal(
-    rate_to_prob(infection_rates[infection_rates!=0]),
+    rate_to_prob(infection_outcome$rates[infection_outcome$rates!=0]),
     c(2.491951e-07, 2.384032e-01, 5.899334e-01),
     tolerance = 1e-3
   )
diff --git a/tests/testthat/test-pev.R b/tests/testthat/test-pev.R
index f3c3d793..6abc67b1 100644
--- a/tests/testthat/test-pev.R
+++ b/tests/testthat/test-pev.R
@@ -152,26 +152,21 @@ test_that('Infection considers pev efficacy', {
 
   infection_outcome <- CompetingOutcome$new(
     targeted_process = function(timestep, target){
-      infection_process_resolved_hazard(timestep, target, variables, renderer, parameters)
+      falciparum_infection_outcome_process(timestep, target, variables, renderer, parameters)
     },
     size = 4
   )
   
   # remove randomness from infection sampling
   bernoulli_mock <- mockery::mock(c(1, 2))
-  mockery::stub(calculate_infections, 'bernoulli_multi_p', bernoulli_mock)
+  mockery::stub(calculate_falciparum_infections, 'bernoulli_multi_p', bernoulli_mock)
 
   # remove randomness from pev parameters
-  mockery::stub(
-    calculate_infections,
-    'sample_pev_param',
-    function(index, profiles, name) {
-      vnapply(index, function(i) profiles[[i]][[name]][[1]]) # return mu
-    },
-    depth = 4
-  )
-
-  infection_rates <- calculate_infections(
+  local_mocked_bindings(sample_pev_param = function(index, profiles, name) {
+    vnapply(index, function(i) profiles[[i]][[name]][[1]]) # return mu
+  })
+  
+  infection_rates <- calculate_falciparum_infections(
     variables = variables,
     bitten_humans = individual::Bitset$new(4)$insert(seq(4)),
     parameters = parameters,
@@ -179,9 +174,9 @@ test_that('Infection considers pev efficacy', {
     timestep = timestep,
     infection_outcome = infection_outcome
   )
-
+  
   expect_equal(
-    rate_to_prob(infection_rates[infection_rates!=0]),
+    rate_to_prob(infection_outcome$rates),
     c(0.590, 0.590, 0.215, 0.244),
     tolerance=1e-3
   )
diff --git a/tests/testthat/test-vivax.R b/tests/testthat/test-vivax.R
index 8fc77956..58ee404b 100644
--- a/tests/testthat/test-vivax.R
+++ b/tests/testthat/test-vivax.R
@@ -255,3 +255,62 @@ 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(vivax_infection_outcome_process, 'bernoulli_multi_p', bernoulli_mock, depth = 2)
+  mockery::stub(vivax_infection_outcome_process, 'calculate_clinical_infections', calc_mock)
+  
+  renderer <- mock_render(1)
+  infected_humans <- individual::Bitset$new(4)$insert(c(1, 2, 3, 4))
+  
+  vivax_infection_outcome_process(
+    timestep = timestep,
+    infected_humans,
+    variables,
+    renderer,
+    parameters,
+    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>