Mass Drug Administriation and Chemoprevention

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

Malaria chemoprevention encompasses a suite of strategies for reducing malaria morbidity, mortality, and transmission through the treatment of at-risk populations with antimalarial drugs. Three such strategies are mass drug administration (MDA), seasonal malaria chemoprevention (SMC), and perennial malaria chemoprevention (PMC). MDA attempts to simultaneously treat an entire population, irrespective of their infection risk or status. SMC targets children aged 3-59 months with a course of an antimalarial through the peak malaria transmission season. PMC is designed for settings with perennially high transmission and seeks to administer a full antimalarial course at defined intervals irrespective of age or disease status.

malariasimulation enables users to simulate the effects of chemoprevention strategies on malaria transmission, and provides the flexibility to design and simulate a range of MDA/SMC/PMC campaigns for a suite of antimalarial drugs. Specifically, users are able to load or specify parameters describing both the characteristics of the antimalarial(s) (e.g. drug efficacy, prophylactic profile) and their distribution (e.g. coverage, distribution timing, and the section of the population to treat). This vignette provides simple, illustrative demonstrations of how MDA and SMC campaigns can be set-up and simulated using malariasimulation. The functionality to simulate PMC campaigns is also built into malariasimulation via the set_pmc() function, but we do not demonstrate the use of PMC in this vignette (for further information on their use, see set_pmc).

Mass Drug Administration

To illustrate how MDA can be simulated using malariasimulation, we will implement a campaign that distributes a single dose of sulphadoxine-pyrimethamine amodiaquine (SP-AQ) to 80% of the entire population once a year over a two year period, initiated at the end of the first year.

Parameterisation

The first step is to generate the list of malariasimulation parameters using the get_parameters() function. Here, we’ll use the default malariasimulation parameter set (run ?get_parameters() to view), but specify the seasonality parameters to simulate a highly seasonal, unimodal rainfall pattern (illustrated in the plot of total adult female mosquito population dynamics below). To do this, we set model_seasonality to TRUE within get_parameters(), and then tune the rainfall pattern using the g0, g, and h parameters (which represent fourier coefficients). In the malariasimulation model, rainfall governs the time-varying mosquito population carrying capacity through density-dependent regulation of the larval population and, therefore, annual patterns of malaria transmission. Next, we use the set_equilibrium() function to tune the initial human and mosquito populations to those required to observe the specified initial EIR.

# Establish the length of time, in daily time steps, over which to simulate:
year <- 365; years <- 3
sim_length <- years * year

# Set the size of the human population and an initial entomological inoculation rate (EIR)
human_population <- 1000
starting_EIR <- 50

# Use the simparams() function to build the list of simulation parameters and then specify
# a seasonal profile with a single peak using the g0, g, and h parameters.
simparams <- get_parameters(
  list(
    human_population = human_population,
    model_seasonality = TRUE,
    g0 = 0.284596,
    g = c(-0.317878,-0.0017527,0.116455),
    h = c(-0.331361,0.293128,-0.0617547)
  )
)

# Use the set_equilibrium() function to update the individual-based model parameters to those 
# required to match the specified initial EIR (starting_EIR) at equilibrium:
simparams <- set_equilibrium(simparams, starting_EIR)

Interventions

Having established a base set of parameters (simparams), the next step is to add the parameters specifying the MDA campaign. We first use the set_drugs() function to update the parameter list (simparams) with the parameter values for the drug(s) we wish to simulate, in this case using the in-built parameter set for SP-AQ (SP_AQ_params). Note that the malariasimulation package also contains in-built parameter sets for dihydroartemisinin-piperaquine (DHA_PQP_params) and artemether lumefantrine (AL_params), and users can also define their own drug parameters.

We then define the time steps on which we want the MDA to occur (stored here in the vector mda_events for legibility), and use the set_mda() function to update the parameter list with our MDA campaign parameters. set_mda() takes as its first two arguments the base parameter list and the time steps on which to schedule MDA events. This function also requires the user to specify the proportion of the population they want to treat using coverages, as well as the age-range of the population we want to treat using min_ages and max_ages, for each timestep on which MDA is scheduled. Note that this enables users to specify individual events to have different drug, coverage, and human-target properties. The drug argument refers to the index of the drug we want to distribute in the MDA in the list of parameters (in this demonstration we have only added a single drug, SP-AQ, using set_drugs(), but we could have specified additional drugs).

# Make a copy of the base simulation parameters to which we can add the MDA campaign parameters:
mdaparams <- simparams

# Update the parameter list with the default parameters for sulphadoxine-pyrimethamine
# amodiaquine (SP-AQ)
mdaparams <- set_drugs(mdaparams, list(SP_AQ_params))

# Specify the days on which to administer the SP-AQ
mda_events <- (c(1, 2) * 365)

# Use set_mda() function to set the proportion of the population that the MDA reaches and the age 
# ranges which can receive treatment.
mdaparams <- set_mda(mdaparams,
                     drug = 1,
                     timesteps = mda_events,
                     coverages = rep(.8, 2),
                     min_ages = rep(0, length(mda_events)),
                     max_ages = rep(200 * 365, length(mda_events)))

