From 3472e76e8800de67f5cdab95cdc0a3b21719811e Mon Sep 17 00:00:00 2001 From: RJSheppard Date: Tue, 17 Sep 2024 18:40:41 +0100 Subject: [PATCH] p.v model allows subpatent infections to occur, so we now calculate lm-detectable and pcr-detectable infections, render these and schedule them. The P.v vignette has been extended to reflect these changes. --- R/human_infection.R | 211 ++++++++++++++------ R/parameters.R | 12 +- R/processes.R | 11 +- man/get_parameters.Rd | 10 +- tests/testthat/test-infection-integration.R | 35 ++-- tests/testthat/test-vivax.R | 118 +++++++++++ vignettes/Plasmodium_vivax.Rmd | 10 +- 7 files changed, 309 insertions(+), 98 deletions(-) diff --git a/R/human_infection.R b/R/human_infection.R index b00e6ec0..33e8614d 100644 --- a/R/human_infection.R +++ b/R/human_infection.R @@ -154,17 +154,19 @@ infection_outcome_process <- function( parameters, prob){ - incidence_renderer( - variables$birth, - renderer, - infected_humans, - 'inc_', - parameters$incidence_rendering_min_ages, - parameters$incidence_rendering_max_ages, - timestep - ) - if (infected_humans$size() > 0) { + + renderer$render('n_infections', infected_humans$size(), timestep) + incidence_renderer( + variables$birth, + renderer, + infected_humans, + 'inc_', + parameters$incidence_rendering_min_ages, + parameters$incidence_rendering_max_ages, + timestep + ) + boost_immunity( variables$ica, infected_humans, @@ -172,6 +174,7 @@ infection_outcome_process <- function( timestep, parameters$uc ) + if(parameters$parasite == "falciparum"){ boost_immunity( variables$id, @@ -180,7 +183,38 @@ infection_outcome_process <- function( timestep, parameters$ud ) + + clinical <- calculate_clinical_infections( + variables, + infected_humans, + parameters, + renderer, + timestep + ) + + treated <- calculate_treated( + variables, + clinical, + parameters, + timestep, + renderer + ) + + update_severe_disease( + timestep, + infected_humans, + variables, + parameters, + renderer + ) + + ## The treated and infected_humans bitsets are re-written so be cautious! + to_D <- treated$not(FALSE)$and(clinical) + to_A <- infected_humans$and(clinical$not(FALSE)) + to_U <- NULL + } else if (parameters$parasite == "vivax"){ + boost_immunity( variables$iaa, infected_humans, @@ -188,45 +222,79 @@ infection_outcome_process <- function( timestep, parameters$ua ) + + ## Only S and U infections are considered in generating lm-det infections + lm_detectable <- calculate_lm_det_infections( + variables, + variables$state$get_index_of(c("S","U"))$and(infected_humans), + parameters, + renderer, + timestep + ) + + # Lm-detectable level infected S and U, and all A infections may receive clinical infections + # There is a different calculation to generate clinical infections, based on current infection level + # LM infections must only pass through the clinical calculation, therefore all "A" infections are included + # "S" and "U" infections must pass through the lm-detectable calculation prior to and in addition to the clinical + # calculation. We therefore consider all "A" infections and only the "S" and "U" infections that are now lm-detectable. + clinical <- calculate_clinical_infections( + variables, + variables$state$get_index_of("A")$and(infected_humans)$or(lm_detectable), + parameters, + renderer, + timestep + ) + + treated <- calculate_treated( + variables, + clinical, + parameters, + timestep, + renderer + ) + + ## The infected_humans,lm_detectable and clinical bitsets are re-written so be cautious! + to_U <- infected_humans$and(lm_detectable$not(F))$and(variables$state$get_index_of(c("S"))) + to_A <- lm_detectable$and(clinical$not(F)) + to_D <- clinical$and(treated$not(F)) } - } - - clinical_infections <- calculate_clinical_infections( - variables, - infected_humans, - parameters, - renderer, - timestep - ) - - if(parameters$parasite == "falciparum"){ - update_severe_disease( - timestep, - infected_humans, - variables, + + schedule_infections( parameters, - renderer + variables, + timestep, + to_D, + to_A, + to_U ) } - - treated <- calculate_treated( - variables, - clinical_infections, - parameters, - timestep, - renderer - ) - - renderer$render('n_infections', infected_humans$size(), timestep) - - schedule_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, - clinical_infections, - treated, - infected_humans, + infections, parameters, + renderer, timestep - ) +) { + + 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)) } #' @title Calculate clinical infections @@ -487,32 +555,27 @@ calculate_successful_treatments <- function( successfully_treated_list } + #' @title Schedule infections #' @description #' Schedule infections in humans after the incubation period -#' @param events a list of all of the model events -#' @param clinical_infections bitset of clinically infected humans -#' @param treated bitset of treated humans -#' @param infections bitset of infected humans #' @param parameters model parameters +#' @param variables a list of all of the model variables +#' @param timestep current timestep +#' @param to_D bitset of humans to move to state D +#' @param to_A bitset of humans to move to state A +#' @param to_U bitset of humans to move to state U #' @noRd schedule_infections <- function( - variables, - clinical_infections, - treated, - infections, parameters, - timestep + variables, + timestep, + to_D, + to_A, + to_U ) { - included <- treated$not(TRUE) - - to_infect <- clinical_infections$and(included) - to_infect_asym <- clinical_infections$copy()$not(TRUE)$and(infections)$and( - included - ) - - if(to_infect$size() > 0) { + if(to_D$size() > 0) { update_infection( variables$state, 'D', @@ -520,18 +583,18 @@ schedule_infections <- function( parameters$cd, variables$progression_rates, 1/parameters$dd, - to_infect + to_D ) } - if(to_infect_asym$size() > 0) { + if(to_A$size() > 0) { if(parameters$parasite == "falciparum"){ # p.f has immunity-determined asymptomatic infectivity update_to_asymptomatic_infection( variables, parameters, timestep, - to_infect_asym + to_A ) } else if (parameters$parasite == "vivax"){ # p.v has constant asymptomatic infectivity @@ -540,9 +603,29 @@ schedule_infections <- function( 'A', variables$infectivity, parameters$ca, - variables$progression_rates, + variables$recovery_rates, 1/parameters$da, - to_infect_asym + to_A + ) + } + } + + if(parameters$parasite == "vivax"){ + # new p.v infections may be subpatent + if(to_U$size() > 0){ + # p.v subpatent recovery rate is immunity dependent + update_infection( + variables$state, + 'U', + variables$infectivity, + parameters$cu, + variables$progression_rates, + 1/anti_parasite_immunity( + parameters$dpcr_min, parameters$dpcr_max, parameters$apcr50, parameters$kpcr, + variables$iaa$get_values(to_U), + variables$iam$get_values(to_U) + ), + to_U ) } } diff --git a/R/parameters.R b/R/parameters.R index 6fe148eb..79ef7629 100644 --- a/R/parameters.R +++ b/R/parameters.R @@ -37,10 +37,10 @@ #' * da - the delay for humans to move from state A to U; default = 195 #' * du - the delay for humans to move from state U to S (p.f only); default = 110 #' -#' duration of pcr detectable infection: du (p.v only): +#' duration of pcr detectable infections: du (p.v only): #' -#' * dpcr_max - Maximum duration of subpatent infection: default = 53.69 -#' * dpcr_min - Minimum duration of subpatent infection: default = 10 +#' * dpcr_max - Maximum duration of PCR-detectable infections: default = 53.69 +#' * dpcr_min - Minimum duration of PCR-detectable infections: default = 10 #' * kpcr - Shape parameter: default = 4.021 #' * apcr50 - Scale parameter: default = 9.8 #' @@ -112,7 +112,7 @@ #' * id0 - scale parameter; default = 1.577533 #' * kd - shape parameter; default = 0.476614 #' -#' probability of patent infection due to anti-parasite immunity (p.v only): +#' probability of light-microscopy detectable infections due to anti-parasite immunity (p.v only): #' #' * philm_max - maximum probability due to no immunity; default = 0.9329 #' * philm_min - maximum reduction due to immunity; default = 0.0131 @@ -141,7 +141,7 @@ #' * cd - infectivity of clinically diseased humans towards mosquitoes; default = 0.068 #' * ca - infectivity of asymptomatic humans towards mosquitoes (p.v only); default = 0.1 #' * gamma1 - parameter for infectivity of asymptomatic humans; default = 1.82425 -#' * cu - infectivity of sub-patent infection; default = 0.0062 +#' * cu - infectivity of subpatent infection; default = 0.0062 #' * ct - infectivity of treated infection; default = 0.021896 #' #' mosquito fixed state transitions (including mortality): @@ -344,7 +344,7 @@ get_parameters <- function(overrides = list(), parasite = "falciparum") { # maternal immunity parameters # probability of pre-erythrocytic infection/blood immunity # probability of asymptomatic detection (p.f only) - # probability of patent infection (due to anti-parasite immunity, p.v only) + # probability of light-microscopy detectable infection (due to anti-parasite immunity, p.v only) # probability of clinical infection # probability of severe infection (p.f only) # infectivity towards mosquitos diff --git a/R/processes.R b/R/processes.R index 72513b0c..b67dc559 100644 --- a/R/processes.R +++ b/R/processes.R @@ -47,7 +47,7 @@ create_processes <- function( immunity_process = create_exponential_decay_process(variables$ica, parameters$rc) ) - + if(parameters$parasite == "falciparum"){ processes <- c( processes, @@ -61,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) ) } diff --git a/man/get_parameters.Rd b/man/get_parameters.Rd index cf58d89a..76e2db6a 100644 --- a/man/get_parameters.Rd +++ b/man/get_parameters.Rd @@ -38,10 +38,10 @@ human fixed state transitions: \item du - the delay for humans to move from state U to S (p.f only); default = 110 } -duration of pcr detectable infection: du (p.v only): +duration of pcr detectable infections: du (p.v only): \itemize{ -\item dpcr_max - Maximum duration of subpatent infection: default = 53.69 -\item dpcr_min - Minimum duration of subpatent infection: default = 10 +\item dpcr_max - Maximum duration of PCR-detectable infections: default = 53.69 +\item dpcr_min - Minimum duration of PCR-detectable infections: default = 10 \item kpcr - Shape parameter: default = 4.021 \item apcr50 - Scale parameter: default = 9.8 } @@ -124,7 +124,7 @@ probability of detection by light-microscopy when asymptomatic (p.f only): \item kd - shape parameter; default = 0.476614 } -probability of patent infection due to anti-parasite immunity (p.v only): +probability of light-microscopy detectable infections due to anti-parasite immunity (p.v only): \itemize{ \item philm_max - maximum probability due to no immunity; default = 0.9329 \item philm_min - maximum reduction due to immunity; default = 0.0131 @@ -156,7 +156,7 @@ infectivity towards mosquitoes: \item cd - infectivity of clinically diseased humans towards mosquitoes; default = 0.068 \item ca - infectivity of asymptomatic humans towards mosquitoes (p.v only); default = 0.1 \item gamma1 - parameter for infectivity of asymptomatic humans; default = 1.82425 -\item cu - infectivity of sub-patent infection; default = 0.0062 +\item cu - infectivity of subpatent infection; default = 0.0062 \item ct - infectivity of treated infection; default = 0.021896 } diff --git a/tests/testthat/test-infection-integration.R b/tests/testthat/test-infection-integration.R index eb4eabb6..ad83c35f 100644 --- a/tests/testthat/test-infection-integration.R +++ b/tests/testthat/test-infection-integration.R @@ -95,8 +95,11 @@ test_that('simulate_infection integrates different types of infection and schedu treated <- individual::Bitset$new(population)$insert(3) treated_mock <- mockery::mock(treated) schedule_mock <- mockery::mock() - - + + to_D <- treated$not(FALSE)$and(clinical) + to_A <- infected$and(clinical$not(FALSE)) + to_U <- NULL + mockery::stub(infection_outcome_process, 'incidence_renderer', mockery::mock()) mockery::stub(infection_outcome_process, 'boost_immunity', boost_immunity_mock) mockery::stub(infection_outcome_process, 'calculate_clinical_infections', clinical_infection_mock) @@ -111,8 +114,7 @@ test_that('simulate_infection integrates different types of infection and schedu renderer, parameters, prob) - - + mockery::expect_args( boost_immunity_mock, 1, @@ -143,6 +145,7 @@ test_that('simulate_infection integrates different types of infection and schedu renderer ) + mockery::mock_args(treated_mock)[[1]][[2]]$to_vector() mockery::expect_args( treated_mock, 1, @@ -156,12 +159,12 @@ test_that('simulate_infection integrates different types of infection and schedu mockery::expect_args( schedule_mock, 1, - variables, - clinical, - treated, - infected, parameters, - timestep + variables, + timestep, + to_D, + to_A, + to_U ) }) @@ -647,16 +650,20 @@ test_that('schedule_infections correctly schedules new infections', { infection_mock <- mockery::mock() asymp_mock <- mockery::mock() + to_D <- treated$not(FALSE)$and(clinical_infections) + to_A <- clinical_infections$not(FALSE)$and(infections) + to_U <- NULL + mockery::stub(schedule_infections, 'update_infection', infection_mock) mockery::stub(schedule_infections, 'update_to_asymptomatic_infection', asymp_mock) schedule_infections( - variables, - clinical_infections, - treated, - infections, parameters, - 42 + variables, + 42, + to_D, + to_A, + to_U ) actual_infected <- mockery::mock_args(infection_mock)[[1]][[7]]$to_vector() diff --git a/tests/testthat/test-vivax.R b/tests/testthat/test-vivax.R index b1310584..8fc77956 100644 --- a/tests/testthat/test-vivax.R +++ b/tests/testthat/test-vivax.R @@ -137,3 +137,121 @@ test_that('that vivax patent prevalence rendering works', { ) }) + +test_that('Test default vivax incidence rendering works', { + + timestep <- 0 + year <- 365 + birth <- individual::IntegerVariable$new( + -c(2, 5, 10, 11) * year + ) + vivax_parameters <- get_parameters( + parasite = "vivax") + + renderer <- mock_render(1) + incidence_renderer( + birth, + renderer, + individual::Bitset$new(4)$insert(c(1, 2, 4)), + 'inc_patent_', + c(0, 2) * year, + c(5, 10) * year, + timestep + ) + + incidence_probability_renderer( + birth, + renderer, + individual::Bitset$new(4)$insert(seq(4)), + c(.1, .2, .3, .4), + 'inc_patent_', + c(0, 2) * year, + c(5, 10) * year, + timestep + ) + + mockery::expect_args( + renderer$render_mock(), + 1, + 'n_inc_patent_0_1825', + 2, + timestep + ) + + mockery::expect_args( + renderer$render_mock(), + 2, + 'n_inc_patent_730_3650', + 2, + timestep + ) + + mockery::expect_args( + renderer$render_mock(), + 3, + 'p_inc_patent_0_1825', + 0.3, + timestep + ) + + mockery::expect_args( + renderer$render_mock(), + 4, + 'p_inc_patent_730_3650', + .6, + timestep + ) +}) + +test_that('vivax schedule_infections correctly schedules new infections', { + parameters <- get_parameters(list(human_population = 20), parasite = "vivax") + variables <- create_variables(parameters) + + variables$state <- individual::CategoricalVariable$new( + c('U', 'A', 'D', 'S', 'Tr'), + rep(c('S','U','A','D','Tr'), times = 4) + ) + + infections <- individual::Bitset$new(20)$insert(1:20) + lm_detectable <- individual::Bitset$new(20)$insert(6:20)$and(variables$state$get_index_of(c("S", "U"))) + clinical <- individual::Bitset$new(20)$insert(11:20)$and(variables$state$get_index_of(c("A"))$or(lm_detectable)) + treated <- individual::Bitset$new(20)$insert(16:20)$and(clinical) + # Only S can be a subpatent infection (1) + # Only S and U can be a patent infection (6, 7) + # S, U and A can be clinical infections (11, 12, 13), but the model re-infects everyone + # Treated only looks at new clinical infections (from SAU, not from D) + to_U <- infections$and(lm_detectable$not(F))$and(variables$state$get_index_of(c("S"))) + to_A <- lm_detectable$and(clinical$not(F)) + to_D <- clinical$and(treated$not(F)) + + infection_mock <- mockery::mock() + mockery::stub(schedule_infections, 'update_infection', infection_mock) + + schedule_infections( + parameters, + variables, + timestep = 42, + to_D, + to_A, + to_U + ) + + actual_infected <- mockery::mock_args(infection_mock)[[1]][[7]]$to_vector() + actual_asymp_infected <- mockery::mock_args(infection_mock)[[2]][[7]]$to_vector() + actual_subpatent_infected <- mockery::mock_args(infection_mock)[[3]][[7]]$to_vector() + + expect_equal( + actual_infected, + c(11:13) + ) + + expect_equal( + actual_asymp_infected, + c(6, 7) + ) + + expect_equal( + actual_subpatent_infected, + c(1) + ) +}) diff --git a/vignettes/Plasmodium_vivax.Rmd b/vignettes/Plasmodium_vivax.Rmd index 07ba7c9d..87c56671 100644 --- a/vignettes/Plasmodium_vivax.Rmd +++ b/vignettes/Plasmodium_vivax.Rmd @@ -53,17 +53,17 @@ The *P. falciparum* model has five human disease compartments: susceptible (S), The *P. vivax* model follows a similar structure to the *P. falciparum* model, and also has five human disease compartments. However, the human disease states modeled explicitly focus on parasite density and detectability, such that we have: susceptible (S), clinical disease (D), **light-microscopy detectable infection (A)**, **PCR detectable infection (U)**, and treated (Tr). +### New Infections + +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 The *P. falciparum* model tracks four kinds of immunity: immunity to blood stage infection (IB), clinical disease (acquired and maternal, ICA and ICM), to severe disease (acquired and maternal, IVA and IVM), and to detectability (ID). 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`). +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