Role of seasonal importation and genetic drift on selection for drug-resistant genotypes of Plasmodium falciparum in high-transmission settings
Abstract
Historically Plasmodium falciparum has followed a pattern of drug resistance first appearing in low-transmission settings before spreading to high-transmission settings. Several features of low-transmission regions are hypothesized as explanations: higher chance of symptoms and treatment seeking, better treatment access, less within-host competition among clones and lower rates of recombination. Here, we test whether importation of drug-resistant parasites is more likely to lead to successful emergence and establishment in low-transmission or high-transmission periods of the same epidemiological setting, using a spatial, individual-based stochastic model of malaria and drug-resistance evolution calibrated for Burkina Faso. Upon controlling for the timing of importation of drug-resistant genotypes and examination of key model variables, we found that drug-resistant genotypes imported during the low-transmission season were (i) more susceptible to stochastic extinction due to the action of genetic drift, and (ii) more likely to lead to establishment of drug resistance when parasites are able to survive early stochastic loss due to drift. This implies that rare importation events are more likely to lead to establishment if they occur during a high-transmission season, but that constant importation (e.g. neighbouring countries with high levels of resistance) may produce a greater risk during low-transmission periods.
1. Introduction
Despite recent advances in malaria control resulting in a reduction of prevalence, Plasmodium falciparum malaria continues to be a major public health concern. The widespread use of artemisinin-based combination therapies (ACTs) has contributed to this reduction in prevalence, but increased usage of ACTs also increases the selective pressure on the parasites to develop drug resistance. Historically, the emergence of drug resistance has followed a pattern of first appearing in low-transmission settings, such as Southeast Asia and South America, followed by later migration to high-transmission settings. This was the case for chloroquine- and sulfadoxine-pyrimethamine-resistant P. falciparum [1–4], and more recently, artemisinin-resistant P. falciparum phenotypes which were identified in western Cambodia in 2007–2008 [5,6]. Since the identification of resistance-associated kelch13 point mutations [5], artemisinin resistance has been identified in other parts of Southeast Asia [7,8], Guyana [9,10], Rwanda [11,12] and Uganda [13–15]. Thus, developing a mechanistic understanding as to the cause of delayed emergence or slower evolution of drug resistance in high-transmission settings is particularly germane in the African context where a reservoir of kelch13 mutations currently exists and has the potential for rapid expansion [16].
Several mechanisms have been proposed to explain this pattern of slower drug-resistance emergence, establishment and/or evolution in higher transmission regions. First, due to higher population-level immunity to malaria, a new infectious mosquito bite is less likely to lead to malaria symptoms in a higher transmission region, resulting in a lower probability that a new infection will be treated by drugs [17–19]. Second, treatment coverage and access are generally lower in high-transmission regions, meaning that new symptomatic infections will also have a lower chance of facing treatment. Third, multi-clonal infections (i.e. infections in which the host is infected with several genetically distinct strains of the parasite) result in within-host competition which may suppress drug-resistant clones due to their cost of resistance or immune-mediated competition [20–23]. Within high-transmission settings, such as sub-Saharan Africa, multi-clonal infections are common [24], thus creating the relevant conditions for drug-resistant and drug-sensitive parasites to be present in the same host. A counteracting factor of this mechanism is that in the presence of the relevant drug therapy, the demise of drug-sensitive parasites in a multi-clonal infection may result in competitive release of the drug-resistant parasite and accelerate its spread [21]. Finally, higher rates of recombination in high-transmission regions may act against multi-genic drug-resistant genotypes by breaking up beneficial combinations of drug-resistance mutations [25] though much work remains to be done on this question.
A growing body of mathematical models suggests that a combination of within-host competition and immune-mediated symptomology are contributors to the cause of the delayed drug resistance emergence in high-transmission settings [18,20,21,26–31]. Of particular note is the work of Bushman et al. [21] who used an individual-based model (IBM) combined with ordinary differential equations to model within-host red blood cells, immune response and parasite dynamics to explore the role of within-host competition. The study found that resistant genotypes initially have a higher risk of extinction in high-transmission settings, but resistance can rapidly spread if extinction is avoided. These findings are supported by Whitlock et al. [20] using a similar IBM approach; however, their model also accounted for variations in the antigenic response to various strains of the parasite. Similarly, Masserey et al. [31] also used an IBM coupled with an emulator based approach to examine the impact that various factors such as drug pharmacokinetics/pharmacodynamics, treatment coverage, parasite biology and environmental factors have on the establishment of drug resistance.
Despite the complexity of the models that have been developed, the effects of spatio-temporal diversity on P. falciparum evolution has not been fully explored [32], and it is unknown in which epidemiological scenarios importation of drug-resistant parasites presents the most risk—a question that we explore here. In the context of the high-transmission regions of sub-Saharan Africa, the malaria burden is not uniformly distributed [33], and a country may contain regions of high and low transmission which may influence the evolutionary environment for resistant genotypes. Another limitation of prior studies is that even high-transmission regions can have significant seasonal variation in transmission patterns with periods of comparatively low transmission occurring outside of the peak transmission season. Accordingly, there has been increasing interest in the role that seasonality plays in malaria transmission, with recent studies suggesting that persistent asymptomatic infections allow for dry season survival [34]. Finally, while importation is a known mechanism through which drug-resistance has been introduced into various countries, the actual risk of establishment or fixation post-importation is unknown. A contributing factor is the inherent complication in surveillance efforts to monitor importation. While data are limited, in a retrospective study of 54 international travellers arriving in Italy from 2014 to 2015 with confirmed cases of P. falciparum malaria, nine genetic markers for drug-resistant genotypes were detected, suggesting that the rate of importation across national borders may be substantial [35], a finding echoed by an earlier study [36].
In this study we explore a straightforward importation mechanism for the introduction of drug resistance in high-transmission regions, through the application of a spatial IBM of malaria that was previously calibrated and validated for Burkina Faso [37]. This simulation also allows us to explore the role of seasonality and how importation of drug-resistant parasites may lead to the emergence and establishment of drug resistance in a realistic high-transmission context with seasonal variation and heterogeneity of malaria transmission. By restricting importation in the simulation to a particular month, we are able to calculate extinction probabilities and follow long-term trajectories to determine what times of year (and what importation rates) pose the most risk for the establishment and spread of drug resistance.
2. Methods
Here we use the term appearance to refer to the period following the importation of one or more artemisinin-resistant genotypes (called 580Y for short, using the most common allele found so far). Not all importations are successful, and an imported genotype may immediately go extinct (i.e. no further transmission), or have a brief period of transmission before going extinct. If the genotype is able to survive the action of genetic drift surrounding its appearance, we say that is has successfully emerged once its allele frequency is greater than 0.001 (10−3) allowing for possible progression to fixation as the dominant strain [38].
2.1. Simulation overview
We used a previously calibrated and validated spatial IBM of malaria and human movement in Burkina Faso [37,39,40]. As the simulation was designed and constructed with the intent of exploring the evolution of drug resistance in P. falciparum [39]¸ it has the appropriate components necessary to explore the mechanisms for delayed emergence without being constructed explicitly for it. The simulation models Burkina Faso as a grid of 10 936 25 km2 cells (approx. 273 400 km2) with 3.6 million simulated individuals whose distribution is consistent with the 2007 population (i.e. 25% of the population, see electronic supplementary material, S1, §2–4). Malaria transmission follows the holoendemic patterns of Burkina Faso with the median P. falciparum prevalence for ages 2 to 10 years of age (PfPR2-10) in a given cell ranging between 7.9% and 67.6% [41]. Transmission also follows a seasonal pattern, with transmission increasing at the start of the rainy season in late-May to early-June, peaking between August and October, and declining to a seasonal low in November [37]. On a regional level, the transmission season may be shorter or longer depending upon the length of the rainy season, with the northern Sahelian climatic region having a heightened transmission of about three months, while the season lasts for about five months in the southern Sudanian climatic region.
Upon model initialization, the simulated landscape initially consists of parasites that are chloroquine resistant, artemisinin and piperaquine sensitive, and either amodiaquine sensitive or resistant with a 50–50 probability. Following model burn-in, mutation by the parasite in the presence of the relevant therapy is enabled (except for the artemisinin resistance locus) based upon previously calibrated mutation rates [42], under the assumption that the large majority of mutation occurs during asexual blood stage replication [43]. We disabled de novo mutation to artemisinin resistance in order to focus on the effects of artemisinin resistance importation alone. When locally transmitted infections occur via new infectious bites, individuals may remain asymptomatic, or progress to clinical symptoms based upon their individual immune response. Upon presenting with clinical symptoms, individuals under 5 years of age follow treatment-seeking rates determined by previous Malaria Indicator Survey findings for their district (52.1% to 87.0%) [44], while individuals over 5 seek treatment at a rate that is 55% lower than the under-5 rate (23.4% to 39.1%). This treatment seeking rate increases at a rate of 3% starting in model year 2019, consistent with the expected expansion of treatment seeking by individuals [45]. Treatment seeking does not ensure that an ACT will be taken, as the private market accounts for 16.8% of treatments in the simulation, consistent with the local treatment landscape [44].
The individual immune response is summarized as follows, with the full scope elucidated in Nguyen et al. [39] and relevant changes included below. Upon being selected by the simulation to be bitten by an infectious mosquito, the individual undergoes a sporozoite challenge during which their immune response may result in sporozoites being cleared before the parasite enters the liver stage [37]. Based upon the individual's immune response the probability of infection ranges from 20% (high immune response) to 80% (low immune response). Successful infections proceed to the blood stage where the total parasitaemia , where R notes the specific P. falciparum clone, is initially set based upon a random uniform draw from the appropriate parasitaemia range (e.g. clinical infections range from 2000 to 200 000 parasites per microlitre of blood).
Clearance of the parasite by the immune system is then based upon the following calculation, fit by Nguyen et al. [39]:
As individuals are infected and clear infections, the individual immune response increases according to rates parametrized in Nguyen et al. [39], and the probability that any new infection will progress to clinical symptoms follows
If individuals go an extended period of time without exposure, decays with a half-life of 400 days [39,47], resulting in an increasing likelihood that an individual will experience malaria symptoms following a new infectious mosquito bite. Note that a half-life of 400 days means that 20% of individuals are expected to lose immunity within 100 days and thus within the duration of the low-transmission malaria season.
Individuals are infected in the population based upon the force of infection (FOI) of parasites present in a given cell, with individual host parasitaemia used to calculate the FOI of P. falciparum clone R, at time t, for all hosts n such that
Equation (2.4) describes the FOI for R where bi is the biting attractiveness of the host, function g describes the saturation of transmission probability with increasing parasite density as described by Ross et al. [48], and β represents a scaling factor used to calibrate the entomological inoculation rate for the given location within the simulation. Thus, is ultimately dependent upon the parasitaemia of R in individuals, resulting in a link to individual immune response and any fitness cost associated with a given parasite. As a result, a high fitness cost incurred by drug-resistant strains will result in a discounting of the FOI, consistent with the competitive advantage of the wild type in absence of drug pressure [49,50].
2.2. Study parametrizations
Using the previously prepared parametrization of Burkina Faso calibrated to the epidemiological situation as of 2017 [37]. While prevalence and case numbers in Burkina Faso are likely to have changed since 2017, the intent of this study is to explore the general behaviour of imported parasites and not exact forecasts of allele frequency. Using this calibration, two simulation studies were conducted, a comprehensive national-scale simulation, and a smaller limited-locality study.
For the first study, the only means for the parasite to acquire artemisinin resistance is through importation. This allows the role of seasonality on importation to be examined by controlling the month of importation, number of importation events per month (1, 3, 6 or 9), and parasite density of the imported infection (symptomatic or asymptomatic infection), yielding a total of 96 combinations. During the month of importation, an importation may occur on any day, with the exact number of imports on a given day determined using a Poisson distribution across the entire month. The location of importation is determined by a weighted draw across the entire population, based upon the total population in each cell within the simulation (electronic supplementary material, S1, §5). Due to the population distribution of Burkina Faso having a lower population in the Sahelian climate zone, this has the effect of biasing importations towards the more densely populated Sudano-Sahelian (containing Ouagadougou, the capital of Burkina Faso) and Sudanian (containing Bobo-Dioulasso, the second largest city in Burkina Faso) climate zones (electronic supplementary material, S1, §3–4).
Upon selection of a cell, the individual to be infected at the location is determined by a uniform random draw from all susceptible individuals at that location. The individual is then infected by an artemisinin-resistant P. falciparum parasite with the pfkelch13 allele 580Y, and the host parasitaemia level is set to the appropriate value for a symptomatic (between 2000 and 200 000 parasites per microlitre of blood) or asymptomatic (less than 1000 parasites per microlitre of blood) infection. In the event of a symptomatic infection, the individual may seek treatment on the basis of their age and the regional treatment seeking rate which ranges from 52.1% to 87.0% for an individual under 5 and 23.4% to 39.1% for over 5. If an individual seeks treatment, they receive either artemether–lumefantrine (68% of treatments), amodiaquine (12.4%), quinine (5.6%), dihydroartemisinin-piperaquine (4.9%), artesunate–amodiaquine (2.5%), artesunate (2.2%), artesunate, sulfadoxine/pyrimethamine (1.6%), chloroquine (1.5%) or mefloquine (1.3%) with the probability of a specific treatment based upon previous survey data and the make-up of the nationally recommended first-line therapies and private-market treatments [37,44,51]. Although importation is implemented by importing a singly infected individual, following importation multi-clonal infections are possible if the individual is infected by another clone. Following model burn-in, the simulation is allowed to run for 20 years. In order to ensure the statistical validity of the results, 50 replicates of each combination were run, for a total of 4800 replicates.
Following evaluation of a drug-resistant genotype establishing under various importation conditions, additional reporting was incorporated in the simulation to capture additional population immunity data and additional replicates (n = 50) in which mutations for artemisinin resistance—pfkelch13 allele 580Y in the simulation—occur using a previously determined mutation rate [42]. These scenarios allowed for an assessment of the mechanisms of delayed drug resistance or establishment to be evaluated, within the constraints of the mechanisms that are implemented within the simulation.
Finally, additional model validation was conducted to ensure conformity of multi-clonal infections and multiplicity of infections (MOI) to field conditions. Within the simulation, the proportion of multi-clonal infections fluctuates on a seasonal basis, and the ranges are in good agreement with previous studies. A sentinel site study in Nanoro, Burkina Faso, located in the Sudano-Sahelian climate zone, was conducted from September 2010 to October 2012 and recorded a mean MOI of 2.732 (±0.056) with a range of 1 to 7 parasite genotypes [52], compared with simulation results of 2.220 ± 0.455 (electronic supplementary material, S1, §1).
2.3. Limited-locality studies
Following completion of the national-scale studies, additional studies were conducted with a limited population or limited geographical scope intended to isolate or reproduce dynamics observed in the national-level model results. Specifically, these more limited models were used to explore if observed seasonal fluctuations in 580Y frequency and treatment seeking behaviour could help to explain observed transmission and infection dynamics. A total of three additional spatial models were prepared, all deriving from the same configuration used for the national-scale model: a single cell with a population of 100 000, a two-by-two grid with a total population of 300 000 individuals, and a three-by-three grid with a total population of 320 000 individuals. All models are based upon four different configurations in which the five-month seasonal pattern of the Sudanian zone is enabled or disabled, and treatment seeking is either balanced (i.e. 50% of under 5, 50% of over 5) or skewed based upon the national upper and lower bounds (i.e. 87% of under 5, 23.4% of over 5).
2.4. Statistical analysis
To examine whether the month of case importation is associated with the frequency of parasites with the 580Y allele, we performed a Kruskal–Wallis test to identify whether differences exist across months, followed by pairwise Wilcoxon rank-sum tests to identify which pairs of months yielded significantly different frequencies following importation. This was repeated, collapsing months into the high (June–October) and low (November–May) transmission seasons to compare these two time periods. This was followed by examining whether different months of importation are associated with a greater probability of emergence of the 580Y allele, defined by having a frequency greater than or equal to 0.001, we used chi-squared tests for proportions. An initial, global chi-squared test identified if any months differed, and subsequent pairwise tests identified months with different probabilities. Within months, we also tested whether symptomatic or asymptomatic importation was associated with probability of establishment. All p-values lower than 10−4 are reported as 10−4.
To show alignment among simulated time trends in total infections, treatment administration, immune response, and frequency of parasites with the 580Y allele, Spearman's correlation coefficients were calculated following visual inspection. For pairs of variables with oscillating trends that did not have aligned peaks and troughs, a lag was applied to align the oscillations, providing insight into how one trend follows another. The use of a lag is appropriate given the inherent delays associated with infection, presentation of symptoms, and transmission of P. falciparum infections. Correlation coefficients were calculated for the three climatic regions for scenarios consisting of de novo mutation, importation during the low-transmission season, and importation during the high-transmission season.
The stochastic model used in these analyses is constructed as a set of related and dependent mechanisms that are fully known at the time of implementation. Mechanistic effects are either present or absent in a simulation like this, and relationships among variables (e.g. correlations between simulation outputs) emerge from a fully known pre-programmed process rather than an unknown process that is under investigation. Therefore, significant p-values showing a difference between two sets of simulation outcomes simply confirm that different pre-programmed mechanisms were used to generate the two sets of results. Nevertheless, as a matter of convention in medical literature using results from stochastic simulations, we have included statistics comparing groups of simulation outcomes.
3. Results
3.1. Role of seasonality on importation
When varying the month of importation of drug resistance, the number of importations and the parasitaemia of the imported individual, there is a clear difference between extinction outcomes and sustained transmission outcomes when comparing low-transmission months with high-transmission months (figure 1). In six of the eight combinations for symptomatic/asymptomatic importation and importation count, drug-resistant genotypes are more likely to establish when imported during low-transmission months (p ≤ 0.0007, Wilcoxon rank-sum, electronic supplementary material, S2, table S2). In the scenarios of three asymptomatic importations or one symptomatic importation, future resistance frequencies appear lower for parasites imported during low-transmission periods (figure 1), but the differences are not statistically significant (p = 0.81 and p = 0.41, respectively) due to the large number of zeros in each set of simulations. This large number of zeros for configurations with one or three asymptomatic importations, or one symptomatic importation per month, underlines the influence of extinction and genetic drift during the importation process. If an imported parasite is unlikely to be sampled by a mosquito, then low-transmission periods would be associated with lower risk of establishment. This is most easily seen in the extinction paths in figure 2 where importation is rare (i.e. one importation event in a given month) and onward transmission occurs with low probability due to the asymptomatic nature of the imported infections; in these scenarios, extinction probabilities are higher for parasites imported during the low-transmission season. However, averaging across all scenarios, pairwise comparisons of months across the low season and the high season indicate that 580Y allele frequency after 10 years is likely to be a median 1.84-fold higher (interquartile range (IQR): 0.90–3.35) if the allele is imported during the low-transmission season (electronic supplementary material, S2, table S4), indicating that if importation events are common (i.e. several per month) low-transmission periods are associated with a higher starting frequency and a higher probability of emergence or establishment for the recently imported parasite.
While any imported parasite has a small chance of surviving past initial appearance, the likelihood of progressing to a frequency of 10−3, suggesting emergence has occurred and the parasite is likely to be observed, in any scenario was generally below 30% (figure 3). As expected, our analysis shows that a higher number of importation events is associated with a higher likelihood of eventual establishment (asymptomatic global χ2 = 109.3, d.f. = 3, p < 0.01; symptomatic global χ2 = 145.9, d.f. = 3, p < 0.01). Median probability of progressing past an allele frequency of 10−3 is 0.08 (IQR: 0.02–0.14; across 56 month-scenario combinations) when importation occurs during the low-transmission season, and 0.02 (IQR: 0.02–0.06; across 40 month-scenario combinations) during the high-transmission season. These results support the finding that when a drug-resistant genotype is imported, assuming it can escape the risk of random extinction, emergence is more probable for imports during the low-transmission season. When importation is common (nine asymptomatic imports per month, figure 4), only 2% to 8% of high-season importations (i.e. between June and October) reach a 580Y allele frequency greater than or equal to 10−3, whereas 4% to 18% of the low season importations reach a 580Y frequency greater than or equal to 10−3. When the months immediately following the high-transmission season are excluded (i.e. November and December), the range is 12% to 18%, suggesting a lagged effect of the seasonal mechanisms and their effects on the success of imported genotypes. Although there is no data source on successful and unsuccessful malaria parasite importations, classical results in population genetics state that newly arrived neutral genotypes in a population generally have a lower than 50% chance of survival past the early stochastic stage when they are rare, consistent with the simulated probabilities above.
3.2. Selection mechanisms during high- and low-transmission seasons
An examination of the short-term changes in malaria dynamics as the transmission season changes suggests three common effects that may influence the selection strength for drug-resistant genotypes: changes in drug coverage, changes in symptoms occurrence and changes in the multi-clonal nature of some infections. In a seasonal malaria setting, the age distribution of cases can change between low and high season, and if treatment coverage depends on age, then selection pressure for drug resistance will vary by season as the total population with clinical infections who seek treatment (i.e. total population treatment coverage) varies by season. A simple demonstration of this can be seen in a small population model (320 000 individuals occupying a 3 × 3 grid) with age-based treatment coverage and seasonal transmission presented as factors in the analysis. When both seasonality (based upon the five-month Sudanian zone in Burkina Faso) and age-based treatment coverage (87% coverage for under 5 and 23.4% coverage for over 5) are present, then treatment coverage varies seasonally, in this example between 54% and 60% (figure 5). The reason for this change is that frequent biting and frequent treatment in the under-5 population means that many bites are occurring while children have residual drug in their blood from a recent treatment. This means that the bite-to-symptoms ratio drops for children under 5 as we move from the low-transmission season to the high-transmission season; this effect is not seen or is weaker in older age groups. The under-5 and over-5 treatment coverages in this example were chosen as representational of the maxima and minima for each group based upon provincial coverages. In our national-scale simulations of Burkina Faso, population treatment coverage changes by 2% to 3% (absolute value) between the seasons resulting in weak to moderate changes in selection pressure (figure 6).
The nature of the simulation allows for the mean level of all individual immune responses to the parasite to be captured, denoted here as . As expected, follows a lagged seasonal cycle, with having the strongest Spearman correlation with number of infections three months prior (electronic supplementary material, S1, §6), which is independent of the number of importations across the six combinations of region and importation timing. This cycle of increasing and decreasing across seasons nominally affects the likelihood of symptoms and treatment seeking at a population level. One consequence of changing between seasons is that the fraction of malaria parasites currently residing in symptomatic patients versus asymptomatic patients changes as well. The quantity φ, defined as the ratio of symptomatic infections to all infections [17,19,53] gives a general description of what proportion of drug-resistant genotypes are currently experiencing positive selection resulting from treatment and what proportion are currently undergoing negative selection imposed by their fitness cost. This symptomatic fraction φ appears to be generally low, ranging from 0.07 to 0.12, with a higher proportion of infections subject to drug pressure at the end of the high season than during the middle of the low season. However, φ is out of phase with the treatment coverage suggesting that the net combined effect of treatment coverage and symptoms presentation may result in negligible changes in evolutionary pressure during the course of the year (figure 6).
The role of the individual immune response also works in tandem with the role of within-host competition occurring in multi-clonal infections. The high MOI—with a median ranging from 1.747 to 2.268 depending upon the scenario and climatic zone (electronic supplementary material, S1, §3)—along with the low frequency of resistant clones during the appearance and emergence of drug-resistant genotypes indicates that within-host competition between drug-sensitive and drug-resistant genotypes is occurring in some small proportion of malaria-positive individuals. The population-level strength of this effect can be examined within the simulation by comparing the proportion of multi-clonal infections carrying a 580Y clone with all multi-clonal infections (figure 7). The model outputs show that this proportion is highest when median MOI is lowest, towards the end of the high-transmission season, again showing that these seasonal forces are acting in opposition preventing the formation of a clear picture of when within-host competition against drug-resistant genotypes should be the strongest.
4. Discussion
Drug resistance has always presented a danger to public health goals. For malaria, however, there is time to prepare as antimalarial drug resistance typically emerges slowly and may take a decade or more to spread geographically. Here, we look at the effects that importation of drug-resistant genotypes has on the national-scale malaria epidemiology (modelled here as the high-transmission settings of Burkina Faso) and we ask whether imported drug-resistant genotypes are more likely to establish if importation occurs in the high-transmission season or in the low-transmission season. We show that genetic drift and starting frequency play an important role in an imported allele's future trajectory, but we are uncertain if change in selection pressure is substantial enough across transmission seasons to alter a drug-resistant genotype's evolutionary path after importation.
The major evolutionary-epidemiological gap identified in our simulation results is that seasonally changing selection pressures for drug resistance are not easily identified as such, and that the direction of change may not always be clear. Three common factors are known to affect drug-resistance evolution across transmission settings—treatment coverage, symptoms presentation and within-host competition in multi-clonal infections—but these factors do not align the same way between seasons in the same epidemiological setting as they do between countries that have different epidemiological settings. For example, in the seasonal setting presented here, treatment coverage goes up seasonally at the same time as symptoms presentation goes down. This is consistent with an expansion of cases during the high-transmission season, resulting in individuals less likely to seek treatment getting infected (figure 5). This expansion of cases also impacts the MOI and the proportion of multi-clonal infections harbouring at least one drug-resistant genotype is highest when median MOI is lowest, leading to an ambiguous picture or perhaps small evolutionary differences in terms of when within-host competition may be acting to reduce the frequency or relative density of drug-resistant genotypes. This appears to be an open question in the evolutionary epidemiology of antimalarial drug resistance evolution, whether selective forces in seasonal malaria settings always arrange themselves in opposing fashion.
An additional open question in the evolution of antimalarial drug resistance is why drug resistance tends to emerge in low-transmission regions, but not low-transmission pockets of high-transmission regions. Part of this question may be answerable if the low-transmission pocket lies in a seasonal malaria setting, in which case the evolutionary forces favouring drug resistance do not all align to point in the same direction at the same time, as is seen in the simulated scenarios here. In non-seasonal settings, the answer may be related to (i) movement between low- and high-transmission pockets, (ii) how long a local pocket has been low transmission, or (iii) general region-wide low access to antimalarial treatment with the low-transmission nature of any pockets determined primarily by environmental factors.
The traditional population-genetic effects of importation and drift do have their expected behaviours in our analysis of drug-resistance importation for malaria. Genetic drift may lead to extinction for imported mutants with higher probabilities of extinction associated with low importation rates and low transmission rates. Imported parasites that are not lost due to drift will progress to emergence and establishment more quickly if the initial importation event occurred in a smaller population, i.e. during low-transmission season.
As in all epidemiological modelling analyses, the model structure itself means that some limitations are present in the analysis and interpretation. First, the way that symptomatic importations are implemented may introduce some bias in favour of the parasite. Specifically, (i) immune response is ignored when an infection is imported, and (ii) treatment seeking behaviour is based upon where the individual resides; however, in practice individuals may be more (or less) likely to seek treatment if symptomatic when passing through a port of entry. Second, genetic background was not included as a factor in the analysis. Imported parasites may be distinct antigenically from resident parasites and may thus be less likely to be controlled by the immune response when transmitted locally in a new setting. Imported parasites that are resistant may also possess compensatory mutations if the history of drug-resistance evolution in the region of origin is old enough. Third, this study focused on importation of pfkelch13 mutants associated with longer clearance half-lives and high rates of treatment failure, but the treatment failure rate of pfkelch13 mutants depends strongly on the presence/absence of certain partner-drug mutations [54,55], as well as the genetic background that these mutations appeared on. Imported parasites with intrinsically high failure rates would probably have an easier time avoiding the effects of drift, spreading and establishing in the population. The success of artemisinin-resistant genotypes is facilitated by and strongly sensitive to the presence of partner-drug resistance mutations [42]. Finally, age, as is well known in malaria, is an important factor associated with malaria history and treatment seeking; and imported parasites in a younger patients will have different parasitaemia levels and different likelihood of seeking treatment depending on the patient's country of origin, recent history in high-transmission settings, and malaria immunity.
While the simulation is a reasonable approximation of Burkina Faso, the complex climatic zones of the country—including Sahelian, Sudano-Sahelian and Sudanian zones—suggest that it may be a challenge to extrapolate patterns observed in the model to the broader African context. In regions of sub-Saharan Africa with highly seasonal malaria transmission, importations may follow the same patterns observed in this model; however, in regions that are only nominally seasonal or have no seasonal transmission the mechanisms driving the establishment of imported drug-resistant genotypes are likely to be different. Spatial heterogeneity of prevalence and treatment access is likely to play a major role in all settings.
The critical public health conversation that this study has implications for is the future of molecular surveillance for specific P. falciparum drug-resistant genotypes that are at risk of being imported from one country to another. In particular, when transmission is highly seasonal, monitoring for known markers of drug resistance may be more important during the low-transmission season, particularly when samples from recrudescent patients are available. This may become increasingly relevant in the African context with the de novo appearance of the drug resistance markers 561H in Rwanda [11], along with 469Y and 675 V in Uganda [14]. In addition to the de novo appearance of these markers on a regional level, at least one instance of the 561H allele has been isolated in Uganda [14,15], suggesting that cross-border migration of drug resistance is already taking place.
A major factor for projecting the future evolution of drug resistance is the increased usage of ACTs within high-transmission settings. While the historical pattern has been for the establishment of imported genotypes following the evolution of drug resistance in low-transmission settings, the recent identification of de novo drug resistance in high-transmission settings [11,13,16] suggests that low-transmission appearance may simply be more likely but not a determinative rule for all drug-resistance emergence events. Given the role that individual immune response plays in creating an environment conducive to the evolution of drug resistance, understanding the possible impact of upcoming vaccinations (i.e. RTS,S/AS01) on selection for—or against—drug resistance by the parasites may play a role in speeding up or slowing down the selection of drug-resistant genotypes [56,57]. Nevertheless, despite some of the known effects of transmission setting on drug-resistance evolution—via differences in drug coverage and symptoms presentation, primarily—these effects do not appear to translate to differences between seasons in the same epidemiological setting. While drift and importation rate do appear to have their traditional effects on the success of recently imported genotypes into a new population, natural selection on drug resistance does not appear to be stronger in one part of the malaria season than another, or these selective differences could not be identified in the Burkina Faso-specific model parametrizations analysed here.
Ethics
This work did not require ethical approval from a human subject or animal welfare committee.
Data accessibility
The source code for the base mathematical model and analysis specific to this manuscript can be found from the GitHub repository: https://github.com/bonilab/malariaibm-spatial-BurkinaFaso-2022 [58]. Within the repository the dataset(s) supporting the conclusions of this article (and its additional files) are accessible from th GitHub repository: https://github.com/bonilab/malariaibm-spatial-BurkinaFaso-2022/tree/main/Data [59].
Supplementary material is available online [60].
Declaration of AI use
We have not used AI-assisted technologies in creating this article.
Authors' contributions
R.J.Z.: conceptualization, data curation, formal analysis, investigation, methodology, software, validation, visualization, writing—original draft, writing—review and editing; J.L.S.: formal analysis, methodology, writing—original draft, writing—review and editing; T.D.N.: methodology, validation, writing—review and editing; T.N.-A.T.: methodology, validation, writing—review and editing; K.T.T.: methodology, validation, writing—review and editing; A.F.S.: validation, writing—review and editing; M.F.B.: conceptualization, funding acquisition, project administration, writing—original draft, writing—review and editing.
All authors gave final approval for publication and agreed to be held accountable for the work performed therein.
Conflict of interest declaration
We declare we have no competing interests.
Funding
This work was supported by National Institutes of Health grants NIAID R01AI153355 (M.F.B., R.J.Z., T.D.N., T.N.-A.T., K.T.T.) and NIAID F32AI167600 (J.L.S.), and the Bill and Melinda Gates Foundation grants and INV-005517 to Pennsylvania State University (M.F.B., R.J.Z., T.D.N., T.N.-A.T., K.T.T.). The funders had no role in study design, data collection and analysis, decision to publish, or preparation of the manuscript.
Acknowledgements
Simulations described in this study were performed on the Pennsylvania State University's Institute for Computational and Data Sciences' Roar supercomputer.