diff --git a/DESCRIPTION b/DESCRIPTION
index b7aaa69a..4250cab2 100644
--- a/DESCRIPTION
+++ b/DESCRIPTION
@@ -68,10 +68,12 @@ Encoding: UTF-8
 LazyData: true
 Remotes:
   mrc-ide/malariaEquilibrium,
+  mrc-ide/malariaEquilibriumVivax,
   mrc-ide/individual
 Imports:
     individual (>= 0.1.17),
     malariaEquilibrium (>= 1.0.1),
+    malariaEquilibriumVivax,
     Rcpp,
     statmod,
     MASS,
diff --git a/NAMESPACE b/NAMESPACE
index 167c0f0e..d3fe01a6 100644
--- a/NAMESPACE
+++ b/NAMESPACE
@@ -10,6 +10,7 @@ export(gamb_params)
 export(get_correlation_parameters)
 export(get_init_carrying_capacity)
 export(get_parameters)
+export(kol_params)
 export(parameterise_mosquito_equilibrium)
 export(parameterise_total_M)
 export(peak_season_offset)
diff --git a/R/biting_process.R b/R/biting_process.R
index 17b75366..b28fe327 100644
--- a/R/biting_process.R
+++ b/R/biting_process.R
@@ -169,14 +169,14 @@ simulate_bites <- function(
       }
     }
 
+    lagged_infectivity$save(sum(human_infectivity * .pi), timestep)
+
     if (is.null(mixing_fn)) {
       infectivity <- lagged_infectivity$get(timestep - parameters$delay_gam)
     } else {
       infectivity <- mixing_fn(timestep=timestep)$inf[[mixing_index]]
     }
 
-    lagged_infectivity$save(sum(human_infectivity * .pi), timestep)
-
     foim <- calculate_foim(a, infectivity)
     renderer$render(paste0('FOIM_', species_name), foim, timestep)
     mu <- death_rate(f, W, Z, s_i, parameters)
diff --git a/R/compatibility.R b/R/compatibility.R
index ec71d3df..6833f695 100644
--- a/R/compatibility.R
+++ b/R/compatibility.R
@@ -4,6 +4,10 @@ inverse_param <- function(name, new_name) {
   function(params) { list(new_name, 1 / params[[name]]) }
 }
 
