The influence of biological, epidemiological, and treatment factors on the establishment and spread of drug-resistant Plasmodium falciparum

1. ThieryMasserey12
2. TamsinLee12
3. MonicaGolumbeanu12
4. AndrewJShattock12
5. SherrieLKelly12
6. IanMHastings3
7. MelissaAPenny[email protected]12

1. Swiss Tropical and Public Health Institute
   AllschwilSwitzerland
2. University of Basel
   BaselSwitzerland
3. Liverpool School of Tropical Medicine
   LiverpoolUnited Kingdom

• Research Article

• Epidemiology and Global Health
• Microbiology and Infectious Disease

• modelling and simulation
• drug resistance
• artemisinin combination therapy
• P. falciparum

• publisher-id77634
• doi10.7554/eLife.77634
• elocation-ide77634

Introduction

Malaria remains a global health priority 103WHO2020. The WHO recommends several artemisinin-based combination therapies (ACTs) to treat uncomplicated Plasmodium falciparum malaria 104WHO2020. ACTs combine a short-acting artemisinin derivative to rapidly reduce parasitaemia during the first 3 days of treatment and a long-acting partner drug to eliminate remaining parasites 104WHO2020. These drug combinations are intended to delay the evolution of drug resistance, which has frequently occurred under monotherapy treatment 27Farooq and Mahajan200499White1999100White2004109Wongsrichanalai et al.2002. However, parasites partially resistant to artemisinin have emerged in the Greater Mekong Subregion (GMS) and, more recently, in Rwanda, Uganda, Guyana, and Papua New Guinea despite the use of ACTs 104WHO202014Chenet et al.201692Uwimana et al.202093Uwimana et al.202168Miotto et al.20204Balikagala et al.2021. Partial artemisinin resistance leads to slower parasite clearance following treatment with ACTs, but not necessarily to treatment failure 104WHO2020. However, high rates of treatment failure have been observed in the GMS due to parasites being less sensitive to artemisinin derivatives and their partner drugs 104WHO2020. To prevent the evolution of drug-resistant parasites and to preserve the efficacy of ACTs or triple artemisinin-based combination therapies (TACT, including a second long-acting drug) now being tested 94Pluijm et al.2020, it is essential to understand which factors drive this process.

The evolution of drug resistance follows a three-step process of mutation, establishment, and spread. First, mutations conferring drug resistance emerge in the population at a rate that depends on multiple factors, such as organism mutation and migration rates 113Wiesch et al.201162Mackinnon2005. Second, establishment is a highly stochastic step as the parasite with the drug-resistant mutation needs to infect other hosts 113Wiesch et al.201162Mackinnon200535Hastings200438Hastings et al.2020. The resistant strain establishes in the population once its frequency is high enough to minimise its risk of stochastic extinction 113Wiesch et al.201162Mackinnon200535Hastings200438Hastings et al.2020. Several forces influence the establishment of mutations. In settings with higher heterogeneity of parasite reproductive success, establishment of mutations is less likely because the effects of stochasticity are more substantial 113Wiesch et al.201135Hastings200438Hastings et al.202033Hastings and Mackinnon1998. This heterogeneity depends on the level of transmission and health system strength 113Wiesch et al.201135Hastings200438Hastings et al.202033Hastings and Mackinnon199851Klein2014. In addition, the more selection favours the resistant strain, the more likely it is to establish 113Wiesch et al.201135Hastings200438Hastings et al.202033Hastings and Mackinnon1998. The strength of selection depends on many factors, such as the parasite and human biology, the transmission setting, drug properties, and health system strength 100White20042Antao and Hastings201141Hughes and Andersson201542Huijben et al.201367Miotto et al.201580Slater et al.201763Mackinnon and Marsh2010. Third, resistance spreads through a region after a resistant mutation has become established. The mutation spreads at a rate that depends on the strength of selection 113Wiesch et al.201138Hastings et al.2020.

It is not fully understood how factors intrinsic to the transmission setting, health system, human and parasite biology, and drug properties interact to influence the establishment and spread of drug-resistant parasites. Mathematical models of infectious disease have not previously been used to systematically assess the joint influence of multiple factors on the establishment and spread of drug-resistant P. falciparum (e.g. 51Klein2014; 80Slater et al.2017; 9Bushman et al.2018; 102Whitlock et al.2021; 34Hastings et al.2002; 97Watson et al.2021; 7Brock et al.2018; 96Watkins and Mosobo1993; 72Pongtavornpinyo et al.2008; 101White et al.2009; 17Chiyaka et al.2009; 26Esteva et al.2009; 53Koella and Antia2003; 56Lee and Penny2019; 57Lee et al.2022; 58Legros and Bonhoeffer2016; 88Tchuenche et al.2011; 91Tumwiine et al.2007) (and to the best of our knowledge other drug-resistant pathogens [virus, bacteria, etc.]). Simple models, based on the Ross and MacDonald model 61Macdonald195774Ross1915, have considered specific components of the epidemiology of resistance and, therefore, are not sophisticated enough to answer questions on how factors have jointly impacted establishment and spread of drug resistance 7Brock et al.201817Chiyaka et al.200926Esteva et al.200953Koella and Antia200356Lee and Penny201988Tchuenche et al.201191Tumwiine et al.2007. Most models have investigated specific transmission scenarios and questions, such as how within-host competition between parasites influences development of drug resistance 9Bushman et al.201872Pongtavornpinyo et al.200858Legros and Bonhoeffer2016, and did not systematically assess the impact of assumptions used on their results. Consequently, previous studies have neither systematically compared the influence of multiple drivers, nor assessed how their influence varies under different transmission settings or health system strengths.

In addition, most models have made simplifications concerning drug action and consequences of partial resistance. They have typically not explicitly modelled the pharmacokinetics and pharmacodynamics of the drugs and have assumed that resistant parasites are fully resistant to the drugs. Parasites partially resistant to artemisinin exhibit an extended ring-stage during which they are not sensitive to artemisinin; however, parasites remain sensitive to artemisinin during other stages 52Klonis et al.201395Wang et al.201778Sá2018107Witkowski et al.2013112Ye et al.2016. In addition, parasites resistant to partner drugs have an increased minimum inhibitory concentration (MIC), meaning that they are not sensitive to low drug concentrations but remain susceptible to high concentrations of partner drugs 11Chaorattanakawee et al.201610Chaorattanakawee et al.201587Tahita et al.2015. Consequently, many models have ignored the residual effect of drugs on resistant parasites and have not investigated the influence of the degree of resistance and drug proprieties on the establishment and spread of drug resistance. Models that have explicitly considered drug action have focused on specific questions such as how half-life impacts the spread of resistance or how resistance to the partner drug influences evolution of artemisinin resistance 34Hastings et al.200297Watson et al.2021. However, they did not investigate how the impact of drug proprieties and the degree of resistance interact with other biological, transmission, and health system factors.