Having generated a base, no-intervention parameter set (simparams) and a version updated to define an MDA campaign (mdaparams), we can now use the run_simulation() function to run simulations for scenarios without interventions and with our MDA campaign.

Simulation

# Run the simulation in the absence of any interventions:
noint_output <- run_simulation(sim_length, simparams)

# Run the simulation with the prescribed MDA campaign
mda_output <- run_simulation(sim_length, mdaparams)

Visualisation

The run_simulation() function returns a dataframe containing time series for a range of variables, including the number of human individuals in each infectious state, and the number of people and detected cases in children aged 2-10. The plots below illustrate how the effect of an MDA campaign on the malaria transmission dynamics can be visualised. The first presents the distribution of people among the infectious states in the MDA simulation. The second compares the prevalence of malaria in children aged 2-10 years old (PfPR2-10) between the no-intervention and MDA scenarios simulated.

# Store the state names, and labels and colours with which to plot them:
states <- c("S_Count", "D_count", "Tr_count", "A_count", "U_count")
state_labels <- c("S", "D", "Tr", "A", "U")

# Open a new plotting window:
plot.new(); par(mar = c(4, 4, 1, 1), new = TRUE)

# Plot the time series
plot(x = mda_output$timestep, y = mda_output$S_count,
     type = "l", col = cols[1], lwd = 2,
     ylim = c(0, 950), ylab = "Individuals", xlab = "Time (days)",
     xaxs = "i", yaxs = "i", cex = 0.8)

# Overlay the time-series for other states
for(i in 1:(length(states)-1)) {
  lines(mda_output[,states[i+1]], col = cols[i+1], lwd = 2)
}

# Add vlines to show when SP-AQ was administered:
abline(v = mda_events, lty = "dashed", lwd = 1.5)
text(x = mda_events+10, y = 900, labels = "MDA int.", adj = 0, cex = 0.8)

# Add gridlines:
grid(lty = 2, col = "grey80", nx = NULL, ny = NULL, lwd = 0.5); box()

# Add a legend:
legend("topleft",
       legend = c(state_labels),
       col = c(cols[1:5]), lty = c(rep(1, 5)),
       lwd = 1, cex = 0.8, box.col = "white", ncol = 2, x.intersp = 0.5)

# Open a new plotting window and add a grid:
plot.new(); par(mar = c(4, 4, 1, 1), new = TRUE)

# Plot malaria prevalence in 2-10 years through time:
plot(x = mda_output$timestep,
     y = mda_output$n_detect_lm_730_3650/mda_output$n_age_730_3650,
     xlab = "Time (days)",
     ylab = expression(paste(italic(Pf),"PR"[2-10])), cex = 0.8,
     ylim = c(0, 1), type = "l", lwd = 2, xaxs = "i", yaxs = "i",
     col = cols[3])

# Add the dynamics for no-intervention simulation
lines(x = noint_output$timestep,
     y = noint_output$n_detect_lm_730_3650/noint_output$n_age_730_3650,
     col = cols[4])

# Add vlines to indicate when SP-AQ were administered:
abline(v = mda_events, lty = "dashed", lwd = 1)
text(x = mda_events+10, y = 0.95, labels = "MDA int.", adj = 0, cex = 0.8)

# Add gridlines:
grid(lty = 2, col = "grey80", nx = NULL, ny = NULL, lwd = 0.5); box()

# Add a legend:
legend(x = 20, y = 0.99, legend = c("MDA", "No-Int"),
       col= c(cols[3:4]), box.col = "white",
       lwd = 1, lty = c(1, 1), cex = 0.8)

Seasonal Malaria Chemoprevention

To demonstrate how to simulate SMC in malariasimulation , in the following section we’ll set-up and simulate a campaign that treats 90% of children aged 3-59 months old with four doses of SP-AQ per year. We will set the SMC campaign to begin in the second year, run for two years, and specify the four SMC events to occur at monthly intervals, starting two months prior to the peak malaria season.

Interventions

While the first step would typically be to establish the base set of parameters, we can use those stored earlier in simparams. We first store a copy of these base parameters and then use set_drugs() to store the preset parameters for SP-AQ in the parameter list. As SMC is conducted during the peak malaria season, we need to use the peak_season_offset() function to determine when the peak malaria season occurs given our specified seasonal profile. This function reads in the parameter list, uses g0, g, h, and rainfall_floor to generate the rainfall profile for a single year, and returns the day on which the maximum rainfall value occurs. We can then take this calculated seasonal peak and time SMC events around it. Here, we’ve specified four monthly drug administration days per year starting two months before the seasonal peak. To illustrate this, we can plot the total adult mosquito population size through time which, in the absence of interventions targeting mosquitoes, will closely match the rainfall pattern, and overlay both the annual seasonal peak and planned SMC events. Once we’re happy with our SMC campaign, we use the set_smc() function to update the parameter list with our drug, timesteps, coverages and target age group (min_ages and max_ages).

# Copy the original simulation parameters
smcparams <- simparams

# Append the parameters for SP-AQ using set_drugs
smcparams <- set_drugs(smcparams, list(SP_AQ_params))

