Skip to content

Commit

Permalink
Merge pull request #127 from mrc-ide/dev
Browse files Browse the repository at this point in the history 
Dev
  • Loading branch information
@giovannic
giovannic authored Jul 29, 2021
2 parents b9b0589 + 32737b5 commit cf3e479
Show file tree
Hide file tree
Showing 18 changed files with 549 additions and 358 deletions.
2 changes: 1 addition & 1 deletion DESCRIPTION
View file
Original file line number Diff line number Diff line change
@@ -1,6 +1,6 @@
Package: malariasimulation
Title: An individual based model for malaria
Version: 1.0.1
Version: 1.0.2
Authors@R: c(
person(
given = "Giovanni",
Expand Down
103 changes: 37 additions & 66 deletions R/biting_process.R
View file
Original file line number Diff line number Diff line change
Expand Up @@ -8,6 +8,8 @@
#' @param variables a list of all of the model variables
#' @param events a list of all of the model events
#' @param parameters model pararmeters
#' @param lagged_infectivity a LaggedValue class with historical sums of infectivity
#' @param lagged_eir a LaggedValue class with historical EIRs
#' @noRd
create_biting_process <- function(
renderer,
Expand All @@ -16,7 +18,7 @@ create_biting_process <- function(
variables,
events,
parameters,
lagged_foim,
lagged_infectivity,
lagged_eir
) {
function(timestep) {
Expand All @@ -32,7 +34,7 @@ create_biting_process <- function(
age,
parameters,
timestep,
lagged_foim,
lagged_infectivity,
lagged_eir
)

Expand All @@ -58,7 +60,7 @@ simulate_bites <- function(
age,
parameters,
timestep,
lagged_foim,
lagged_infectivity,
lagged_eir
) {
bitten_humans <- individual::Bitset$new(parameters$human_population)
Expand Down Expand Up @@ -95,30 +97,8 @@ simulate_bites <- function(
W <- average_p_successful(p_bitten$prob_bitten_survives, .pi, Q0)
Z <- average_p_repelled(p_bitten$prob_repelled, .pi, Q0)
f <- blood_meal_rate(s_i, Z, parameters)

lambda <- .effective_biting_rates(
s_i,
psi,
.pi,
zeta,
p_bitten,
W,
Z,
f,
parameters
)

renderer$render(
paste0('lambda_', parameters$species[[s_i]]),
sum(lambda),
timestep
)

renderer$render(
paste0('normal_lambda_', parameters$species[[s_i]]),
sum(lambda) / mean(psi),
timestep
)
a <- .human_blood_meal_rate(f, s_i, W, parameters)
lambda <- effective_biting_rates(a, .pi, p_bitten)

if (parameters$individual_mosquitoes) {
species_index <- variables$species$get_index_of(
Expand All @@ -136,18 +116,19 @@ simulate_bites <- function(
n_infectious <- calculate_infectious_compartmental(solver_states)
}

lagged_eir[[s_i]]$save(n_infectious * sum(lambda), timestep)
lagged_eir[[s_i]]$save(n_infectious * a, timestep)
species_eir <- lagged_eir[[s_i]]$get(timestep - parameters$de)
EIR <- EIR + species_eir / sum(psi)
n_bites <- rpois(1, species_eir)
EIR <- EIR + species_eir
n_bites <- rpois(1, species_eir * mean(psi))
if (n_bites > 0) {
bitten_humans$insert(
fast_weighted_sample(n_bites, lambda)
)
}

foim <- calculate_foim(human_infectivity, lambda, psi)
lagged_foim$save(foim, timestep)
infectivity <- lagged_infectivity$get(timestep - parameters$delay_gam)
lagged_infectivity$save(sum(human_infectivity * .pi), timestep)
foim <- calculate_foim(a, infectivity)
renderer$render(paste0('FOIM_', s_i), foim, timestep)
mu <- death_rate(f, W, Z, s_i, parameters)
renderer$render(paste0('mu_', s_i), mu, timestep)
Expand All @@ -161,7 +142,7 @@ simulate_bites <- function(

biting_effects_individual(
variables,
lagged_foim$get(timestep - parameters$delay_gam),
foim,
events,
s_i,
susceptible_species_index,
Expand All @@ -174,7 +155,7 @@ simulate_bites <- function(
adult_mosquito_model_update(
models[[s_i]],
mu,
lagged_foim$get(timestep - parameters$delay_gam),
foim,
solver_states[[ADULT_ODE_INDICES['Sm']]],
f
)
Expand All @@ -192,37 +173,15 @@ simulate_bites <- function(
# =================

calculate_eir <- function(species, solvers, variables, parameters, timestep) {
lambda <- effective_biting_rates(species, variables, parameters, timestep)
a <- human_blood_meal_rate(species, variables, parameters, timestep)
infectious <- calculate_infectious(species, solvers, variables, parameters)
infectious * sum(lambda)
}

effective_biting_rates <- function(species, variables, parameters, timestep) {
age <- get_age(variables$birth$get_values(), timestep)
psi <- unique_biting_rate(age, parameters)
zeta <- variables$zeta$get_values()
p_bitten <- prob_bitten(timestep, variables, species, parameters)
.pi <- human_pi(zeta, psi)
Q0 <- parameters$Q0[[species]]
W <- average_p_successful(p_bitten$prob_bitten_survives, .pi, Q0)
Z <- average_p_repelled(p_bitten$prob_repelled, .pi, Q0)
f <- blood_meal_rate(species, Z, parameters)
.effective_biting_rates(species, psi, .pi, zeta, p_bitten, W, Z, f, parameters)
infectious * a * mean(psi)
}

.effective_biting_rates <- function(
species,
psi,
.pi,
zeta,
p_bitten,
W,
Z,
f,
parameters
) {
a <- human_blood_meal_rate(f, species, W, parameters)
a * intervention_coefficient(p_bitten) * zeta * psi
effective_biting_rates <- function(a, .pi, p_bitten) {
a * .pi * p_bitten$prob_bitten / sum(.pi * p_bitten$prob_bitten_survives)
}

calculate_infectious <- function(species, solvers, variables, parameters) {
Expand Down Expand Up @@ -274,7 +233,20 @@ blood_meal_rate <- function(v, z, parameters) {
1 / (interrupted_foraging_time + gonotrophic_cycle)
}

human_blood_meal_rate <- function(f, v, W, parameters) {
human_blood_meal_rate <- function(species, variables, parameters, timestep) {
age <- get_age(variables$birth$get_values(), timestep)
psi <- unique_biting_rate(age, parameters)
zeta <- variables$zeta$get_values()
p_bitten <- prob_bitten(timestep, variables, species, parameters)
.pi <- human_pi(zeta, psi)
Q0 <- parameters$Q0[[species]]
W <- average_p_successful(p_bitten$prob_bitten_survives, .pi, Q0)
Z <- average_p_repelled(p_bitten$prob_repelled, .pi, Q0)
f <- blood_meal_rate(species, Z, parameters)
.human_blood_meal_rate(f, species, W, parameters)
}

.human_blood_meal_rate <- function(f, v, W, parameters) {
Q <- 1 - (1 - parameters$Q0[[v]]) / W
Q * f
}
Expand All @@ -294,10 +266,9 @@ unique_biting_rate <- function(age, parameters) {

#' @title Calculate the force of infection towards mosquitoes
#'
#' @param human_infectivity a vector of infectivities for each human
#' @param lambda a vector of biting rates for each human
#' @param psi a vector of age-based relative biting rates for each human
#' @param a human blood meal rate
#' @param infectivity_sum a sum of infectivity weighted by relative biting rate
#' @noRd
calculate_foim <- function(human_infectivity, lambda, psi) {
sum(human_infectivity * lambda) / mean(psi)
calculate_foim <- function(a, infectivity_sum) {
a * infectivity_sum
}
25 changes: 6 additions & 19 deletions R/parameters.R
View file
Original file line number Diff line number Diff line change
Expand Up @@ -146,17 +146,10 @@
#' please set vector control strategies using `set_betnets` and `set_spraying`
#'
#' * bednets - boolean for if bednets are enabled; default = FALSE
#' * rn - probability mosquito is repelled by the bednet; default = 0.56
#' * rnm - minimum probability mosquito is repelled by the bednet ; default = 0.24
#' * dn0 - probability killed by the bednet; default = 0.533
#' * phi_bednets - proportion of bites taken in bed; default = 0.85
#' * k0 - proportion of females bloodfed with no net; default = 0.699
#' * spraying - boolean for if indoor spraying is enabled; default = FALSE
#' * rs - probability repelled by indoor spraying; default = 0.2
#' * phi_indoors - proportion of bites taken indoors; default = 0.97
#' * phi_bednets - proportion of bites taken in bed; default = 0.89
#' * endophily - proportion of mosquitoes resting indoors after feeding with no
#' intervention; default = 0.813
#' * gammas - the half-life of spraying efficacy (timesteps); default = 91.25
#' * gamman - the half-life of bednet efficacy (timesteps); default = 963.6
#' * phi_indoors - proportion of bites taken indoors; default = 0.90
#'
#' treatment parameters:
#' please set treatment parameters with the convenience functions in
Expand Down Expand Up @@ -337,20 +330,14 @@ get_parameters <- function(overrides = list()) {
species_proportions = 1,
blood_meal_rates = 1/3,
Q0 = .92,
endophily = .813,
foraging_time = .69,
# bed nets
bednets = FALSE,
rn = .56,
rnm = .24,
dn0 = .533,
phi_bednets = .89,
gamman = 2.64 * 365,
phi_bednets = .85,
k0 = .699,
# indoor spraying
spraying = FALSE,
rs = .2,
phi_indoors = .97,
gammas = .25 * 365,
phi_indoors = .90,
# treatment
drug_efficacy = numeric(0),
drug_rel_c = numeric(0),
Expand Down
22 changes: 14 additions & 8 deletions R/processes.R
View file
Original file line number Diff line number Diff line change
Expand Up @@ -54,12 +54,12 @@ create_processes <- function(
# schedule infections for humans and set last_boosted_*
# move mosquitoes into incubating state
# kill mosquitoes caught in vector control
init_eir <- lapply(
lagged_eir <- lapply(
seq_along(parameters$species),
function(species) {
LaggedValue$new(
parameters$de + 1,
calculate_eir(
max_lag = parameters$de + 1,
default = calculate_eir(
species,
solvers,
variables,
Expand All @@ -70,6 +70,15 @@ create_processes <- function(
}
)

age <- get_age(variables$birth$get_values(), 0)
psi <- unique_biting_rate(age, parameters)
.pi <- human_pi(psi, variables$zeta$get_values())
init_infectivity <- sum(.pi * variables$infectivity$get_values())
lagged_infectivity <- LaggedValue$new(
max_lag = parameters$delay_gam + 2,
default = init_infectivity
)

processes <- c(
processes,
create_biting_process(
Expand All @@ -79,11 +88,8 @@ create_processes <- function(
variables,
events,
parameters,
LaggedValue$new(
parameters$delay_gam + 2,
parameters$init_foim
), # lagged FOIM
init_eir # lagged EIR
lagged_infectivity,
lagged_eir
),
create_mortality_process(variables, events, renderer, parameters),
create_progression_process(
Expand Down
6 changes: 2 additions & 4 deletions R/variables.R
View file
Original file line number Diff line number Diff line change
Expand Up @@ -48,8 +48,6 @@ create_variables <- function(parameters) {
initial_age <- round(rexp(size, rate=1/parameters$average_age))

if (parameters$enable_heterogeneity) {
#zeta_norm <- rnorm(size)
#het_nodes <- gaussian_quadrature(parameters$n_heterogeneity_groups)
quads <- statmod::gauss.quad.prob(
parameters$n_heterogeneity_groups,
dist='normal'
Expand Down Expand Up @@ -209,8 +207,8 @@ create_variables <- function(parameters) {
tbv_vaccinated <- individual::DoubleVariable$new(rep(-1, size))

# Init vector controls
net_time <- individual::DoubleVariable$new(rep(-1, size))
spray_time <- individual::DoubleVariable$new(rep(-1, size))
net_time <- individual::IntegerVariable$new(rep(-1, size))
spray_time <- individual::IntegerVariable$new(rep(-1, size))

variables <- list(
state = state,
Expand Down
Loading

0 comments on commit cf3e479

Please sign in to comment.