Bait for mosquitos. Increases daily mortality in mosquitos. Modelled by @KeithJF82

Mosquitoes can also be in a "latent infection" state, see here.

Steps:

• Add the extra state an processes for mosquitoes
• See the effect this has on the equilibrium of the model

TODO:

• Change average population age to 30
• Run model with +1k humans
• Output demography over time
• Check immunity over time disaggregated by age
• Check EIR vs prevalence/incidence
• Check states for over 200 years

Alternative vaccine model. Ask @JDChallenger and Thomas.

Correlations are hard to understand. We should provide outputs so that researchers can easily inspect the populations that are being drawn using the correlation parameters.

What would these outputs look like? Perhaps @pwinskill has a basic idea?

• Graphs (x and y axes, groups...)

NOTE: it's been recommended to keep graphs out of our core R code. So we'll probably make a vignette to demonstrate how researchers can plot.

Initial convo and example

Steps to recreate:

1. install dev version
2. render clinical incidence

params <- get_parameters(list(
  clinical_incidence_rendering_min_ages = 0,
  clinical_incidence_rendering_max_ages = 5 * 365
))
output <- malariasimulation::run_simulation(10, params)

3. Print output names

Expected:

Output should include n_0_1825

Actual:

There is no denominator

Related to bitset$not

Affects mosquito hatch rate. Modelled by @EllieSherrardSmith

Each state and parameter should have a nice description, like Bob's description of the "I" state here

We want to be able to correlate MDA/RTS,S and TBV.

This is documented in Griffin et al. 2010 Protocol S2 "2.4 Correlations between Interventions"

RTS,S vaccinations are administered to a proportion of the population within an age range. That is primary vaccination coverage. This is how you would vaccinate 80% of 5 to 17 month olds:

year <- 365
month <- 30

rtssparams <- set_rtss(
  rtssparams,
  start = start_timestep,
  end = end_timestep,
  frequency = year,
  min_ages = 5 * month,
  max_ages = 17 * month,
  boosters = 18 * month,
  coverage = 0.8
)

Currently, all vaccinated individuals will receive booster vaccinations, boosters timesteps after their primary vaccination.

We would like to specify a coverage for each booster dose, which will select a random subset of the vaccinated individuals to receive a booster. It could be parameterised something like this:

year <- 365
month <- 30

rtssparams <- set_rtss(
  rtssparams,
  start = start_timestep,
  end = end_timestep,
  frequency = year,
  min_ages = 5 * month,
  max_ages = 17 * month,
  boosters = 18 * month,
  coverage = 0.8,
  boosters_coverage = 0.9
)

Notes to developer:

• General parameterisation is in R/parameters.R. Vaccine strategy parameterisation is in R.vaccine_parameters.R
• Vaccination code is in R/rtss_processes.R
• Test code is in tests/testthat/test-rtss.R
• The main implementation changes will probably be in create_rtss_booster_listener

Vaccine modelling feedback

• Model LLIN
• Configure Jamie's model to apply LLIN individually
• Compare outputs to Jamie's model

Currently bednet paramters rn and dn are are constant for the whole of the simulation:

# ...
# set bednet profile
bednetparams <- get_parameters(list(
    rn = .56,
    rnm = .24,
    dn0 = .533,
    gamman = 2.64 * 365
))

# set bednet strategy
bednetparams <- set_bednets(
  bednetparams,
  timesteps = 4:7 * year,
  coverages = rep(.8, 4),
  retention = 5 * year
)
# ...

But in actuality mosquitoes become resistant to bednets. We would like to model changes to dn/rn. Questions:

• are the dn/rn changes normally modelled as a curve or discrete changes?
• And are these changes independent of the type of bednet? See #54

• Model IRS
• Configure Jamie's model to apply IRS individually
• Compare outputs to Jamie's model

The model currently models drug prophylaxis using a weibull distribution. But we have published some models using the gamma distribution. We would like to be able to use either. This would involve:

• Changing set_drugs to allow the parameterisation of different distributions for prophylaxis
• Changing calculate_treated to incorporate the correct model and outputs into the overall model

Example paper: https://bmcmedicine.biomedcentral.com/articles/10.1186/s12916-020-1494-3

Building vignettes is a typical part of the R building and dev process, so it should be tractable in a few mins max.

Todo:

• move the validation vignette to a separate repository
• set up triggers for releases to build and report the validation vignette
• make intervention vignettes much smaller so that they execute in < 5 mins

Currently we can only model one type of bednet:

# ...
# set bednet profile
bednetparams <- get_parameters(list(
    rn = .56,
    rnm = .24,
    dn0 = .533,
    gamman = 2.64 * 365
))

# set bednet strategy
bednetparams <- set_bednets(
  bednetparams,
  timesteps = 4:7 * year,
  coverages = rep(.8, 4),
  retention = 5 * year
)
# ...

But we would like to model different kinds of nets at different times:

# ...
# set bednet profiles
bednetparams <- get_parameters()
bednetparams <- set_bednet_profiles(bednetparams, c(BEDNET_PROFILE, NEW_IMPROVED_BEDNET_PROFILE))

# set bednet strategy
bednetparams <- set_bednets(
  bednetparams,
  timesteps = 4:7 * year,
  coverages = rep(.8, 4),
  bednet_profiles = c(1, 1, 2, 2),
  retention = 5 * year
)
# ...

Currently we can only define one cohort for each intervention. We would like to model many...

Interface TBD...

The parameters dn0, rn, rs, gammas, gamman need to be time varying. (Are ll of these required?)

Todo:

• Update set_bednets and set_spraying to update and validate these time varying parameters (vector_control_parameters.R)
• Update prob_bitten to update repelling and survival probabilities based on timestep
• Update the vector control vignette to demonstrate its use

Hi Giovanni,
I've run the vaccines vignette, implementation of the transmission-blocking vaccine isn't affecting prevalence or incidence of malaria in the population (easier to ascertain by increasing population size or turning off seasonality). The RTS,S vaccine seems fine. I'm not quite sure what's causing the issue, but the version of set_tbv() we were using last year contained the arguments 'start', 'end' & 'frequency', instead of 'timesteps'. So I wondered whether this is causing the issue. Happy to discuss in more detail if useful, and there's no huge rush to sort out. Thanks!

• SMC
• Configure Jamie's model to apply SMC individually
• Compare outputs to Jamie's model

See here

• Model LLIN
• Configure Jamie's model to apply LLIN individually
• Compare outputs to Jamie's model

Tx is not needed for this milestone and introduces unnecessary complexity

As a user, I'd like to view the species of each mosquito at each timestep in the simulation.

Related to individual's issue

You'd want to add something like:

parameters <- get_parameters(list(render_individual_species = TRUE))
output <- run_simulation(timesteps = 100, parameters = parameters)
output$species
# > character matrix of size parameters$mosquito_limit x timesteps with the species of each mosquito

Speed up execution with a compartmental model for mosquitos

Currently we can only model one type of indoor spraying

# ...
# set spraying profile
sprayingparams <- get_parameters(list(
    rs = .2,
    gammas = .5 * 365
))

# set spraying strategy
sprayingparams <- set_spraying(
  sprayingparams,
  timesteps = 4:7 * year,
  coverages = rep(.8, 4)
)
# ...

We would like to add more parameters to the spraying profile. And parameterise strategies based on the different profiles:

# ...
# set spraying profiles
sprayingparams <- get_parameters()
sprayingparams <- set_spraying_profiles(sprayingparams, c(SPRAYING_PROFILE_1, SPRAYING_PROFILE_2))

# set bednet strategy
sprayingparams <- set_spraying(
  sprayingparams,
  timesteps = 4:7 * year,
  coverages = rep(.8, 4),
  spraying_profiles = c(1, 1, 2, 2)
)
# ...

