In mrc-ide/malariasimulation: An individual based model for malaria

knitr::opts_chunk$set(
  collapse = TRUE,
  comment = "#>",
  dpi=300,
  fig.width=7
)

# Load the requisite packages:
library(malariasimulation)
# Set colour palette:
cols <- c("#E69F00", "#56B4E9", "#009E73", "#F0E442", "#0072B2", "#D55E00", "#CC79A7")

As alluded to in the Model Structure, Vector Control: IRS, and Vector Control: Bed net vignettes, it is possible to account for varying proportions of mosquito species in the setting you are modelling using the set_species() function. IRS and bed nets could be expected to have different impacts depending on the proportion of each mosquito species because of variations in indoor resting, insecticide resistance, and proportion of bites on humans versus animals by species. If you have specified more than 1 species, then the arguments for set_spraying() and/or set_bednets() must be populated with values for each species at each timestep that the intervention is implemented.

There are preset parameters for An. gambiae, An. arabiensis, An. funestus and An. stephensi that can be set by the helper objects gamb_params, arab_params, fun_params and steph_params, respectively. The default values for each species in these helper functions are from Sherrard-Smith et al., 2018. The parameters are:

species: the mosquito species name or signifier
blood_meal_rates: the blood meal rates for each species
foraging_time: time spent taking blood meals
Q0: proportion of blood meals taken on humans
phi_bednets: proportion of bites taken in bed
phi_indoors: proportion of bites taken indoors

We will demonstrate how to specify different mosquito species and how this could alter intervention impact using an example with IRS.

We will create a plotting function to visualise the output.

# Plotting functions
plot_prev <- function() {
  plot(x = output_endophilic$timestep, y = output_endophilic$n_detect_lm_730_3650 / output_endophilic$n_age_730_3650, 
       type = "l", col = cols[3], lwd = 1,
       xlab = "Time (days)", ylab = expression(paste(italic(Pf),"PR"[2-10])),
       xaxs = "i", yaxs = "i", ylim = c(0,1))
  lines(x = output_exophilic$timestep, y = output_exophilic$n_detect_lm_730_3650 / output_exophilic$n_age_730_3650,
        col = cols[5], lwd = 1)
  abline(v = sprayingtimesteps, lty = 2, lwd = 2.5, col = "black")
  text(x = sprayingtimesteps + 10, y = 0.9, labels = "Spraying\nint.", adj = 0, cex = 0.8)
  grid(lty = 2, col = "grey80", lwd = 0.5)
  legend("bottomleft", box.lty = 0, 
         legend = c("Prevalence for endophilic mosquitoes", "Prevalence for exophilic mosquitoes"),
         col = c(cols[3], cols[5]), lty = c(1,1), lwd = 2, cex = 0.8, y.intersp = 1.3)
}

species_plot <- function(Mos_pops){
  plot(x = Mos_pops[,1], y = Mos_pops[,2], type = "l", col = cols[2],
       ylim = c(0,max(Mos_pops[,-1]*1.25)), ylab = "Mosquito population size", xlab = "Days",
       xaxs = "i", yaxs = "i", lwd = 2)
  grid(lty = 2, col = "grey80", lwd = 0.5)
  sapply(3:4, function(x){
    points(x = Mos_pops[,1], y = Mos_pops[,x], type = "l", col = cols[x])})
  legend("topright", legend = c("A. arab","A. fun","A. gamb"),
         col = cols[-1], lty = 1, lwd = 2, ncol = 1, cex = 0.8, bty = "n")
  }

Setting mosquito species parameters

Single endophilic mosquito species

Use the get_parameters() function to generate a list of parameters, accepting most of the default values, but modifying seasonality values to model a seasonal setting. We will use set_species to model mosquitoes similar to An. funestus but with a higher propensity to bite indoors (which we will name "endophilic"). Then, we use the set_equilibrium() function to to initialise the model at a given entomological inoculation rate (EIR).

We used the set_spraying() function to set an IRS intervention. This function takes as arguments the parameter list, timesteps of spraying, coverage of IRS in the population and a series of parameters related to the insecticide used in the IRS. The proportion of mosquitoes dying following entering a hut is dependent on the parameters ls_theta, the initial efficacy, and ls_gamma, how it changes over time. The proportion of mosquitoes successfully feeding is dependent on ks_theta, the initial impact of the insecticide in IRS, and ks_gamma, how the impact changes over time. Finally, the proportion of mosquitoes being deterred away from a sprayed hut depends on ms_theta, the initial impact of IRS, and ms_gamma, the change in impact over time. See a more comprehensive explanation in the Supplementary Information of Sherrard-Smith et al., 2018.

year <- 365
month <- 30
sim_length <- 3 * year
human_population <- 1000
starting_EIR <- 50

simparams <- get_parameters(
  list(
    human_population = human_population,
    # seasonality parameters
    model_seasonality = TRUE, 
    g0 = 0.285277,
    g = c(-0.0248801, -0.0529426, -0.0168910),
    h = c(-0.0216681, -0.0242904, -0.0073646)
  )
)

peak <- peak_season_offset(simparams)

# Create an example mosquito species (named endophilic) with a high value for `phi_indoors`
endophilic_mosquito_params <- fun_params
endophilic_mosquito_params$phi_indoors <- 0.9
endophilic_mosquito_params$species <- 'endophilic'

