From a75a67122830492d7cf1710864ec7644acd74375 Mon Sep 17 00:00:00 2001 From: RJSheppard Date: Wed, 18 Sep 2024 13:24:08 +0100 Subject: [PATCH] Added in functionality to start at equilibrium, based on Michael White's equilibrium solution. This is accessed during the initial equilibrium calculation in compatability and during the variable creation. The equilibrium solution is found in malariaEquilibriumVivax. Note that we draw individuals from a joint probability distribution over age and het group to generate human state, immunities and number of hypnozoite batches. tests included to check correct equilibrium is generated and maintained. I have also added the An. koliensis parameters as found in Michael White's repo in order to compare equilibrium solutions. A few errors were identified. Total infection rates were not being calculated correctly in human_infection.R. Lagged infectivity needs to be saved before the current infectivity is extracted during the biting_process.R. Also, when immunity is boosted through infection, this should incorporate the decrease in immunity from decay (there may be better ways of achieving this). --- DESCRIPTION | 2 + NAMESPACE | 1 + R/biting_process.R | 4 +- R/compatibility.R | 133 ++++++++++++-- R/human_infection.R | 37 ++-- R/processes.R | 21 +++ R/variables.R | 124 ++++++++++--- R/vector_parameters.R | 24 +++ man/kol_params.Rd | 30 +++ man/set_equilibrium.Rd | 2 +- tests/testthat/test-vivax-biology.R | 214 ++++++++++++++++++++++ tests/testthat/test-vivax-compartmental.R | 161 ++++++++++++++++ tests/testthat/test-vivax-eq.R | 86 +++++++++ vignettes/Plasmodium_vivax.Rmd | 4 + 14 files changed, 778 insertions(+), 65 deletions(-) create mode 100644 man/kol_params.Rd create mode 100644 tests/testthat/test-vivax-biology.R create mode 100644 tests/testthat/test-vivax-compartmental.R create mode 100644 tests/testthat/test-vivax-eq.R 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.