# Use the peak_season_offset() to calculate the yearly offset (in daily timesteps) for the peak mosq.
# season
peak <- peak_season_offset(smcparams)

# Create a variable for total mosquito population through time:
noint_output$mosq_total = noint_output$Sm_gamb_count +
  noint_output$Im_gamb_count +
  noint_output$Pm_gamb_count

# Schedule drug administration times (in daily time steps) before, during and after the seasonal peak:
admin_days <- c(-60, -30, 0, 30)

# Use the peak-offset, number of simulation years and number of drug admin. days to calculate the days
# on which to administer drugs in each year
smc_days <- rep((365 * seq(1, years-1, by = 1)), each = length(admin_days)) + peak + rep(admin_days, 2)

# Turn off scientific notation for the plot:
options(scipen = 666)

# Open a new plotting window,set the margins and plot the total adult female mosquito population
# through time:
plot.new(); par(new = TRUE, mar = c(4, 4, 1, 1))
plot(x = noint_output$timestep, noint_output$mosq_total,
     xlab = "Time (days)", ylab = "Adult Female Mosquitos",
     ylim = c(0, max(noint_output$mosq_total)*1.1),
     type = "l", lwd = 2, cex = 0.8,
     xaxs = "i", yaxs = "i",
     col = cols[3])

# Add gridlines
grid(lty = 2, col = "grey80", nx = 11, ny = NULL, lwd = 0.5); box()

# Overlay the SMC days and the calculated seasonal peak
abline(v = smc_days, lty = "dashed", lwd = 2)
abline(v =  c(0, 1, 2) * 365 + peak, lty = 2, lwd = 2, col = cols[4])

# Add a legend:
legend("topleft",
       legend = c(expression("M"["Tot"]), "SMC", "Peak"),
       col= c(cols[3], "black", cols[4]), box.col = "white",
       lty = c(1, 2, 2), lwd = 2, cex = 0.9)


# Update the SMC parameters list:
smcparams <- set_smc(
  smcparams,
  drug = 1,
  timesteps = smc_days,
  coverages = rep(0.9, length(smc_days)),
  min_ages = rep(3 * 30, length(smc_days)),
  max_ages = rep(59 * 30, length(smc_days))
)

Simulations

Having generated our SMC parameter list, we can run the simulation using the run_simulation() function.

# Run the SMC simulation
smc_output <- run_simulation(sim_length, smcparams)

Visualisation

We can now plot the simulation output to visualise the effect of our SMC campaign on malaria transmission dynamics, again focusing on the distribution of people in each of the infectious states and on malaria prevalence in children aged 2-10 (PfPR2-10).

# Open a new plotting window and add a grid:
plot.new(); par(new = TRUE, mar = c(4, 4, 1, 1))

# Plot the time series of human infectious states through time under the SMC campaign
plot(x = smc_output$timestep, y = smc_output$S_count,
     type = "l", col = cols[1], lwd = 2,
     ylim = c(0, 950), ylab = "Individuals", xlab = "Time (days)",
     xaxs = "i", yaxs = "i", cex = 0.8)

# Overlay the time-series for other human infection states
for(i in 1:(length(states)-1)) {
  lines(smc_output[,states[i+1]], col = cols[i+1], lwd = 2)
}

# Add vlines to indicate when SMC drugs were administered:
abline(v = smc_days, lty = "dashed", lwd = 1)
text(x = smc_days[c(4,8)]+10, y = 900, labels = "SMC\nint.", adj = 0, cex = 0.8)

# Add gridlines:
grid(lty = 2, col = "grey80", nx = NULL, ny = NULL, lwd = 0.5)

# Add a legend:
legend("topleft",
       legend = c(state_labels),
       col = c(cols[1:5]), lty = c(rep(1, 5)),
       box.col = "white", lwd = 2, cex = 0.8, ncol = 2, x.intersp = 0.5) 

# Open a new plotting window and add a grid:
plot.new(); par(new = TRUE, mar = c(4, 4, 1, 1))

# Plot malaria prevalence in 2-10 years through time:
plot(x = smc_output$timestep, y = smc_output$n_detect_lm_730_3650/smc_output$n_age_730_3650,
     xlab = "Time (days)", 
     ylab = expression(paste(italic(Pf),"PR"[2-10])), cex = 0.8,
     ylim = c(0, 1), type = "l", lwd = 2, xaxs = "i", yaxs = "i",
     col = cols[3])

# Add the dynamics for no-intervention simulation
lines(x = noint_output$timestep,
     y = noint_output$n_detect_lm_730_3650/noint_output$n_age_730_3650,
     col = cols[4])

# Add lines to indicate SMC events:
abline(v = smc_days, lty = "dashed", lwd = 1)
text(x = smc_days[c(4,8)]+10, y = 0.95, labels = "SMC\nint.", adj = 0, cex = 0.8)

# Add gridlines:
grid(lty = 2, col = "grey80", nx = 11, ny = 10, lwd = 0.5)

# Add a legend:
legend("topleft",
       legend = c("SMC", "No-Int"),
       col= c(cols[3:4]), box.col = "white",
       lwd = 1, lty = c(1, 1), cex = 0.8)