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..fc5a98ba 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$copy()$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
+      prob <- 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 <- 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)
   }
-
-  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..11426cdb 100644
--- a/R/processes.R
+++ b/R/processes.R
@@ -47,26 +47,6 @@ create_processes <- function(
     immunity_process = create_exponential_decay_process(variables$ica,
                                                         parameters$rc)
   )
-  
-  if(parameters$parasite == "falciparum"){
-    processes <- c(
-      processes,
-      # Blood immunity
-      create_exponential_decay_process(variables$ib, parameters$rb),
-      # Immunity to severe disease
-      create_exponential_decay_process(variables$ivm, parameters$rvm),
-      create_exponential_decay_process(variables$iva, parameters$rva),
-      # Immunity to detectability
-      create_exponential_decay_process(variables$id, parameters$rid)
-    )
-  } else if (parameters$parasite == "vivax"){
-    processes <- c(
-      processes,
-      # Anti-parasite immunity
-      create_exponential_decay_process(variables$iam, parameters$rm),
-      create_exponential_decay_process(variables$iaa, parameters$ra)
-    )
-  }
 
   if(parameters$parasite == "falciparum"){
     processes <- c(
@@ -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>