New spraying profiles will be modelled roughly like this (to be confirmed!):

vector_control_decay_sprays_rs_ss <- function(since_spray, 
                                              irs_decay_mort1, irs_decay_mort2,
                                              irs_decay_succ1, irs_decay_succ2,
                                              irs_k0, # the maximum successful feeding in the absence of indoor interventions
                                              irs_decay_det1, irs_decay_det2) {
  
  mort = 1 / (1 + exp(-(irs_decay_mort1 + irs_decay_mort2 * since_spray)))
  succ = irs_k0 / (1 + exp(-(irs_decay_succ1 + irs_decay_succ2 * since_spray)))
  det =  1 / (1 + exp(-(irs_decay_det1  +  irs_decay_det2 * since_spray)))
  rep = 1 - mort - succ
  
  kp_s = (1 - det) * succ
  jp_s = (1 - det) * rep + det
  lp_s =  (1 - det) * mort
  
  rs = (1 - kp_s/irs_k0)*(jp_s/(lp_s + jp_s))
  ds = (1 - kp_s/irs_k0)*(lp_s/(lp_s + jp_s))  
  ss = 1 - rs - ds
}

This is a follow up from #48

The current model does this by:

• running a simulation
• executing this function, (with side-effects on the simulation)
• and finding the root of the function outputs here

We would like to create behaviour like:

simparams <- calibrate_model_for_prevalence(simparams, pfpr2_10)

Which would run the simulation in a similar way and return simparams with an updated total_M to recreate the baseline prevalence.

Developer notes:

• Other paramterisation/calibration functions are currently in R/parameters.R.
• It would be appropriate to add a test in tests/testthat/test-biology.R to make sure that it indeed produces a pfpr2-10 that is close to the desired one.

sample_mothers and which are some of the biggest bottlenecks to simulation at the moment (profile attached). They appear heavily in create_mortality_process. This process kills humans, replaces them with newborns and assigns a maternal immunity from the same heterogeneity group.

Translating this process to C++ should give a big speedup. See Individual documentation for how to do this.

What's a good paper for RTS,S?

• Model RTS,S
• Configure Jamie's model to apply RTS,S individually
• Compare outputs to Jamie's model

This could bring easy speedups

Add treatment to the model (ACTs)

• model T and Prophylaxis states for humans
• compare to equilibrium solution

Here's an example drat repo, maintained with gh actions

Currently, set_clinical_treatment() allows you to specify the probability of seeking treatment ft, the drugs drugs and their realtive coverages coverages.

It would be advantageous for these inputs to be allowed to vary over time. This could be implemented in a similar way to other interventions by allowing the user to set a timesteps parameter to specify when changes occur. Perhaps, ft, and coverages could then simply be list inputs each of length == length(timesteps) (and drugs the full compliment of drugs over all timesteps?).

Which paper is good for MDA?

• Model MDA
• Configure Jamie's model to apply MDA individually
• Compare outputs to Jamie's model

Model won't run with multiple species parameters plugged into set_species(). But it works with just one species!

Test case below; gets this error:"Error in state_values[species_range] <- rep(c("Sm", "Pm", "Im"), times = mosquito_counts[ADULT_ODE_INDICES]) : replacement has length zero"

Coming from here? https://github.com/mrc-ide/malariasimulation/blob/master/R/variables.R#L270

testparams <- malariasimulation::get_parameters()
testparams <- set_species(testparams,species=list(gamb_params, arab_params, fun_params),proportions=c(1,0,0))
testoutput <- run_simulation(100, testparams)

Make a vignette entry describing the basic assumptions which differ, including:

• Individual mosquitoes
• Prophylaxis compartment vs curve
• Reinfection of Asymptomatics
• Maternal immunity parent age (15 - 35 vs 20)
• Heterogeneity (continuous vs quadrature)
• Lag of FOI vs lag of EIR

