diff --git a/R/disease_progression.R b/R/disease_progression.R index 436a52c0..3347bedc 100644 --- a/R/disease_progression.R +++ b/R/disease_progression.R @@ -6,18 +6,34 @@ #' @param progression_outcome competing hazards object for disease progression rates #' @noRd create_progression_rates_process <- function( + parameters, variables, progression_outcome ) { function(timestep){ target <- variables$state$get_index_of("S")$not() + progression_rates <- variables$progression_rates$get_values(target) + if (parameters$parasite == "vivax"){ + # p.v subpatent recovery is immunity-dependent + progression_rates <- variables$progression_rates$get_values(target) + u_index <- variables$state$get_index_of("U") + target_u <- bitset_index(target, u_index) + progression_rates[target_u] <- + 1 / anti_parasite_immunity( + min = parameters$dpcr_min, + max = parameters$dpcr_max, + a50 = parameters$apcr50, + k = parameters$kpcr, + iaa = variables$iaa$get_values(index = u_index), + iam = variables$iam$get_values(index = u_index) + ) + } progression_outcome$set_rates( target, - variables$progression_rates$get_values(target)) + progression_rates) } } - #' @title Disease progression outcomes #' @description Following resolution of competing hazards, update state and #' infectivity of sampled individuals diff --git a/R/human_infection.R b/R/human_infection.R index 09597149..3977c0ba 100644 --- a/R/human_infection.R +++ b/R/human_infection.R @@ -172,13 +172,23 @@ infection_outcome_process <- function( timestep, parameters$uc ) - boost_immunity( - variables$id, - infected_humans, - variables$last_boosted_id, - timestep, - parameters$ud - ) + if(parameters$parasite == "falciparum"){ + boost_immunity( + variables$id, + infected_humans, + variables$last_boosted_id, + timestep, + parameters$ud + ) + } else if (parameters$parasite == "vivax"){ + boost_immunity( + variables$iaa, + infected_humans, + variables$last_boosted_iaa, + timestep, + parameters$ua + ) + } } clinical_infections <- calculate_clinical_infections( @@ -624,3 +634,9 @@ blood_immunity <- function(ib, parameters) { (1 + (ib / parameters$ib0) ** parameters$kb) ) } + +# Implemented from White et al., 2018 - Supplementary Information +anti_parasite_immunity <- function(min, max, a50, k, iaa, iam){ + min + (max - min) / ( + 1 + ((iaa + iam) / a50) ** k) +} \ No newline at end of file diff --git a/R/mortality_processes.R b/R/mortality_processes.R index 9e527e2a..6396414e 100644 --- a/R/mortality_processes.R +++ b/R/mortality_processes.R @@ -69,6 +69,10 @@ sample_maternal_immunity <- function(variables, target, timestep, parameters) { birth_ivm <- variables$iva$get_values(mothers) * parameters$pvm variables$icm$queue_update(birth_icm, target_group) variables$ivm$queue_update(birth_ivm, target_group) + if(parameters$parasite == "vivax"){ + birth_iam <- variables$iaa$get_values(mothers) * parameters$pcm + variables$iam$queue_update(birth_iam, target_group) + } } } } @@ -93,15 +97,19 @@ reset_target <- function(variables, events, target, state, parameters, timestep) # non-maternal immunity variables$last_boosted_ica$queue_update(-1, target) variables$last_boosted_iva$queue_update(-1, target) - variables$last_boosted_id$queue_update(-1, target) variables$ica$queue_update(0, target) variables$iva$queue_update(0, target) - variables$id$queue_update(0, target) variables$state$queue_update(state, target) if(parameters$parasite == "falciparum"){ variables$last_boosted_ib$queue_update(-1, target) + variables$last_boosted_id$queue_update(-1, target) variables$ib$queue_update(0, target) + variables$id$queue_update(0, target) + + } else if (parameters$parasite == "vivax"){ + variables$last_boosted_iaa$queue_update(-1, target) + variables$iaa$queue_update(0, target) } # treatment diff --git a/R/processes.R b/R/processes.R index 39385ac3..b9b2ee82 100644 --- a/R/processes.R +++ b/R/processes.R @@ -49,9 +49,7 @@ create_processes <- function( immunity_process = create_exponential_decay_process(variables$ica, parameters$rc), immunity_process = create_exponential_decay_process(variables$iva, - parameters$rva), - immunity_process = create_exponential_decay_process(variables$id, - parameters$rid) + parameters$rva) ) if(parameters$parasite == "falciparum"){ @@ -59,7 +57,16 @@ create_processes <- function( processes, # Blood immunity immunity_process = create_exponential_decay_process(variables$ib, - parameters$rb) + parameters$rb), + # Immunity to detectability + 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) ) } @@ -126,6 +133,7 @@ create_processes <- function( processes <- c( processes, progression_process = create_progression_rates_process( + parameters, variables, progression_outcome ), @@ -185,9 +193,11 @@ create_processes <- function( # Rendering # ========= - imm_var_names <- c('ica', 'icm', 'id', 'iva', 'ivm') + imm_var_names <- c('ica', 'icm', 'iva', 'ivm') if(parameters$parasite == "falciparum"){ - imm_var_names <- c(imm_var_names, 'ib') + imm_var_names <- c(imm_var_names, 'ib', 'id') + } else if (parameters$parasite == "vivax"){ + imm_var_names <- c(imm_var_names, 'iaa', 'iam') } processes <- c( diff --git a/R/render.R b/R/render.R index 5ba9d79d..34f4bd99 100644 --- a/R/render.R +++ b/R/render.R @@ -6,9 +6,10 @@ in_age_range <- function(birth, timestep, lower, upper) { #' #' @description renders prevalence numerators and denominators for individuals #' detected by light microscopy and pcr, where those infected asymptomatically by -#' P. falciparum have reduced probability of infection due to detectability +#' p.f has reduced probability of infection due to detectability #' immunity (reported as an integer sample: n_detect_lm, and summing over -#' detection probabilities: p_detect_lm) +#' detection probabilities: p_detect_lm), whearas p.v human states are defined +#' explicitly by lm and pcr detectability. #' #' @param state human infection state #' @param birth variable for birth of the individual @@ -26,12 +27,17 @@ create_prevalence_renderer <- function( ) { function(timestep) { asymptomatic <- state$get_index_of('A') - prob <- probability_of_detection( - get_age(birth$get_values(asymptomatic), timestep), - immunity$get_values(asymptomatic), - parameters - ) - asymptomatic_detected <- bitset_at(asymptomatic, bernoulli_multi_p(prob)) + + if(parameters$parasite == "falciparum"){ + prob <- probability_of_detection( + get_age(birth$get_values(asymptomatic), timestep), + immunity$get_values(asymptomatic), + parameters + ) + asymptomatic_detected <- bitset_at(asymptomatic, bernoulli_multi_p(prob)) + } else if (parameters$parasite == "vivax") { + asymptomatic_detected <- asymptomatic + } clinically_detected <- state$get_index_of(c('Tr', 'D')) detected <- clinically_detected$copy()$or(asymptomatic_detected) @@ -46,13 +52,15 @@ create_prevalence_renderer <- function( in_age$copy()$and(detected)$size(), timestep ) - renderer$render( - paste0('p_detect_lm_', lower, '_', upper), - in_age$copy()$and(clinically_detected)$size() + sum( - prob[bitset_index(asymptomatic, in_age)] - ), - timestep - ) + if(parameters$parasite == "falciparum"){ + renderer$render( + paste0('p_detect_lm_', lower, '_', upper), + in_age$copy()$and(clinically_detected)$size() + sum( + prob[bitset_index(asymptomatic, in_age)] + ), + timestep + ) + } renderer$render( paste0('n_detect_pcr_', lower, '_', upper), pcr_detected$copy()$and(in_age)$size(), diff --git a/R/variables.R b/R/variables.R index db635ba8..8567998c 100644 --- a/R/variables.R +++ b/R/variables.R @@ -10,12 +10,14 @@ #' * birth - an integer representing the timestep when this individual was born #' * last_boosted_* - the last timestep at which this individual's immunity was #' boosted for tracking grace periods in the boost of immunity +#' * IAM - Maternal anti-parasite immunity (p.v only) #' * ICM - Maternal immunity to clinical disease -#' * IVM - Maternal immunity to severe disease +#' * IVM - Maternal immunity to severe disease (p.f only) #' * IB - Pre-erythrocytic immunity (p.f only) +#' * IAA - Acquired anti-parasite immunity (p.v only) #' * ICA - Acquired immunity to clinical disease -#' * IVA - Acquired immunity to severe disease -#' * ID - Acquired immunity to detectability +#' * IVA - Acquired immunity to severe disease (p.f only) +#' * ID - Acquired immunity to detectability (p.f 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 @@ -96,9 +98,9 @@ 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) + last_boosted_ica <- individual::DoubleVariable$new(rep(-1, size)) last_boosted_iva <- individual::DoubleVariable$new(rep(-1, size)) - last_boosted_id <- individual::DoubleVariable$new(rep(-1, size)) # Maternal immunity icm <- individual::DoubleVariable$new( @@ -136,8 +138,47 @@ create_variables <- function(parameters) { 'IB' ) ) + + # Acquired immunity to detectability + last_boosted_id <- individual::DoubleVariable$new(rep(-1, size)) + id <- individual::DoubleVariable$new( + initial_immunity( + parameters$init_id, + initial_age, + groups, + eq, + parameters, + 'ID' + ) + ) + + } else if (parameters$parasite == "vivax"){ + # Acquired anti-parasite immunity + last_boosted_iaa <- individual::DoubleVariable$new(rep(-1, size)) + iaa <- individual::DoubleVariable$new( + initial_immunity( + parameters$init_iaa, + initial_age, + groups, + eq, + parameters, + 'IAA' + ) + ) + + # Maternal anti-parasite immunity + iam <- individual::DoubleVariable$new( + initial_immunity( + parameters$init_iam, + initial_age, + groups, + eq, + parameters, + 'IAM' + ) + ) } - + # Acquired immunity to clinical disease ica <- individual::DoubleVariable$new( initial_immunity( @@ -149,6 +190,7 @@ create_variables <- function(parameters) { 'ICA' ) ) + # Acquired immunity to severe disease iva <- individual::DoubleVariable$new( initial_immunity( @@ -160,18 +202,7 @@ create_variables <- function(parameters) { 'IVA' ) ) - # Acquired immunity to detectability - id <- individual::DoubleVariable$new( - initial_immunity( - parameters$init_id, - initial_age, - groups, - eq, - parameters, - 'ID' - ) - ) - + # Initialise infectiousness of humans -> mosquitoes # NOTE: not yet supporting initialisation of infectiousness of Treated individuals infectivity_values <- rep(0, get_human_population(parameters, 0)) @@ -205,7 +236,17 @@ create_variables <- function(parameters) { progression_rate_values <- rep(0, get_human_population(parameters, 0)) progression_rate_values[diseased] <- 1/parameters$dd progression_rate_values[asymptomatic] <- 1/parameters$da - progression_rate_values[subpatent] <- 1/parameters$du + if(parameters$parasite == "falciparum"){ + # p.f subpatent recovery rate is constant + progression_rate_values[subpatent] <- 1/parameters$du + } else if (parameters$parasite == "vivax"){ + # p.v subpatent recovery rate is immunity-dependent + progression_rate_values[subpatent] <- 1/anti_parasite_immunity( + parameters$dpcr_min, parameters$dpcr_max, parameters$apcr50, parameters$kpcr, + iaa$get_values(subpatent), + iam$get_values(subpatent) + ) + } progression_rate_values[treated] <- 1/parameters$dt # Initialise the disease progression rate variable @@ -229,12 +270,10 @@ create_variables <- function(parameters) { birth = birth, last_boosted_ica = last_boosted_ica, last_boosted_iva = last_boosted_iva, - last_boosted_id = last_boosted_id, icm = icm, ivm = ivm, ica = ica, iva = iva, - id = id, zeta = zeta, zeta_group = zeta_group, infectivity = infectivity, @@ -252,7 +291,15 @@ create_variables <- function(parameters) { if(parameters$parasite == "falciparum"){ variables <- c(variables, last_boosted_ib = last_boosted_ib, - ib = ib + last_boosted_id = last_boosted_id, + ib = ib, + id = id + ) + } else if (parameters$parasite == "vivax"){ + variables <- c(variables, + last_boosted_iaa = last_boosted_iaa, + iaa = iaa, + iam = iam ) } diff --git a/data-raw/parasite_parameters.csv b/data-raw/parasite_parameters.csv index 255254c0..073e5a37 100644 --- a/data-raw/parasite_parameters.csv +++ b/data-raw/parasite_parameters.csv @@ -86,16 +86,12 @@ vivax,f,0.02439024, White (2018); doi: 10.1186/s12936-018-2318-1 vivax,gammal,0.002610966, White (2018): doi: 10.1186/s12936-018-2318-1 vivax,init_hyp,0, vivax,kmax,10, White (2018); doi: 10.1038/s41467-018-05860-8 -vivax,du,110,to_be_removed vivax,init_iva,0,to_be_removed vivax,init_ivm,0,to_be_removed -vivax,init_id,0,to_be_removed vivax,uv,11.4321,to_be_removed -vivax,ud,9.44512,to_be_removed vivax,pvm,0.195768,to_be_removed vivax,rvm,76.8365,to_be_removed vivax,rva,10950,to_be_removed -vivax,rid,3650,to_be_removed vivax,theta0,0.0749886,to_be_removed vivax,theta1,0.0001191,to_be_removed vivax,iv0,1.09629,to_be_removed @@ -103,9 +99,3 @@ vivax,kv,2.00048,to_be_removed vivax,fv0,0.141195,to_be_removed vivax,av,2493.41,to_be_removed vivax,gammav,2.91282,to_be_removed -vivax,fd0,0.007055,to_be_removed -vivax,ad,7993.5,to_be_removed -vivax,gammad,4.8183,to_be_removed -vivax,d1,0.160527,to_be_removed -vivax,id0,1.577533,to_be_removed -vivax,kd,0.476614,to_be_removed diff --git a/data/parasite_parameters.rda b/data/parasite_parameters.rda index 13717f51..542045c7 100644 Binary files a/data/parasite_parameters.rda and b/data/parasite_parameters.rda differ diff --git a/tests/testthat/test-vivax.R b/tests/testthat/test-vivax.R index acd5f089..4b4a150f 100644 --- a/tests/testthat/test-vivax.R +++ b/tests/testthat/test-vivax.R @@ -24,10 +24,15 @@ test_that('Test difference between falciparum and vivax parameter lists', { expect_identical( in_falciparum_not_vivax, - c("init_ib", "rb", "ub", "b0", "b1", "ib0", "kb", # blood immunity parameters - "gamma1") # asymptomatic infected infectivity towards mosquitos parameter + c("du", # human delays + "init_ib", "init_id", ## initial immunities + "rb", "rid", # rates of immune loss + "ub", "ud", # immunity non-boosting periods + "b0", "b1", "ib0", "kb", # blood immunity parameters + "fd0", "ad", "gammad", "d1", "id0", "kd", # immunity to detectability parameters + "gamma1") # asymptomatic infected infectivity towards mosquitos parameter ) - + expect_identical( in_vivax_not_falciparum, c("dpcr_max", "dpcr_min", "kpcr", "apcr50", # human sub-patent state delay @@ -39,3 +44,93 @@ test_that('Test difference between falciparum and vivax parameter lists', { ) ) }) + +## Test anti-parasite immunity subpatent progression functions +test_that('Test anti-parasite immunity function', { + vivax_parameters <- get_parameters(parasite = "vivax", + overrides = list(s_proportion = 0, + d_proportion = 0, + a_proportion = 0, + u_proportion = 1, + t_proportion = 0)) + variables <- create_variables(vivax_parameters) + index <- variables$state$get_index_of("U") + dpcr_min <- vivax_parameters$dpcr_min + dpcr_max <- vivax_parameters$dpcr_max + apcr50 <- vivax_parameters$apcr50 + kpcr <- vivax_parameters$kpcr + + expect_equal(object = anti_parasite_immunity( + dpcr_min, dpcr_max, apcr50, kpcr, + variables$iaa$get_values(index), + variables$iam$get_values(index)), + expected = rep(dpcr_max,100)) + + ## Change initial values of ID, and IDM, check they are the same + variables$iaa <- individual::DoubleVariable$new(1:100) + UAA_durations <- anti_parasite_immunity( + dpcr_min, dpcr_max, apcr50, kpcr, + variables$iaa$get_values(index), + variables$iam$get_values(index)) + + variables$iaa <- individual::DoubleVariable$new(rep(0,100)) + variables$iam <- individual::DoubleVariable$new(1:100) + UAM_durations <- anti_parasite_immunity( + dpcr_min, dpcr_max, apcr50, kpcr, + variables$iaa$get_values(index), + variables$iam$get_values(index)) + + expect_equal(object = UAA_durations, expected = UAM_durations) + + ## Check convergence to min_du at high immunity + variables$iaa <- individual::DoubleVariable$new(rep(1E7,100)) + variables$iam <- individual::DoubleVariable$new(rep(0,100)) + expect_equal(object = anti_parasite_immunity( + dpcr_min, dpcr_max, apcr50, kpcr, + variables$iaa$get_values(index), + variables$iam$get_values(index)), + expected = rep(dpcr_min,100), + tolerance = 1E-2) +}) + +test_that('that vivax patent prevalence rendering works', { + + timestep <- 1 + state <- individual::CategoricalVariable$new( + c('U', 'A', 'D', 'S', 'Tr'), + c('U', 'A', 'D', 'D', 'D', 'S') + ) + birth <- individual::IntegerVariable$new( + -c(3, 4, 5, 1, 11, 6) * 365 + ) + immunity <- individual::DoubleVariable$new(rep(1, 6)) + vivax_parameters <- get_parameters(parasite = "vivax") + renderer <- mock_render(1) + + process <- create_prevalence_renderer( + state, + birth, + immunity, + vivax_parameters, + renderer + ) + + process(timestep) + + mockery::expect_args( + renderer$render_mock(), + 1, + 'n_detect_lm_730_3650', + 2, + timestep + ) + + mockery::expect_args( + renderer$render_mock(), + 2, + 'n_detect_pcr_730_3650', + 3, + timestep + ) + +}) diff --git a/vignettes/Plasmodium_vivax.Rmd b/vignettes/Plasmodium_vivax.Rmd index 5ed608d4..07ba7c9d 100644 --- a/vignettes/Plasmodium_vivax.Rmd +++ b/vignettes/Plasmodium_vivax.Rmd @@ -59,6 +59,12 @@ The *P. falciparum* model tracks four kinds of immunity: immunity to blood stage The *P. vivax* model tracks two kinds of immunity: immunity to clinical infection (acquired and maternal, ICA and ICM) and anti-parasite immunity (acquired and maternal, IAA and IAM). We do not track immunity to blood stage infections, severe immunity or immunity to detectability. +### Anti-parasite immunity + +Anti-parasite immunity has effects in two places. The first is in the separation of PCR-detectable infections from LM-detectable infections (where greater immunity reduces the change of an infection with higher parasitaemia). The second is in the calculation of a PCR-detectable infection (where greater immunity results in a shorter duration before recovery). + +The *P. falciparum* model does not model infections to the sub-patent compartment (U) and has a constant sub-patent infection duration (`du = 110`). + ### Infectivity of LM-detectable infections 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`).