In this study, we developed a disease model with an emulator-based approach to quantify the influence of factors intrinsic to the biology of the parasite and human, the transmission setting, the health system strength, and the drug properties on the establishment and spread of drug-resistant parasites. Our approach is based on a detailed individual-based malaria model, OpenMalaria (https://github.com/SwissTPH/openmalaria/wiki), that includes a mechanistic within-host model (based on 70Molineaux et al.2002). We first adapted our model, OpenMalaria, to explicitly include mechanistic drug action models at the individual level (as a one, two, or three-compartment pharmacokinetic model with a pharmacodynamics component of parasite killing 5Bertrand and Mentré200846Kay et al.201343Johnston et al.2014106Winter and Hastings2011) and to track multiple parasite genotypes to which we could assign fitness costs and drug susceptibility (i.e. pharmacodynamics) properties. We then built an emulator-based workflow to quantify, through a series of global sensitivity analyses, the influence of multiple factors on the establishment and spread of parasites having different degrees of resistance to artemisinin derivatives and/or their partner drugs when used in monotherapy and combination (as ACTs). Emulators are predictive models that can approximate the relationship between input and output parameters of complex models and can run much faster than complex models to perform global sensitivity analyses more efficiently 30Grow and Hilton2018. OpenMalaria is a mechanistic model, so the observed dynamics at the population level (e.g. the spread of resistant genotypes) emerges from the relationship between the different model components and their input parameters. These dynamics can only be understood and tested through extensive analyses as undertaken here. Identifying which factors (e.g. drug properties and/or setting characteristics) favour the evolution of resistance, enables us to identify drug properties or strategies to slow or mitigate resistance and guides the development and implementation of more sustainable therapies.

Results

Development of drug resistance

We investigated the establishment and spread of drug-resistant genotypes by varying the degrees of resistance for three different treatment profiles. The first treatment profile considered was a monotherapy using a short-acting drug. The short-acting drug had a short half-life and a high killing efficacy typical of artemisinin derivatives (Figure 1A and B). Patients received a daily dose of the short-acting drug for 6 days (see Materials and methods). To mimic the mechanism of resistance to artemisinin derivatives, we assumed that genotypes resistant to the short-acting drug had lower maximum killing rates (Emax) than sensitive ones (Figure 1B) (see Materials and methods). We defined the degree of resistance to the short-acting drug as the relative decrease of the Emax of the resistant genotype compared with the sensitive one. The second treatment profile was also a monotherapy but with a long-acting drug. The long-acting drug had a longer half-life and a lower Emax than the short-acting drug, typical of partner drugs used for ACTs (such as mefloquine, piperaquine, and lumefantrine) (Figure 1A and B). Patients received a daily dose of the long-acting drug for 3 days (see Materials and methods). We assumed that genotypes resistant to the long-acting drug had higher half-maximal effective concentrations (EC50) than sensitive ones (Figure 1B) (see Materials and methods). We defined the degree of resistance to the long-acting drug as the relative increase of the EC50 of the resistant genotype compared with the sensitive genotype. Note that monotherapies for malaria are no longer recommended, but we investigated drivers of resistance under monotherapy to identify the determinants specific to each drug profile. The last treatment profile was a daily dose of a combination of the short-acting and the long-acting drugs for 3 days, simulating ACTs. In this case, we focused on resistance to the short-acting drug, as artemisinin is the shared compound of all ACTs and is of greater concern. Thus, the resistant and sensitive genotypes refer to sensitivity to the short-acting drug. We measured selection for resistance to the short-acting drug against a background of differing sensitivity to the long-acting partner drug, whose effectiveness was varied as described in Table 1. We assumed that the genotypes sensitive and resistant to the short-acting drug had identical sensitivities to the partner drug (i.e. there was no cross-resistance).

Figure 1.

Our analysis had two steps. First, we quantified the impact of factors listed in Table 1 on the spread of drug-resistant parasites through global sensitivity analyses using an emulator trained on our model simulations (Figure 1C, purple area [right side], see Materials and methods). For each simulation, we tracked a drug-sensitive genotype and a drug-resistant genotype, and we estimated the rate of spread using the selection coefficient, which measures the rate at which the logit of the resistant genotype frequency increases each parasite generation (see Materials and methods, note that a selection coefficient below zero implies that resistance does not spread in the population) 38Hastings et al.2020. Then, we assessed the probability of establishment for a subset of resistant genotypes with known and positive selection coefficients to observe the relationship between selection coefficient and the probability of establishment in different settings (Figure 1C, orange area [left side], see Materials and methods). We could then extrapolate the probability of establishing any mutations with a known selection coefficient, which made the process more efficient since estimating the probability of establishment requires running many more stochastic realisations than estimating the selection coefficient due to the stochasticity of this step.

Key drivers of the spread of drug-resistant parasites

Under monotherapy, access to treatment (the probability of symptomatic cases to receive treatment within 2 weeks from the onset of symptoms) and degree of resistance of a monotherapy were the main drivers of the spread of resistance (Figure 2A). For the short-acting and the long-acting drugs used as monotherapy, the selection coefficient increased with increasing access to treatment (Figure 2—figure supplement 1). In addition, higher degrees of resistance of the resistant genotype to the short-acting drug (relative decrease in the resistant genotype Emax compared with the sensitive one) and the long-acting drug (relative increase in the resistant genotype EC50 compared with the sensitive one) promoted the spread of parasites resistant to the short-acting and the long-acting drugs, respectively (Figure 2—figure supplement 1).

#library("ggh4x")
library("ggplot2")
library("gridExtra")
library("cowplot")

### Figure 2A ####
# load the data
data <- read.csv(file = "data/Figure2-Sourcedata1.csv", header  = TRUE)

# Order the level of the factors (in the order of the highest sobol indices to lower one)
data$Factor <- factor(data$Factor, levels = c(
  "Diangostic",
  "Maximum killing rate of drug A",
  "Maximum killing rate of drug B",
  "EIR",
  "Cmax/IC50 of drug A",
  "Cmax/IC50 of drug B",
  "Half-life of drug A",
  "Half-life of drug B",
  "Fitness",
  "Resistance level of drug A",
  "Resistance level of drug B",
  "Access"
))

# Make drug variable as a factor for visualisation
data$drug <- factor(data$drug, levels = c("Drug A", "Drug B", "Drug A + Drug B"))

# Select only the non-seasonal setting
data_2 <- data[data$Setting == "Spread_sesonality1", ]

# define the break for the sobol indices on the y values
break_y <- c(0, 0.25, 0.5, 0.75, 1)
Label_yy <- c(0, 0.25, 0.5, 0.75, 1)

# Define the label for each drug archetype
Dr.labs <- c("Drug A", "Drug B", "Drug A + B")
names(Dr.labs) <- c("A", "B", "A+B")

# select only the first order indicies
data_3 <- data_2[data_2$Effect == "First", ]

constant<-2
# visualise as a columns
PA <- ggplot(data_3, aes(x = drug, y = First, fill = Factor)) +
  geom_col(color = "black", width = 0.6) +
  xlab("") +
  ylab("First-order indices") +
  scale_fill_manual(
    values = c(
      "#126429",
      "#5087C1", "#273871",
      "#999933",
      "#E6959F", "#AA4499",
      "#88CCEE", "#009E73",
      "#DDCC77",
      "#882255", "#661100",
      "#888888"),
    name = "Factors:",
    breaks = c("Access", "Resistance level of drug A", "Resistance level of drug B", "Fitness", "Half-life of drug A", "Half-life ofdrug B", "Cmax/IC50 of drug A", "Cmax/IC50 of drug B", "EIR", "Maximum killing rate of drug A", "Maximum killing rate of drug B", "Diangostic"),
    labels = c("Access to treatment (%)", 
               "Degree of resistance to\nthe short-acting drug", 
               "Degree of resistance to\nthe long-acting drug", 
               "Fitness cost", 
               "Half-life of the\nshort-acting drug (days)", 
               "Half-life of the\nlong-acting drug (days)", 
               "Cmax/EC50 of the\nshort-acting drug", 
               "Cmax/EC50 of the\nlong-acting drug", 
               "EIR (inoculations per person per year)", 
               "Emax of the\nshort-acting drug (per day)", 
               "Emax of the\nlong-acting drug (per day)", 
               "Diagnostic detection limit\n(parasites/microliter)")) +
  scale_x_discrete(labels = c("Drug A" = "Short-acting\ndrug", "Drug B" = "Long-acting\ndrug", "Drug A + Drug B" = "Short-acting + \nLong-acting drugs")) +
  theme_bw() +
  theme(
    axis.text.x = element_text(size = 15/constant, colour = "black"),
    axis.text.y = element_text(size = 15/constant),
    axis.title.x = element_text(size = 15/constant, face = "bold"),
    axis.title.y = element_text(size = 15/constant, face = "bold"),
    plot.title = element_text(size = 20/constant, hjust = 0.5, face = "bold")) +
  theme(legend.text = element_text(size = 15/constant)) +
  theme(legend.title = element_text(size = 15/constant, face = "bold")) +
  ggtitle(label = "")+
  theme(legend.key.size = unit(0.3, "cm"))

### Figure 2B ####
# Load the data
Quantil_final_final <- read.csv(file = "data/Figure2-Sourcedata2.csv", header  = TRUE)

Quantil_final_final <- Quantil_final_final[, c(1:5, 10)]

# add data about the drug archeytpe
Quantil_final_final$Drug <- "Short-acting + Long-acting drugs"

# Define the break for the y axis
break_y <- c(1, 10, 20, 30, 39)
Label_yy <- c("Min", "", "", "", "Max")

#---- visualise three most important factor in the non seasonal setting ----

# select the data
Quantil_final_2 <- Quantil_final_final[Quantil_final_final$Seasonality == "sesonality1", ]
Quantil_final_2 <- Quantil_final_2[Quantil_final_2$G == "Access" | Quantil_final_2$G == "Resistance_Level_long" | Quantil_final_2$G == "Resistance_Level", ]


constant<-2
# visualise
PB <- ggplot(data = Quantil_final_2) +
  geom_line(aes(x = x, y = M, color = G), size = 2/constant) +
  geom_ribbon(aes(x = x, ymin = L, ymax = U, fill = G), alpha = 0.1) +
  facet_grid(. ~ Drug) +
  theme_bw() +
  scale_y_continuous(name = "Selection coefficient\n(resistance to the short-acting drug)") +
  scale_x_continuous(name = "Factor values", breaks = break_y, labels = Label_yy) +
  scale_color_manual(
    name = "Factors:",
    values = c(
      "#126429",
      "#709FCD", "#273871"),
    breaks = c("Access", "Resistance_Level", "Resistance_Level_long", "half_life_long", "Fitness", "eir", "half_life_short", "C_max_IC50"), 
    labels = c("Access to treatment (%)", "Degree of resistance to the short-acting drug", "Degree of resistance to the long-acting drug", "\nHalf-life\nof drug B", "Fitness cost", "EIR", "\nHalf-life\nof drug A", "Cmax/IC50\nof drug B")) +
  scale_fill_manual(
    name = "Factors:",
    values = c(
      "#117733",
      "#6699CC", "#332288"),
    breaks = c("Access", "Resistance_Level", "Resistance_Level_long", "half_life_long", "Fitness", "eir", "half_life_short", "C_max_IC50"), 
    labels = c("Access to treatment (%)", "Degree of resistance to the short-acting drug", "Degree of resistance to the long-acting drug", "\nHalf-life\nof drug B", "Fitness cost", "EIR", "\nHalf-life\nof drug A", "Cmax/IC50\nof drug B")) +
  theme(
    axis.text.x = element_text(size = 15/2),
    axis.text.y = element_text(size = 15/2),
    axis.title.x = element_text(size = 15/2, face = "bold"),
    axis.title.y = element_text(size = 15/2, face = "bold"),
    plot.title = element_text(size = 20/2, hjust = 0.5, face = "bold")) +
  theme(legend.text = element_text(size = 15/2)) +
  theme(legend.title = element_text(size = 15/2, face = "bold")) +
  theme(
    strip.text.x = element_text(size = 15/2, color = "black", face = "bold"),
    strip.text.y = element_text(size = 15/2, color = "black", face = "bold")) +
  # theme(legend.position = "bottom", legend.direction = "vertical")+
  theme(legend.key.size = unit(0.3, "cm"))

grid.arrange(PA,PB)

# Load the data
Quantil_final_final <- read.csv(file = "data/Figure2-figuresupplement1-Sourcedata1.csv", sep = ",", header = T)

# select non seasonal setting
Quantil_final_2 <- Quantil_final_final[Quantil_final_final$Seasonality == "sesonality1",]

# select factor of interest
Quantil_final_2 <-  Quantil_final_2[Quantil_final_2$G == "Access" |
                                      Quantil_final_2$G == "Resistance level of drug B" |
                                      Quantil_final_2$G == "Resistance level of drug A",]

# Define a label for lower value of each factor and one for highest value
break_y <- c(1, 10, 20, 30, 39)
Label_yy <- c("Min", "", "", "", "Max")

# Define label for each drug type
Dr.labs <- c("Short-acting drug", "Long-acting drug")
names(Dr.labs) <- c("Drug A", "Drug B")

constant<-2
# Visualise the results
  ggplot(data = Quantil_final_2) +
  geom_line(aes(x = x, y = M, color = G), size = 2/constant) +
  geom_ribbon(aes(
    x = x,
    ymin = L,
    ymax = U,
    fill = G)
    , alpha = 0.15) +
  facet_grid(. ~ drug,  labeller = labeller(drug = Dr.labs)) +
  theme_bw() +
  ggtitle("") +
  scale_y_continuous(name = "Selection coefficient") +
  scale_x_continuous(name = "Factor values",
                     breaks = break_y,
                     labels = Label_yy) +
  theme(
    axis.text.x = element_text(size = 16/constant),
    axis.text.y = element_text(size = 16/constant),
    axis.title.x = element_text(size = 18/constant, face = "bold"),
    axis.title.y = element_text(size = 18/constant, face = "bold"),
    plot.title = element_text(size = 18/constant, hjust = 0.5, face = "bold")) +
  theme(legend.text = element_text(size = 18/constant)) +
  theme(legend.title = element_text(size = 18/constant, face = "bold")) +
  scale_color_manual(
    name = "Factors:",
    values = c("#126429", "#709FCD", "#273871"),
    breaks = c(
      "Access",
      "Resistance level of drug A",
      "Resistance level of drug B",
      "half_life",
      "Fitness",
      "eir",
      "C_max_IC50"),
    labels = c(
      "Access to treatment (%)",
      "Degree of resistance to the short-acting drug",
      "Degree of resistance to the long-acting drug",
      "Half-life",
      "Fitness cost",
      "EIR",
      "Cmax/IC50")) +
  scale_fill_manual(
    name = "Factors:",
    values = c("#117733", "#6699CC", "#332288"),
    breaks = c(
      "Access",
      "Resistance level of drug A",
      "Resistance level of drug B",
      "half_life",
      "Fitness",
      "eir",
      "C_max_IC50"),
    labels = c(
      "Access to treatment (%)",
      "Degree of resistance to the short-acting drug",
      "Degree of resistance to the long-acting drug",
      "Half-life",
      "Fitness cost",
      "EIR",
      "Cmax/IC50")) +
  ggtitle(label = "") +
  theme(
    strip.text.x = element_text(size = 18/constant, color = "black", face = "bold"),
    strip.text.y = element_text(size = 18/constant, color = "black", face = "bold"))+
  theme(legend.position = "bottom", legend.direction = "vertical")+
  theme(legend.key.size = unit(1/constant, "cm")) 

When the short-acting and the long-acting drugs were used in combination in our model, we referred to the resistant and sensitive genotypes as the genotypes resistant and sensitive to the short-acting drug, respectively. However, both genotypes could have some degree of resistance to the long-acting drug. In this case, the most important driver of spread was the degree of resistance of both genotypes to the long-acting drug (Figure 2A). The median selection coefficient was below zero when both genotypes were susceptible to the long-acting drug (the minimum degree of resistance to the long-acting drug) (Figure 2B), indicating that using an efficient partner drug can limit the spread of artemisinin resistance. The spread of parasites resistant to the short-acting drug was accelerated when parasites were also resistant to the long-acting drug, highlighting that resistance to the long-acting drug can facilitate the spread of artemisinin resistance. We further illustrated with concrete examples (Appendix: section 1.1) how the spread of partial resistance to the short-acting drug accelerates with higher degrees of resistance to the long-acting drug. These results further confirmed that resistance to partner drugs facilitates the spread of resistance to artemisinin, highlighting the importance of combining artemisinin derivatives with an efficient partner drug.

Variation in the influence of factors across settings and degrees of resistance

We compared the effects of drug properties and levels of fitness cost on estimated selection coefficients for a fixed set of degrees of resistance, levels of access to treatment, transmission intensities, seasonality patterns, and levels of adherence to treatment (percentage of treatment doses adhered by patients). Figure 3 summarises the impact of key factors influencing estimated selection coefficients in seasonal transmission settings with a population fully adherent to treatment (the impact of factors was similar across seasonality pattern and levels of adherence to treatment Figure 3—figure supplements 1–!number(2)). The impact of all factors in each setting is shown in Figure 3—figure supplements 1–!number(2).

# load the data and add the drug name
Quantil_final_final<-read.csv(file = "data/Figure3-Sourcedata1.csv", sep = ",", header = T)


# Transform constrained variable into a factor
Quantil_final_final$drug <- factor(Quantil_final_final$drug, levels = c("A", "B", "A+B"))
Quantil_final_final$Dosage <- factor(Quantil_final_final$Dosage, levels = c("1", "0"))
Quantil_final_final$EIR <- factor(Quantil_final_final$EIR, levels = c("5", "10", "500"))

# Creat a label for each constrained variable
S.labs <- c("No seasonality", "Seasonality")
names(S.labs) <- c("sesonality1", "sesonality2")
T.labs <- c("Low access to treatment", "High access to treatment")
names(T.labs) <- c("0.04", "0.5")
D.labs <- c("High adherence\n to treatment", "Low adherence\n to treatment")
names(D.labs) <- c("1", "0")
Dr.labs <- c("Short-acting\ndrug", "Long-acting\ndrug", "Short-acting +\n Long-acting drugs")
names(Dr.labs) <- c("A", "B", "A+B")
Dr.labs_2 <- c("", "", "")
names(Dr.labs_2) <- c("A", "B", "A+B")
R.labs <- c("Low degree\n of resistance", "Low degree\n of resistance", "High degree\n of resistance", "High degree\n of resistance")
names(R.labs) <- c("7", "2.5", "18", "10")
R.labs_2 <- c("", "", "", "")
names(R.labs_2) <- c("7", "2.5", "18", "10")

# creat the label for the y axis
break_y <- c(1, 10, 20, 30, 39)
Label_yy <- c("Min", "", "", "", "Max")

# select the data
Quantil_final_2 <- Quantil_final_final[Quantil_final_final$G == "Fitness" | Quantil_final_final$G == "half_life_short" | Quantil_final_final$G == "C_max_IC50_short" | Quantil_final_final$G == "half_life_long" | Quantil_final_final$G == "C_max_IC50_long", ]
Quantil_final_3 <- Quantil_final_2[Quantil_final_2$EIR == 500 | Quantil_final_2$EIR == 5, ]
Quantil_final_3 <- Quantil_final_3
#Quantil_final_3 <- Quantil_final_3[Quantil_final_3$Access == 0.5 & Quantil_final_3$Seasonality == "sesonality1" | Quantil_final_3$Access == 0.04 & Quantil_final_3$Dosage == 1 & Quantil_final_3$Seasonality == "sesonality1" | Quantil_final_3$Access == 0.5 & Quantil_final_3$Dosage == 1 & Quantil_final_3$Seasonality == "sesonality2", ]
#Quantil_final_3 <- Quantil_final_3[Quantil_final_3$Resistance_level == 2.5 | Quantil_final_3$Resistance_level == 7, ]

Quantil_final_3$x[Quantil_final_3$G == "Fitness"] <- (Quantil_final_3$x[Quantil_final_3$G == "Fitness"] - 40) * -1
Quantil_final_3$x[Quantil_final_3$G == "C_max_IC50_short" & Quantil_final_3$drug == "A"] <- (Quantil_final_3$x[Quantil_final_3$G == "C_max_IC50_short" & Quantil_final_3$drug == "A"] - 40) * -1
Quantil_final_3$x[Quantil_final_3$G == "C_max_IC50_short" & Quantil_final_3$drug == "A+B"] <- (Quantil_final_3$x[Quantil_final_3$G == "C_max_IC50_short" & Quantil_final_3$drug == "A+B"] - 40) * -1



Quantil_final_4 <- Quantil_final_3[Quantil_final_3$Seasonality == "sesonality2", ]
Quantil_final_4 <- Quantil_final_4[Quantil_final_4$Dosage == 1, ]
Quantil_final_4<-Quantil_final_4[Quantil_final_4$G=="half_life_short" | Quantil_final_4$G=="Fitness" |Quantil_final_4$G=="C_max_IC50_long" |Quantil_final_4$G=="half_life_long",]
Quantil_final_4$G<-as.factor(Quantil_final_4$G)

# Do the plot for drug A
Quantil_final_4_1 <- Quantil_final_4[Quantil_final_4$drug == "A", ]
Quantil_final_4_1<-Quantil_final_4_1[Quantil_final_4_1$G=="half_life_short" & Quantil_final_4_1$Access==0.5 | Quantil_final_4_1$G=="Fitness" & Quantil_final_4_1$Access==0.04,]

break_y <- c(1, 10, 20, 30, 39)
Label_y <- c("Min", "", "", "", "Max")

Label_yy_1 <- c("", "", "", "", "")
constant<-2.5
pl11 <- ggplot(data = Quantil_final_4_1) +
  geom_line(aes(x = x, y = M, color = G, linetype = EIR), size = 1.9/constant, alpha = 1) +
  # facet_nested(drug~Access+Seasonality+Dosage,labeller=labeller(drug= Dr.labs,Access=T.labs, Dosage=D.labs, Seasonality=S.labs,Resistance_level=R.labs)) +
  facet_nested(drug ~ Access + Resistance_level, labeller = labeller(drug = Dr.labs, Access = T.labs, Dosage = D.labs, Seasonality = S.labs, Resistance_level = R.labs), scales = "free", independent = "x") +
  theme_bw() +
  ggtitle("") +
  scale_y_continuous(name = "Selection coefficient",limits=c(-0.125,0.6)) +
  scale_x_continuous(name = "",breaks = break_y, labels = Label_y) +
  theme(
    axis.text.x = element_text(size = 15/constant),
    axis.text.y = element_text(size = 15/constant),
    axis.title.x =element_text(size = 16/constant, face = "bold"),
    axis.title.y = element_text(size = 16/constant, face = "bold"),
    plot.title = element_text(size = 18/constant, hjust = 0.5, face = "bold")) +
  theme(legend.text = element_text(size = 12/constant)) +
  theme(legend.title = element_text(size = 12/constant, face = "bold")) +
  theme(legend.margin = margin(0,0,0,0, unit="cm"))+
  guides(
    linetype = guide_legend(keywidth = 5/constant, keyheight = 0.05/constant, override.aes = list(size = 2.5/constant)),
    colour = guide_legend(keywidth = 2/constant, keyheight = 0.05/constant, override.aes = list(size = 5/constant)),
    legend.spacing.y = unit(-1500, "cm")) +
  scale_color_manual(
    name = "Factors:",
    values = c(
      "#999933",
      "#CC6677", "#AA4499",
      "#88CCEE", "#44AA99","#44AA99","#44AA99","#44AA99"),
    breaks = c("Fitness", "half_life_short", "half_life_long", "C_max_IC50_long"),
    labels = c("Fitness cost", "Half-life (days) of the\nshort-acting drug", "Half-life (days) of the\nlong-acting drug", "Cmax/EC50 of the\nlong-acting drug"), 
    drop = FALSE) +
  scale_linetype_manual(
    values = c("solid", "solid", "dashed"),
    name = "EIR:",
    breaks = c("5", "10", "500"),
    labels = c("5", "10", "500")) +
  ggtitle(label = "") +
  theme(
    strip.text.x = element_text(size = 16/constant, color = "black", face = "bold"),
    strip.text.y = element_text(size = 16/constant, color = "black", face = "bold")) +
  theme(plot.margin = unit(c(-1, 0.2, 1, 0.2)/constant, "cm")) +
  geom_hline(yintercept = 0, linetype = "dotted", color = "#999999", size = 1.5) +
  theme(legend.position="top", legend.box="vertical")  +
  guides(fill=guide_legend(nrow=2))

#pl11

# do the plot for drug B
Quantil_final_4_2 <- Quantil_final_4[Quantil_final_4$drug == "B", ]
Quantil_final_4_2<-Quantil_final_4_2[Quantil_final_4_2$G=="half_life_long" & Quantil_final_4_2$Access==0.5  & Quantil_final_4_2$Resistance_level==2.5 | Quantil_final_4_2$G=="C_max_IC50_long" & Quantil_final_4_2$Access==0.5  & Quantil_final_4_2$Resistance_level==10| Quantil_final_4_2$G=="Fitness" & Quantil_final_4_2$Access==0.04,]

# adjust the label
Label_yy_1 <- c("", "", "", "", "")
S.labs2 <- c("", "")
names(S.labs2) <- c("sesonality1", "sesonality2")
T.labs2 <- c("", "")
names(T.labs2) <- c("0.04", "0.5")
D.labs2 <- c("", "")
names(D.labs2) <- c("1", "0")

pl12 <- ggplot(data = Quantil_final_4_2) +
  geom_line(aes(x = x, y = M, color = G, linetype = EIR), size = 1.9/constant, alpha = 1) +
  # facet_nested(drug~Access+Seasonality+Dosage,labeller=labeller(drug= Dr.labs,Access=T.labs, Dosage=D.labs, Seasonality=S.labs,Resistance_level=R.labs)) +
  facet_nested(drug  ~ Access  + Resistance_level, labeller = labeller(drug = Dr.labs, Access = T.labs2, Dosage = D.labs, Seasonality = S.labs2, Resistance_level = R.labs)) +
  theme_bw() +
  ggtitle("") +
  scale_y_continuous(name = "Selection coefficient", limits=c(-0.125,0.6)) +
  scale_x_continuous(name = "", breaks = break_y, labels = Label_yy) +
  theme(
    axis.text.x =  element_text(size = 15/constant),
    axis.text.y = element_text(size = 15/constant),
    axis.title.x = element_text(size = 16/constant, face = "bold"),
    axis.title.y = element_text(size = 16/constant, face = "bold"),
    plot.title = element_text(size = 18/constant, hjust = 0.5, face = "bold"),) +
  scale_color_manual(
    name = "Factors:",
    values = c(
      "#999933",
      "#CC6677", "#AA4499",
      "#88CCEE", "#44AA99"),
    breaks = c("Fitness", "half_life_short", "half_life_long", "C_max_IC50_short", "C_max_IC50_long"),
    labels = c("Fitness cost", "Half-life drug A (days)", "Half-life drug B (day)", "Cmax/EC50 drug A", "Cmax/EC50 drug B")) +
  scale_linetype_manual(
    values = c("solid", "solid", "dashed"),
    name = "EIR:",
    breaks = c("5", "10", "500"),
    labels = c("5", "10", "500")) +
  ggtitle(label = "") +
  theme(
    strip.text.x = element_blank(),
    strip.text.y = element_text(size = 16/constant, color = "black", face = "bold")) +
  theme(plot.margin = unit(c(-1, 0.2, 1, 0.2)/constant, "cm")) +
  theme(strip.background.x = element_rect(fill = "white", colour = "white")) +
  geom_hline(yintercept = 0, linetype = "dotted", color = "#999999", size = 1.5) +
  theme(legend.position = "none")

#pl12

# do the plot for drug A and B
Quantil_final_4_3 <- Quantil_final_4[Quantil_final_4$drug == "A+B", ]
Quantil_final_4_3<-Quantil_final_4_3[Quantil_final_4_3$G=="C_max_IC50_long" & Quantil_final_4_3$Access==0.5 | Quantil_final_4_3$G=="Fitness" & Quantil_final_4_3$Access==0.04,]

# adjust the label
Label_yy_1 <- c("", "", "", "", "")
S.labs2 <- c("", "")
names(S.labs2) <- c("sesonality1", "sesonality2")
T.labs2 <- c("", "")
names(T.labs2) <- c("0.04", "0.5")
D.labs2 <- c("", "")
names(D.labs2) <- c("1", "0")

pl13 <- ggplot(data = Quantil_final_4_3) +
  geom_line(aes(x = x, y = M, color = G, linetype = EIR), size = 1.9/constant, alpha = 1) +
  # facet_nested(drug~Access+Seasonality+Dosage,labeller=labeller(drug= Dr.labs,Access=T.labs, Dosage=D.labs, Seasonality=S.labs,Resistance_level=R.labs)) +
  facet_nested(drug  ~ Access + Seasonality + Resistance_level, labeller = labeller(drug = Dr.labs, Access = T.labs2, Dosage = D.labs, Seasonality = S.labs2, Resistance_level = R.labs)) +
  theme_bw() +
  ggtitle("") +
  scale_y_continuous(name = "Selection coefficient",limits=c(-0.125,0.6)) +
  scale_x_continuous(name = "", breaks = break_y, labels = Label_yy) +
  theme(
    axis.text.x = element_text(size = 15/constant),
    axis.text.y = element_text(size = 15/constant),
    axis.title.x = element_text(size = 16/constant, face = "bold", hjust = 0.3),
    axis.title.y = element_text(size = 16/constant, face = "bold"),
    plot.title = element_text(size = 18, hjust = 0.75, face = "bold")) +
  scale_color_manual(
    name = "Factors:",
    values = c(
      "#999933",
      "#CC6677", "#AA4499",
      "#88CCEE", "#44AA99"),
    breaks = c("Fitness", "half_life_short", "half_life_long", "C_max_IC50_short", "C_max_IC50_long"),
    labels = c("Fitness cost", "Half-life drug A (day)", "Half-life drug B (day)", "Cmax/EC50 drug A", "Cmax/EC50 drug B")) +
  scale_linetype_manual(
    values = c("solid", "solid", "dashed"),
    name = "EIR:",
    breaks = c("5", "10", "500"),
    labels = c("5", "10", "500")) +
  ggtitle(label = "") +
  theme(
    strip.text.x = element_blank(),
    strip.text.y = element_text(size = 16/constant, color = "black", face = "bold")) +
  theme(plot.margin = unit(c(-2, 0.2, 1, 0.2)/constant, "cm")) +
  theme(strip.background.x = element_rect(fill = "white", colour = "white")) +
  geom_hline(yintercept = 0, linetype = "dotted", color = "#999999", size = 1.5) +
  theme(legend.position = "none")

#pl13

# merge all the plots
Test_1 <- plot_grid(pl11, pl12, pl13,
                    ncol = 1, nrow = 3, rel_heights = c(1.55, 1, 1), scale = c(1, 1, 1))

Test_1<-Test_1+ 
  draw_label("Fitness cost", x = 0.12, y = 0.0308,size=16/2.5, hjust = 0, vjust = 0, fontface ="bold")+
  draw_label("Fitness cost", x = 0.12, y = 0.3125,size=16/2.5, hjust = 0, vjust = 0, fontface ="bold")+
  draw_label("Fitness cost", x = 0.12, y = 0.5942,size=16/2.5, hjust = 0, vjust = 0, fontface ="bold")+
  
  draw_label("Fitness cost", x = 0.3475, y = 0.0308,size=16/2.5, hjust = 0, vjust = 0, fontface ="bold")+
  draw_label("Fitness cost", x = 0.3475, y = 0.3125,size=16/2.5, hjust = 0, vjust = 0, fontface ="bold")+
  draw_label("Fitness cost", x = 0.3475, y = 0.5942,size=16/2.5, hjust = 0, vjust = 0, fontface ="bold")+
  
  draw_label("Cmax/EC50 of the\nlong-acting drug", x = 0.554, y = 0.0165,size=16/2.5, hjust = 0, vjust = 0, fontface ="bold")+
  draw_label("Half-life (days) of the\nlong-acting drug", x = 0.541, y =  0.2982,size=16/2.5, hjust = 0, vjust = 0, fontface ="bold")+
  draw_label("Half-life (days) of the\nshort-acting drug", x = 0.541, y = 0.5799,size=16/2.5, hjust = 0, vjust = 0, fontface ="bold")+
  
  draw_label("Cmax/EC50 of the\nlong-acting drug", x = 0.7825, y = 0.0165,size=16/2.5, hjust = 0, vjust = 0, fontface ="bold")+
  draw_label("Cmax/EC50 of the\nlong-acting drug", x = 0.7825, y =  0.2982,size=16/2.5, hjust = 0, vjust = 0, fontface ="bold")+
  draw_label("Half-life (days) of the\nshort-acting drug", x = 0.769, y = 0.5799,size=16/2.5, hjust = 0, vjust = 0, fontface ="bold")

Test_1 <- plot_grid(pl11, pl12, pl13,
                    ncol = 1, nrow = 3, rel_heights = c(1.8, 1, 1), scale = c(1, 1, 1))
Test_1<-Test_1+ 
  draw_label("Fitness cost", x = 0.118, y = 0.04275,size=16/2.5, hjust = 0, vjust = 0, fontface ="bold")+
  draw_label("Fitness cost", x = 0.118, y = 0.305,size=16/2.5, hjust = 0, vjust = 0, fontface ="bold")+
  draw_label("Fitness cost", x = 0.118, y = 0.5673,size=16/2.5, hjust = 0, vjust = 0, fontface ="bold")+
  
  draw_label("Fitness cost", x = 0.336, y = 0.04275,size=16/2.5, hjust = 0, vjust = 0, fontface ="bold")+
  draw_label("Fitness cost", x = 0.336, y = 0.305,size=16/2.5, hjust = 0, vjust = 0, fontface ="bold")+
  draw_label("Fitness cost", x = 0.336, y = 0.5673,size=16/2.5, hjust = 0, vjust = 0, fontface ="bold")+
  
  draw_label("Cmax/EC50 of the\nlong-acting drug", x = 0.528, y = 0.024,size=16/2.5, hjust = 0, vjust = 0, fontface ="bold")+
  draw_label("Half-life (days) of the\nlong-acting drug", x = 0.5125, y =  0.28625,size=16/2.5, hjust = 0, vjust = 0, fontface ="bold")+
  draw_label("Half-life (days) of the\nshort-acting drug", x = 0.5125, y = 0.549,size=16/2.5, hjust = 0, vjust = 0, fontface ="bold")+
  
  draw_label("Cmax/EC50 of the\nlong-acting drug", x = 0.7465, y = 0.024,size=16/2.5, hjust = 0, vjust = 0, fontface ="bold")+
  draw_label("Cmax/EC50 of the\nlong-acting drug", x = 0.7465, y =  0.28625,size=16/2.5, hjust = 0, vjust = 0, fontface ="bold")+
  draw_label("Half-life (days) of the\nshort-acting drug", x = 0.73, y = 0.549,size=16/2.5, hjust = 0, vjust = 0, fontface ="bold")

Test_1

# Transform constrained variable into a factor
Quantil_final_final$drug <- factor(Quantil_final_final$drug, levels = c("A", "B", "A+B"))
Quantil_final_final$Dosage <- factor(Quantil_final_final$Dosage, levels = c("1", "0"))
Quantil_final_final$EIR <- factor(Quantil_final_final$EIR, levels = c("5", "10", "500"))

# Creat a label for each constrained variable
S.labs <- c("No seasonality", "Seasonality")
names(S.labs) <- c("sesonality1", "sesonality2")
T.labs <- c("Low access to treatment", "High access to treatment")
names(T.labs) <- c("0.04", "0.5")
D.labs <- c("High adherence\n to treatment", "Low adherence\n to treatment")
names(D.labs) <- c("1", "0")
Dr.labs <- c("Short-acting drug", "Long-acting drug", "Short-acting +\nLong-acting drugs")
names(Dr.labs) <- c("A", "B", "A+B")
Dr.labs_2 <- c("", "", "")
names(Dr.labs_2) <- c("A", "B", "A+B")
R.labs <- c("Low degree\n of resistance", "Low degree\n of resistance", "High degree\n of resistance", "High degree\n of resistance")
names(R.labs) <- c("7", "2.5", "18", "10")
R.labs_2 <- c("", "", "", "")
names(R.labs_2) <- c("7", "2.5", "18", "10")

# creat the label for the y axis
break_y <- c(1, 10, 20, 30, 39)
Label_yy <- c("Min", "", "", "", "Max")

# Select the data
Quantil_final_2 <- Quantil_final_final[Quantil_final_final$G == "Fitness" | Quantil_final_final$G == "half_life_short" | Quantil_final_final$G == "C_max_IC50_short" | Quantil_final_final$G == "half_life_long" | Quantil_final_final$G == "C_max_IC50_long", ]
Quantil_final_3 <- Quantil_final_2[Quantil_final_2$EIR == 500 | Quantil_final_2$EIR == 5, ]
Quantil_final_3 <- Quantil_final_3

# Transform IC50 and Fitness in the good direction
Quantil_final_3$x[Quantil_final_3$G == "Fitness"] <- (Quantil_final_3$x[Quantil_final_3$G == "Fitness"] - 40) * -1
Quantil_final_3$x[Quantil_final_3$G == "C_max_IC50_short" & Quantil_final_3$drug == "A"] <- (Quantil_final_3$x[Quantil_final_3$G == "C_max_IC50_short" & Quantil_final_3$drug == "A"] - 40) * -1
Quantil_final_3$x[Quantil_final_3$G == "C_max_IC50_short" & Quantil_final_3$drug == "A+B"] <- (Quantil_final_3$x[Quantil_final_3$G == "C_max_IC50_short" & Quantil_final_3$drug == "A+B"] - 40) * -1

#----- Plot for high level of access to treatment (Figure 3) ----
# plot all high acess to treatment
Quantil_final_4 <- Quantil_final_3[Quantil_final_3$Access == 0.5, ]
Quantil_final_4 <- Quantil_final_4[!(Quantil_final_4$drug == "A+B" & Quantil_final_4$G == "C_max_IC50_short"), ]
Quantil_final_4 <- Quantil_final_4[!(Quantil_final_4$drug == "A+B" & Quantil_final_4$G == "half_life_short"), ]

constant<-2
ggplot(data = Quantil_final_4) +
  geom_line(aes(x = x, y = M, color = G, linetype = EIR), size = 1.9/constant, alpha = 1) +
  # facet_nested(drug~Access+Seasonality+Dosage,labeller=labeller(drug= Dr.labs,Access=T.labs, Dosage=D.labs, Seasonality=S.labs,Resistance_level=R.labs)) +
  facet_nested(drug + Resistance_level ~  Seasonality + Dosage, labeller = labeller(drug = Dr.labs, Dosage = D.labs, Seasonality = S.labs, Resistance_level = R.labs)) +
  theme_bw() +
  ggtitle("") +
  scale_y_continuous(name = "Selection coefficient") +
  scale_x_continuous(name = "Factor values", breaks = break_y, labels = Label_yy) +
  theme(
    axis.text.x = element_text(size = 15/constant),
    axis.text.y = element_text(size = 15/constant),
    axis.title.x = element_text(size = 16/constant, face = "bold"),
    axis.title.y = element_text(size = 16/constant, face = "bold"),
    plot.title = element_text(size = 20/constant, hjust = 0.5, face = "bold")) +
  theme(legend.text = element_text(size = 12/constant)) +
  theme(legend.title = element_text(size = 12/constant, face = "bold")) +
  scale_color_manual(
    name = "Factors:",
    values = c(
      "#999933",
      "#E6959F", "#AA4499",
      "#88CCEE", "#009E73"),
    breaks = c("Fitness", "half_life_short", "half_life_long", "C_max_IC50_short", "C_max_IC50_long"),
    labels = c("Fitness cost", "Half-life  (days) of the\nshort-acting drug", "Half-life (days) of the\nlong-acting drug", "Cmax/EC50 of the\nshort-acting drug", "Cmax/EC50 of the\nlong-acting drug")) +
  scale_linetype_manual(
    values = c("solid", "solid", "dashed"),
    name = "EIR (inoculations per person per year):",
    breaks = c("5", "10", "500"),
    labels = c("5", "10", "500")) +
  guides(
    linetype = guide_legend(keywidth = 5/constant, keyheight = 0.05/constant, override.aes = list(size = 2.5/constant)),
    colour = guide_legend(keywidth = 2/constant, keyheight = 0.05/constant, override.aes = list(size = 5/constant)),
    legend.spacing.y = unit(150, "cm")) +
  ggtitle(label = "") +
  theme(
    strip.text.x = element_text(size = 16/constant, color = "black", face = "bold"),
    strip.text.y = element_text(size = 16/constant, color = "black", face = "bold")) +
  theme(plot.margin = unit(c(0.05, 0.05,0.05 , 0.05), "cm")) +
  geom_hline(yintercept = 0, linetype = "dotted", color = "#999999", size = 1.5/constant) +
  theme(legend.position="top", legend.box="vertical")  +
  guides(fill=guide_legend(nrow=2))

#----- Plot for high level of access to treatment (Figure 3) ----

# plot all high acess to treatment
Quantil_final_4 <- Quantil_final_3[Quantil_final_3$Access == 0.04, ]
Quantil_final_4 <- Quantil_final_4[!(Quantil_final_4$drug == "A+B" & Quantil_final_4$G == "C_max_IC50_short"), ]
Quantil_final_4 <- Quantil_final_4[!(Quantil_final_4$drug == "A+B" & Quantil_final_4$G == "half_life_short"), ]

constant<-2
ggplot(data = Quantil_final_4) +
  geom_line(aes(x = x, y = M, color = G, linetype = EIR), size = 1.9/constant, alpha = 1) +
  # facet_nested(drug~Access+Seasonality+Dosage,labeller=labeller(drug= Dr.labs,Access=T.labs, Dosage=D.labs, Seasonality=S.labs,Resistance_level=R.labs)) +
  facet_nested(drug + Resistance_level ~  Seasonality + Dosage, labeller = labeller(drug = Dr.labs, Dosage = D.labs, Seasonality = S.labs, Resistance_level = R.labs)) +
  theme_bw() +
  ggtitle("") +
  scale_y_continuous(name = "Selection coefficient") +
  scale_x_continuous(name = "Factor values", breaks = break_y, labels = Label_yy) +
  theme(
    axis.text.x = element_text(size = 15/constant),
    axis.text.y = element_text(size = 15/constant),
    axis.title.x = element_text(size = 16/constant, face = "bold"),
    axis.title.y = element_text(size = 16/constant, face = "bold"),
    plot.title = element_text(size = 20/constant, hjust = 0.5, face = "bold")) +
  theme(legend.text = element_text(size = 12/constant)) +
  theme(legend.title = element_text(size = 12/constant, face = "bold")) +
  scale_color_manual(
    name = "Factors:",
    values = c(
      "#999933",
      "#E6959F", "#AA4499",
      "#88CCEE", "#009E73"),
    breaks = c("Fitness", "half_life_short", "half_life_long", "C_max_IC50_short", "C_max_IC50_long"),
    labels = c("Fitness cost", "Half-life  (days) of the\nshort-acting drug", "Half-life (days) of the\nlong-acting drug", "Cmax/EC50 of the\nshort-acting drug", "Cmax/EC50 of the\nlong-acting drug")) +
  scale_linetype_manual(
    values = c("solid", "solid", "dashed"),
    name = "EIR (inoculations per person per year):",
    breaks = c("5", "10", "500"),
    labels = c("5", "10", "500")) +
  guides(
    linetype = guide_legend(keywidth = 5/constant, keyheight = 0.05/constant, override.aes = list(size = 2.5/constant)),
    colour = guide_legend(keywidth = 2/constant, keyheight = 0.05/constant, override.aes = list(size = 5/constant)),
    legend.spacing.y = unit(150, "cm")) +
  ggtitle(label = "") +
  theme(
    strip.text.x = element_text(size = 16/constant, color = "black", face = "bold"),
    strip.text.y = element_text(size = 16/constant, color = "black", face = "bold")) +
  theme(plot.margin = unit(c(0.05, 0.05,0.05, 0.05), "cm")) +
  geom_hline(yintercept = 0, linetype = "dotted", color = "#999999", size = 1.5/constant) +
  theme(legend.position="top", legend.box="vertical")  +
  guides(fill=guide_legend(nrow=2))

# load the data
data <- read.csv(file = "data/Figure3-figuresupplement3-Sourcedata1.csv", header  = TRUE)

# Transform the constrained variable into factor
data$Resistance_level <- factor(data$Resistance_level, levels = c("7", "18"))
data$Dosage <- factor(data$Dosage, levels = c("4", "0"))
data$EIR <- factor(data$EIR, levels = c("5", "10", "500"))
data$Treatment_access <- factor(data$Treatment_access, levels = c("0.5", "0.04"))

# Creat a label for each constrained variable
T.labs <- c("High access to treatment", "Low access to treatment")
names(T.labs) <- c("0.5", "0.04")
R.labs <- c("High degree\n of resistance", "Low degree\n of resistance")
names(R.labs) <- c("18", "7") 
D.labs <- c("High adherence\n to treatment", "Low adherence\n to treatment")
names(D.labs) <- c("4", "0")
S.labs <- c("No seasonality", "Seasonality")
names(S.labs) <- c("sesonality1", "sesonality2")

# define the break for the y axis
break_y <- c(0, 0.25, 0.5, 0.75, 1)
Label_yy <- c(0, 0.25, 0.5, 0.75, 1)


# select the data
data_2 <- data[data$Effect == "First", ]

# order the level of each factors
data_2$Factor <- factor(data_2$Factor, levels = c("Diagnostic", "MKR", "IC50_S", "half_life", "Fitness"))

constant<-2
# visualise
pd <- position_dodge
  ggplot(data_2, aes(x = EIR, y = First, fill = Factor)) +
  geom_col(color = "black", width = 0.6) +
  facet_nested(Treatment_access + Resistance_level ~ Seasonality+ Dosage, labeller = labeller(Resistance_level = R.labs, Treatment_access = T.labs, Dosage = D.labs, Seasonality = S.labs)) +
  scale_fill_manual(
    name = "Factors:",
    values = c(
      "#999933",
      "#E6959F",
      "#88CCEE",
      "#882255",
      "#888888"),
    breaks = c("Fitness", "half_life", "IC50_S", "MKR", "Diagnostic"),
    labels = c("\nFitness cost\n", "\nHalf-life (days) of the\nshort-acting drug", "\nCmax/EC50 of the\nshort-acting drug", "\nEmax (per day) of the\nshort-acting drug", "\nDiagnostic detection\nlimit (parasites/ul)")) +
  theme(axis.title = element_text(face = "bold")) +
  ylab("First-order indices") +
  xlab("EIR (inoculations per person per year)") +
  theme_bw() +
  theme(
    axis.text.x = element_text(size = 15/constant),
    axis.text.y = element_text(size = 15/constant),
    axis.title.x = element_text(size = 16/constant, face = "bold"),
    axis.title.y = element_text(size = 16/constant, face = "bold"),
    plot.title = element_text(size = 20/constant, hjust = 0.5, face = "bold")) +
  theme(legend.text = element_text(size = 15/constant)) +
  theme(legend.title = element_text(size = 15/constant, face = "bold")) +
  ggtitle(label = "Short-acting drug") +
  theme(
    strip.text.x = element_text(size = 16/constant, color = "black", face = "bold"),
    strip.text.y = element_text(size = 16/constant, color = "black", face = "bold")) +
  scale_y_continuous(breaks = break_y, labels = Label_yy) +
  theme(legend.position="top", legend.box="vertical")  +
  guides(fill=guide_legend(nrow=2))

# load the data
data <- read.csv(file = "data/Figure3-figuresupplement4-Sourcedata1.csv", header  = TRUE)

# transform the constrained variable into factor
data$Resistance_level <- factor(data$Resistance_level, levels = c("2.5", "10"))
data$Dosage <- factor(data$Dosage, levels = c("1", "0"))
data$EIR <- factor(data$EIR, levels = c("5", "10", "500"))
data$Treatment_access <- factor(data$Treatment_access, levels = c("0.5","0.04"))

# creat a lavel for each constrained variable
T.labs <- c("High access to treatment", "Low access to treatment")
names(T.labs) <- c("0.5", "0.04")
R.labs <- c("High degree\n of resistance", "Low degree\n of resistance")
names(R.labs) <- c("10", "2.5") 
D.labs <- c("High adherence\n to treatment", "Low adherence\n to treatment")
names(D.labs) <- c("1", "0")
S.labs <- c("No seasonality", "Seasonality")
names(S.labs) <- c("sesonality1", "sesonality2")

# define the break for the y axis
break_y <- c(0, 0.25, 0.5, 0.75, 1)
Label_yy <- c(0, 0.25, 0.5, 0.75, 1)

# ---- visualize the first order indices of each factors -----

# select the data
data_2 <- data[data$Effect == "First", ]

# define the level of the variable factors
data_2$Factor <- factor(data_2$Factor, levels = c("MKR", "Diangostic", "Cmax/IC50", "half-life", "Fitness"))

constant<-2
# Visualize

  ggplot(data_2, aes(x = EIR, y = First, fill = Factor)) +
  geom_col(color = "black", width = 0.6) +
  facet_nested(Treatment_access + Resistance_level ~ Seasonality+ Dosage, labeller = labeller(Resistance_level = R.labs, Treatment_access = T.labs, Dosage = D.labs, Seasonality = S.labs)) +
  scale_fill_manual(
    name = "Factors:",
    values = c(
      "#999933",
      "#AA4499",
      "#009E73",
      "#661100",
      "#888888"),
    breaks = c("Fitness", "half-life", "Cmax/IC50", "MKR", "Diangostic"),
    labels = c("\nFitness cost\n", "\nHalf-life (days) of the\nlong-acting drug", "\nCmax/EC50 of the\nlong-acting drug", "\nEmax (per day) of the\nlong-acting drug", "\nDiagnostic detection\nlimit (parasites/ul)")) +
  theme(axis.title = element_text(face = "bold")) +
  ylab("First-order indices") +
  xlab("EIR (inoculations per person per year)") +
  theme_bw() +
  theme(
    axis.text.x = element_text(size = 15/constant),
    axis.text.y = element_text(size = 15/constant),
    axis.title.x = element_text(size = 16/constant, face = "bold"),
    axis.title.y = element_text(size = 16/constant, face = "bold"),
    plot.title = element_text(size = 20/constant, hjust = 0.5, face = "bold")) +
  theme(legend.text = element_text(size = 16/constant)) +
  theme(legend.title = element_text(size = 16/constant, face = "bold")) +
  ggtitle(label = "Long-acting drug") +
  theme(
    strip.text.x = element_text(size = 16/constant, color = "black", face = "bold"),
    strip.text.y = element_text(size = 16/constant, color = "black", face = "bold")) +
  scale_y_continuous(breaks = break_y, labels = Label_yy)+
  theme(legend.position="top", legend.box="vertical")  +
  guides(fill=guide_legend(nrow=2))

# load the data
data <- read.csv(file = "data/Figure3-figure-supplement5-Sourcedata1.csv", header  = TRUE)

# transform the constrained variable into factor
data$Resistance_level <- factor(data$Resistance_level, levels = c("7", "18"))
data$Dosage <- factor(data$Dosage, levels = c("1", "0"))
data$EIR <- factor(data$EIR, levels = c("5", "10", "500"))
data$Treatment_access <- factor(data$Treatment_access, levels = c("0.5","0.04"))

# creat a lavel for each constrained variable
T.labs <- c("High access to treatment", "Low access to treatment")
names(T.labs) <- c("0.5", "0.04")
R.labs <- c("High degree\n of resistance", "Low degree\n of resistance")
names(R.labs) <- c("18", "7") 
D.labs <- c("High adherence\n to treatment", "Low adherence\n to treatment")
names(D.labs) <- c("1", "0")
S.labs <- c("No seasonality", "Seasonality")
names(S.labs) <- c("sesonality1", "sesonality2")

# define the break for the y axis
break_y <- c(0, 0.25, 0.5, 0.75, 1)
Label_yy <- c(0, 0.25, 0.5, 0.75, 1)

# ---- visualise the first order indices of each factors -----

# select the data
data_2 <- data[data$Effect == "First", ]

# order the level of each factors
data_2$Factor <- factor(data_2$Factor, levels = c(
  "Diangostic",
  "MKR_short",
  "MKR_long",
  "IC50_S_short",
  "half-life_short",
  "half-life_long",
  "Cmax/IC50",
  "Fitness"))

constant<-2
# visualise
ggplot(data_2, aes(x = EIR, y = First, fill = Factor)) +
  geom_col(color = "black", width = 0.6) +
  facet_nested(Treatment_access + Resistance_level ~ Seasonality + Dosage, labeller = labeller(Resistance_level = R.labs, Treatment_access = T.labs, Dosage = D.labs, Seasonality = S.labs)) +
  scale_fill_manual(
    values = c(
      "#999933",
      "#E6959F", "#AA4499",
      "#88CCEE", "#009E73",
      "#882255", "#661100",
      "#888888"),
    name = "Factors:",
    breaks = c("Fitness", "half-life_short", "half-life_long", "IC50_S_short", "Cmax/IC50", "MKR_short", "MKR_long", "Diangostic"),
    labels = c("\nFitness cost\n", "\nHalf-life (days) of the\nshort-acting drug", "\nHalf-life  of the\nlong-acting drug", "\nCmax/EC50 of the\nshort-acting drug\n", "\nCmax/EC50  of the\nlong-acting drug\n", "\nEmax (per day) of the\nshort-acting drug", "\nEmax (per day) of the\nshort-acting drug", "\nDiagnostic detection\nlimit (parasites/ul)")) +
  theme(axis.title = element_text(face = "bold")) +
  ylab("First-order indices") +
  xlab("EIR (inoculations per person per year)") +
  theme_bw() +
  theme(
    axis.text.x = element_text(size = 15/constant),
    axis.text.y = element_text(size = 15/constant),
    axis.title.x = element_text(size = 16/constant, face = "bold"),
    axis.title.y = element_text(size = 16/constant, face = "bold"),
    plot.title = element_text(size = 20/constant, hjust = 0.5, face = "bold")) +
  theme(legend.text = element_text(size = 16/constant)) +
  theme(legend.title = element_text(size = 16/constant, face = "bold")) +
  ggtitle(label = "Short-acting + Long-acting drugs") +
  theme(
    strip.text.x = element_text(size = 16/constant, color = "black", face = "bold"),
    strip.text.y = element_text(size = 16/constant, color = "black", face = "bold")) +
  scale_y_continuous(breaks = break_y, labels = Label_yy) +
  theme(legend.position="top", legend.box="vertical")  +
  guides(fill=guide_legend(nrow=2))

::: {warning=FALSE}

::: {#fig3s5}

# Load the data
Scenario_liste <- read.csv(file = "data/Figure3-figuresupplement679-Sourcedata1.csv", header  = TRUE)

# Change arm information into factor
Scenario_liste$drug <- factor(Scenario_liste$drug,
                              levels = c("Drug A", "Drug B", "Drug A + Drug B"))

Scenario_liste$eir <- factor(Scenario_liste$eir,
                             levels = c("5", "10", "500"))

Scenario_liste$Access <- factor(Scenario_liste$Access,
                                levels = c("0.04", "0.5"))

Scenario_liste$Resistance_Level <- factor(Scenario_liste$Resistance_Level,
                                          levels = c("2.5", "7", "10", "18"))

Scenario_liste$Dosage <- factor(Scenario_liste$Dosage,
                                levels = c("1", "0"))

Scenario_liste$seasonality <- factor(Scenario_liste$seasonality,
                                     levels = c("sesonality1", "sesonality2"))


#---- Plot at high level of access to treatment (Figure_3_supplement_5) ----

# Select arm that have a high level of treatment access
Scenario_liste_2 <- Scenario_liste[Scenario_liste$Access == 0.5, ]
Scenario_liste_2$Resistance_Level[Scenario_liste_2$Resistance_Level == 2.5] <- 7
Scenario_liste_2$Resistance_Level[Scenario_liste_2$Resistance_Level == 10] <- 18
Scenario_liste_2 <- Scenario_liste_2[Scenario_liste_2$Resistance_Level == 7, ]

# Define the label
R.labs <-  c(
  "Low resistance\n level",
  "Low resistance\n level",
  "High resistance\n level",
  "High resistance\n level")
names(R.labs) <- c("7", "2.5", "18", "10")

Dr.labs <- c("Short-acting drug", "Long-acting drug", "Short-acting +\nLong-acting drugs")
names(Dr.labs) <- c("Drug A", "Drug B", "Drug A + Drug B")


constant<-2
# Plot by level of adherence to treatment
pd <- position_dodge(0.9)
PP1 <- ggplot(Scenario_liste_2, aes(x = Dosage, y = Indicator, fill = Dosage)) +
  geom_boxplot(position = pd) +
  labs(title = "", x = "Adherence to treatment (%)", y = "Selection coefficient") +
  scale_x_discrete(labels = c("1" = "100%", "0" = "60%")) +
  facet_grid( ~ drug,labeller = labeller(drug = Dr.labs)) +
  scale_fill_manual(
    name = "Adherence to\n  treatment",
    values = c("lightblue", "blue"),
    breaks = c(1, 0),
    labels = c("100 %", "60 %")) +
  scale_color_manual(
    name = "Resistance level",
    values = c("black", "grey"),
    breaks = c(2.5, 10),
    labels = c("Low (2.5)", "High (10)")) +
  theme_bw() +
  theme(
    axis.text.x = element_text(size = 15/constant),
    axis.text.y = element_text(size = 15/constant),
    axis.title.y = element_text(size = 16/constant, face = "bold"),
    axis.title.x = element_text(size = 16/constant, face = "bold"),
    plot.title = element_text(size = 18/constant, hjust = 0.5, face = "bold"),
    legend.text = element_text(size = 16/constant),
    legend.title = element_text(size = 16/constant)) +
  theme(
    strip.text.x = element_text(size = 16/constant, color = "black", face = "bold"),
    strip.text.y = element_text(size = 16/constant, color = "black", face = "bold")) +
  theme(legend.position = "none")

# Plot by level of seasonality
pd <- position_dodge(0.9)
PP2 <-ggplot(Scenario_liste_2,
             aes(x = seasonality, y = Indicator,  fill = seasonality)) +
  geom_boxplot(position = pd) +
  facet_grid( ~ drug,labeller = labeller(drug = Dr.labs)) +
  labs(title = "", x = "Seasonality", y = "Selection coefficient") +
  scale_x_discrete(labels = c("sesonality1" = "No", "sesonality2" = "Yes")) +
  scale_fill_manual(
    name = "Seasonality",
    values = c("lightblue", "blue"),
    breaks = c("sesonality1", "sesonality2"),
    labels = c("No", "Yes")) +
  scale_color_manual(
    name = "Resistance level",
    values = c("black", "grey"),
    breaks = c(2.5, 10),
    labels = c("Low (2.5)", "High (10)")) +
  theme_bw() +
  theme(axis.text.x = element_text(size = 15/constant),
        axis.text.y = element_text(size = 15/constant),
        axis.title.y = element_text(size = 16/constant, face = "bold"),
        axis.title.x = element_text(size = 16/constant, face = "bold"),
        plot.title = element_text(size = 18/constant, hjust = 0.5, face = "bold"),
        legend.text = element_text(size = 16/constant),
        legend.title = element_text(size = 16/constant)) +
  theme(
    strip.text.x = element_text(size = 16/constant, color = "black", face = "bold"),
    strip.text.y = element_text(size = 16/constant, color = "black", face = "bold")) +
  theme(legend.position = "none")


# Plot by level of transmission
pd <- position_dodge(0.9)
PP3 <- ggplot(Scenario_liste_2, aes(x = eir, y = Indicator, fill = eir)) +
  geom_boxplot(position = pd) +
  labs(title = "", x = "EIR (inoculations per person per year)", y = "Selection coefficient") +
  facet_grid( ~ drug,labeller = labeller(drug = Dr.labs)) +
  scale_x_discrete(labels = c("5" = "5", "10" = "10", "500" = "500")) +
  scale_fill_manual(name = "EIR",
                    values = c("white", "lightblue", "blue")) +
  scale_color_manual(
    name = "Resistance level",
    values = c("black", "grey"),
    breaks = c(2.5, 10),
    labels = c("Low (2.5)", "High (10)")) +
  theme_bw() +
  theme(
    axis.text.x = element_text(size = 15/constant),
    axis.text.y = element_text(size = 15/constant),
    axis.title.y = element_text(size = 16/constant, face = "bold"),
    axis.title.x = element_text(size = 16/constant, face = "bold"),
    plot.title = element_text(size = 18/constant, hjust = 0.5, face = "bold"),
    legend.text = element_text(size = 16/constant),
    legend.title = element_text(size = 16/constant)) +
  theme(
    strip.text.x = element_text(size = 16/constant, color = "black", face = "bold"),
    strip.text.y = element_text(size = 16/constant, color = "black", face = "bold")) +
  theme(legend.position = "none")


# Merge all plot
grid.arrange(PP3,PP2,PP1,ncol=1)

# Select arm that have a high level of treatment access
Scenario_liste_2 <- Scenario_liste[Scenario_liste$Access == 0.5, ]
Scenario_liste_2$Resistance_Level[Scenario_liste_2$Resistance_Level == 2.5] <- 7
Scenario_liste_2$Resistance_Level[Scenario_liste_2$Resistance_Level == 10] <- 18
Scenario_liste_2 <- Scenario_liste_2[Scenario_liste_2$Resistance_Level == 18, ]

# Define the label
R.labs <-  c(
  "Low resistance\n level",
  "Low resistance\n level",
  "High resistance\n level",
  "High resistance\n level")
names(R.labs) <- c("7", "2.5", "18", "10")

Dr.labs <- c("Short-acting drug", "Long-acting drug", "Short-acting +\nLong-acting drugs")
names(Dr.labs) <- c("Drug A", "Drug B", "Drug A + Drug B")


constant<-2
# Plot by level of adherence to treatment
pd <- position_dodge(0.9)
PP1 <- ggplot(Scenario_liste_2, aes(x = Dosage, y = Indicator, fill = Dosage)) +
  geom_boxplot(position = pd) +
  labs(title = "", x = "Adherence to treatment (%)", y = "Selection coefficient") +
  scale_x_discrete(labels = c("1" = "100%", "0" = "60%")) +
  facet_grid( ~ drug,labeller = labeller(drug = Dr.labs)) +
  scale_fill_manual(
    name = "Adherence to\n  treatment",
    values = c("lightblue", "blue"),
    breaks = c(1, 0),
    labels = c("100 %", "60 %")) +
  scale_color_manual(
    name = "Resistance level",
    values = c("black", "grey"),
    breaks = c(2.5, 10),
    labels = c("Low (2.5)", "High (10)")) +
  theme_bw() +
  theme(
    axis.text.x = element_text(size = 15/constant),
    axis.text.y = element_text(size = 15/constant),
    axis.title.y = element_text(size = 16/constant, face = "bold"),
    axis.title.x = element_text(size = 16/constant, face = "bold"),
    plot.title = element_text(size = 18/constant, hjust = 0.5, face = "bold"),
    legend.text = element_text(size = 16/constant),
    legend.title = element_text(size = 16/constant)) +
  theme(
    strip.text.x = element_text(size = 16/constant, color = "black", face = "bold"),
    strip.text.y = element_text(size = 16/constant, color = "black", face = "bold")) +
  theme(legend.position = "none")

# Plot by level of seasonality
pd <- position_dodge(0.9)
PP2 <-ggplot(Scenario_liste_2,
             aes(x = seasonality, y = Indicator,  fill = seasonality)) +
  geom_boxplot(position = pd) +
  facet_grid( ~ drug,labeller = labeller(drug = Dr.labs)) +
  labs(title = "", x = "Seasonality", y = "Selection coefficient") +
  scale_x_discrete(labels = c("sesonality1" = "No", "sesonality2" = "Yes")) +
  scale_fill_manual(
    name = "Seasonality",
    values = c("lightblue", "blue"),
    breaks = c("sesonality1", "sesonality2"),
    labels = c("No", "Yes")) +
  scale_color_manual(
    name = "Resistance level",
    values = c("black", "grey"),
    breaks = c(2.5, 10),
    labels = c("Low (2.5)", "High (10)")) +
  theme_bw() +
  theme(axis.text.x = element_text(size = 15/constant),
        axis.text.y = element_text(size = 15/constant),
        axis.title.y = element_text(size = 16/constant, face = "bold"),
        axis.title.x = element_text(size = 16/constant, face = "bold"),
        plot.title = element_text(size = 18/constant, hjust = 0.5, face = "bold"),
        legend.text = element_text(size = 16/constant),
        legend.title = element_text(size = 16/constant)) +
  theme(
    strip.text.x = element_text(size = 16/constant, color = "black", face = "bold"),
    strip.text.y = element_text(size = 16/constant, color = "black", face = "bold")) +
  theme(legend.position = "none")


# Plot by level of transmission
pd <- position_dodge(0.9)
PP3 <- ggplot(Scenario_liste_2, aes(x = eir, y = Indicator, fill = eir)) +
  geom_boxplot(position = pd) +
  labs(title = "", x = "EIR (inoculations per person per year)", y = "Selection coefficient") +
  facet_grid( ~ drug,labeller = labeller(drug = Dr.labs)) +
  scale_x_discrete(labels = c("5" = "5", "10" = "10", "500" = "500")) +
  scale_fill_manual(name = "EIR",
                    values = c("white", "lightblue", "blue")) +
  scale_color_manual(
    name = "Resistance level",
    values = c("black", "grey"),
    breaks = c(2.5, 10),
    labels = c("Low (2.5)", "High (10)")) +
  theme_bw() +
  theme(
    axis.text.x = element_text(size = 15/constant),
    axis.text.y = element_text(size = 15/constant),
    axis.title.y = element_text(size = 16/constant, face = "bold"),
    axis.title.x = element_text(size = 16/constant, face = "bold"),
    plot.title = element_text(size = 18/constant, hjust = 0.5, face = "bold"),
    legend.text = element_text(size = 16/constant),
    legend.title = element_text(size = 16/constant)) +
  theme(
    strip.text.x = element_text(size = 16/constant, color = "black", face = "bold"),
    strip.text.y = element_text(size = 16/constant, color = "black", face = "bold")) +
  theme(legend.position = "none")


# Merge all plot
grid.arrange(PP3,PP2,PP1,ncol=1)

# Load the data
Scenario_liste<-read.csv(file = "data/Figure3-figuresupplement8-Sourcedata1.csv", header=T)
Scenario_liste_2<-Scenario_liste[Scenario_liste$Access==0.5 & Scenario_liste$seed==2 & Scenario_liste$eir<=400,]

Scenario_liste_3<-NULL
Scenario_liste_3$eir<-c(Scenario_liste_2$eir,Scenario_liste_2$eir)
Scenario_liste_3$Percentage<-c(Scenario_liste_2$Treated_Total,Scenario_liste_2$Treated_infectTotal)
Scenario_liste_3$Total<-c(rep("Person", length(Scenario_liste_2$eir)),rep("infection", length(Scenario_liste_2$eir)))
Scenario_liste_3<-as.data.frame(Scenario_liste_3)

#Plot
constant<-2
ggplot(Scenario_liste_3,aes(x=eir, y=Percentage*100, col=Total))+
  geom_point(size=3/constant)+
  ylim(0,30)+
  ylab("Percentage (%)")+
  xlab("EIR (inoculation per person per year)")+
  scale_color_manual(name  ="Percentage of:",
                     values=c(1,2),
                     #values=c(SteppedSequential5Steps[9],SteppedSequential5Steps[17]),
                     breaks=c("person", "infection"),
                     labels=c("People that received treatment", "Infected people that received treatment")
  )+
  theme_bw() +
  theme(
    axis.text.x = element_text(size = 16/constant),
    axis.text.y = element_text(size = 16/constant),
    axis.title.x = element_text(size = 18/constant, face = "bold"),
    axis.title.y = element_text(size = 18/constant, face = "bold")) +
  theme(legend.text = element_text(size = 18/constant)) +
  theme(legend.title = element_text(size = 18/constant, face = "bold")) +
  theme(legend.key.size = unit(0.9, "cm")) +
  theme(legend.key.width = unit(1.0, "cm"))+ 
  theme(legend.position = c(0.5, 0.8))

# Scenario_liste_2 <- Scenario_liste[Scenario_liste$Access == 0.04, ]
# Scenario_liste_2$Resistance_Level[Scenario_liste_2$Resistance_Level == 2.5] <- 7
# Scenario_liste_2$Resistance_Level[Scenario_liste_2$Resistance_Level == 10] <- 18
# Scenario_liste_2 <- Scenario_liste_2[Scenario_liste_2$Resistance_Level == 7, ]
# 
# # define the label
# R.labs <- c(
#   "Low resistance\n level",
#   "Low resistance\n level",
#   "High resistance\n level",
#   "High resistance\n level")
# 
# names(R.labs) <- c("7", "2.5", "18", "10")
# 
# constant<-2
# # plot by level of treatment adherence
# pd <- position_dodge(0.9)
# PP1 <- ggplot(Scenario_liste_2, aes(x = Dosage, y = Indicator, fill = Dosage)) +
#   geom_boxplot(position = pd) +
#   labs(title = "", x = "Adherence to treatment (%)", y = "Selection coefficient") +
#   scale_x_discrete(labels = c("1" = "100%", "0" = "60%")) +
#   facet_grid( ~ drug,labeller = labeller(drug = Dr.labs)) +
#   scale_fill_manual(
#     name = "Adherence to\n  treatment",
#     values = c("lightblue", "blue"),
#     breaks = c(1, 0),
#     labels = c("100 %", "60 %")) +
#   scale_color_manual(
#     name = "Resistance level",
#     values = c("black", "grey"),
#     breaks = c(2.5, 10),
#     labels = c("Low (2.5)", "High (10)")) +
#   theme_bw() +
#   theme(
#     axis.text.x = element_text(size = 15/constant),
#     axis.text.y = element_text(size = 15/constant),
#     axis.title.y = element_text(size = 16/constant, face = "bold"),
#     axis.title.x = element_text(size = 16/constant, face = "bold"),
#     plot.title = element_text(size = 18/constant, hjust = 0.5, face = "bold"),
#     legend.text = element_text(size = 16/constant),
#     legend.title = element_text(size = 16/constant)) +
#   theme(
#     strip.text.x = element_text(size = 16/constant, color = "black", face = "bold"),
#     strip.text.y = element_text(size = 16/constant, color = "black", face = "bold")) +
#   theme(legend.position = "none")
# 
# # plot by level of seasonality
# pd <- position_dodge(0.9)
# PP2 <- ggplot(Scenario_liste_2,
#               aes(x = seasonality, y = Indicator,  fill = seasonality)) +
#   geom_boxplot(position = pd) +
#   facet_grid( ~ drug,labeller = labeller(drug = Dr.labs)) +
#   labs(title = "", x = "Seasonality", y = "Selection coefficient") +
#   scale_x_discrete(labels = c("sesonality1" = "No", "sesonality2" = "Yes")) +
#   scale_fill_manual(
#     name = "Seasonality",
#     values = c("lightblue", "blue"),
#     breaks = c("sesonality1", "sesonality2"),
#     labels = c("No", "Yes")) +
#   scale_color_manual(
#     name = "Resistance level",
#     values = c("black", "grey"),
#     breaks = c(2.5, 10),
#     labels = c("Low (2.5)", "High (10)")) +
#   theme_bw() +
#   theme(
#     axis.text.x = element_text(size = 15/constant),
#     axis.text.y = element_text(size = 15/constant),
#     axis.title.y = element_text(size = 16/constant, face = "bold"),
#     axis.title.x = element_text(size = 16/constant, face = "bold"),
#     plot.title = element_text(size = 18/constant, hjust = 0.5, face = "bold"),
#     legend.text = element_text(size = 16/constant),
#     legend.title = element_text(size = 16/constant)) +
#   theme(
#     strip.text.x = element_text(size = 16/constant, color = "black", face = "bold"),
#     strip.text.y = element_text(size = 16/constant, color = "black", face = "bold")) +
#   theme(legend.position = "none")
# 
# # plot by level of transmission
# pd <- position_dodge(0.9)
# PP3 <- ggplot(Scenario_liste_2, aes(x = eir, y = Indicator, fill = eir)) +
#   geom_boxplot(position = pd) +
#   labs(title = "", x = "EIR (inoculations per person per year)", y = "Selection coefficient") +
#   facet_grid( ~ drug,labeller = labeller(drug = Dr.labs)) +
#   scale_x_discrete(labels = c("5" = "5", "10" = "10", "500" = "500")) +
#   scale_fill_manual(name = "EIR",
#                     values = c("white", "lightblue", "blue")) +
#   scale_color_manual(
#     name = "Resistance level",
#     values = c("black", "grey"),
#     breaks = c(2.5, 10),
#     labels = c("Low (2.5)", "High (10)")) +
#   theme_bw() +
#   theme(
#     axis.text.x = element_text(size = 15/constant),
#     axis.text.y = element_text(size = 15/constant),
#     axis.title.y = element_text(size = 16/constant, face = "bold"),
#     axis.title.x = element_text(size = 16/constant, face = "bold"),
#     plot.title = element_text(size = 18/constant, hjust = 0.5, face = "bold"),
#     legend.text = element_text(size = 16/constant),
#     legend.title = element_text(size = 16/constant)) +
#   theme(
#     strip.text.x = element_text(size = 16/constant, color = "black", face = "bold"),
#     strip.text.y = element_text(size = 16/constant, color = "black", face = "bold")) +
#   theme(legend.position = "none")

# Merge all plot
# plot_grid(
#   PP3,
#   PP2,
#   PP1,
#   ncol = 1,
#   nrow = 3,
#   rel_widths = c(1, 1),
#   scale = 1,
#   labels = c('A', 'B', "C"),
#   label_size = 18/constant)

Across settings with a low access to treatment, we found that fitness cost had the largest influence on the selection coefficient (Figure 3—figure supplements 3–!number(5)). The fitness cost of a resistant genotype was defined as the relative decrease in the resistant genotype multiplication rate within an untreated human host compared with the sensitive genotype. Consequently, high fitness costs prevented the spread of resistance (Figure 3, Figure 3—figure supplement 2). At a high level of access to treatment, the effect of fitness cost was reduced, and drug properties played a critical role in the spread of drug resistance, and their influence varied for each treatment profile as described below.

For the short-acting drug used as monotherapy, the half-life had an important influence on the rate of spread (Figure 3—figure supplement 3). A long half-life reduced the spread of resistant parasites by extending the period during which the drug killed partially resistant parasites (Figure 3, Figure 3—figure supplement 1). The half-life effect was stronger for parasites with a lower degree of resistance (Figure 3, Figure 3—figure supplement 3), as highly resistant parasites were less affected by the drug. Furthermore, the spread of the resistant genotype was faster in populations with low adherence to treatment (Figure 3—figure supplement 6) because with fewer treatment doses, the parasite was exposed to the drug for a shorter time, leading to higher parasite survival. Overall, these results highlight that the time during which the parasite is exposed to artemisinin is a critical driver of the spread of partial artemisinin resistance.

For parasites with a low degree of resistance to the long-acting drug used as monotherapy, the drug half-life also had a key influence on the selection coefficient (Figure 3—figure supplement 4). However, long half-lives were associated with large selection coefficients (Figure 3, Figure 3—figure supplement 1). Drugs with a long half-life have an extended period of low drug concentration in treated patients during which only resistant parasites can infect the host. This period of low drug concentration is called the selection window 34Hastings et al.200247Kay and Hastings2015. These results confirm that the selection window plays a crucial role in the spread of resistance to long-acting drugs.

The spread of parasites with a high degree of resistance to the long-acting drug used as monotherapy was also accelerated by longer drug half-life (Figure 3—figure supplement 4). For these resistant parasites, the ratio Cmax/EC50 also had an important influence on the rate of spread (Figure 3—figure supplement 4). This ratio captured the duration of the drug killing effect on the sensitive genotype by assessing the proximity between the EC50 of the sensitive genotype and Cmax. A lower Cmax/EC50 ratio captures a shorter duration of the drug killing effect for the sensitive genotype and, consequently, also a lower duration of the drug killing effect against the resistant genotype (higher EC50). Thus, when the drug had a low Cmax/EC50 ratio, the duration of the drug killing effect was not long enough to ensure the successful clearances of parasites with a higher degree of resistance (higher EC50), favouring their spread (Figure 3, Figure 3—figure supplement 1). Furthermore, for parasites with a low degree of resistance (lower EC50), the ratio Cmax/EC50 also influenced the rate of spread in settings with a low level of treatment adherence, since low adherence reduces a Cmax leading to treatment failure (Figure 3—figure supplements 1 and !number(4)). These results highlight the importance of treatment adherence to assure that the drug concentration is high enough to eliminate partially resistant genotypes and limit their spread.

When the genotype was resistant to the short-acting drug in a population that used the short-acting and the long-acting drugs in combination, factors related to the long-acting drug had the most influence on the selection coefficient (Figure 3, Figure 3—figure supplement 5). When the Cmax/EC50 ratio of the long-acting drug was large, the duration of the killing effect of the long-acting drug on parasites resistant to the short-acting drug was higher, reducing their spread (Figure 3, Figure 3—figure supplement 1). In addition, the rate of spread rose when the level of adherence to treatment was low (Figure 3—figure supplement 6). These results highlight that the spread of partial resistance to artemisinin strongly depends on the capacity of the partner drug to kill them.

The influence of the transmission intensity (represented by entomological inoculation rate [EIR]) and its seasonality on the selection coefficient varied by treatment profiles and degrees of resistance. When the parasite was resistant to the short-acting drug (when used in monotherapy or combination), selection coefficients were higher in settings with lower EIR (Figure 3—figure supplements 6 and !number(7)). We observed a similar trend for parasites with a high degree of resistance to the long-acting drug used in monotherapy (Figure 3—figure supplement 7). Two factors account for this trend. First, the selection of parasites resistant to the short-acting drug (low and high degrees of resistance) and parasites highly resistant to the long-acting drug depends on the proportion of treated infections. A higher portion of treated infections can lead to a higher proportion of delayed parasite clearance or treatment failure of drug-resistant infections allowing these resistant genotypes to spread. The proportion of treated infections is higher at lower EIR due to the lower level of immunity (Figure 3—figure supplement 8). Furthermore, lower immunity levels for individuals living in low transmission settings may also increase the risk of treatment failure and favour the spread of resistance. Second, there is a higher proportion of individuals coinfected by both genotypes at higher EIR. This higher level of co-infection enhanced within-host competition between genotypes, which inhibits the multiplication of resistant parasites within hosts due to their fitness cost and thus limits their spread. Similarly, the spread of resistant parasites was higher in the seasonal settings than in non-seasonal settings (Figure 3—figure supplements 6 and !number(7)) due to the reduction of immunity levels and a decline in within-host competition between genotypes during the low transmission season of the seasonal settings. Overall, these results indicate that the spread of partial artemisinin resistance is faster in seasonal settings with low transmission levels.

However, for parasites with a low degree of resistance to the long-acting drug used in monotherapy, selection coefficients were higher in settings with a large EIR (Figure 3—figure supplement 6). This arises because the proportion of patients with low drug concentrations persisting from previous treatments increases at higher EIR where higher infection rates increase the overall usage of treatment (Figure 3—figure supplement 8). These low drug concentrations may fall within the selective window and hence drive the spread of parasites partially resistant to the long-acting drug. Note that this trend was only observed for settings with high access to treatment. In settings with low access to treatment, we observe similar trends as for parasites resistant to the short-acting drug (Figure 3—figure supplement 9) since here, the impact of the selection window was more negligible. These results highlight that the selection window of the long-acting drug can change the interplay between the transmission setting and the spread of drug resistance.

Probability of establishment of drug resistance and its key drivers

Population genetic theory has shown that the probability of establishment of a mutation depends on two factors: (i) the size of its selection coefficient (i.e. establishment becomes more likely as the mutation becomes more advantageous) and (ii) the degree of heterogeneity in the number of parasite offspring. This occurs because higher heterogeneity increases stochastic fluctuations of allele number, so increases the chance that the mutation is lost despite its advantage 113Wiesch et al.201135Hastings200438Hastings et al.202033Hastings and Mackinnon1998. The probability of establishment can also be altered by temporal fluctuation in the population size or magnitude of the selection coefficient 98Waxman2011. Both effects are likely to be present in seasonal settings of malaria transmission where population size fluctuates, and selection intensity may also change if the level of drug use fluctuates in response to the seasonality of transmission. We avoid these complications by investigating only non-seasonal settings, from which we selected 10 different resistant genotypes having a known selection coefficient and quantified their probability of establishment (see Materials and methods). By doing so, we evaluated the relationship between the selection coefficient and probability of establishment and assessed how this relationship varies across settings due to variation in the heterogeneity of parasite reproductive success.

As expected, the establishment of a mutation was more probable when its selection coefficient was high (Figure 4). For each treatment profile, the probability of establishment of mutations with similar selection coefficients was higher at low EIR than at high EIR (Figure 4), especially for mutations with a high selection coefficient. The lower probability of establishment in higher transmission settings suggests that higher transmission levels increase the heterogeneity in parasite reproductive success, reducing the chance to transmit advantageous mutations. Two factors increase the heterogeneity of parasite reproductive success in settings with a high EIR. First, in higher transmission settings, there is higher variability in the number of parasites with distinct genotypes carried by a host and which are competing for reproductive success. Thus, the greater variability of within-host competition between hosts leads to greater heterogeneity of parasite reproductive success. Second, settings with a high EIR have a large variation in the level of individual immunity. Host immunity influences parasite reproductive success by reducing parasite growth within the human host. Therefore, in high transmission settings, the greater variation of immunity leads to higher heterogeneity of parasite reproductive success and reduces the chance that the emerging mutation will be transmitted despite its advantage.

# Load the data
Scenario_liste <- read.csv(file = "data/Figure4-Sourcedata1.csv", header  = TRUE)

# Define setting variable to be a factor
Scenario_liste$eir <- factor(Scenario_liste$eir,
                             levels = c("5", "10", "500"))

Scenario_liste$Dosage <- factor(Scenario_liste$Dosage,
                                levels = c("1", "0"))


Scenario_liste$drug <- factor(Scenario_liste$drug,
                              levels = c("A", "B", "A+B"))

# Select the data from the setting that we want show on the plot
Scenario_liste_2 <- Scenario_liste[Scenario_liste$eir == 5 | Scenario_liste$eir == 500,]
Scenario_liste_2 <- Scenario_liste_2[Scenario_liste_2$Dosage == 1,]

# Define label for variable that are a factor
D.labs <- c("Adherence = 100 %", "Adherence = 60 %")
names(D.labs) <- c("1", "0")

E.labs <- c("EIR = 5", "EIR = 10", "EIR = 500")
names(E.labs) <- c("5", "10", "500")

Dr.labs <- c("Short-acting drug", "Long-acting drug", "Short-acting +\nLong-acting drugs")
names(Dr.labs) <- c("A", "B", "A+B")

# Define a constant use to adjust size of plot 
constant<-2.5

# Estimate the 95CI
Scenario_liste_2$Pe_L<-Scenario_liste_2$Pe-1.96*sqrt((Scenario_liste_2$Pe*(1-Scenario_liste_2$Pe))/Scenario_liste_2$Number_mutation)
Scenario_liste_2$Pe_U<-Scenario_liste_2$Pe+1.96*sqrt((Scenario_liste_2$Pe*(1-Scenario_liste_2$Pe))/Scenario_liste_2$Number_mutation)


# Plot the data
ggplot(data = Scenario_liste_2) +
  geom_line(aes(x = Indicator, y = Pe, linetype = eir,color = eir), size = 2 / constant) +
  geom_ribbon(aes(x = Indicator, ymin = Pe_L, ymax = Pe_U, fill = eir), alpha = 0.2) +
  facet_grid( ~ drug,
              scale = "free",
              labeller = labeller(eir = E.labs, drug = Dr.labs)) +
  theme_bw() +
  scale_linetype_manual(
    values = c("solid", "solid", "twodash"),
    name = "EIR (inoculations per person per year):",
    breaks = c("5", "10", "500"),
    labels = c("5", "10", "500")  ) +
  scale_color_manual(
    values = c("#CC8E51","#0072B2", "#2C85B2"), # "#009E73","#0072B2", "#D55E00","#56B4E9","#CC79A7"
    name = "EIR (inoculations per person per year):",
    breaks = c("5", "10", "500"),
    labels = c("5", "10", "500")
  ) +
  scale_fill_manual(
    values = c("#CC8E51","#0072B2", "#2C85B2"), # "#009E73","#0072B2", "#D55E00","#56B4E9","#CC79A7"
    name = "EIR (inoculations per person per year):",
    breaks = c("5", "10", "500"),
    labels = c("5", "10", "500")
  ) +
  scale_x_continuous(name = "Selection coefficient") +
  scale_y_continuous(name = "Probability of establishment") +
  theme(
    axis.text.x = element_text(size = 16 / constant),
    axis.text.y = element_text(size = 16 / constant),
    axis.title.x = element_text(size = 18 / constant, face = "bold"),
    axis.title.y = element_text(size = 18 / constant, face = "bold"),
    plot.title = element_text(
      size = 20 / constant,
      hjust = 0.5,
      face = "bold")) +
  theme(legend.text = element_text(size = 18 / constant)) +
  theme(legend.title = element_text(size = 18 / constant, face = "bold")) +
  ggtitle(label = "") +
  theme(
    strip.text.x = element_text(
      size = 18 / constant,
      color = "black",
      face = "bold"),
    strip.text.y = element_text(
      size = 18 / constant,
      color = "black",
      face = "bold")) +
  theme(plot.margin = unit(c(0.1, 0.1, 0.1, 0.1), "cm")) +
  guides(linetype = guide_legend(keywidth = 4 / constant, keyheight = 1 /constant)) +
  theme(
    legend.position = "top",
    legend.direction = "horizontal",
    legend.spacing.y = unit(-5, "cm"),
    legend.margin = margin(-0.3, -0.3, -0.3, -0.3, unit = "cm"))

Discussion

Understanding which disease, transmission, epidemiological, health system, and drug factors systemically drive the evolution of drug resistance is challenging. A full understanding requires vast observational data or clinical trials on a scale that is not possible or mathematical models that are sufficiently detailed to capture all these factors while remaining computationally feasible to simultaneously assess the impact of these factors. In response to this need, we updated a detailed individual-based model of malaria dynamics to include a full pharmacological (i.e. PK/PD) description of antimalarial treatments. We introduced a global sensitivity analysis approach based on emulators for computationally intensive models to systematically assess which factors jointly drive the evolution of drug-resistant parasites. As discussed below, our approach allowed us to understand the guiding principles of the evolution of drug resistance against ACTs and to explain the difference in trends observed in the GMS and in malaria endemic Africa. Improving our understanding of the factors that lead to drug resistance establishment and spread allows us to identify strategies to mitigate these dynamics and guides initial considerations for developing more sustainable malaria treatment.

Our results support the belief that evolution of resistance to ACTs begins with the establishment and spread of parasites resistant to the partner drug and once the protective effect of the partner drug is reduced, drug selection falls on the artemisinin component, and parasites then start to acquire resistance to artemisinin derivatives (e.g. 97Watson et al.2021; 37Hastings et al.2016). The fact that resistance to the partner drug appears before resistance to artemisinin derivatives was supported by two points elucidated in our study. First, resistance to the partner drug depends on the period of low concentration of this drug during which only resistant parasites can multiply within the host (known as the selection window). As artemisinin derivatives are short-acting, they cannot prevent patients from being reinfected by parasites resistant to the partner drug during this selection window. Second, resistance to the partner drug was the most critical factor that enhanced establishment and spread of partial artemisinin resistance. Without resistance to the partner drug, parasites partially resistant to artemisinin could only spread at a low rate as the partner drug could still eliminate them, thereby removing their selective advantage. Our results are in line with recent molecular data which show that parasites resistant to partner drugs (piperaquine and mefloquine) were already present in the GMS before partial artemisinin resistance emerged and that the spread of resistance to artemisinin accelerated when it became linked to resistance to the partner drugs 1Amato et al.201831Hamilton et al.2019110Wongsrichanalai and Meshnick2008. Thus, the presence of partner drug resistance has probably facilitated the spread of resistance to artemisinin in the GMS. In contrast, in Africa, to date, only a low degree of resistance to the most commonly used partner drugs (lumefantrine and amodiaquine) are present 104WHO202024Ehrlich et al.2021, which has likely limited establishment of resistance to artemisinin derivatives. We additionally note that the evolution of drug resistance in the GMS may have been favoured by the low transmission intensity (annual EIR range approximate from less than 1 to 25 inoculations per person per year 22Edwards et al.201913Chaumeau et al.201823Edwards et al.2019) compare to Africa where the transmission intensity is overall higher (annual EIR range from less than 1 to more than 500 inoculations per person per year 39Hay et al.2000111Yamba et al.2020). Similar to previous studies 99White19999Bushman et al.201872Pongtavornpinyo et al.200857Lee et al.202232Hastings1997, establishment of drug resistance in our model was more likely in low transmission settings due to the reduced level of within-host competition between genotypes, as well as population immunity.

Our results suggest that a key strategy to mitigate the evolution of partial artemisinin resistance is to ensure that the partner drug efficiently kills the partially resistant parasite. Therefore, to delay the establishment of artemisinin resistance in Africa and to mitigate the spread of partial artemisinin resistance in regions where it is already established, we should ensure that limited or no genotypes are resistant to the partner drug for first-line ACT. One approach to ensure this is to implement robust molecular surveillance of resistance markers and to specify more sustainable treatment policies, such as changing first-line ACTs upon detection of resistance or when the frequency of resistant parasites reach a threshold as recommended by the WHO 104WHO2020. Furthermore, consistent with our results, adherence should continue to be promoted, as lower treatment compliance can lead to treatment failure even in the absence of resistance to the partner drug 104WHO202079Siddiqui et al.20158Bruxvoort et al.2014.

Our results suggest that future antimalarial therapies should shorten the selection windows of long-acting partner drugs. We show that resistance to long-acting drugs is the first step in the evolution of resistance to ACTs, and it depends mainly on the length of the selection window. We confirm that the selection window strongly depends on the drug half-life, also consistent with previous studies 99White199980Slater et al.201734Hastings et al.200296Watkins and Mosobo199372Pongtavornpinyo et al.200847Kay and Hastings2015. Consequently, reducing the half-life of the partner drug in an ACT regimen could reduce the spread of resistance. However, unless selection windows are substantially minimised or completely eliminated, the evolution of resistance would not totally be prevented 47Kay and Hastings2015. Thus, a more sustainable option for ACTs would be to use TACTs. TACTs involve combining an artemisinin derivative with two long-acting drugs 54Krishna2019.

If or when TACTs are to be widely used, our results emphasise that the two long-acting drugs should have matching half-lives to ensure that parasites are not exposed to residual drug concentrations of only one of the two partner drugs (noting that this is simple in principle, but more difficult in practice 36Hastings and Hodel2014). In addition, the parasite population should be devoid of parasites resistant to either of the two long-acting drugs. If resistance to one partner drug already exists in the population, the second partner drug would not be protected, and mutations conferring resistance to this second drug could be selected. However, additional forces will play a role in the evolution of resistance to drugs used in TACTs, such as if the drugs combined lead to opposite selection pressure. Additional analyses should assess which factors promote drug resistance under TACTs to guide their development. Note that the development of new partner drugs for TACTs may be challenging because combining three drugs is likely to increase the risk of toxicity and the treatment price, and future antimalarial medicines must remain tolerated by patients and affordable 54Krishna2019.

Another approach to delay the evolution of partial artemisinin resistance could focus on extending the period of action of artemisinin derivatives. In our monotherapy analysis on the spread of a genotype partially resistant to artemisinin, we found that the spread of partially resistant genotypes decreased when the drug was present in patients for a longer time, such as if it had a long half-life and there was proper treatment adherence. This result arises because partially resistant parasites are still affected by the drug 52Klonis et al.201378Sá2018107Witkowski et al.2013112Ye et al.2016. Thus, increasing their exposure to the drug leads to higher killing and reduced spread. Increasing the exposure to artemisinin derivatives can be achieved by using the artemisinin derivative having the longest half-life and, as highlighted in other studies 48Kay et al.201520Dogovski et al.201549Khoury et al.2020, can be done by increasing the number of doses and days that patients receive treatment. However, it is worth noting that extending the dosage regimen will be efficient only with adequate adherence to treatment, which may be challenging to achieve in practice. Also, as artemisinin derivatives are co-administrated with at least one long-acting drug, increasing the number of doses of this combination therapy would require reducing the concentration of the partner drug to prevent the partner drug from reaching toxic concentrations.

The evolution of drug resistance is a three-step process consisting of mutation, establishment, and spread. Mutation rates in malaria can be measured. Spread, quantified by the selection coefficient, is also easy to measure. However, the probability of establishment and its relation to the selection coefficient constituted a significant knowledge gap. Standard population genetic models assume that the number of secondary infections follows a Poisson distribution 35Hastings200418Crow and Motoo2017. Under this assumption, for selection coefficients lower than 0.2 (according to an informal literature review in 38Hastings et al.2020, most selection coefficient estimates for malaria drug resistance mutations from the field fall between 0.02 and 0.12), the probability of establishment is approximately equal to twice the selection coefficient 35Hastings200418Crow and Motoo2017. However, the number of secondary malaria infections more likely follows a negative binomial distribution due to the high heterogeneity of transmission, which may substantially reduce the probability of establishment (Box 2 of 35Hastings2004). In this modelling study, we were uniquely able to quantify the link between selection coefficients and the probability of establishment of mutations. On average, we estimated that, for selection coefficients lower than 0.2, the probability of establishment was equal to 0.87 times the selection coefficient. Therefore, our findings suggest that the variation in the number of secondary infections of P. falciparum must be much greater than the Poisson distribution assumed by standard population genetics models, and this higher variation reduces the probability of establishment of emerging mutations. Note that higher heterogeneity in parasite reproductive success may exist in the real-world than as simulated in our model due to factors not captured by our model (such as geographical heterogeneity of exposure to mosquito bites). These factors may further decrease the probability of establishment 51Klein2014.

As with all modelling studies, our approach has several limitations, primarily arising from constraints imposed by the model. First, our drug action model does not capture stage-specific killing effects, so we could not model parasites partially resistant to artemisinin being insensitive to the drug only during extended ring-stage 52Klonis et al.201378Sá2018107Witkowski et al.2013112Ye et al.2016, although previous analyses suggested this would be captured by our variation in the maximum killing rate 40Hodel et al.2016. Nevertheless, if we modelled a reduction of the drug effect restricted to the ring-stage, we expect to obtain similar results. That is, a long half-life and high treatment adherence would increase the likelihood that the drug is present within patients during any stage other than the ring-stage, and thus the drug would kill more resistant parasites.

Second, our model did not capture the impact of artemisinin resistance on gametocytes. Previous studies have highlighted that artemisinin kills gametocytes, and patients infected with parasites partially resistant to artemisinin exhibit higher gametocyte densities than patients infected with sensitive parasites 3Ashley et al.2014108Witmer et al.2020. We did not model the impact of artemisinin and resistance on gametocytes. This effect is likely to accelerate the spread of partial resistance. However, the relationship between the different factors reported in this study should be unchanged.

Third, our model, OpenMalaria, does not capture the recombination of P. falciparum parasites in mosquitoes (it does not track different genotypes in mosquitoes, and the genotype of a new infections is based on the genotype frequency in humans). In practice, this means we can only investigate the spread of resistance at one locus at a time (because if there is no genetic variability at other loci, then the lack of recombination has no impact). Our results, therefore, apply to the case when resistance is already fixed for one drug before resistance starts to spread to the second; we cannot model the simultaneous spread of resistance to two or more drugs. Moreover, the resistant genotype had a fixed degree of resistance across the simulation and could not acquire additional mutations that provide higher degrees of resistance. Nevertheless, by varying the degree of resistance in our analysis, we were able to assess the changing pattern of selection that occurred with increasing degrees of resistance. The impact of recombination when genetic variability does exist at more than one locus involved in resistance has been investigated previously by simpler genetic models whose main results are as follows. When multiple mutations are needed to confer drug resistance, recombination can slow the spread of drug resistance by separating these mutations 32Hastings199721Dye and Williams1997. Resistance to some partner drugs requires multiple mutations 104WHO2020. Partial artemisinin resistance is caused by a mutation in a single gene, but recombination may still impact its spread by separating this mutation from mutations that can minimise the fitness cost associated with resistance 86Stokes et al.2021. Recombination is more likely to impact the spread of resistant parasites in high transmission settings where recombination between different parasite genotypes is more likely. In addition, the impact of recombination depends on the frequencies of mutations involved in the resistant phenotype 21Dye and Williams1997. When their frequencies are low, recombination will have a stronger effect as resistant parasites are more likely to recombine with sensitive parasites leading to the separation of these mutations. When their frequencies are high, the impact of recombination is reduced as resistant parasites are more likely to recombine with resistant parasites. A consequence of not including recombination is that in high transmission settings, we have probably overestimated the probability of establishment of resistant parasites that have multiple mutations involved in the drug-resistant phenotype. This means that the difference between the probabilities of establishment in low and high transmission settings is likely greater than reported here. In addition, we may have overestimated the spread of these resistant parasites when these mutations are present in low frequencies.

Lastly, to investigate the establishment of drug-resistant parasites, we modelled the emergence of mutations through importation. Consequently, our estimations represent the establishment of mutations imported into a population or mutations emerging in mosquitoes (assuming that the mosquito has only transmitted the mutated genotype and not the wild type genotype to the individual). A mutation emerging during the blood-stage within the human host may have a lower probability of establishment because sensitive parasites would be present in the host, leading to competition between them. It is still unclear whether mutations conferring drug resistance arise during the blood-stage (due to the high parasite numbers) or during the sexual stage in mosquitoes (because recombination generates many genetic variations). Nevertheless, the probabilities of establishment estimated in this study are consistent with the probabilities of establishment estimated by a previous study 35Hastings2004.

In summary, our results confirm that mutations conferring malaria drug resistance are more likely to establish in low transmission settings. Our results demonstrate that the establishment and spread of resistance to artemisinin derivatives have likely been facilitated by pre-existing resistance to partner drugs. Thus, it is essential to prioritise monitoring and to limit the spread of resistance to partner drugs in current or future ACT regimens. If resistance to the partner drug is confirmed, response strategies should prioritise monitoring molecular markers and treatment failure and switching to an ACT with an effective partner drug should be considered. In addition, our results show that drug properties play an essential role in the evolution of drug resistance. Thus, the ongoing development of new antimalarial combinations should limit selection windows of partner drugs by matching half-lives, hopefully leading to longer lasting combination treatments against malaria. In the medium-term, for existing ACTs, it would be advantageous to increase the time of parasite exposure to the short-acting artemisinin derivate and/or to include a second long-acting partner drug with a matching half-life to the other long-acting partner drug (triple ACTs 54Krishna2019) and for which limited or no parasite resistance exists in the target population.

Materials and methods

Simulation model and the parameterisation of treatment profiles and resistant genotypes

Overview of our OpenMalaria model

Our individual-based model, OpenMalaria, simulates the dynamics of P. falciparum in humans and links it to a periodically forced deterministic model of P. falciparum in mosquitoes 16Chitnis et al.201281Smith et al.200683Smith et al.2008. The model structure and fitting are described in detail elsewhere 81Smith et al.200683Smith et al.2008, including open-access code (https://github.com/SwissTPH/openmalaria) and documentation (https://github.com/SwissTPH/openmalaria/wiki), and a recently published manuscript provides a new calibration 73Reiker et al.2021. Here, we have summarised the main components of OpenMalaria and its latest developments in version 40.1, which enabled us to model the establishment and spread of drug-resistant parasites.

OpenMalaria includes multiple sub-models in which mosquito and infection events, parasite, and human attributes are updated every 5 days. A demography component maintains a constant human population size and age structure across the simulation. Multiple parasite genotypes and their initial frequency can be defined in more recent model versions. For each infection, a mechanistic model simulates the parasite dynamics within the host and incorporates innate, variant, and acquired immunity 70Molineaux et al.2002. The within-host model allows for concurrent infection of multiple parasite genotypes within the same host and captures indirect competition between genotypes based on host immunity, which regulates the overall parasite load. The user can specify a reduction of the within-host multiplication factors of each genotype to model a fitness cost associated with the mutation. The host’s parasite density determines the symptoms and mortality of patients and diagnostic test results. The occurrence and severity of patient symptoms depend on their pyrogenic threshold, which increases (until saturation) with recent parasite exposure and decays over time 82Smith et al.2006. Severe episodes of malaria occur due to a high parasite density or due to co-morbidities 76Ross et al.2006. Malaria mortality can be a consequence of a severe episode or an uncomplicated episode with co-morbidity 76Ross et al.200677Ross and Smith2006. The model also takes into account neonatal deaths 76Ross et al.200677Ross and Smith2006. Immunity to asexual parasites prevents severe cases by decreasing the parasite multiplication rate within the host. Individual immunity depends on the cumulative parasite and infection exposure frequency, as well as maternal immunity in infants for several months 65Maire et al.2006.

The case management component of OpenMalaria describes the use of treatment for uncomplicated and severe cases and depends on access to health services and whether patients have previously been treated for the same episode 89Tediosi et al.2006. The disease model includes explicit PK/PD models that capture the process whereby drugs reduce the parasite multiplication rate in treated hosts 5Bertrand and Mentré2008106Winter and Hastings2011. Pharmacodynamics parameters are parameterised individually for each genotype to allow different degrees of drug susceptibility to be modelled.

The entomological component of OpenMalaria simulates the mosquito vector feeding behaviours and tracks the infectious status of mosquitoes 16Chitnis et al.2012. The periodicity of this model allows seasonal patterns of transmission to be captured. The probability that a feeding mosquito becomes infected depends on the parasite density within bitten individuals 75Ross et al.2006. No recombination is modelled between the different genotypes in the mosquitoes. The number of newly infected hosts depends on the simulated EIR of the vector model 16Chitnis et al.2012. The genotype of new infections is based on the genotype frequencies in humans from the previous five time steps 75Ross et al.2006.

Parameterisation of the treatment profiles

This study investigated factors influencing the establishment and spread of parasites resistant to three different treatment profiles.

The first treatment profile modelled was a short-acting drug administered as monotherapy. The short-acting drug has a short half-life and a high killing efficacy, simulating artemisinin derivatives (Figure 1A and B). We modelled the pharmacokinetics of the short-acting drug using a one-compartment model, which is considered sufficient when modelling short-acting antimalarials 46Kay et al.2013106Winter and Hastings2011. We varied key PK/PD parameters (half-life, EC50, Emax) in the global sensitivity analysis to assess their influence on the rate of spread of resistance. The EC50 ranged from 0.0016 to 0.009 mg/l to include the EC50 of artemether, artesunate, and dihydroartemisinin 46Kay et al.2013106Winter and Hastings2011. The half-life parameter ranges represented the values for artemether, artesunate, and dihydroartemisinin used by 46Kay et al.2013; 106Winter and Hastings2011 (Table 1). Note that in 46Kay et al.2013, the Emax of all short-acting drugs was equal to 27.6 per day. However, we varied the killing rate and included higher values to investigate its effects on the rate of spread (Table 1). To ensure that the short-acting drug killed the sensitive parasites efficiently for any combination of parameters, we extended the treatment course from a daily drug dose for 3 days to a daily drug dose for 6 days. Moreover, we parameterised the dosage and constant parameter values to that for dihydroartemisinin (Appendix 1—table 1), as it is the artemisinin derivate with the shortest elimination half-life and highest EC50 46Kay et al.2013106Winter and Hastings2011. By doing so, we also ensured that the short-acting drug had the typical profile of an artemisinin derivative.

The second treatment profile modelled was a long-acting drug administered as monotherapy. The long-acting drug had a long half-life and a lower Emax than the short-acting drug (Figure 1A and B), typical of partner drugs used for ACTs. We modelled the PK of the long-acting drug with a two-compartment model, which is more typical of the clinical PK of partner drugs 5Bertrand and Mentré2008. As for the short-acting drug, key PK/PD parameters (half-life, EC50, Emax, and dosage) were varied in the global sensitivity analysis. The EC50 ranged from 0.01 to 0.03 mg/l to include the EC50 of mefloquine, piperaquine, and lumefantrine used by 46Kay et al.2013; 106Winter and Hastings2011. The half-life range corresponded to the value reported for mefloquine, piperaquine, and lumefantrine in 12Charles et al.2007; 85Staehli Hodel et al.2013; 44Jullien et al.2014; 45Karunajeewa et al.2008; 64Maganda et al.2015 (Table 1). We increased the Emax range from 3.45 per day (as reported in 106Winter and Hastings2011) to 5.00 per day to investigate the effect on the rate of spread (Table 1). We also assessed the impact of Cmax on the rate of spread for the long-acting drug because the Cmax varies between ACTs partner drugs and has a strong influence on the post-treatment killing effect of the long-acting drug 36Hastings and Hodel2014. We varied drug dosage from 30mg/kg to 40 mg/kg to examine the influence of variation of Cmax on the spread rate for the long-acting drug. The lower limit of 30 mg/kg was fixed to ensure that the long-acting drug killed the sensitive genotype efficiently for any parameter combination. The treatment course involved a daily drug dose for three consecutive days. To ensure that the long-acting drug had the profile of typical partner drugs, the values of the constant parameters were parameterised to the values of piperaquine reported in 106Winter and Hastings2011; 85Staehli Hodel et al.2013 (Appendix 1—table 2).

The last treatment profile was a combination of short- and long-acting drugs, simulating ACT. We tracked the concentration of each drug independently. We used the same models, parameter values and ranges for the two drugs as when both drugs were used as monotherapy. However, the treatment course involved a daily dose of both drugs for 3 days, as recommended by the WHO for most ACTs 105WHO2021. In OpenMalaria, the killing effects of the two drugs were calculated independently and acted simultaneously on the parasites.

Parameterisation of the drug-resistant genotypes

For each simulation, we tracked two genotypes, one drug-resistant and one drug-sensitive. We investigated the spread of resistant parasites with different degrees of resistance (Table 1). We modelled the phenotype of drug resistance and the degree of resistance differently for each drug profile.

Previous studies have shown that parasites partially resistant to artemisinin exhibit an extended ring-stage during which they are not sensitive to artemisinin (even at high drug concentrations) but remain sensitive to the drug during other stages of the blood replication cycle 52Klonis et al.201395Wang et al.201778Sá2018107Witkowski et al.2013112Ye et al.2016. OpenMalaria does not model the specific drug-killing effect for the different steps of the blood-stage. As in 59Lohy Das et al.2017; 60Lohy Das et al.2018, we assumed that parasites resistant to the short-acting drug had a reduced Emax compared with sensitive ones (Figure 1B). This assumption captured the fact that, overall, the short-acting drug killed fewer resistant parasites than sensitive ones at any drug concentration because they are not sensitive to artemisinin during the ring-stage and that this stage-specific effect is best incorporated into PK/PD modelling by variation in Emax 40Hodel et al.2016.

Previous studies reported that parasites resistant to long-acting drugs typically have an increased EC50 11Chaorattanakawee et al.201610Chaorattanakawee et al.201587Tahita et al.2015. Thus, as in other models, we defined parasites resistant to the long-acting drug to have a higher EC50 than the sensitive ones (Figure 1B; 46Kay et al.2013; 106Winter and Hastings2011). With an increased EC50, the resistant parasites were less susceptible to the drug at low drug concentrations. Thus, these resistant genotypes were more likely to survive drug treatment and are more likely to successfully infect new hosts with higher residual drug concentrations 47Kay and Hastings2015.

Considering short- and long-acting drugs in combination, the resistant genotype was resistant to the short-acting drug. But in the global sensitivity analysis, both the sensitive and resistant genotypes could have some degree of resistance to the long-acting drug. The decreased susceptibility to the long-acting drug was the same for both sensitive and resistant genotypes, meaning that we assumed the two genotypes differed only in one mutation, which conferred resistance to the short-acting drug. This assumption allowed us to ignore the effect of recombination in the mosquitoes. In effect, this assumed that the allele defining the degree of resistance to the long-acting drug was fixed in the population.

Approach to identify the key drivers of the spread of drug-resistant parasites

Through global sensitivity analyses, we quantified how the factors in Table 1 influenced the spread of drug-resistant parasites for each treatment profile. First, we estimated the effect of each factor in a non-seasonal setting with a population fully adherent to treatment. Based on these results, we identified specific settings for further analysis, which used constrained sensitivity analyses to investigate the impact of varying drug properties and fitness costs in a fixed set of settings (i.e. in low and high transmission settings, with low and high treatment levels of monotherapy or combination therapy) and with a fixed degree of resistance. In these additional constrained sensitivity analyses, we also investigated the effect of drivers in seasonal transmission settings (based on the seasonality pattern of a setting in Tanzania [66Maire et al.2006, Appendix 1—figure 2]) and where populations adhere to either 100 or 67% of treatment doses.

Due to the number of factors investigated, each global sensitivity analysis required a large number of simulations (see details below) that is computationally infeasible for detailed individual-based models. Therefore, we trained an HGP 6Binois and Gramacy2021 on a limited set of OpenMalaria simulations (3500–11,500 simulations). We then used the trained emulator to predict the output of OpenMalaria for a large number of simulations and used these outputs to perform the global sensitivity analysis (Figure 1C), adapting a similar approach to 28Golumbeanu et al.2022 and 73Reiker et al.2021. Our approach involved: (i) randomly sampling combinations of parameters; (ii) simulating and estimating the rate of spread of the resistant genotype for each parameter combination in OpenMalaria; (iii) training an HGP to learn the relationship between the input (for the different drivers) and output (the rate of spread) with iterative improvements to fitting through adaptive sampling; and (iv) performing a global sensitivity analysis based on the Sobol variance decomposition using the trained emulator 50Kilian et al.2000. Each step of the workflow is detailed below.

Random sample combinations of parameters

We randomly sampled 250 different parameter combinations from the parameter space shown in Table 1 using an LHS algorithm 29Gramacy2007. The parameter ranges were defined as follows. We defined the ranges for the properties of the short-acting drug and the long-acting drug to include the typical parameter values of artemisinin derivatives and long-acting partner drugs, respectively 46Kay et al.2013106Winter and Hastings201112Charles et al.200785Staehli Hodel et al.201344Jullien et al.201445Karunajeewa et al.200864Maganda et al.2015. The range of the degree of resistance captured the spread of drug-resistant parasites, which vary from fully sensitive to having almost no drug sensitivity. The fitness costs were extracted from studies investigating the decline of chloroquine-resistant parasites after the drug pressure was removed 55Kublin et al.200369Mita et al.2003. The variation in annual EIR captured settings with low transmission to those with high transmission. The range of access to treatment captured settings with low to high level of access to treatment. The variation in the diagnostic detection limit captured the range of sensitivity of typical diagnostics used for malaria (such as rapid diagnostic test, microscopy, and PCR) 50Kilian et al.200071Murray et al.2008.

Simulate and estimate the rate of spread of the drug-resistant genotype

We quantified the rate of spread through the selection coefficient, a measure widely used in population genetics to assess the strength of selection on a genotype 38Hastings et al.2020. The selection coefficient is the rate at which the logit of the resistant genotype frequency increases each parasite generation and should be linear throughout the spread 38Hastings et al.2020. Population genetics theory often assumes an infinite population size to remove stochastic fluctuation of the allele frequency also called genetic drift 38Hastings et al.2020. However, in our model the parasite population size is finite, so stochastic fluctuations are present. Thus, we should avoid estimating the selection coefficient when there is a low frequency of the resistant genotype (from a small human population size, a low EIR, and a small initial frequency of the resistant genotype) because the resistant genotype may become extinct due to the stochastic fluctuation. In addition, the effects of genetic drift that occurs when a genotype is present at a low frequency may cause non-linearity during resistance spread which may obscure the estimation of the selection coefficient 38Hastings et al.2020.

Following the approach described in 38Hastings et al.2020, we assumed an initial percentage of infected humans carrying the resistant genotype of 50%. A high initial percentage minimises the impact of random fluctuation on our estimation, and the subsequent risk of extinction, without affecting our estimate because the selection coefficient was not frequency-dependent (Appendix 1—figure 3). We simulated the spread of resistant parasites in a human population of 100,000 individuals with an age structure typical of some countries in Africa (17.7% of people under 5 years of age) 25Ekström et al.2016. We ran each parameter combination on five stochastic realisations. The simulation started with a burn-in period of 100 years to reach the expected level of immunity in the population and an additional 30 years to reach EIR equilibrium (Appendix 1—figure 4). Both genotypes were sensitive to the drug during this period, so the percentage of infected humans carrying the resistant genotype remained stable. After the burn-in period, we introduced the fitness cost and the drug for which the resistant genotype had reduced sensitivity. We then estimated the selection coefficient, s, as,

where p(t) is the frequency of the resistant genotype in inoculations (the number of inoculations carrying the resistant genotype divided by the total number of inoculations resistant and sensitive genotypes), t is the number of parasite generations after introducing the new drug at t=0. We assumed that a parasite generation is 2 months (60 days) as in 38Hastings et al.2020. We started the regression at one parasite generation after introducing the new drug (at 60 days). We stopped the regression 12 generations later, at 720 days, because, as shown in 38Hastings et al.2020, it was computationally convenient and returned stable selection coefficient estimates. The regression was stopped sooner if the frequency of inoculations carrying the resistant genotype was higher than 90% or lower than 30% to prevent tracking a small number of a single genotype for which genetic drift is strong. In seasonal settings, the rate of spread of the resistant genotype varied throughout the year. Consequently, we estimated the selection coefficient using a moving average of the frequency of the resistant genotype in inoculations (Appendix 1—figure 5). This method prevented biasing the selection coefficient according to the period included in the regression.

Once the selection coefficient was estimated, it could be converted to the number of parasite generations needed for the frequency of the resistant genotype in inoculations to increase from p(1) to p(t),

We could then convert the number of parasite generations to time in years, a more relevant public health measure than the selection coefficient itself.

Train the emulator and improve its accuracy

We randomly split our data into a training dataset containing 80% of simulations and a test dataset containing 20% of simulations. We trained the HGP on the training dataset using the function mleHetGPfrom the R package ‘hetGP’ 6Binois and Gramacy2021. We chose to use HGP as it was successfully used in two previous studies that performed global sensitivity analyses of OpenMalaria 73Reiker et al.202128Golumbeanu et al.2022. In addition, 73Reiker et al.2021 tested different emulators and found that HGP provided the best fit with a limited number of simulations (analysis not shown in the published study). To assess the accuracy of the emulator, for the test dataset we assessed the correlation coefficient and root mean squared error between selection coefficients estimated with the emulator and selection coefficients estimated using OpenMalaria. We iteratively improved the accuracy of our emulator through adaptive sampling. Adaptive sampling involved resampling 100 parameter combinations in the parameter space where we were less confident (higher variation) in the HGP prediction and repeating the entire process until the emulator had a satisfactory level of accuracy. The satisfactory level of accuracy was defined based on the correlation coefficient and the root means squared error between the estimated selection coefficient and expected selection coefficient for the test dataset (Appendix 1—figures 6–!number(12)).

Global sensitivity analysis

Using the emulator, we undertook global sensitivity analyses using Sobol’s method 84Sobol2001. This method attributed fractions of the selection coefficient variance to each input 84Sobol2001. To do this, we first generated two random datasets with a sample size of 100,000 using an LHS algorithm 29Gramacy2007 that sampled within the parameter ranges of Table 1. When then estimated selection coefficients for these datasets with the trained emulators. Note that without emulators, we would have to run these simulations in OpenMalaria, which would not have been feasible due to computational requirements. We then used the function soboljansen from the R package ‘sensitivity’ to perform the global sensitivity analysis with 150,000 bootstrap replicates and the two datasets 15Cheng et al.2021. With this function, we estimated first-order and total Sobol' indices simultaneously. The first-order indices represent contributions of each parameter’s main effect to the model output variance. The total effect represents the contribution of each parameter to the model output variance considering their interactions with other factors. We report only the first-order indices in the Results section because we did not observe many interactions between these factors. Some parameters supported the spread of resistance (increased the selection coefficient), whilst others hindered the spread (decreased the selection coefficient). To visualise the direction of the effect of each parameter, we calculated the 25th, 50th, and 75th quantiles of the estimated selection coefficient of the two random datasets over the corresponding parameter ranges.

Establishment of drug resistance

As explained in the Introduction, the establishment of resistant mutations is a stochastic process that depends on the selection coefficient of the mutation and the heterogeneity of parasites reproductive success in the setting, which in turn depends on the transmission level and the health system strength 113Wiesch et al.201135Hastings200438Hastings et al.202033Hastings and Mackinnon199851Klein2014. Estimating the probability of establishment requires running many stochastic realisations due to the stochasticity of this step. To be more computationally efficient, we assessed the probability of establishment of a subset of 10 resistant genotypes with a known selection coefficient per setting and treatment profile. Based on the observed relationships between the selection coefficient and the probability of establishment for each treatment profile and setting, we could then extrapolate the probability of establishment of any mutations having a known selection coefficient.

To estimate the probability of establishment, we modelled the emergence of resistant genotypes in a fully susceptible parasite population. We used the approach described in 38Hastings et al.2020, in which resistant infections were imported into the population at a low rate. In OpenMalaria, imported infections have the same frequencies of genotypes as in initialisation, thus we cannot import only resistant infections. Therefore, to import resistant infections in a population infected only by sensitive parasites, we followed the step described below (Appendix 1—figure 13). We first defined a 50% frequency of resistant parasites in infected humans. The simulation started with a burn-in phase of 100 years, during which both genotypes were sensitive to treatment. This meant that the frequency of the resistant parasites was stable (at 50%). In the second phase, we introduced a drug to which resistant parasites were hypersensitive (the drug EC50 was 100 times lower in the resistant genotype than the sensitive one). The second phase ran for 100 years, and once complete, the parasite population was fully susceptible. In the third phase, we imported new infections at a rate low enough to ensure that the previously imported resistant genotype either established or went extinct before a new resistant infection was imported (Appendix 1: section 5.1). The third phase ran until the resistant genotype established (frequency of the resistant genotype in infected humans is equal to 50%).

The probability of establishment, P_e, can be estimated based on the average number of resistant infections that are imported until the resistant genotype establishes, N_e, as follows (the probability of a successful event can be estimated as one divided by the mean number of independent trials required to achieve the first success 19Dekking et al.2005),

We simulated 300 stochastic realisations, R, and estimated P_e, as,

where N_m,j is the number of imported resistant infections until the resistant genotype established in run j. Re-arranging the formula shows that Pe is equal to the number of resistant genotypes established in all stochastic realisations (this number is equal to R as only one resistant genotype established per stochastic realisation) divided by the total number of resistant infections imported into all stochastic realisations (includes resistant genotypes that became extinct and established). We estimated the 95% confined intervals of Pe (Wilson methods 19Dekking et al.2005), as,

Note that in each stochastic realisation, we estimated N_m, as,

where t_e is defined as the last time that the number of infections with a resistant genotype was equal to zero, that is the time (in years) until the arrival of the imported resistant infection that led to the successful establishment of the resistant genotype. N_i is the number of imported resistant infections per year. Note that OpenMalaria specifies the number of imported infections, V, in numbers of imported infections per 1000 people per year, and half of the imported infections were sensitive. Thus, the number of imported resistant infections that occurred until one established can be estimated as,

where N is the human population size. We simulated a population size of 10,000 individuals to increase computational feasibility for the large number of simulations required for our extensive global sensitivity analyses. A larger population was unnecessary, as the population size does not influence the probability of establishment unless it is extremely small 98Waxman2011. This was not the case in our simulation, which had a minimum of 3018 infections in the low transmission setting.

Data and software availability

We did not use individual participant-level data. Parameters values used in the model were informed from the literature as referred to in the main text or the Appendix. The source code for OpenMalaria was developed using the C++language and is available at https://github.com/SwissTPH/openmalaria090Thüring et al.2022and a documentation is available at https://github.com/SwissTPH/openmalaria/wiki. The analysis script was developed using the R software and is available at https://zenodo.org/badge/latestdoi/458217287. All data and codes used to produce the figures are available at https://zenodo.org/badge/latestdoi/458226427. In addition, the data used to produce the figure are included in the manuscript.