Role of seasonal importation and genetic drift on selection for drug-resistant genotypes of Plasmodium falciparum in high-transmission settings

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 km² cells (approx. 273 400 km²) 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 (PfPR_2-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 $D_{\mathrm{R}}$, 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]:

where $D_{R,t}$ represents the parasite density at time t, $M_t$ represents the host immune response using a scale from zero (no immunity) to one (full immunity) with the real value between these two bounds as individual is exposed to the parasite and loses immunity over time, and C_R represents the fitness cost of the given strain with zero representing the wild-type (no fitness cost) to a maximum daily fitness cost of 0.006. The faster clearance of parasites when $C_R>0$ is the only mechanism through which within-host competition is realized in the simulation. The parasite density is updated asynchronously every 7 days ($K=7)$ and prior validation ensured that this did not differ from daily updating [39]. In the event of a new infection (resulting in a multi-clonal infection), or clearance of a prior infection, the parasite density is updated prior to the next scheduled 7-day update interval. However, if an individual receives an antimalarial, parasites are cleared based upon the daily killing rate following the parametrization in Nguyen et al. [46] through a standard concentration-effect curve. The standard therapeutic dose of artemether, artesunate and dihydroartemisinin have a daily killing rate of 99.9% on the wild-type parasite; however, resistant parasites require a higher therapeutic dose to achieve the same results, and accordingly have a lower daily killing rate. A total of 64 genotypes are included in this parametrization, with wild-type and mutant genotypes at the pfcrt 76 locus, pfmdr1 86 and 184 loci, pfkelch13 580 locus, copy number of pfmdr1 and copy number of plasmepsin genes.

As individuals are infected and clear infections, the individual immune response $M_t$ increases according to rates parametrized in Nguyen et al. [39], and the probability that any new infection will progress to clinical symptoms follows

where ${\theta }_{\mathrm{mid}}=0.15$ describes the point at which immunity confers a 50% chance of developing symptoms, and z describes the convexity of the relationship between the level of immunity and the likelihood of developing symptoms [47]. The probability of symptom occurrence can span the full range of zero to one (see figure S2 in Nguyen et al. [39]), and under this parametrization when $M_t\ge 0.245$ the probability of progression to clinical symptoms is less then or equal to 10%, and the same can be inferred for the mean population-level immunity variable, ϴ_pop.

If individuals go an extended period of time without exposure, $M_t$ 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

andwhere equation (2.3) describes the total parasitaemia $D_i$ of an individual, where the quantity of the parasite density of clone j in the host is given by ${\delta }_{\hspace{0.17em}j.i}$, with ${\gamma }_{\hspace{0.17em}j,i}$ representing the presence/absence of gametocyte production of the clone j, where one represents normal production and zero represents none. Due to the nature of the simulation, the residual gametocytaemia following a cured infection (i.e. the asexual parasitaemia is zero) is not simulated. To compensate for this the gametocytaemic period is shifted earlier so that infectious hosts have the same number of infectious days as real infections, and overall, this does not impact the simulation since gametocytocidal drugs are not included in the simulation.

Equation (2.4) describes the FOI for R where b_i 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, ${\mathrm{\Lambda }}_{t,R}$ 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⁻⁴ are reported as 10⁻⁴.

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.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.

• Download figure
• Open in new tab
• Download PowerPoint

• Download figure
• Open in new tab
• Download PowerPoint

While any imported parasite has a small chance of surviving past initial appearance, the likelihood of progressing to a frequency of 10⁻³, 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 χ² = 109.3, d.f. = 3, p < 0.01; symptomatic global χ² = 145.9, d.f. = 3, p < 0.01). Median probability of progressing past an allele frequency of 10⁻³ 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⁻³, whereas 4% to 18% of the low season importations reach a 580Y frequency greater than or equal to 10⁻³. 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.

• Download figure
• Open in new tab
• Download PowerPoint

• Download figure
• Open in new tab
• Download PowerPoint

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).

• Download figure
• Open in new tab
• Download PowerPoint

• Download figure
• Open in new tab
• Download PowerPoint

The nature of the simulation allows for the mean level of all individual immune responses to the parasite to be captured, denoted here as ${\theta }_{\mathrm{pop}}$. As expected, follows a lagged seasonal cycle, with ${\theta }_{\mathrm{pop}}$ 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 ${\theta }_{\mathrm{pop}}$ increasing and decreasing across seasons nominally affects the likelihood of symptoms and treatment seeking at a population level. One consequence of ${\theta }_{\mathrm{pop}}$ 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.

• Download figure
• Open in new tab
• Download PowerPoint