# Load the requisite packages:
library(malariasimulation)
# Set colour palette:
cols <- c("#E69F00", "#56B4E9", "#009E73", "#CC79A7","#F0E442", "#0072B2", "#D55E00")
One of the major threats to the continued success of efforts to
reduce the burden of malaria is the evolution and spread of resistance
to the antimalarial drugs used to treat uncomplicated cases of malaria.
The most effective frontline antimalarials are a class of drugs called
artemisinin combination therapies (ACTs). ACTs combine a fast-acting,
short-lived artemisinin derivative (e.g. artemether), with a
slower-acting, longer-lasting partner drug (e.g. lumefantrine) that
clears the remaining parasites from the patient’s system. Efforts to
understand the effect of resistance to ACTs on malaria morbidity and
mortality, and to develop strategies to control or mitigate the spread
of resistance, would benefit from insights derived from mathematical
modelling. Building on the model developed by Slater et al. (2016),
malariasimulation
provides the functionality to simulate
the effects of resistance to the artemisinin and/or partner drug
components to multiple, independent ACTs, on malaria transmission
dynamics.
Resistance to the artemisinin component of an ACT can result either
in slow parasite clearance (SPC), in which treatment with an ACT takes
longer than 3 days to fully clear patients with resistant parasites, or
early treatment failure (ETF), in which the ACT fails to clear the
infection and the individual develops a clinical infection. Resistance
to the partner drug, where the partner drug fails to clear the parasite
after the artemisinin derivative is depleted, results in infections
recrudescing to either clinical (D) or asymptomatic infections (A).
Resistance to the partner drug can also result in individuals developing
a novel, resistant infection following treatment, as the prophylaxis
provided by the ACT fails to protect the individual against reinfection
by a resistant strain. In the following vignette, we illustrate how to
parameterise and run malariasimulation
simulations with
resistance to ACTs deployed as a clinical treatment
Simulations capturing the effects of resistance to clinical treatment
using antimalarial drugs are parameterised using the
set_antimalarial_resistance()
function. This function
appends user-defined resistance parameters to a
malariasimulation
parameter list and accepts ten inputs.
The first is a list of malariasimulation
parameters to
append the resistance parameters to, and the second the index of the
drug
for which resistance is being parameterised, as set
using the set_drugs()
function. The
set_antimalarial_resistance()
function requires the
timesteps
, artemisinin_resistance_proportion
,
partner_drug_resistance_proportion_proportion
,
slow_parasite_clearance_probability
,
early_treatment_failure_probability
,
late_clinical_failure_probability
,
late_parasitological_failure_prob
, and
reinfection_during_prophylaxis_probability
inputs to be of
equal length so that, for each time step in which an update occurs, a
value is available for each parameter. Finally, the
slow_parasite_clearance_time
parameter represents the mean
residence time, in days, for artemisinin-resistant individuals
experiencing slow parasite clearance (SPC) in the Treated compartment,
and must be input as a single, positive value.
To illustrate how to parameterise resistance to an ACT using the
set_antimalarial_resistance()
function, we’ll set-up and
run three simulations. The first simulates malaria transmission in the
absence of interventions or resistance. The second simulates a simple
regime of clinical treatment in which 80% of clinical cases are treated
with artemether lumefantrine (AL), beginning after one year, in the
absence of antimalarial resistance. The third simulates the same
clinical treatment programme but with resistance to the artemisinin
component of AL emerging after two years. For illustrative purposes, we
assume that the proportion of infections resistant to the artemisinin
component of AL increases from 0% to 80%, and that these infections have
a 90% chance of resulting in early treatment failure.
First, we load the default malariasimulation
model
parameters, using the overrides
argument to increase the
human population. The human and mosquito population parameters are then
calibrated to a user-specified initial EIR using the
set_equilibrium()
function. Next, we load the in-built
parameters for the antimalarial drug AL and append them to the parameter
list using set_drugs()
. We can then use
set_clinical_treatment()
to specify a clinical treatment
regime, beginning after one year, that treats, on average, 80% of the
clinical cases of malaria with AL (AL_params
). The
set_antimalarial_resistance()
function is then used to
specify a scenario in which resistance is initially absent from the
population, but after two years the proportion of malaria infections
that are resistant to the artemisinin component of AL rises to 80%. We
also set the probability that artemisinin-resistant infections result in
early treatment failure to 0.9. In the current instance, we’ve set the
proportion of infections resistant to the AL partner drug to 0% and the
probabilities of other resistant infection outcomes to zero for
simplicity.
# Specify the time steps over which to simulate:
timesteps <- 365 * 3
# Specify an initial EIR to calibrate to:
init_EIR <- 8
# Specify a time step for treatment to begin:
treatment_start <- (1 * 365) + 1
# Specify a time step for resistance to emerge:
resistance_start <- (2 * 365) + 1
# Load the base malariasimulation parameter set:
simparams <- get_parameters(
overrides = list(
human_population = 10000)
)
# Calibrate to the initial EIR:
simparams <- set_equilibrium(parameters = simparams, init_EIR = init_EIR)
# Append the parameters for artemether lumefantrine (AL) to the parameter list:
simparams_clin_treatment <- set_drugs(parameters = simparams, drugs = list(AL_params))
# Use the set_clinical_treatment() function to specify a treatment regime for AL:
simparams_clin_treatment <- set_clinical_treatment(parameters = simparams_clin_treatment,
drug = 1,
timesteps = treatment_start,
coverages = 0.8)
# Use the set_antimalarial_resistance() function to specify resistance to AL:
simparams_resistance <- set_antimalarial_resistance(parameters = simparams_clin_treatment,
drug = 1,
timesteps = c(0, resistance_start),
artemisinin_resistance_proportion = c(0, 0.8),
partner_drug_resistance_proportion = rep(0, 2),
slow_parasite_clearance_probability = rep(0, 2),
early_treatment_failure_probability = c(0, 0.9),
late_clinical_failure_probability = rep(0, 2),
late_parasitological_failure_probability = rep(0, 2),
reinfection_during_prophylaxis_probability = rep(0, 2),
slow_parasite_clearance_time = 10)
We can now use the run_simulation()
to simulate the
three scenarios for which we have established parameter lists above.
# Baseline: No-intervention, no resistance simulation:
sim_out_baseline <- run_simulation(timesteps = timesteps, parameters = simparams)
# Clinical treatment with no antimalarial drug resistance:
sim_out_clin_treatment <- run_simulation(timesteps = timesteps, parameters = simparams_clin_treatment)
# Clinical treatment with antimalarial drug resistance:
sim_out_resistance <- run_simulation(timesteps = timesteps, parameters = simparams_resistance)
With the outputs from the run_simulation()
function, we
can visualise the effect of resistance on malaria transmission by
plotting the prevalence of malaria in children aged 2-10 years old
(PfPR2-10) over time.
# In each timestep, calculate PfPR_2-10 and add it to as a new column for each simulation output:
sim_out_baseline$pfpr210 = sim_out_baseline$n_detect_lm_730_3650/sim_out_baseline$n_age_730_3650
sim_out_clin_treatment$pfpr210 = sim_out_clin_treatment$n_detect_lm_730_3650/sim_out_clin_treatment$n_age_730_3650
sim_out_resistance$pfpr210 = sim_out_resistance$n_detect_lm_730_3650/sim_out_resistance$n_age_730_3650
# Plot the prevalence through time in each of the three simulated scenarios:
plot.new(); par(mar = c(4, 4, 1, 1), new = TRUE)
plot(x = sim_out_baseline$timestep, y = sim_out_baseline$pfpr210,
xlab = "Time (days)",
ylab = expression(paste(italic(Pf),"PR"[2-10])), cex = 0.7,
ylim = c(0, 1), type = "l", lwd = 2, xaxs = "i", yaxs = "i",
col = cols[3])
# Add a line for the clinical treatment scenario:
lines(x = sim_out_clin_treatment$timestep,
y = sim_out_clin_treatment$pfpr210,
col = cols[4])
# Add a line for the clinical treatment with resistance scenario:
lines(x = sim_out_resistance$timestep,
y = sim_out_resistance$pfpr210,
col = cols[7])
# Add lines to indicate the initiation of treatment and resistance:
abline(v = treatment_start, lty = "dashed")
abline(v = resistance_start, lty = "dashed")
# Annotate the added vlines:
text(x = treatment_start + 10, y = 0.9, labels = "Start of\nTreatment", adj = 0, cex = 0.7)
text(x = resistance_start + 10, y = 0.9, labels = "Start of\nResistance", adj = 0, cex = 0.7)
# 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("Baseline", "Treatment", "Resistance"),
col= c(cols[3:4], cols[7]), box.col = "white",
lwd = 1, lty = c(1, 1), cex = 0.7)
In the absence of clinical treatment or resistance, prevalence is comparable between all three scenarios for the first year. Following the initiation of clinical treatment at the beginning of the second year, PfPR2-10 approximately halves relative to the no-intervention baseline. However, following the introduction of artemisinin resistance at the beginning of the third year, early treatment failure causes the PfPR2-10 to increase to an intermediate level in the resistance scenario.
We can also capture scenarios in which resistance to a drug changes through time. To illustrate, we’ll establish and simulate a scenario in which resistance to sulfadoxine-pyrimethamine amodiaquine (SP-AQ) is absent from the population in the first year, but emerges in the third year and doubles in proportion each year thereafter until 100% of infections are artemisinin resistant. For simplicity, we’ll assume only artemisinin resistance is present in the population, and resistance to artemisinin results only, and always, in early treatment failure.
First, we store in vectors the artemisinin resistance proportions and
the time steps on which they will be updated in the simulation. We also
create a vector of early treatment failure probabilities which, for
simplicity, we assume remain at 1 for each update. Next, we load the
default malariasimulation
parameter set, specifying a
larger population size and seasonal transmission, and append the
parameters for SP-AQ (SP_AQ_params
) using the
set_drugs()
function. We’ll specify a simple treatment
regimen using set_clinical_treatment()
where, on average,
80% of clinical cases are treated with SP-AQ, beginning after one year.
We then specify a resistance schedule in which artemisinin resistance is
introduced at a proportion of 0.2 after 3 years, and doubles each year
thereafter until all infections are artemisinin resistant. Finally, we
calibrate the human and mosquito population parameters to a defined
entomological inoculation rate (EIR) and are ready to run the
simulation.
# Specify the time steps over which to simulate:
timesteps <- 365 * 8
# Set an initial EIR value:
initial_eir <- 8
# Specify the proportion of infections that are artemisinin resistant:
resistance_updates <- c(0, 0.2, 0.4, 0.8, 1)
# Specify the time steps in which to update the resistance parameters:
resistance_update_timesteps <- c(0, seq(3*365, 6*365, by = 365))
# Specify the probability artemisinin resistance infections result in early treatment failure:
early_treatment_failure_updates <- rep(1, length(resistance_update_timesteps))
# Load the base malariasimulation parameter set, with seasonal transmission:
simparams <- get_parameters(
overrides = list(
human_population = 1000,
model_seasonality = TRUE,
g0 = 0.284596,
g = c(-0.317878,-0.0017527,0.116455),
h = c(-0.331361,0.293128,-0.0617547)
))
# Append the parameters for sulfadomamine pyrimethaine (SP-AQ) to the parameter list:
simparams <- set_drugs(parameters = simparams, drugs = list(SP_AQ_params))
# Use the set_clinical_treatment() function to specify a treatment regime for SP-AQ:
simparams <- set_clinical_treatment(parameters = simparams,
drug = 1,
timesteps = 365 * 1,
coverages = 0.8)
# Use the set_antimalarial_resistance() function to specify resistance to SP-AQ:
simparams <- set_antimalarial_resistance(parameters = simparams,
drug = 1,
timesteps = resistance_update_timesteps,
artemisinin_resistance_proportion = resistance_updates,
partner_drug_resistance_proportion = rep(0, length(resistance_update_timesteps)),
slow_parasite_clearance_probability = rep(0, length(resistance_update_timesteps)),
early_treatment_failure_probability = early_treatment_failure_updates,
late_clinical_failure_probability = rep(0, length(resistance_update_timesteps)),
late_parasitological_failure_probability = rep(0, length(resistance_update_timesteps)),
reinfection_during_prophylaxis_probability = rep(0, length(resistance_update_timesteps)),
slow_parasite_clearance_time = 10)
# Calibrate the parameters to an initial EIR:
simparams <- set_equilibrium(parameters = simparams, init_EIR = initial_eir)
We can now use our parameter list to run the simulation using the the
run_simulation()
function.
We can visualise the effect of increasing resistance through time by plotting the PfPR2-10. We’ve added vertical lines to indicate when clinical treatment begins, and when the proportion of infections resistant to artemisinin is updated.
# Calculate the prevalence (PfPR210) through time:
dynamic_resistance_output$pfpr210 <- dynamic_resistance_output$n_detect_lm_730_3650/dynamic_resistance_output$n_age_730_3650
# Open a new plotting window and add a grid:
plot.new(); par(mar = c(4, 4, 1, 1), new = TRUE)
plot(x = dynamic_resistance_output$timestep,
y = dynamic_resistance_output$pfpr210,
xaxt = "n",
xlab = "Time (years)",
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])
# Specify x-axis ticks and labels:
axis(1, at = seq(0, 8 * 365, by = 365), labels = seq(0, 8 * 365, by = 365)/365)
# Add a line indicating the start of the clinical treatment:
abline(v = 365, lty = "dotted")
# Add lines indicating when resistance is updated:
abline(v = resistance_update_timesteps, lty = "dashed")
# Add a line highlighting the maximum PfPR_2-10 value prior to treatment or resistance:
abline(h = max(dynamic_resistance_output$pfpr210[1:365]), col = "red")
# Add annotations for the vlines:
text(x = 365 + 30, y = 0.6, labels = "Clin. treat. begins", adj = 0, cex = 0.8, srt = 90)
text(x = resistance_update_timesteps[2:5] + 30, y = 0.6, labels = paste0("Art. Res. = ", resistance_updates[2:5]),
adj = 0, cex = 0.8, srt = 90)
Looking at the figure, we can see that the PfPR2-10 decreases over the two years following the onset of clinical treatment in the absence of artemisinin resistance. However, as resistance is introduced and increases through time, the PfPR2-10 increases towards the pre-intervention seasonal peak as SP-AQ becomes increasingly ineffective in the treatment of clinical cases of malaria.
As we’ve discussed, resistance to an ACT can manifest in multiple ways. For instance, resistance to the artemisinin component of an ACT can result in either early treatment failure or slow parasite clearance.
Using malariasimulation
, we can simulate the effects of
multiple potential outcomes of resistance on malaria transmission
dynamics. To illustrate, we’ll parameterise and run a series of
simulations in is i) no clinical treatment, ii) no resistance, iii)
resistance with early treatment failure, iv) resistance with slow
parasite clearance, and v) resistance with early treatment failure and
slow parasite clearance.
# Determine the number of timesteps to run for:
timesteps <- 365 * 6
# Set up a list to store the simulation parameter lists in:
simulation_parameters <- list()
# Establish a list of the base parameters with no clinical treatment or antimalarial resistance:
get_parameters(overrides = list(human_population = 1000)) -> simulation_parameters$base
# Establish a parameter list with clinical treatment starting after one year:
simulation_parameters$base |>
set_drugs(drugs = list(AL_params)) |>
set_clinical_treatment(drug = 1, timesteps = (365 * 1) + 1, coverages = c(0.8)) |>
set_equilibrium(init_EIR = 16) -> simulation_parameters$treatment
# Set the equilibrium for the base parameters:
simulation_parameters$base |>
set_equilibrium(init_EIR = 16) -> simulation_parameters$base
# Establish a parameter list with clinical treatment and early treatment failure
simulation_parameters$treatment |>
set_antimalarial_resistance(drug = 1,
timesteps = c(0, (365 * 3) + 1),
artemisinin_resistance_proportion = c(0, 0.8),
partner_drug_resistance_proportion = c(0, 0),
slow_parasite_clearance_probability = c(0, 0),
early_treatment_failure_probability = c(0, 0.8),
late_clinical_failure_probability = c(0, 0),
late_parasitological_failure_probability = c(0, 0),
reinfection_during_prophylaxis_probability = c(0, 0),
slow_parasite_clearance_time = 10) -> simulation_parameters$etf
# Establish a parameter list with clinical treatment and slow parasite clearance
simulation_parameters$treatment |>
set_antimalarial_resistance(drug = 1,
timesteps = c(0, (365 * 3) + 1),
artemisinin_resistance_proportion = c(0, 0.8),
partner_drug_resistance_proportion = c(0, 0),
slow_parasite_clearance_probability = c(0, 0.8),
early_treatment_failure_probability = c(0, 0),
late_clinical_failure_probability = c(0, 0),
late_parasitological_failure_probability = c(0, 0),
reinfection_during_prophylaxis_probability = c(0, 0),
slow_parasite_clearance_time = 10) -> simulation_parameters$spc
# Establish a parameter list with clinical treatment, early treatment failure and slow parasite clearance:
simulation_parameters$treatment |>
set_antimalarial_resistance(drug = 1,
timesteps = c(0, (365 * 3) + 1),
artemisinin_resistance_proportion = c(0, 0.8),
partner_drug_resistance_proportion = c(0, 0),
slow_parasite_clearance_probability = c(0, 0.8),
early_treatment_failure_probability = c(0, 0.8),
late_clinical_failure_probability = c(0, 0),
late_parasitological_failure_probability = c(0, 0),
reinfection_during_prophylaxis_probability = c(0, 0),
slow_parasite_clearance_time = 10) -> simulation_parameters$etf_spc
We can now use our lists of malariasimulation
parameters
to run the simulations.
# Open a list to store the simulation outputs in:
simulation_outputs <- list()
# Run the simulations
for(i in seq(length(simulation_parameters))) {
# Run the i-th simulation
simulation_temp <- run_simulation(timesteps = timesteps,
parameters = simulation_parameters[[i]])
# Append the simulation identifier:
simulation_temp$identifier <- names(simulation_parameters)[i]
# Append the ith simulation outputs to the combined simulation dataframe:
simulation_outputs[[names(simulation_parameters)[i]]] <- simulation_temp
# Print the number of columns in the i-th simulation outputs dataframe:
print(ncol(simulation_temp))
}
#> [1] 31
#> [1] 35
#> [1] 37
#> [1] 37
#> [1] 37
We can compare the effects of independent resistance outcomes with combined resistance outcomes visually. In the following plot, we compare the PfPR2-10 between a baseline without any clinical treatment or antimalarial resistance, a clinical-treatment only run, a clinical treatment with early treatment failure run, a clinical treatment with slow parasite clearance run, and a clinical treatment with both early treatment failure and slow parasite clearance run.
# 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 = simulation_outputs$base$timestep,
y = simulation_outputs$base$n_detect_lm_730_3650/simulation_outputs$base$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 = simulation_outputs$treatment$timestep,
y = simulation_outputs$treatment$n_detect_lm_730_3650/simulation_outputs$treatment$n_age_730_3650,
col = cols[4])
lines(x = simulation_outputs$etf$timestep,
y = simulation_outputs$etf$n_detect_lm_730_3650/simulation_outputs$etf$n_age_730_3650,
col = cols[5])
lines(x = simulation_outputs$spc$timestep,
y = simulation_outputs$spc$n_detect_lm_730_3650/simulation_outputs$spc$n_age_730_3650,
col = cols[6])
lines(x = simulation_outputs$etf_spc$timestep,
y = simulation_outputs$etf_spc$n_detect_lm_730_3650/simulation_outputs$etf_spc$n_age_730_3650,
col = cols[7])
# Add vlines to indicate when SP-AQ were administered:
abline(v = 365 + 1, lty = "dashed", lwd = 1)
text(x = (365 * 1) - 40, y = 0.6, labels = "Treatment Introduced", adj = 0, cex = 0.8, srt = 90)
abline(v = (365 * 3) + 1, lty = "dashed", lwd = 1)
text(x = (365 * 3) - 40, y = 0.6, labels = "Resistance Introduced", adj = 0, cex = 0.8, srt = 90)
# Add gridlines:
grid(lty = 2, col = "grey80", nx = NULL, ny = NULL, lwd = 0.5); box()
# Add a legend:
legend(x = 3000, y = 0.99, legend = c("Baseline", "Treatment", "ETF-only", "SPC-only", "ETF and SPC"),
col= c(cols[3:7]), box.col = "white",
lwd = 1, lty = c(1, 1), cex = 0.8)
Slater, H.C., Griffin, J.T., Ghani, A.C. and Okell, L.C., 2016. Assessing the potential impact of artemisinin and partner drug resistance in sub-Saharan Africa. Malaria journal, 15(1), pp.1-11.