We currently render n_2_10 and prev_2_10 - which confusingly is the total number of individuals between 2 and 10 years of age in the D state.

Research questions require the prevalence and/or incidence of clinical and/or severe disease within arbitrary age ranges. They could be parameterised like this:

simparams <- get_parameters(
  prevalence_rendering_min_ages = c(0, 2) * year,
  prevalence_rendering_max_ages = c(5, 10) * year,
  severe_incidence_rendering_min_ages = 10 * year,
  severe_incidence_rendering_max_ages = 15 * year
)

That would produce extra columns in the output dataframe with names like pv_0_1825, pv_730_3650, inc_3650_5400.

Notes to the developer:

• There is a prevalence rendering process in R/render.R
• To easily render incidence you will have to write an infection listener in R/processes.R -> create_event_based_processes

parameterise_mosquito_equilibrium currently asks for a limit_grace parameter.

parameters <- malariasimulation::get_parameters()
parameters <- malariasimulation::parameterise_mosquito_equilibrium(parameters, EIR=50, limit_grace=10)
# we write limit_grace=10 to allocate memory for 10 * total_M mosquitoes in the simulation

Researchers shouldn't have to worry about limit_grace, mosquito_limit or any memory allocation issues. We can use the seasonality parameters to figure that out

TBC: Requires steps to recreate

Error in if (n_bites > 0) { : missing value where TRUE/FALSE needed
Calls: lapply ... <Anonymous> -> execute_any_process -> p -> simulate_bites
In addition: Warning message:
In rpois(1, species_eir * mean(psi)) : NAs produced
Execution halted

Mosquitoes currently have a fixed incubation period of parameters$dem. But it should be drawn from a log_uniform distribution

Any longer and checking equilibrium solutions will be prohibitively slow. This will involve:

• Profiling the longest running processes
• Changing simple processes into primitives contingent on this
• Translate long process functions into C++
• Optimise individual_index_t

Their interactions need to be a bit clearer.

• Rename MosquitoModel -> AquaticMosquitoModel
• Add docs for MosquitoModel, AdultMosquitoModel and Solver classes
• Rename create_ode -> create_equations. They're not technically ODEs
• Add docs for create_equations

To make it easy to check the model against the equilibrium solution we should have the following:

• easy paramterisation of model
• create a switch for seasonality
• translate (rename) new parameters to legacy ones
• estimate the EIR for a particular M
• create output plots for 100 years
• mark equilibrium on the plots
• create an RMarkdown for recreation

You can create the following output from VectorControl.Rmd, with the following parameter changes:

sim_length <- 20 * year
human_population <- 30000
starting_EIR <- 0.1
repetitions <- 1

Does this stabilise after some time or is something inconsistent?

We were discussing parameterising the prophylaxis curves for the different drug types in the model. We noticed that in Okell et al. 2014 Nat Comms, the form of the prophylaxis curve is a Weibull survival function: P=exp((-t/lambda)^w), where lambda is the scale parameter and w is the shape parameter. In the malariasimulation model (human_infection.R), the prophylaxis curve is the Weibull pdf (dweibull) instead of the survival curve. This seems to imply a period of low (zero) prophylaxis at times immediately after treatment, with prophylaxis taking effect as the mean of the distribution is approached. We wanted to check this, and suggest that the survival curve should be used instead of the weibull pdf? Thus, prophylaxis is initially at its maximum level (directly after treatment) and then declines monotonically.

Lucy has given us the following weibull survival curve parameterisations (for the two ACT drugs):
time<-0:50

scale_al<-10.6
shape_al<-11.3

scale_dp<-28.1
shape_dp<-4.4

p_protect_al<- exp(-((time/scale_al)^shape_al))
p_protect_dp<- exp(-((time/scale_dp)^shape_dp))

plot(time,p_protect_dp,type="l", ylab="probability of protection")
lines(time,p_protect_al,col="blue")

Thanks!

We want to see the values for the default parameters at a distance on get_parameters