## Set mosquito species  with a high propensity for indoor biting
simparams <- set_species(
  simparams, 
  species = list(endophilic_mosquito_params),
  proportions = c(1)
)

sprayingtimesteps <- c(1, 2) * year + peak - 3 * month # There is a round of IRS in the 1st and second year 3 months prior to peak transmission.

simparams <- set_spraying(
  simparams,
  timesteps = sprayingtimesteps,
  coverages = rep(.8, 2), # # Each round covers 80% of the population 
  # nrows=length(timesteps), ncols=length(species) 
  ls_theta = matrix(2.025, nrow=length(sprayingtimesteps), ncol=1), # Matrix of mortality parameters per round of IRS and per species 
  ls_gamma = matrix(-0.009, nrow=length(sprayingtimesteps), ncol=1), # Matrix of mortality parameters per round of IRS and per species
  ks_theta = matrix(-2.222, nrow=length(sprayingtimesteps), ncol=1), # Matrix of feeding success parameters per round of IRS and per species
  ks_gamma = matrix(0.008, nrow=length(sprayingtimesteps), ncol=1), # Matrix of feeding success parameters per round of IRS and per species
  ms_theta = matrix(-1.232, nrow=length(sprayingtimesteps), ncol=1), # Matrix of deterrence parameters per round of IRS and per species
  ms_gamma = matrix(-0.009, nrow=length(sprayingtimesteps), ncol=1) # Matrix of deterrence parameters per round of IRS and per species
)

simparams <- set_equilibrium(simparams, starting_EIR)

# Running simulation with IRS
output_endophilic <- run_simulation(timesteps = sim_length, parameters = simparams)

We can see below that only the endophilic species is modelled.

simparams$species
simparams$species_proportions

Single exophilic mosquito species

We will run the same model with IRS as above, but this time with an example mosquito species similar to An. funestus, but with a lower propensity to bite indoors (which we will name "exophilic"). Note that now there are two rows for the ls_theta, ls_gamma, ks_theta, ks_gamma, ms_theta, and ms_gamma arguments (rows represent timesteps where changes occur, columns represent additional species). See the Vector Control: IRS vignette for more information about setting IRS with different mosquito species.

# Create an example mosquito species (named exophilic) with a low value for `phi_indoors`
exophilic_mosquito_params <- fun_params
exophilic_mosquito_params$phi_indoors <- 0.2
exophilic_mosquito_params$species <- 'exophilic'

## Set mosquito species  with a low propensity for indoor biting
simparams <- set_species(
  simparams, 
  species = list(exophilic_mosquito_params),
  proportions = c(1)
)

peak <- peak_season_offset(simparams)

sprayingtimesteps <- c(1, 2) * year + peak - 3 * month # There is a round of IRS in the 1st and second year 3 months prior to peak transmission.

simparams <- set_spraying(
  simparams,
  timesteps = sprayingtimesteps,
  coverages = rep(.8, 2), # # Each round covers 80% of the population 
  # nrows=length(timesteps), ncols=length(species) 
  ls_theta = matrix(2.025, nrow=length(sprayingtimesteps), ncol=1), # Matrix of mortality parameters
  ls_gamma = matrix(-0.009, nrow=length(sprayingtimesteps), ncol=1), # Matrix of mortality parameters per round of IRS and per species
  ks_theta = matrix(-2.222, nrow=length(sprayingtimesteps), ncol=1), # Matrix of feeding success parameters per round of IRS and per species
  ks_gamma = matrix(0.008, nrow=length(sprayingtimesteps), ncol=1), # Matrix of feeding success parameters per round of IRS and per species
  ms_theta = matrix(-1.232, nrow=length(sprayingtimesteps), ncol=1), # Matrix of deterrence parameters per round of IRS and per species
  ms_gamma = matrix(-0.009, nrow=length(sprayingtimesteps), ncol=1) # Matrix of deterrence parameters per round of IRS and per species
)

output_exophilic <- run_simulation(timesteps = sim_length, parameters = simparams)

Plot adult female infectious mosquitoes by species over time

In the plot below, we can see that IRS is much more effective when the endophilic mosquito species is modelled compared to the scenario where an exophilic species is modelled. In this case, IRS will not be as effective because a larger proportion of bites take place outside of the home.

plot_prev()

Setting multiple mosquito species

Finally, we give an example of how to set multiple mosquito species.

# Update parameter list with species distributions
simparams <- get_parameters(
  list(
    human_population = human_population,
    # seasonality parameters
    model_seasonality = TRUE, 
    g0 = 0.285277,
    g = c(-0.0248801, -0.0529426, -0.0168910),
    h = c(-0.0216681, -0.0242904, -0.0073646)
  )
)

params_species <- set_species(parameters = simparams,
                              species = list(arab_params, fun_params, gamb_params),
                              proportions = c(0.1,0.3,0.6))

# Run simulation
species_simulation <- run_simulation(timesteps = sim_length, parameters = params_species)

## Plot species distributions
Mos_sp_dist_sim <- species_simulation[,c("timestep", "total_M_arab", "total_M_fun", "total_M_gamb")]
species_plot(Mos_sp_dist_sim)