+product_param <- function(new_name, name_1, name_2) {
+  function(params) { list(new_name, params[[name_1]] * params[[name_2]]) }
+}
+
 mean_param <- function(new_name, name, weights) {
   function(params) {
     list(new_name, weighted.mean(params[[name]], params[[weights]]))
@@ -120,6 +124,56 @@ back_translations = list(
   kv = 'kv'
 )
 
+vivax_translations = list(
+  
+  mean_age  = 'average_age',
+  rho_age   = 'rho',
+  age_0     = 'a0',
+  N_het     = 'n_heterogeneity_groups',
+  
+  bb  = 'b',      ## mosquito -> human transmission probability
+  
+  c_PCR = 'cu',   ## human -> mosquito transmission probability (PCR)
+  c_LM  = 'ca',   ## human -> mosquito transmission probability (LM-detectable)
+  c_D   = 'cd',   ## human -> mosquito transmission probability (disease state)
+  c_T   = 'ct',   ## human -> mosquito transmission probability (treatment)
+  
+  d_E   = 'de',   ## duration of liver-stage latency
+  r_D   = inverse_param('dd', 'r_D'),   ## duraton of disease = 1/rate
+  r_T   = inverse_param('dt', 'r_T'),   ## duraton of prophylaxis = 1/rate
+  
+  r_par = inverse_param('ra', 'r_par'),    ## rate of decay of anti-parasite immunity
+  r_clin= inverse_param('rc', 'r_clin'),    ## rate of decay of clinical immunity
+  
+  mu_M  = mean_param('mum', 'mum', 'species_proportions'),  ## mosquito death rate = 1/(mosquito life expectancy)
+  Q0    = mean_param('Q0', 'Q0', 'species_proportions'),
+  blood_meal_rates = mean_param('blood_meal_rates', 'blood_meal_rates', 'species_proportions'),
+  tau_M = 'dem', 	## duration of sporogony
+  
+  ff      = 'f',      ## relapse rate
+  gamma_L = 'gammal', ## duration of liver-stage carriage
+  K_max   = 'kmax', ## maximum hypnozoite batches
+  
+  u_par        = 'ua',        ## refractory period for anti-parasite immune boosting
+  phi_LM_max   = 'philm_max', ## probability of LM_detectable infection with no immunity
+  phi_LM_min   = 'philm_min', ## probability of LM_detectable infection with maximum immunity
+  A_LM_50pc    = 'alm50',     ## blood-stage immunity scale parameter
+  K_LM         = 'klm',       ## blood-stage immunity shape parameter
+  u_clin       = 'uc',        ## refractory period for clinical immune boosting
+  phi_D_max    = 'phi0',      ## probability of clinical episode with no immunity
+  phi_D_min    = product_param("phi_D_min", "phi0", "phi1"),    ## probability of clinical episode with maximum immunity
+  A_D_50pc     = 'ic0',       ## clinical immunity scale parameter
+  K_D          = 'kc',        ## clinical immunity shape parameter
+  A_d_PCR_50pc = 'apcr50',    ## scale parameter for effect of anti-parasite immunity on PCR-detectable infection
+  K_d_PCR      = 'kpcr',      ## shape parameter for effect of anti-parasite immunity on PCR-detectable infection
+  d_PCR_max    = 'dpcr_max',  ## maximum duration on PCR-detectable infection
+  d_PCR_min    = 'dpcr_min',  ## maximum duration of PCR-detectable infection
+  d_LM         = 'da',        ## duration of LM-detectable infection
+  P_MI         = 'pcm',       ## Proportion of immunity acquired maternally
+  d_MI         = 'rm'         ## Rate of waning of maternal immunity
+  
+)
+
 #' @description translate parameter keys from the malariaEquilibrium format
 #' to ones compatible with this IBM 
 #' @param params with keys in the malariaEquilibrium format
@@ -162,6 +216,23 @@ translate_parameters <- function(params) {
   translated
 }
 
+#' @description translate parameter keys from the malariaVivaxEquilibrium format
+#' to ones compatible with this IBM
+#' @param params with keys in the malariaVivaxEquilibrium format
+#' @noRd
+translate_vivax_parameters <- function(params) {
+  translated <- params
+  for (i in 1:length(vivax_translations)) {
+    if (is.character(vivax_translations[[i]])) {
+      translated[[names(vivax_translations)[i]]] <- params[[vivax_translations[[i]]]]
+    }
+    if (is.function(vivax_translations[[i]])) {
+      translated[[names(vivax_translations)[i]]] <- vivax_translations[[i]](params)[[2]]
+    }
+  }
+  translated
+}
+
 #' @title remove parameter keys from the malariaEquilibrium format that are not used
 #' in this IBM 
 #' @param params with keys in the malariaEquilibrium format
@@ -182,7 +253,7 @@ remove_unused_equilibrium <- function(params) {
 #' @title Set equilibrium
 #' @description This will update the IBM parameters to match the
 #' equilibrium parameters and set up the initial human and mosquito population
-#' to acheive init_EIR
+#' to achieve init_EIR
 #' @param parameters model parameters to update
 #' @param init_EIR the desired initial EIR (infectious bites per person per day over the entire human
 #' population)
@@ -190,28 +261,50 @@ remove_unused_equilibrium <- function(params) {
 #' The default malariaEquilibrium parameters will be used
 #' @export
 set_equilibrium <- function(parameters, init_EIR, eq_params = NULL) {
-  if (is.null(eq_params)) {
-    eq_params <- translate_parameters(parameters)
-  } else {
+  if(parameters$parasite == "falciparum"){
+    if (is.null(eq_params)) {
+      eq_params <- translate_parameters(parameters)
+    } else {
+      parameters <- c(
+        translate_equilibrium(remove_unused_equilibrium(eq_params)),
+        parameters
+      )
+    }
+    eq <- malariaEquilibrium::human_equilibrium(
+      EIR = init_EIR,
+      ft = sum(get_treatment_coverages(parameters, 1)),
+      p = eq_params,
+      age = EQUILIBRIUM_AGES,
+      h = malariaEquilibrium::gq_normal(parameters$n_heterogeneity_groups)
+    )
+    parameters <- c(
+      list(
+        init_foim = eq$FOIM,
+        init_EIR = init_EIR,
+        eq_params = eq_params
+      ),
+      parameters
+    )
+  } else if (parameters$parasite == "vivax"){
+    
+    if (!(is.null(eq_params))) {
+      stop("Importing MalariaEquilibriumVivax parameters is not supported")
+    }
+    
+    eq <- malariaEquilibriumVivax::vivax_equilibrium(
+      EIR = init_EIR,
+      ft = sum(get_treatment_coverages(parameters, 1)),
+      p = translate_vivax_parameters(parameters),
+      age = EQUILIBRIUM_AGES
+    )
+    
     parameters <- c(
-      translate_equilibrium(remove_unused_equilibrium(eq_params)),
+      list(
+        init_foim = eq$FOIM,
+        init_EIR = init_EIR
+      ),
       parameters
     )
   }
-  eq <- malariaEquilibrium::human_equilibrium(
-    EIR = init_EIR,
-    ft = sum(get_treatment_coverages(parameters, 1)),
-    p = eq_params,
-    age = EQUILIBRIUM_AGES,
-    h = malariaEquilibrium::gq_normal(parameters$n_heterogeneity_groups)
-  )
-  parameters <- c(
-    list(
-      init_foim = eq$FOIM,
-      init_EIR = init_EIR,
-      eq_params = eq_params
-    ),
-    parameters
-  )
   parameterise_mosquito_equilibrium(parameters, init_EIR)
 }
diff --git a/R/human_infection.R b/R/human_infection.R
index fc5a98ba..c8342de0 100644
--- a/R/human_infection.R
+++ b/R/human_infection.R
@@ -119,7 +119,6 @@ calculate_infections <- function(
       )
     }
     
-    # infection_rates <- rep(0, length = parameters$human_population)
     if(parameters$parasite == "falciparum"){
       
       ## p.f models blood immunity
@@ -205,15 +204,16 @@ infection_outcome_process <- function(
       timestep
     )
     
-    boost_immunity(
-      variables$ica,
-      infected_humans,
-      variables$last_boosted_ica,
-      timestep,
-      parameters$uc
-    )
-
     if(parameters$parasite == "falciparum"){
+      
+      boost_immunity(
+        variables$ica,
+        infected_humans,
+        variables$last_boosted_ica,
+        timestep,
+        parameters$uc
+      )
+      
       boost_immunity(
         variables$id,
         infected_humans,
@@ -261,13 +261,23 @@ infection_outcome_process <- function(
         renderer,
         timestep
       )
-
+      
+      boost_immunity(
+        variables$ica,
+        infected_humans,
+        variables$last_boosted_ica,
+        timestep,
+        parameters$uc,
+        1/parameters$rc
+      )
+      
       boost_immunity(
         variables$iaa,
         infected_humans,
         variables$last_boosted_iaa,
         timestep,
-        parameters$ua
+        parameters$ua,
+        1/parameters$ra
       )
       
       ## Only S and U infections are considered in generating lm-det infections
@@ -744,7 +754,8 @@ boost_immunity <- function(
     exposed_index,
     last_boosted_variable,
     timestep,
-    delay
+    delay,
+    decay_rate = 0
 ) {
   # record who can be boosted
   exposed_index_vector <- exposed_index$to_vector()
@@ -754,7 +765,7 @@ boost_immunity <- function(
   if (sum(to_boost) > 0) {
     # boost the variable
     immunity_variable$queue_update(
-      immunity_variable$get_values(exposed_to_boost) + 1,
+      immunity_variable$get_values(exposed_to_boost) * (1-decay_rate) + 1,
       exposed_to_boost
     )
     # record last boosted
diff --git a/R/processes.R b/R/processes.R
index 11426cdb..a05a2b6b 100644
--- a/R/processes.R
+++ b/R/processes.R
@@ -47,6 +47,27 @@ 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),
+      create_hypnozoite_batch_decay_process(variables$hypnozoites, parameters$gammal)
+    )
+  }
 
   if(parameters$parasite == "falciparum"){
     processes <- c(
diff --git a/R/variables.R b/R/variables.R
index b6bd1f80..9745d026 100644
--- a/R/variables.R
+++ b/R/variables.R
@@ -82,26 +82,49 @@ create_variables <- function(parameters) {
       to_char_vector(groups)
     )
     if (!is.null(parameters$init_EIR)) {
-      eq <-	list(
-        malariaEquilibrium::human_equilibrium_no_het(
-          parameters$init_EIR,
-          sum(get_treatment_coverages(parameters, 0)),
-          parameters$eq_params,
-          EQUILIBRIUM_AGES
+      if(parameters$parasite == "falciparum"){
+        eq <- list(
+          malariaEquilibrium::human_equilibrium_no_het(
+            parameters$init_EIR,
+            sum(get_treatment_coverages(parameters, 0)),
+            parameters$eq_params,
+            EQUILIBRIUM_AGES
+          )
         )
-      )
+      } else if (parameters$parasite == "vivax"){
+        eq <- malariaEquilibriumVivax::vivax_equilibrium(
+          EIR = parameters$init_EIR,
+          ft = sum(get_treatment_coverages(parameters, 0)),
+          p = translate_vivax_parameters(parameters),
+          age = EQUILIBRIUM_AGES
+        )$states
+      }
     } else {
       eq <- NULL
     }
   }
 
   states <- c('S', 'D', 'A', 'U', 'Tr')
-  initial_states <- initial_state(parameters, initial_age, groups, eq, states)
+  
+  if(parameters$parasite == "falciparum"){
+    initial_states <- initial_state(parameters, initial_age, groups, eq, states)
+    hypnozoite_v <- NULL
+    
+  } else if (parameters$parasite == "vivax"){
+    
+    eq_v_output <- initial_state_vivax(parameters, initial_age, groups, eq, states)
+    
+    # Human states
+    initial_states <- eq_v_output$human_states
+    hypnozoite_v <- eq_v_output$hypnozoites_v
+    
+    ## Initial hypnozoites
+    hypnozoites <- individual::IntegerVariable$new(hypnozoite_v)
+    
+  }
+  
   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(
@@ -111,7 +134,8 @@ create_variables <- function(parameters) {
       groups,
       eq,
       parameters,
-      'ICM'
+      'ICM',
+      hypnozoite_v
     )
   )
   
@@ -124,7 +148,8 @@ create_variables <- function(parameters) {
       groups,
       eq,
       parameters,
-      'ICA'
+      'ICA',
+      hypnozoite_v
     )
   )
 
@@ -189,7 +214,8 @@ create_variables <- function(parameters) {
         groups,
         eq,
         parameters,
-        'IAM'
+        'IAM',
+        hypnozoite_v
       )
     )
     
@@ -202,11 +228,12 @@ create_variables <- function(parameters) {
         groups,
         eq,
         parameters,
-        'IAA'
+        'IAA',
+        hypnozoite_v
       )
     )
   }
-
+  
   # Initialise infectiousness of humans -> mosquitoes
   # NOTE: not yet supporting initialisation of infectiousness of Treated individuals
   infectivity_values <- rep(0, get_human_population(parameters, 0))
@@ -386,7 +413,8 @@ initial_immunity <- function(
   groups = NULL,
   eq = NULL,
   parameters = NULL,
-  eq_name = NULL
+  eq_name = NULL,
+  hyp = NULL
   ) {
   if (!is.null(eq)) {
     age <- age / 365
@@ -395,7 +423,12 @@ initial_immunity <- function(
       function(i) {
         g <- groups[[i]]
         a <- age[[i]]
-        eq[[g]][which.max(a < eq[[g]][, 'age']), eq_name]
+        if(parameters$parasite == "falciparum"){
+          eq[[g]][which.max(a < eq[[g]][, 'age']), eq_name]
+        } else if (parameters$parasite == "vivax"){
+          h <- hyp[[i]]
+          eq[[eq_name]][which.max(a < eq$Age[-1]), g, h+1]
+        }
       }
     ))
   }
@@ -423,6 +456,30 @@ initial_state <- function(parameters, age, groups, eq, states) {
   rep(ibm_states, times = calculate_initial_counts(parameters))
 }
 
+initial_state_vivax <- function(parameters, age, groups, eq, states) {
+  # vivax human states and hypnozoites must be calculated over a combined probability distribution
+  ibm_states <- states
+  if (!is.null(eq)) {
+    eq_states <- c('S', 'D', 'A', 'U', 'T')
+    age <- age / 365
+    human_states_pop <- sapply(
+      seq_along(age),
+      function(i) {
+        g <- groups[[i]]
+        a <- age[[i]]
+        human_states <- c(expand.grid("hyp" = 0:parameters$kmax, "human_state" = eq_states))
+        probs <- c(sapply(eq_states, function(state){eq[[state]][which.max(a < eq$Age[-1]), g, ]}))
+        smp <- sample(x = 1:length(probs), size = 1, replace = T, prob = probs)
+        return(c(human_states$human_state[smp],human_states$hyp[smp]))
+      }
+    )
+    return(list(human_states = states[human_states_pop[1,]],
+                hypnozoites_v = human_states_pop[2,]))
+  }
+  return(list(human_states = rep(ibm_states, calculate_initial_counts(parameters)),
+              hypnozoites_v = rep(parameters$init_hyp, parameters$human_population)))
+}
+
 calculate_initial_counts <- function(parameters) {
   pop <- get_human_population(parameters, 0)
   initial_counts <- round(
@@ -441,17 +498,26 @@ calculate_initial_counts <- function(parameters) {
 
 calculate_eq <- function(het_nodes, parameters) {
   ft <- sum(get_treatment_coverages(parameters, 0))
-	lapply(
-		het_nodes,
-		function(n) {
-			malariaEquilibrium::human_equilibrium_no_het(
-				parameters$init_EIR * calculate_zeta(n, parameters),
-				ft,
-				parameters$eq_params,
-        EQUILIBRIUM_AGES
-			)
-		}
-	)
+  if(parameters$parasite == "falciparum"){
+    lapply(
+      het_nodes,
+      function(n) {
+        malariaEquilibrium::human_equilibrium_no_het(
+          parameters$init_EIR * calculate_zeta(n, parameters),
+          ft,
+          parameters$eq_params,
+          EQUILIBRIUM_AGES
+        )
+      }
+    )
+  } else if (parameters$parasite == "vivax"){
+    eq <- malariaEquilibriumVivax::vivax_equilibrium(
+      EIR = parameters$init_EIR,
+      ft = sum(get_treatment_coverages(parameters, 0)),
+      age = EQUILIBRIUM_AGES,
+      p = translate_vivax_parameters(parameters)
+    )$states
+  }
 }
 
 calculate_zeta <- function(zeta_norm, parameters) {
diff --git a/R/vector_parameters.R b/R/vector_parameters.R
index 2024d253..11040d82 100644
--- a/R/vector_parameters.R
+++ b/R/vector_parameters.R
@@ -94,6 +94,30 @@ steph_params <- list(
   mum = .112
 )
 
+#' @title Preset parameters for the An. koliensis vector
+#' @details Default parameters:
+#' species: "kol"
+#' blood_meal_rates: 0.3333333
+#' foraging_time: .68
+#' Q0: 0.5
+#' phi_bednets: 0.9
+#' phi_indoors: 0.5
+#' mum: 0.1666667
+#'
+#' parameters reference:
+#' https://github.com/MWhite-InstitutPasteur/Pvivax_IBM/blob/master/koliensis_parameters.txt
+#' Value for blood meal rate is assumed in absence of  phi are from:
+#' @export
+kol_params <- list(
+  species = 'kol',
+  blood_meal_rates = 1/3,
+  foraging_time = 0.68,
+  Q0 = 0.5,
+  phi_bednets = 0.9,
+  phi_indoors = 0.5,
+  mum = 0.1666667
+)
+
 #' @title Parameterise the mosquito species to use in the model
 #'
 #' @param parameters the model parameters
diff --git a/man/kol_params.Rd b/man/kol_params.Rd
new file mode 100644
index 00000000..3047a46c
--- /dev/null
+++ b/man/kol_params.Rd
@@ -0,0 +1,30 @@
+% Generated by roxygen2: do not edit by hand
+% Please edit documentation in R/vector_parameters.R
+\docType{data}
+\name{kol_params}
+\alias{kol_params}
+\title{Preset parameters for the An. koliensis vector}
+\format{
+An object of class \code{list} of length 7.
+}
+\usage{
+kol_params
+}
+\description{
+Preset parameters for the An. koliensis vector
+}
+\details{
+Default parameters:
+species: "kol"
+blood_meal_rates: 0.3333333
+foraging_time: .68
+Q0: 0.5
+phi_bednets: 0.9
+phi_indoors: 0.5
+mum: 0.1666667
+
+parameters reference:
+https://github.com/MWhite-InstitutPasteur/Pvivax_IBM/blob/master/koliensis_parameters.txt
+Value for blood meal rate is assumed in absence of  phi are from:
+}
+\keyword{datasets}
diff --git a/man/set_equilibrium.Rd b/man/set_equilibrium.Rd
index 0ed691b8..1268cee9 100644
--- a/man/set_equilibrium.Rd
+++ b/man/set_equilibrium.Rd
@@ -18,5 +18,5 @@ The default malariaEquilibrium parameters will be used}
 \description{
 This will update the IBM parameters to match the
 equilibrium parameters and set up the initial human and mosquito population
-to acheive init_EIR
+to achieve init_EIR
 }
diff --git a/tests/testthat/test-vivax-biology.R b/tests/testthat/test-vivax-biology.R
new file mode 100644
index 00000000..df471279
--- /dev/null
+++ b/tests/testthat/test-vivax-biology.R
@@ -0,0 +1,214 @@
+test_that('total_M and EIR functions are consistent with equilibrium EIR', {
+  skip_on_ci()
+  expected_EIR <- c(1, 5, 10, 100)
+  population <- 1000
+  parameters <- get_parameters(parasite = "vivax",
+                               list(
+                                 human_population = population,
+                                 enable_heterogeneity = FALSE
+                               )
+  )
+  actual_EIR <- vnapply(expected_EIR, function(EIR) {
+    parameters <- set_equilibrium(parameters, EIR)
+    
+    #set up arguments for EIR calculation
+    variables <- create_variables(parameters)
+    vector_models <- parameterise_mosquito_models(parameters, 1)
+    solvers <- parameterise_solvers(vector_models, parameters)
+    estimated_eir <- calculate_eir(1, solvers, variables, parameters, 0)
+    mean(estimated_eir) / population * 365
+  })
+  
+  expect_equal(
+    actual_EIR,
+    expected_EIR,
+    tolerance = 1e-4
+  )
+})
+
+test_that('total_M and EIR functions are consistent with equilibrium EIR (with het)', {
+  skip_on_ci()
+  population <- 1000
+  
+  expected_EIR <- c(1, 5, 10, 100)
+  
+  parameters <- get_parameters(parasite = "vivax",
+                               list(human_population = population))
+  
+  actual_EIR <- vnapply(expected_EIR, function(EIR) {
+    parameters <- set_equilibrium(parameters, EIR)
+    
+    variables <- create_variables(parameters)
+    vector_models <- parameterise_mosquito_models(parameters, 1)
+    solvers <- parameterise_solvers(vector_models, parameters)
+    estimated_eir <- calculate_eir(1, solvers, variables, parameters, 0)
+    mean(estimated_eir) / population * 365
+  })
+  
+  expect_equal(
+    actual_EIR,
+    expected_EIR,
+    tolerance = 1e-4
+  )
+})
+
+test_that('FOIM is consistent with equilibrium', {
+  skip_on_ci()
+  population <- 10000
+  
+  EIRs <- c(1, 5, 10, 100)
+  
+  eq_params <- get_parameters(parasite = "vivax")
+  
+  ages <- 0:800 / 10
+  expected_foim <- vnapply(
+    EIRs,
+    function(EIR) {
+      malariaEquilibriumVivax::vivax_equilibrium(
+        age = EQUILIBRIUM_AGES,
+        ft = 0,
+        EIR = EIR,
+        p =  translate_vivax_parameters(eq_params))$FOIM
+    }
+  )
+  
+  actual_foim <- vnapply(
+    EIRs,
+    function(EIR) {
+      parameters <- get_parameters(parasite = "vivax", list(human_population = population))
+      parameters <- set_equilibrium(parameters, EIR)
+      variables <- create_variables(parameters)
+      a <- human_blood_meal_rate(1, variables, parameters, 0)
+      age <- get_age(variables$birth$get_values(), 0)
+      psi <- unique_biting_rate(age, parameters)
+      zeta <- variables$zeta$get_values()
+      .pi <- human_pi(psi, zeta)
+      calculate_foim(a, sum(.pi * variables$infectivity$get_values()))
+    }
+  )
+  expect_equal(
+    expected_foim,
+    actual_foim,
+    tolerance = 5e-3
+  )
+})
+
+test_that('phi is consistent with equilibrium at high EIR (no het)', {
+  skip_on_ci()
+  population <- 1e5
+  
+  EIR <- 100
+  ages <- 0:800 / 10
+  
+  parameters <- get_parameters(parasite = "vivax", list(
+    human_population = population,
+    n_heterogeneity_groups = 1,
+    enable_heterogeneity = FALSE
+  ))
+  parameters <- set_equilibrium(parameters, EIR)
+  
+  eq <- malariaEquilibriumVivax::vivax_equilibrium(
+    EIR = EIR,
+    ft = 0,
+    p = translate_vivax_parameters(parameters),
+    age = EQUILIBRIUM_AGES)$states
+  
+  variables <- create_variables(parameters)
+  expect_equal(
+    mean(
+      clinical_immunity(
+        variables$ica$get_values(),
+        variables$icm$get_values(),
+        parameters
+      )
+    ),
+    sum(c(eq$phi_clin * eq$HH)),
+    tolerance=1e-2
+  )
+  
+  # Check phi_patent
+  expect_equal(
+    mean(
+      anti_parasite_immunity(min = parameters$philm_min, max = parameters$philm_max,
+                             a50 = parameters$alm50, k = parameters$klm,
+                             iaa = variables$iaa$get_values(),
+                             iam = variables$iam$get_values()
+      )
+    ),
+    sum(c(eq$phi_patent * eq$HH)),
+    tolerance=1e-2
+  )
+  
+  # Check r_PCR
+  expect_equal(
+    mean(
+      1/anti_parasite_immunity(min = parameters$dpcr_min, max = parameters$dpcr_max,
+                               a50 = parameters$apcr50, k = parameters$kpcr,
+                               iaa = variables$iaa$get_values(),
+                               iam = variables$iam$get_values()
+      )
+    ),
+    sum(c(eq$r_PCR * eq$HH)),
+    tolerance=1e-2
+  )
+  
+})
+
+test_that('phi is consistent with equilibrium at high EIR', {
+  skip_on_ci()
+  population <- 1e5
+  
+  EIR <- 100
+  ages <- 0:999 / 10
+  
+  parameters <- get_parameters(parasite = "vivax",
+                               list(human_population = population))
+  parameters <- set_equilibrium(parameters, EIR)
+  
+  eq <- malariaEquilibriumVivax::vivax_equilibrium(
+    age = EQUILIBRIUM_AGES,
+    ft = 0,
+    EIR = EIR,
+    p = translate_vivax_parameters(parameters))$states
+  
+  variables <- create_variables(parameters)
+  het <- statmod::gauss.quad.prob(parameters$n_heterogeneity_groups, dist = "normal")
+  expect_equal(
+    mean(
+      clinical_immunity(
+        variables$ica$get_values(),
+        variables$icm$get_values(),
+        parameters
+      )
+    ),
+    sum(c(eq$phi_clin * eq$HH)),
+    tolerance=1e-2
+  )
+  
+  # Check phi_patent
+  expect_equal(
+    mean(
+      anti_parasite_immunity(min = parameters$philm_min, max = parameters$philm_max,
+                             a50 = parameters$alm50, k = parameters$klm,
+                             iaa = variables$iaa$get_values(),
+                             iam = variables$iam$get_values()
+      )
+    ),
+    sum(c(eq$phi_patent * eq$HH)),
+    tolerance=1e-2
+  )
+  
+  # Check r_PCR
+  expect_equal(
+    mean(
+      1/anti_parasite_immunity(min = parameters$dpcr_min, max = parameters$dpcr_max,
+                               a50 = parameters$apcr50, k = parameters$kpcr,
+                               iaa = variables$iaa$get_values(),
+                               iam = variables$iam$get_values()
+      )
+    ),
+    sum(c(eq$r_PCR * eq$HH)),
+    tolerance=1e-2
+  )
+  
+})
diff --git a/tests/testthat/test-vivax-compartmental.R b/tests/testthat/test-vivax-compartmental.R
new file mode 100644
index 00000000..06a69ed0
--- /dev/null
+++ b/tests/testthat/test-vivax-compartmental.R
@@ -0,0 +1,161 @@
+test_that('vivax ODE stays at equilibrium with a constant total_M', {
+  parameters <- get_parameters(list(
+    individual_mosquitoes = TRUE
+  ), parasite = "vivax")
+  total_M <- 1000
+  timesteps <- 365 * 10
+  f <- parameters$blood_meal_rates
+  models <- parameterise_mosquito_models(parameters, timesteps)
+  solvers <- parameterise_solvers(models, parameters)
+  
+  counts <- c()
+  
+  for (t in seq(timesteps)) {
+    counts <- rbind(counts, c(t, solvers[[1]]$get_states()))
+    aquatic_mosquito_model_update(
+      models[[1]]$.model,
+      total_M,
+      f,
+      parameters$mum
+    )
+    solvers[[1]]$step()
+  }
+  
+  expected <- c()
+  equilibrium <- initial_mosquito_counts(
+    parameters,
+    1,
+    parameters$init_foim,
+    parameters$total_M
+  )[ODE_INDICES]
+  for (t in seq(timesteps)) {
+    expected <- rbind(expected, c(t, equilibrium))
+  }
+  
+  expect_equal(counts, expected, tolerance=1e-4)
+})
+
+test_that('vivax Adult ODE stays at equilibrium with a constant foim and mu', {
+  parameters <- get_parameters(parasite = "vivax")
+  parameters <- set_equilibrium(parameters, 100.)
+  f <- parameters$blood_meal_rates
+  timesteps <- 365 * 10
+  models <- parameterise_mosquito_models(parameters, timesteps)
+  solvers <- parameterise_solvers(models, parameters)
+  
+  counts <- c()
+  
+  for (t in seq(timesteps)) {
+    states <- solvers[[1]]$get_states()
+    counts <- rbind(counts, c(t, states))
+    adult_mosquito_model_update(
+      models[[1]]$.model,
+      parameters$mum,
+      parameters$init_foim,
+      states[ADULT_ODE_INDICES['Sm']],
+      f
+    )
+    solvers[[1]]$step()
+  }
+  
+  expected <- c()
+  equilibrium <- initial_mosquito_counts(
+    parameters,
+    1,
+    parameters$init_foim,
+    parameters$total_M
+  )
+  
+  for (t in seq(timesteps)) {
+    expected <- rbind(expected, c(t, equilibrium))
+  }
+  expect_equal(counts, expected, tolerance=1e-3)
+})
+
+test_that('vivax ODE stays at equilibrium with low total_M', {
+  total_M <- 10
+  parameters <- get_parameters(list(
+    total_M = total_M,
+    individual_mosquitoes = TRUE
+  ), parasite = "vivax")
+  timesteps <- 365 * 10
+  f <- parameters$blood_meal_rates
+  models <- parameterise_mosquito_models(parameters, timesteps)
+  solvers <- parameterise_solvers(models, parameters)
+  
+  
+  counts <- c()
+  
+  for (t in seq(timesteps)) {
+    counts <- rbind(counts, c(t, solvers[[1]]$get_states()))
+    aquatic_mosquito_model_update(
+      models[[1]]$.model,
+      total_M,
+      f,
+      parameters$mum
+    )
+    solvers[[1]]$step()
+  }
+  
+  expected <- c()
+  equilibrium <- initial_mosquito_counts(
+    parameters,
+    1,
+    parameters$init_foim,
+    parameters$total_M
+  )[ODE_INDICES]
+  
+  for (t in seq(timesteps)) {
+    expected <- rbind(expected, c(t, equilibrium))
+  }
+  
+  expect_equal(counts, expected, tolerance=1e-4)
+})
+
+
+test_that('vivax changing total_M stabilises', {
+  total_M_0 <- 500
+  total_M_1 <- 400
+  parameters <- get_parameters(list(
+    total_M = total_M_0,
+    individual_mosquitoes = TRUE
+  ), parasite = "vivax")
+  f <- parameters$blood_meal_rates
+  timesteps <- 365 * 10
+  models <- parameterise_mosquito_models(parameters, timesteps)
+  solvers <- parameterise_solvers(models, parameters)
+  
+  change <- 50
+  burn_in <- 365 * 5
+  
+  counts <- c()
+  
+  for (t in seq(timesteps)) {
+    counts <- rbind(counts, c(t, solvers[[1]]$get_states()))
+    if (t < change) {
+      total_M <- total_M_0
+    } else {
+      total_M <- total_M_1
+    }
+    aquatic_mosquito_model_update(
+      models[[1]]$.model,
+      total_M,
+      f,
+      parameters$mum
+    )
+    solvers[[1]]$step()
+  }
+  
+  initial_eq <- initial_mosquito_counts(
+    parameters,
+    1,
+    parameters$init_foim,
+    parameters$total_M
+  )[ODE_INDICES]
+  final_eq <- counts[burn_in, ODE_INDICES + 1]
+  
+  expect_equal(counts[1,], c(1, initial_eq), tolerance=1e-4)
+  expect_equal(counts[timesteps,], c(timesteps, final_eq), tolerance=1e-4)
+  expect_false(isTRUE(all.equal(initial_eq, final_eq)))
+  expect_false(any(is.na(counts)))
+})
diff --git a/tests/testthat/test-vivax-eq.R b/tests/testthat/test-vivax-eq.R
new file mode 100644
index 00000000..f008d889
--- /dev/null
+++ b/tests/testthat/test-vivax-eq.R
@@ -0,0 +1,86 @@
+test_that('Initial states are consistent with equilibrium', {
+  skip_on_ci()
+  population <- 10000
+  
+  EIRs <- c(1, 5, 10, 100, 0.1*365)
+  
+  eq_params <- get_parameters(parasite = "vivax", overrides = list(human_population = population))
+  
+  expected_states <- sapply(
+    EIRs,
+    function(EIR) {
+      eq <- malariaEquilibriumVivax::vivax_equilibrium(
+        age = EQUILIBRIUM_AGES,
+        ft = 0,
+        EIR = EIR,
+        p = translate_vivax_parameters(eq_params))$states
+      return(sapply(c("S","D","A","U","T"), function(state){sum(eq[[state]])}))
+    })
+  
+  actual_states <- sapply(
+    EIRs,
+    function(EIR) {
+      parms <- set_equilibrium(eq_params, init_EIR = EIR)
+      vars <- create_variables(parameters = parms)
+      
+      return(c(vars$state$get_size_of("S"),
+               vars$state$get_size_of("D"),
+               vars$state$get_size_of("A"),
+               vars$state$get_size_of("U"),
+               vars$state$get_size_of("Tr")))
+    })/population
+  
+  expect_equal(object = c(expected_states), expected =  c(actual_states), tolerance = 1E-2)
+  
+})
+
+test_that('Initial immunities are consistent with equilibrium', {
+  skip_on_ci()
+  population <- 10000
+  
+  EIRs <- c(1, 5, 10, 100, 0.1*365)
+  
+  eq_params <- get_parameters(parasite = "vivax", overrides = list(human_population = population))
+  eq_params <- set_species(parameters = eq_params, species = list(arab_params), proportions = 1)
+  
+  expected_averages <- sapply(
+    EIRs,
+    function(EIR) {
+      het <- malariaEquilibrium::gq_normal(eq_params$n_heterogeneity_groups)
+      eq <- malariaEquilibriumVivax::vivax_equilibrium(
+        age = EQUILIBRIUM_AGES,
+        ft = 0,
+        EIR = EIR,
+        p = translate_vivax_parameters(eq_params))$states
+      
+      return(c(sapply(c("ICA","ICM","IAA","IAM"), function(x){sum(eq$HH * eq[[x]])}),
+               sum(colSums(eq$HH, dims = 2) * c(1:dim(eq$HH)[3]-1))))
+    })
+  
+  
+  actual_averages <- sapply(
+    EIRs,
+    function(EIR) {
+      
+      parms <- set_equilibrium(eq_params, init_EIR = EIR)
+      vars <- create_variables(parameters = parms)
+      
+      return(c(mean(vars$ica$get_values()),
+               mean(vars$icm$get_values()),
+               mean(vars$iaa$get_values()),
+               mean(vars$iam$get_values()),
+               mean(vars$hypnozoites$get_values())))})
+  
+  expect_equal(object = c(expected_averages), expected =  c(actual_averages), tolerance = 1E-1)
+  
+})
+
+test_that('vivax equilibrium works with multiple species', {
+  skip_on_ci()
+  
+  simparams_pv <- get_parameters(parasite="vivax")
+  params_species_pv <- set_species(parameters = simparams_pv,
+                                   species = list(arab_params, fun_params, gamb_params),
+                                   proportions = c(0.1,0.3,0.6))
+  expect_no_error(set_equilibrium(params_species_pv, 6))
+})
diff --git a/vignettes/Plasmodium_vivax.Rmd b/vignettes/Plasmodium_vivax.Rmd
index 67c0227e..a9a5b22e 100644
--- a/vignettes/Plasmodium_vivax.Rmd
+++ b/vignettes/Plasmodium_vivax.Rmd
@@ -75,6 +75,10 @@ While the *P. falciparum* model calculates the onward infectivity of asymptomati
 
 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.
 
+### Equilibrium
+
+The malariaEquilibriumVivax package has been designed to calculate the equilibrium solution for the *P. vivax* model. This equilibrium is calculated by assigning an initial EIR value via the `set_equilibrium` function, as  in the *P. falciparum* model. Based on the age and heterogeneity group of an individual, it assigns a human disease state, clinical and anti-parasite immunities and number of hypnozoite batches.
+
 ### 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>