Immune selection suppresses the emergence of drug resistance in malaria parasites but facilitates its spread

Abstract

Although drug resistance in Plasmodium falciparum typically evolves in regions of low transmission, resistance spreads readily following introduction to regions with a heavier disease burden. This suggests that the origin and the spread of resistance are governed by different processes, and that high transmission intensity specifically impedes the origin. Factors associated with high transmission, such as highly immune hosts and competition within genetically diverse infections, are associated with suppression of resistant lineages within hosts. However, interactions between these factors have rarely been investigated and the specific relationship between adaptive immunity and selection for resistance has not been explored. Here, we developed a multiscale, agent-based model of Plasmodium parasites, hosts, and vectors to examine how host and parasite dynamics shape the evolution of resistance in populations with different transmission intensities. We found that selection for antigenic novelty (“immune selection”) suppressed the evolution of resistance in high transmission settings. We show that high levels of population immunity increased the strength of immune selection relative to selection for resistance. As a result, immune selection delayed the evolution of resistance in high transmission populations by allowing novel, sensitive lineages to remain in circulation at the expense of the spread of a resistant lineage.

In contrast, in low transmission settings, we observed that resistant strains were able to sweep to high population prevalence without interference. Additionally, we found that the relationship between immune selection and resistance changed when resistance was widespread. Once resistance was common enough to be found on many antigenic backgrounds, immune selection stably maintained resistant parasites in the population by allowing them to proliferate, even in untreated hosts, when resistance was linked to a novel epitope. Our results suggest that immune selection plays a role in the global pattern of resistance evolution.

Introduction

Although widespread drug resistance in Plasmodium falciparum malaria has been an ongoing public health challenge for decades, several observations about the evolution of resistance seem counterintuitive. For one, despite hundreds of millions of cases of malaria per year [1], resistance mutations arise surprisingly rarely. Resistance mutations occur readily in the lab [2], and the parasite population size and mutation rate ensure that every infection should contain multiple resistance mutations [3, 4]. However, widely circulating resistant lineages can typically be traced to just a handful of origins [4–9]. Furthermore, the points of origin are not, geographically, where one might expect. Sub-Saharan Africa suffers 90% of the global burden of malaria [1], but antimalarial resistance rarely originates there. Instead, resistant lineages typically originate in regions with relatively low transmission intensity, such as South America or Southeast Asia. Once introduced to Africa, however, these resistant lineages spread readily. For example, resistance to chloroquine originated five times in South America and Southeast Asia [10] but only appeared in Africa once introduced from Asia [8] in 1978. Within ten years, chloroquine resistance had been detected in every country in tropical Africa [11]. Together, these observations suggest that there is a dynamic that suppresses the spread of a resistance mutation immediately after its origin and, further, that this dynamic is intensified in high transmission settings. Unfortunately, the spread of an established resistance mutation is demonstrably less constrained. One by one, P. falciparum has evolved resistance to every drug introduced. The current front-line treatment, artemisinin combination therapy (ACT), is still generally effective in most countries, but artemisinin resistance is rising. The situation is especially dire in Southeast Asia, where the prevalence of a multidrug-resistant lineage exceeds 50% throughout most of the region and reaches as high as 80% at some sites [12]. Further, the average rate of treatment failure in response to one artemisinin combination exceeds 50% across sites in Cambodia and Vietnam [13]. With millions of lives at stake, fully understanding the evolution of resistance is urgent.

Transmission intensity shapes both the genetic make-up of the circulating parasite population and the within-host environment those parasites experience. Most malaria infections are genetically diverse [14–17], and this diversity increases in high transmission regions where hosts experience frequent, overlapping infections [18–21]. In a multiply-infected host, different lineages compete for the same resources while subject to both specific and nonspecific immune regulation. When drug-resistant and drug-sensitive parasites are present in the same host, sensitive lineages may suppress the growth and transmission of resistant lineages [22–24]. Frequent infections also result in high levels of adaptive immunity. Adaptive immunity is not sterilizing; instead, infections in a highly immune host remain asymptomatic [25]. For example, individuals in holoendemic regions of Africa receive upwards of 500 infectious bites per year [26] and symptomatic infections become increasingly rare after age five [25]. Conversely, in regions with low transmission intensity, hosts may never develop broad immunity. Within hosts, high levels of adaptive immunity are associated with increased clearance of resistant parasites [27]. While the influence of transmission intensity on diversity, competition, and immunity are relatively well-understood, the interactions between these forces may shape the evolution of resistance in complex ways.

Competitive suppression of drug resistant parasites by sensitive ones is often considered a consequence of a fitness cost of resistance. Estimates of the selection coefficient of resistance have been based on observational epidemiological data as well as direct competition assays in vitro [28–31]. Measuring the intrinsic fitness of any resistance mutation may be complicated by the presence of compensatory mutations [32–34]. Furthermore, the consequence of any potential fitness cost on selection will depend on additional factors, including the complexity of infection, host immunity, and treatment rate [35]. Competitive suppression can occur even when resistance is cost-free. For example, nonspecific host immunity can suppress resistant lineages. In high transmission settings, low-frequency resistant mutants will often arrive in a currently infected host whose active immune responses are likely to lead to the indiscriminate loss of rare variants [36]. Even in naive hosts, numerical dominance of a sensitive strain can lead to the loss of a rare resistant strain as density-dependent immune regulation of the entire parasite population begins [37]. The power of host immunity to limit the within-host parasite population size also decreases the likelihood that a resistance mutation occurs in the first place [38].

These ecological constraints occur independently of selection on resistant variants, but high levels of population immunity alter the adaptive landscape as well. Due to its strain-specificity, adaptive immunity generates selective pressure which favours the proliferation of lineages with novel epitopes in a partially immune host. All else equal, hosts gain immunity to a subset of circulating strains with each exposure. As the number of exposures increases for a particular host, so does the breadth of that host’s immune portfolio. Consequently, the proportion of circulating strains that are able to infect that host decreases, shaping the composition of future infections. This relationship has been observed within malaria patients, with less diverse infections found in older, more immune, individuals [39]. As a result, in highly immune host populations, parasites are subject to an additional dimension of selection in the form of selection for antigenic novelty.

While mathematical and computational models have been extraordinarily useful for investigating the processes driving patterns of resistance evolution in malaria parasites [36, 37, 40–49], as far as we are aware, none have explored the potential for interactions between selection for antigenic novelty and selection for drug resistance. Additionally, many have necessarily made simplifying assumptions, such as simulating homogeneous host populations, introducing mutations at a relatively high frequency, or using deterministic models of evolutionary dynamics. However, host heterogeneity, such as an individual’s exposure history, will affect the environment of a new mutation [41], a rare mutation will be subject to different selective pressure than one at higher frequency, and Plasmodium life history is characterized by severe bottlenecks alternating with massive population expansion, intensifying the contribution of both selection and drift [45, 48, 50, 51]. As a result, it is difficult to draw conclusions about the origin of a resistance mutation using these prior approaches.

Here we develop a multi-scale, agent-based, stochastic model of Plasmodium infection, transmission, and evolution in a host (and vector) population. We explicitly track parasites and their genomes; hosts, their immune systems, and their red blood cell resources; mosquito vectors; and the interactions between all three organisms. Importantly, both resistance mutations and antigenic mutations are random and recurrent. Because resistance is not inherently linked to any specific antigenic background, selection works independently on both. While our model makes some explicit assumptions (e.g., host and vector population sizes, the number of antigenically-unique parasite strains), evolutionary drivers like transmission intensity, multiplicity of infection, genetic diversity in the parasite population, and characteristics of host immunity are emergent properties. We show that immune selection interferes with selection for resistance in high transmission settings, where broad population immunity increases the strength of selection for antigenic novelty. As a result, sensitive lineages are maintained in the population, antagonizing the early spread of resistance. However, once resistance is common, immune selection promotes its maintenance regardless of transmission intensity. Additionally, we find that immune novelty is necessary to maintain resistance at a stable, high prevalence. Our results provide a plausible mechanism to explain the observed pattern of resistance evolution, in which high transmission intensity impedes the origin, but not the spread, of resistance.

Results

Endemicity and infection diversity increase with transmission intensity and number of antigenic strains

To briefly summarize the simulations, each vector, parasite, and host was individually modelled. Each parasite had a discrete, three locus genome. Two of the loci were associated with both resistance and growth rate, such that a mutation at either or both produced full drug resistance as well as a cost in the form of a growth rate reduction (though we also explore the effects of removing this cost, below). The allele at the third locus determined interactions with the host immune responses. We refer to all parasites that share an allele here as a “strain”, which is equivalent to a serotype. Mutations at each locus were random and recurrent. At the time of population initiation, every parasite was genetically identical. Over time, genetic diversity evolved through the forces of mutation, selection, and drift. The maximum number of possible strain variants was capped at either one strain, representing no antigenic diversity, or 30 strains. Because strain identity evolved independently of resistance, a single strain may contain both sensitive and resistant lineages. Parasites reproduce asexually within the host and undergo sexual recombination within the vector. After a within-vector maturation period of ten days, parasites will be transmitted the next time the vector feeds. An average of 12 parasites are transmitted to a host per feeding. During an infection, hosts developed strain-specific immunity, which shaped the density and duration of current and future infections (S1 and S2 Figs). Populations were founded with two different vector population sizes, producing low and high transmission intensity, and simulations were allowed to run for 2000 time steps, by which point genetic diversity in the parasite population and measures of endemicity in the host population had reached an equilibrium state.

Differences in transmission intensity in our model, driven by the number of vectors or strains, produced epidemiological patterns which reflected those observed in real populations. Using untreated, replicate populations at equilibrium, we measured aspects of endemicity and population dynamics. With 300 vectors and one strain, each host received an average of nine bites per year from infectious mosquitoes, akin to the entomological inoculation rate, or EIR. With higher numbers of vectors and strains—1200 and 30, respectively—mean EIR was over 200 (Fig 1A). EIR was lower with one strain because, without antigenic novelty, transmission was rare from any infection that followed a host’s primary infection (S3 Fig), due to suppression of the parasite population by host immunity.

[Figure omitted. See PDF.]

Fig 1. Baseline values in equilibrium populations.

Ten populations for each condition were allowed to evolve to equilibrium, and mean values were measured over 100 time steps. Figures show the density distribution of the data. (A) Annual entomological inoculation rate (EIR) (B) Percentage of hosts exposed (C) Percentage of hosts with active infections (D) Complexity of infection, or average number of strains per infection. (E) Average number of strains to which hosts had been exposed.

https://doi.org/10.1371/journal.pcbi.1008577.g001

At equilibrium, the proportion of hosts with active infections and the proportion of hosts which had been exposed increased with number of strains and vectors (Fig 1B and 1C). With 30 strains, almost all hosts had been exposed and most were actively infected throughout the sampling window, indicating that hosts experienced repeated, overlapping infections, beginning early in life. With one strain, the proportion of infected hosts was less than the proportion of exposed hosts, indicating that hosts recovered from infections before potentially being infected again (Fig 1B and 1C). In accord with empirical measurements, high transmission intensity produced more complex infections (Fig 1D) and hosts were exposed to more strains earlier in life (S4 and S5 Figs).

Because adaptive immunity is strain-specific, looking only at the fraction of hosts exposed to infection misses subtleties of the immune (and therefore fitness) landscape experienced by parasites. For example, there were only modest differences in the proportion of hosts exposed between the one strain, 1200 vector condition and the 30 strain conditions. However, the proportion of circulating strains to which those hosts were exposed varied substantially. In the former case, there was only one strain, so no antigenic novelty remained; when there were 30 possible strains, hosts on average were exposed to half and two-thirds of those strains in the low (300) and high (1200) vector conditions, respectively (Fig 1E). As a result, despite similar overall levels of exposure in the host population, differences in strain-specific exposures gave rise to different levels of population standing immunity, which we define as the likelihood that an arbitrary host has immunity to an arbitrary parasite. Operationally, we defined this metric as the average percent activation of immunity against a particular strain across all hosts. Thus, the conditions could be ranked from highest population immunity to lowest as: one strain and 1200 vectors, one strain and 300 vectors, 30 strains and 1200 vectors, and 30 strains and 300 vectors.

Increasing the vector number changed the patterns of strain-specific immunity in the host population over time (Fig 2). In populations with 300 vectors, average immunity to individual strains fluctuated over time, frequently to the point that immunity to a particular strain had almost completely waned in the population. With 1200 vectors, mean immunity was higher, more stable, and more evenly distributed across strains, suggesting that frequent exposure maintained strain diversity. High levels of population immunity kept strains circulating at intermediate frequencies by strengthening selection for rare strains and against common strains. This form of selection, known as negative frequency-dependent selection, is consistent with field observations of antigenic diversity [52–55].

[Figure omitted. See PDF.]

Fig 2. Population-wide mean strain-specific immunity over time in untreated equilibrium populations.

Immunity to each strain is averaged over all hosts in the population at each time point and is represented as percent activation, where 1 represents the maximum possible strain-specific immunity and 0 represents no immunity. Left: 300 vectors. Right: 1200 vectors. Data shows one representative simulation for each condition.

https://doi.org/10.1371/journal.pcbi.1008577.g002

Low population immunity facilitated the spread of resistance

After populations reached equilibrium, we introduced drug treatment at a set rate, such that at any given time, 30% of hosts with symptomatic infections were treated. An infection was considered to be symptomatic if parasites exceeded a threshold abundance within the host. Population endemicity dropped immediately after treatment introduction, but rebounded over time as resistance evolved and spread. We allowed evolution to proceed until the prevalence of treatment-resistant infections in the host population reached an equilibrium.

We found that antigenic diversity was the most important factor modulating the evolution of resistance, but there were additional, complex interactions driven by the number of vectors and feedback from the changing prevalence of resistance. We defined two metrics to summarize evolutionary dynamics. First, Tfail is the time at which 10% of infections were composed of at least 10% resistant parasites. Resistance readily evolved and spread to this low prevalence among infected hosts by 200 time steps in all conditions (Fig 3A). The 30 strain and 300 vector condition reached Tfail most quickly (p <.001 and A = 0.1864 (large effect size) for the most similar pairwise comparison), but the decrease was relatively small, and there was no significant difference between the other conditions (p = .541 for the least similar pairwise comparison). The similarity in Tfail between conditions demonstrates that later results were not due to differences between mutation supply or a lack of selection for resistance.

[Figure omitted. See PDF.]

Fig 3. Time to resistance.

Mean time to resistance was measured over 25 replicate populations. A: 10% prevalence of treatment failure (Tfail). B: 75% prevalence of treatment failure (Tubiq). Treatment failure is assumed to occur when 10% of the parasites within an infection are genetically resistant. Figures show the density distribution of the data. Note different Y axes.

https://doi.org/10.1371/journal.pcbi.1008577.g003

Second, we defined Tubiq as the time at which 75% of infections were composed of at least 10% resistant parasites. In simulations with one strain, populations never reached Tubiq (Fig 3B). Instead, equilibrium prevalence of treatment-resistant infections was relatively low, with large fluctuations over time (Fig 4). The prevalence of resistance was lower with 1200 vectors because resistant lineages could only reach densities sufficient for transmission in naive hosts, and there were fewer naive hosts in high transmission populations (p <.001, A = 0.122 (large effect size)). At higher treatment rates, resistance evolved more quickly and to a higher prevalence. Regardless of treatment rate, one strain conditions had a lower equilibrium prevalence of resistance and a slower rate of resistance evolution relative to the 30 strain conditions (S6 and S7 Figs). Within hosts, mixed infections—those harboring both sensitive and resistant parasites—were rare (Fig 5A and 5C); in most infections, the parasite population was either almost entirely sensitive or almost entirely resistant. This is because the relative fitness of a resistant lineage in an infection was entirely dependent on the host’s treatment status. Sensitive lineages were rapidly eliminated in treated hosts, but in an untreated host, sensitive lineages outcompeted resistant lineages. As a result, sensitive lineages continued to proliferate in the population.

[Figure omitted. See PDF.]

Fig 4. Equilibrium resistance prevalence.

After treated populations had reached equilibrium, the average prevalence of treatment failure (10% within-host frequency of resistance) was measured over 100 time steps in 25 replicate populations. Figures show the density distribution of the data.

https://doi.org/10.1371/journal.pcbi.1008577.g004

[Figure omitted. See PDF.]

Fig 5. Evolutionary trajectory of resistance.

Following introduction of treatment at time 0, the prevalence of within-host frequencies of resistance was monitored over time. The top (purple) line indicates the prevalence of all treatment-resistant infections, defined as infections consisting of at least 10% resistant parasites. Moving down from the top, the prevalence of infections consisting of between 10% and 50% resistant parasites is shown in blue. Below that, the teal section shows the prevalence of infections in which the frequency of resistant parasites is between 50% and 90%. Finally, the bottom section in green shows the prevalence of infections in which resistance is at near fixation within the host, with a frequency exceeding 99.9%. Top row (A and B): 300 vectors. Bottom row (C and D): 1200 vectors. Left column (A and C): one strain. Right column (B and D): 30 strains. The minimal separation between lines in the one strain conditions demonstrates that mixed infections are rare, compared to 30 strain conditions. Data shows one representative simulation for each condition.

https://doi.org/10.1371/journal.pcbi.1008577.g005

With 30 strains, treatment-resistant infections reached a high, stable prevalence that exceeded our threshold for ubiquity in both vector conditions although, as in the one strain conditions, prevalence was lower with 1200 vectors (p <.001, A = 0.399 (small effect size)) (Fig 4). Higher levels of exposure (i.e., more vectors) also slowed the spread of resistance (Fig 3B). Tubiq was approximately 250 time steps with 300 vectors, compared to over 1000 time steps for 1200 vectors. Although increased strain number increased the prevalence of treatment-resistant infections in the population, it also suppressed the frequency of resistance within individual hosts. A substantial and stable proportion of infections contained both sensitive and resistant parasites at equilibrium (Fig 5B and 5D). With 1200 vectors, mixed infections were common before equilibrium as well, but they were rare with 300 vectors.

The differences can be understood in terms of an adaptive landscape, in which the parasite traits under selection were growth rate, resistance, and immune novelty. Because all reproduction within the host was clonal, selection on the resistance trait depended critically on its antigenic background. In the case of a rare strain carrying the drug-sensitive allele, or a common strain carrying a resistance mutation, selection at one locus interfered with selection on linked loci. This decreased the efficiency of selection and produced clonal interference, in which both variants were maintained in the population, preventing resistance from going to fixation. This was responsible for the ongoing presence of sensitive strains and the delayed Tubiq in high transmission populations, where immune novelty was more rare. Linkage could also facilitate selection. Resistance arising in an untreated host on a strain to which that host had previously been exposed was quickly purged due to the combined effects of immunity and the cost of resistance. In contrast, a resistance mutation linked to a novel strain in a treated host had multiple within-host advantages, and thus was under the strongest positive selection. Even in an untreated host, a resistance mutation linked to a novel strain could hitchhike to high prevalence. In this manner, immune selection partially compensated for the cost of resistance and stabilized resistance in the population. The breadth of strain-specific immunity in the host population determined the probability of each of these scenarios. In general, the lower the population immunity, the higher the likelihood that resistance mutations originated in a host environment in which they were highly beneficial.

Removing cost decreased the impact of population immunity on resistance evolution

The reduction in growth rate caused by resistance could have several consequences, such as altering immune dynamics or producing less dense infections which affect transmission. High transmission also increases the likelihood that hosts are co-infected with both sensitive and resistant lineages, introducing the possibility that within-host competition suppresses the evolution of resistance. To eliminate these differences, we allowed populations to evolve as before, but removed the cost of resistance.

As expected, the absence of cost accelerated the rate of resistance evolution (i.e., decreased Tfail and Tubiq) in all conditions, but the effect was most significant with one strain, where populations were able to evolve stable resistance and reach ubiquity (Fig 6B), which did not occur when there was a cost to resistance (recall Fig 3B). With only one strain, a higher vector number allowed resistance to reach ubiquity more rapidly (p = .005, A = 0.268 (medium effect size)), demonstrating that high transmission intensity dispersed mutations more rapidly throughout the host population. While we cannot rule out contributions from within-host ecological processes, such as prior residency or the loss of rare variants due to nonspecific host immunity, these factors are expected to have the greatest effect in high transmission settings. The fact that high transmission sped up the spread of resistance shows that the advantage conferred by increased dispersal rate outweighed any other disadvantages.

[Figure omitted. See PDF.]

Fig 6. Time to resistance with costless resistance.

Mean time to resistance was measured over 25 replicate populations with no cost to resistance. A: 10% prevalence of treatment failure (Tfail). B: 75% prevalence of treatment failure (Tubiq). Figures show the density distribution of the data. Note different Y axes.

https://doi.org/10.1371/journal.pcbi.1008577.g006

In both of the 30 strain conditions, removing cost sped up the evolution of resistance (Fig 6). This effect was most significant with 1200 vectors. However, unlike in the one strain conditions, high transmission intensity did not accelerate the spread of resistance relative to the low transmission condition (p = .449) because immune selection dynamics counterbalanced the advantage conferred by increased transmission opportunities. In other words, immune selection alone slowed the spread of a beneficial mutation.

When resistance was costly, a resistance mutation could only proliferate in untreated hosts if it was genetically linked to a novel strain. Without cost, and without multiple strains, resistance mutations were only selectively neutral or beneficial, depending on the treatment status of the host. Therefore, transmission rate was the major factor determining their rate of spread. With multiple strains, a lineage carrying a resistance mutation may still be selected against within the host environment if host immunity recognized the linked strain, but selection against it would have been weaker than when resistance was costly. Similarly, a resistance mutation linked to a novel strain would have been under positive selection, even in untreated hosts, increasing its ability to hitchhike to a high prevalence. This hastened the rate of resistance evolution.

Patterns of within-host frequency were qualitatively similar to simulations in which resistance is costly (S12 Fig). With 30 strains, mixed infections were present at equilibrium, regardless of the number of vectors, as well as before equilibrium with high transmission. Mixed (sensitive and resistant) infections were rare with one strain. This supports the notion that the suppression of within-host frequency observed with costly resistance in 30 strain conditions can be attributed to immune selection, rather than solely to cost.

Immune selection for novel strains antagonized selection for resistance

In order to understand how immune interactions drive the evolution and maintenance of resistance at the parasite population level, we monitored strain dynamics over time, including strain-specific immunity averaged across all hosts, the frequency of each strain within the entire circulating parasite population, and the proportion each strain contributed to the total level of drug resistance in the parasite population.

Shortly after treatment was introduced, the number of circulating strains dropped precipitously (Fig 7), as did the host population infection rate and the total size of the parasite metapopulation (i.e., all parasites, summed across all infections). These losses were especially severe with 300 vectors simulations, in which strain diversity was lower prior to treatment introduction.

[Figure omitted. See PDF.]

Fig 7. Strain-specific patterns of resistance.

Following introduction of treatment at time 0, immunity and resistance were monitored for each strain over time. Top row: 300 vectors. Bottom row: 1200 vectors. Left column: Mean population immunity to each strain, as in Fig 2. Middle column: the frequency of each strain in the population. Right column: The relative contribution of each strain to the total amount of resistance in the entire parasite population. The first line, in blue, denotes Tfail, and the second line, in red, denotes Tubiq. Data shows one representative simulation.

https://doi.org/10.1371/journal.pcbi.1008577.g007

With 300 vectors (Fig 7, top row), resistance rapidly evolved on the few strains which persisted after treatment introduction. In a phenomenon resembling a selective sweep, these resistant strains steadily increased in prevalence until resistance reached ubiquity. Most strain diversity was transient, but the frequency of a highly resistant strain steadily increased (strain 30 in this representative simulation; middle column), despite high levels of immunity against it. Strain 30 only declined after ubiquity, when resistance was found on many strains (Fig 7, right column, and S13 Fig). This demonstrates that there was selection for immune novelty, but also that it failed to outweigh selection for resistance until resistance was ubiquitous.

With increasing vector number, more strains carried resistance mutations, producing more co-circulation and widespread standing adaptive immunity. Novelty in one host was less likely to translate to novelty in another host due to broader immune portfolios, so none of the strains were able to sweep to a high prevalence. Unlike with lower vector numbers, strain diversity remained high, even though most circulating strains were sensitive. Resistance only became ubiquitous once it was found on multiple strain backgrounds, indicating that balancing selection for strain diversity remained strong relative to selection for resistance (compare bottom row, center and right column in Fig 7). Additionally, changes in strain frequency tracked closely with changes in strain immunity, allowing rare sensitive strains to increase in prevalence. For example, at time 500, there was very little immunity to strains 1 and 2. Subsequently, they rose in frequency despite not carrying resistance mutations (S13 Fig). This demonstrates that, under high transmission intensity, immune selection outweighed selection for resistance.

The effect of vector number on pre-equilibrium strain diversity explains the difference in pre-equilibrium within-host frequency. With 300 vectors, only a few, highly resistant strains were circulating, and selection for resistance was much stronger than selection for novelty. This allowed quasi-fixation of resistance within hosts and promoted the rapid rise to ubiquity (Fig 5B). As the frequency of resistance rose, the strength of selection for resistance weakened. Immune selection only became strong enough, relative to selection for resistance, to suppress within-host resistance frequency after Tubiq. This was not the case with 1200 vectors, in which within-host frequency was suppressed throughout (Fig 5D). Immune selection was strong enough to maintain sensitive strains in the population, even when selection for resistance was at its strongest. While a higher transmission rate could disperse resistance more quickly, it also meant that sensitive lineages continued to be transmitted. The ongoing transmission of sensitive lineages delayed Tubiq.

At both transmission intensities, the strain background of the resistance trait had the greatest consequence shortly after the origin of a new resistance mutation, while it was linked to a single strain. If immunity against the strain background was high, resistance could be lost due to selective interference between the traits. However, as resistant lineages rose in frequency, most strain mutations happened on a resistant background, increasing the antigenic diversity of resistant parasites and reducing interference. The suppressive impact of immune selection was weakened because resistance was no longer associated with a minority of strains. As a result, the availability of antigenic novelty maintained a high, stable prevalence of resistance.

Discussion

Here we have demonstrated that there is a complex relationship between immune selection within hosts and the spread of drug resistant parasites between hosts, and that the net relationship changes based on the prevalence of the resistance trait in the parasite population. Immune selection was antagonistic during the initial spread of resistance. As the frequency of a rare resistant lineage increased, so did population immunity against that lineage, and therefore the strength of immune selection against it. This effect was strongest when transmission intensity was high and resulted in a delayed time to ubiquitous treatment failure. However, once resistance was common, immune selection maintained costly resistance at a high prevalence within the population, in part because linkage to a novel strain could counteract the cost of resistance. Although the antagonistic effect of immune selection was intensified by cost, immune selection suppressed within-host frequency of resistance even when resistance was not costly because sensitive lineages on a novel strain background could outcompete resistant lineages linked to common strains.

These results expand on other studies that found that within-host ecological processes were sufficient to suppress resistance [36, 37]. A significant difference in this investigation is that recurrent strain mutation meant that immune selection alone was likely to produce an antigenically-diverse infection, and that a resistance mutation could be found on any strain. This increased the complexity of within-host interactions and introduced negative frequency-dependent selection on strain background. Competitive suppression and immune selection are not mutually exclusive and may work in unison in natural populations.

Immune selection provides a plausible explanation for the observation that resistance does not usually originate in high-transmission regions, but spreads rapidly once introduced. African-origin resistance mutations may be suppressed by population immunity, but Asian Plasmodium spp. lineages have different common epitopes than those in Africa [56], and it is likely that Asian antigens are less cross-reactive with native African antigens. This would put a newly-introduced Asian strain under positive immune selection, allowing it to spread widely, along with any associated drug resistance traits.

There are many other factors we did not explore here. Due to limitations of computational feasibility and the complexity of the task, this work was constrained to a small region of the parameter space, including a steep cost of resistance and a fixed cross reactivity that assumed equal relatedness between strains. Wide ranges of values have been measured for many of these traits, and further, these values may change over time [32]. Computational feasibility also limited the size of the host population which, in turn, required us to use short host lifespans to prevent populations from becoming uniformly saturated with asymptomatic infections, regardless of transmission intensity. Many avenues for future investigation remain. For example, full resistance was produced by a point mutation, but resistance to some treatments, including ACTs and sulfadoxine-pyrimethamine, is multigenic, while other resistant phenotypes are associated with a host of background mutations [57]. Epistatic interactions could influence the evolution of multilocus resistance [41], as would the presence of a tolerant intermediate [41, 45] or compensatory mutations [34]. Although we found no effect from recombination (S9 Fig), this may be attributable to the simple genetic architecture and short genome, both of which decreased linkage disequilibrium, weakening any potential impact from recombination [58]. Additionally, vector dynamics were highly simplified. Plasmodium spp. have many genes that are under selection within the mosquito [59, 60], meaning that the malaria parasite is subject to selection in multiple environments over the course of its life cycle. This complicates predictions about the spread of resistance based on within-host dynamics. Finally, the environment was static, but malaria is often seasonal in low transmission regions, which has implications for selection and complexity of infection [19]. This model is intended to be versatile and adaptable, and could be expanded to include any of these factors.

With rates of treatment failure rising and few alternatives available, a further understanding of immune selection could provide another tool in the fight to preserve the efficacy of existing antimalarials. This is more critical than ever now that artemisinin resistance has been detected in Sub-Saharan Africa [61, 62]. One potential strategy is drug cycling, in which drugs are withdrawn from use as the rate of treatment failure rises in a population and then reintroduced once resistance wanes [63]. Other investigations have supported the simultaneous use of multiple first-line therapies [46, 64]. Quantifying the relationship between population immunity, resistance prevalence, and strength of selection for resistance could be used to identify frequencies of resistance which act as tipping points, in order to most effectively determine when to use individual drugs and to identify optimal treatment combinations. Immune selection may also produce an unforeseen benefit of vaccination. To date, attempts to make a widely effective vaccine against P. falciparum have been challenged by its antigenic diversity and its extraordinary ability to evade the immune system [65, 66]. The benefits of the most effective vaccine to date, RTS,S/AS01, are modest and short-lived. It prevents approximately 36% of infections [67] in the first year after vaccination, but efficacy drops nearly to zero within four years [68]. However, while the partial immunity produced by an imperfect vaccine may not totally prevent infection, it might alter immune selection dynamics, increasing the strength of selection for antigenic novelty and potentially interfering with selection for resistance. Similarly, bed nets and vector control campaigns have been highly effective at reducing malaria transmission, but it is not clear how the reduction in transmission will affect the evolution of drug resistance in the parasite. By manipulating host immunity, drug treatment regimes, or vector population over time, this model could be used to evaluate the consequences of malaria control measures. These findings emphasize the importance of considering the potential for immune selection in investigations of Plasmodium spp evolution.

Model and methods

Overview

The purpose of our model is to explore the role and relative importance of selection driven by host immune responses and by drug treatment on the evolution of malaria parasites. We seek to further understand the observed pattern of drug resistance, in which it rarely arises in high transmission settings, but spreading easily once introduced from elsewhere [8]. Our model is individual-based at the levels of parasites, vectors, and hosts. Parasites are nested within hosts or vectors, parasite genotypes are explicitly tracked, and the state of the within-host environment (e.g., immunity and red blood cells) is defined at the individual host level. Simulations are conducted over discrete time steps, equivalent to days. Detailed methods and parameters can be found in S1 Text and S1 and S2 Tables. Here, we describe the building blocks (i.e., agents) of the model, the processes that the model captures, and the simulations that form the basis of our analysis.

Individuals, state variables, and assumptions

Mosquitoes.

Individual mosquito vectors are assumed to live for 20 days. Upon death, they are immediately replaced by a new mosquito of age zero. There is no background mortality nor age-related mortality increase. Malaria infection is assumed to have no impact on mosquito survival. The mosquito population size remains constant over the course of a given simulation, representing a carrying capacity in a given environment.

In the absence of infection, mosquitoes feed every three days [69], on days at which their age is divisible by three, so that population feeding is asynchronous. If the mosquito becomes infected, a one day feeding delay is added to represent the decreased feeding behaviour that has been observed in pre-infectious mosquitoes infected with Plasmodium spp. [70–72]. We assume that mosquitoes ingest two microliters of blood from a single host at each feeding.

Hosts.

In our simulations, hosts live for a maximum of 500 time steps. While this is short, given our interpretation of time steps as days, we do this for two reasons. First, a high rate of transmission was required to prevent eradication when drug therapy was introduced to the population. As a result, endemicity in the low transmission conditions substantially exceeded the endemicity of the natural resistance hotspots in Southeast Asia and South America. At older host ages, just as in holoendemic regions, constant, asymptomatic infection became the rule under all conditions. Because the population size was limited by computational constraints, increasing the host lifespan caused the population to be saturated with infection. Excluding older hosts increased the difference between the high and low transmission conditions because the rate of exposure was slower in low transmission conditions. This allowed us to focus on the part of the lifespan where the populations diverge the most.

Hosts die from random background mortality at a fixed rate. Hosts can also die from malaria infection if their within-host count of uninfected red blood cells falls below a minimum threshold or if the percent of cells infected (parasitemia) exceeds a maximum threshold. When a host dies due to background mortality, infection, or reaching the maximum age, they are are immediately replaced by a new host of age zero. Thus, the host population size, N, remains constant over the course of the simulation.

We prioritized parameter choices that created realistic variation in infection outcomes because infection duration, immunity levels, and asymptomatic transmission all had potential to affect the within-host environment and the dynamics of selection. In human populations, the duration of infections is highly variable and difficult to measure, and both acute and chronic infections are common [73–77]. Variation in host response to infection was generated in two ways. First, each host had a randomly selected hyperparasitemia mortality threshold drawn from a Poisson distribution with mean of 11%, a value based on the WHO threshold for severe P. falciparum malaria [78]. Second, the rate of growth of the adaptive immune response varied between hosts. These factors produced a range of outcomes and courses of disease that is comparable to empirical observations.

Within a naive host, our default parameters give rise to peak parasitemia between day 12–14 post-infection. Approximately 25% of hosts die at this time during a primary infection and an additional 25% of infections are cleared (S1 Fig). Early clearance is highly time-sensitive and occurs when adaptive immunity increases quickly enough to produce a substantial effect while innate immunity is still at peak (S2 Fig). The rest of the infections last at least 150 days and are characterized by dampened oscillations produced by antigenic escape [77]. Mortality is rare during later stages of infection. Infections are asymptomatic during most of that time, though occasionally oscillations will cause transient symptoms. Asymptomatic infections are still transmissible. Subsequent infections are typically low density and truncated. Reinfections with a previously exposed strain usually remain subpatent and are unable to produce gametocytes (i.e., parasites transmissible to mosquitoes), but a novel strain may produce a symptomatic, transmissible infection (S3 Fig).

Parasites.

Individual parasites are tracked within hosts and within vectors. Each parasite has a unique genotype consisting of three unlinked loci, two of which encode resistance (and therefore determine growth rate). The allele at the third locus determines interactions with the host immune responses. We refer to all parasites that share an allele here as a “strain”. Strains are equivalent to serotypes, and all members of the same strain have identical epitopes and elicit the same immune response. The genome uses the standard nucleotides A, G, T, and C. Each base is subject to mutation at a given rate at each replication cycle. The genomic mutation rate—applying to the loci determining resistance and growth rate—is 2.5 × 10-5 per residue per replication. This value makes it likely that a small number of resistant mutants will be generated in each infection, especially during primary infections, in accord with numerical expectations [3]. We assume a finite number of possible strains. At each replication cycle, there is a probability of mutation to a different random strain, chosen with equal probability.

The strain mutation rate is 1 × 10-5 per replication, allowing variation in growth rate/resistance to increase within a single antigenic strain as a consequence of selection driven by competition and/or drug treatment. Antigenic diversity within the population is primarily due to immune selection. During peak infection in a naive host, 100,000–200,000 parasites are produced per time step. This will generate approximately 1–2 strain mutants per time step. Without immune selection, new strain mutants are almost immediately lost due to genetic drift (S14 Fig). At lower strain mutation rates, population strain diversity was lower, was more stable over time, and was similar between transmission intensities (S11 Fig).

A parasite genotype is drug resistant if either or both resistance-associated loci are mutated to a value of “C”. This approximates single-locus resistance, such as resistance to chloroquine. There is no epistasis between loci, so any genotype with at least one resistance mutation produces a fully resistant phenotype, effectively doubling the resistance mutation rate. The cost of resistance is modelled as an intrinsic growth rate reduction. In untreated naive hosts, an infection containing only resistant mutants will grow at 40% of the rate of an infection containing only sensitive parasites. In an untreated mixed infection, expansion of a resistant lineage will additionally be limited by immunity elicited against the faster-growing sensitive population.

Processes

Transmission to hosts and within-host infection dynamics.

An infection in a host begins when an infectious mosquito takes a blood meal. We assume this occurs with a probability of one. The number of motile sporozoites (parasites transmissible to hosts) injected into the host are drawn from a Poisson distribution with a mean of twelve, which ensures a blood-stage infection that leads to gradual build up of immunity within the host. In reality, these parasites travel to the liver and undergo several rounds of replication; in our model, sporozoites spend seven days dormant, during which time they are not recognized by the host immune system. At the end of this period, each sporozoite produces five merozoites—the red blood cell (RBC) infecting parasite stage.

We assume hosts have a maximum of 8.5 × 106 RBCs and produce up to 3.7 × 105 new RBCs per day. We choose these numbers for computational practicality and note that they essentially correspond to modeling dynamics in one μL of blood in a mouse. In an uninfected host, RBCs are subject to background mortality. In an infected host, uninfected RBC are subject to bystander killing, in which infected cells increase the mortality rate of uninfected cells [79–81]. Infected RBCs have a higher background mortality rate due to fragility [82] and clearance by the spleen [78].

The number of erythrocytes, a, infected in each host at each time point is calculated aswhere U represents the number of uninfected erythrocytes per microliter, wi represents the growth rate of an individual merozoite, of which there are m, and ϑ is a constant representing infection probability, which can be considered as the probability that a merozoite encounters and infects an RBC [36]. At each time step, a merozoites are chosen randomly with replacement to infect RBCs, but the probability that any individual merozoite will infect an RBC is weighted by its growth rate. The number of RBCs infected cannot exceed the summed growth rate of all merozoites. In other words, ten merozoites with a growth rate of 1 will infect no more than ten RBCs; ten merozoites with a growth rate of 0.5 will infect no more than five RBCs.

Once infecting an RBC, parasites take one of two developmental pathways. In 95% of infected cells, parasites replicate asexually, releasing ten progeny parasites one day after infection. In 5% of infected cells, parasites instead go on to produce gametocytes, the stage of the parasite that is transmissible to mosquitoes. While this value is constant and does not reflect predictions about parasite phenotypic plasticity [83, 84], it corresponds to overall values observed in mice [85]. We assume gametocytes take two days to mature; thus, the within-host life cycle roughly captures the biological details of P. chabaudi infection in mice. Again, these choices were for computational practicality. While altering the cell cycle duration of parasites or blood density of hosts may quantitatively alter the timing of evolutionary dynamics, we were concerned with the relative timing across different contexts.

Merozoites have an instantaneous background mortality drawn from Poisson distribution with a mean of 20%. Any merozoites which do not infect an erythrocyte are cleared in each time step. Gametocytes are removed via daily background mortality. We do not separately track different gametocyte sexes.

The host immune response is composed of both innate and adaptive responses. We model immunity phenomenologically in terms of killing rate. Both arms of immunity are dependent on the density of infected erythrocytes, and all blood-stage parasites are susceptible to their effects. Total immune activity against a given parasite strain is the sum of four processes: innate activation, adaptive activation, cross-reactive activation, and immune saturation (which reduces total activity). An overview of each is described below, with a full description in S1 Text.

Innate immunity is strictly density dependent, peaking early and rapidly declining. Its growth and its effects are strain-independent. Adaptive immunity is specific to each antigenic strain and the immune response to a given strain is not decreased by co-infection with additional strains. At each time step in an infected host, adaptive immunity grows as a function of an individual host’s adaptive growth rate and infected erythrocyte density, and is decreased by antigenic escape (see below). Once a host is no longer exposed to a strain, adaptive immunity to that strain decays slowly, even if the host is infected by other strains.

While adaptive immunity grows and decays independently for each strain, strains are cross-reactive, and total adaptive immune activity against a given strain is increased by adaptive immunity to other strains. We assume that each strain is identically cross-reactive to the others.

Antigenic escape is conceptually based on PfEMP1 proteins encoded by var genes [86, 87]. PfEMP1 variation is independent of serotype. The rate of antigenic escape shapes the duration and course of infections [88]. Early in an infection, there are many possible alternative configurations. Initially, adaptive immunity fluctuates with changes in density of infected erythrocytes. As the duration of exposure increases, the rate of escape slows as fewer novel variants remain and eventually halts once all possible configurations have been exhausted, resulting in sustained growth of adaptive immunity even with a low density of infected erythrocytes.

Exposure is measured independently for each strain, such that if strain two is introduced into a host with a long duration of exposure to strain one, the immune response against strain two will be bolstered by cross-reactivity against one, but the growth of the strain two-specific immune response will be slow initially because strain two retains full potential for antigenic escape.

Finally, the immune system’s capacity to clear parasites saturates at high parasite densities. There is a sigmoidal relationship between the efficacy of immune response and parasite density [89].

Transmission to vectors and within-vector infection dynamics.

Mosquitoes can become infected with malaria when feeding on an infected host. The probability of infection, p, is a sigmoidal function of within-host gametocyte density, given as L, per microlitre. Following [29], we used the relationship determined in [90] with shape values based on data from [91]. The probability of infection is thus given as,

If the vector is infected, it will ingest 2L gametocytes, assuming a 2 μL bloodmeal size. The number of gametocytes is assumed to not be a limiting factor, so gametocytes which are ingested by a vector are not removed from the host’s pool of gametocytes. From the ingested gametocytes, 1% are randomly chosen to sexually reproduce with another gametocyte from the same bloodmeal, with free recombination between all three loci. Vector immunity dynamics are not explicitly modelled, but the other 99% of the ingested parasites are assumed to die before becoming oocysts. Remarkably, this is less severe than the average 2,754-fold decrease observed in experimental infections [92]. Ten days after the bloodmeal [59], each oocyst bursts to produce 12 sporozoites, each of which undergoes mutation independently with the same parameters as within-host mutation. Once sporozoites are produced, the mosquito is infectious.

Drug treatment.

By default, 30% of symptomatic hosts within a population will be treated. Pyrogenic threshold–the parasite density at which a patient will exhibit fever as a symptom–varies widely between individuals and may be affected by host age and regional transmission intensity [93, 94]. Here, hosts with merozoite density exceeding a threshold of 5000/μL are considered to be symptomatic and thus become eligible for treatment. Symptomatic hosts are selected randomly to be treated for three days with a drug that kills 99.9% of sensitive parasites per day, which is sufficient to clear sensitive infections and is in line with expectations for artemisinin combination therapies [95]. Treatment persists even if the host recovers before the end of the three days. Drug concentration is binary with no half-life.

Simulations

Initialization and population burn-ins.

The host population size is set at 100. Populations are initiated without vectors, mutation, or host mortality by infecting one host every five time steps for 250 time steps, producing a base population with a wide range of adaptive immunity and infection age. Parasite genomes begin with all bases set to A. After 250 time steps, host mortality is restored, with uniformly-distributed ages between 0 and 500, and mosquitoes are introduced, with uniformly-distributed ages between zero and 20. Given the other model assumptions, the mosquito age-structure will be retained throughout the simulations. The distribution of host ages, however, will shift over time due to background and infection-induced mortality. Populations are allowed to run for 2000 time steps. By the end of that time, measures of endemicity and strain diversity have reached equilibrium. At that time, drug treatment is introduced.

Emergent properties.

The size of the mosquito population, set at either 300 or 1200 individual vectors, is the source of top-down regulation of transmission intensity. However, the realized transmission intensity is the product of interactions between treatment, frequency of resistance, immunity and genetic diversity, each of which is, in turn, affected by transmission intensity as well as interactions with one another. As a result, complex patterns can emerge from a few rules, and the nature of the interactions is difficult to intuit.

Resistance evolution dynamics.

Because we model resistant parasites at both the within-host and the host population levels, we track the evolution of resistance by measuring prevalence of infections in the host population which are composed of at least 10% resistant parasites, a frequency which was sufficient to cause treatment failure in our simulations. Therefore, time to widespread treatment failure, given as Tfail, is defined as the time point at which 10% of infections are composed of at least 10% resistant parasites. Time to ubiquity, at which treatment failure predominates in a host population, is defined as the time point at which 75% of infections are composed of at least 10% resistance parasites. Tubiq was only measured for conditions in which the average equilibrium prevalence of treatment failure exceeded 75%. Simulations are run until they reach an equilibrium resistance prevalence.

Statistical analysis.

Analysis was conducted using the R statistical package [96]. Treatment outcomes were compared using the Wilcoxon Rank Sum test, and normal approximated p-values are reported. For conditions in which the p-value was significant at the 95% confidence interval, treatment effect size was quantified using Vargha and Delaney’s A [97], implemented in the R library effsize.

Supporting information

S1 Text. Supplemental methods and data.

https://doi.org/10.1371/journal.pcbi.1008577.s001

(PDF)

S1 Table. Table of symbols used.

https://doi.org/10.1371/journal.pcbi.1008577.s002

(PDF)

S2 Table. List of parameters and default values.

https://doi.org/10.1371/journal.pcbi.1008577.s003

(PDF)

S1 Fig. Duration of untreated primary infections.

Infection clearance time was measured for 500 infections in untreated hosts. Approximately 50% of infections become chronic. Approximately 25% of infections have a rapid clearance time, represented in the early peak. Infections which result in mortality are not represented on this graph, but account for approximately 25% of primary infections.

https://doi.org/10.1371/journal.pcbi.1008577.s004

(TIF)

S2 Fig. Parasite and immune dynamics during primary infections.

Two naive hosts are infected at time 0 with a 50 merozoites of a single strain. No mutation, treatment, or reinfection occurs. (A) Merozoite density/μL over time. Host 1 (teal) clears the infection rapidly. Host 2 (purple) shows the more common course of infection, under our default parameters. (B) Innate immunity over time. (C) Adaptive immunity over time.

https://doi.org/10.1371/journal.pcbi.1008577.s005

(TIF)

S3 Fig. Infection trajectory during secondary exposure.

The hosts from S2 Fig were inoculated with 50 merozoites at time 300. (A) Reinfection with the same strain. (B) Infection with a novel strain.

https://doi.org/10.1371/journal.pcbi.1008577.s006

(TIF)

S4 Fig. Strains exposed by age.

High transmission increased both the rate of exposure to new strains as well as the total number of strains to which hosts are ever exposed. Age is given in days. Values were calculated from 10 replicate untreated equilibrium populations. Lines represent means and shaded regions are standard deviations.

https://doi.org/10.1371/journal.pcbi.1008577.s007

(TIF)

S5 Fig. Complexity of infection by age.

High transmission produced more complex infections at every age, reaching a peak at approximately 100 time steps. Even though older hosts had been exposed to more strains S4 Fig), complexity of infection did not increase after this peak, and in fact, slightly decreased with high transmission. Values were calculated from 10 untreated replicate equilibrium populations. Lines represent means and shaded regions are standard deviations.

https://doi.org/10.1371/journal.pcbi.1008577.s008

(TIF)

S6 Fig. The equilibrium prevalence of resistance by treatment rate.

The equilibrium prevalence of treatment failure was measured in 20 replicate simulations, excluding replicates in which malaria was eradicated. Treatment rate for each simulation is indicated by the title of the graph.

https://doi.org/10.1371/journal.pcbi.1008577.s009

(TIF)

S7 Fig. Time to resistance with a 70% treatment rate.

Mean time to resistance was measured over 20 replicate populations, excluding replicates in which malaria was eradicated. A: 10% prevalence of treatment failure (Tfail). B: 75% prevalence of treatment failure (Tubiq). Note different Y axes.

https://doi.org/10.1371/journal.pcbi.1008577.s010

(TIF)

S8 Fig. Effective treatment rate.

Effective treatment rates were measured for 20 replicate equilibrium populations over 250 time steps with means and standard deviations shown. Top row: The proportion of all infections that were treated. A: One strain, B: 30 strains. Bottom row: The proportion of total parasite population gametocytes that are in a treated host. C: One strain. D: 30 strains.

https://doi.org/10.1371/journal.pcbi.1008577.s011

(TIF)

S9 Fig. The effect of recombination on time to resistance.

In populations with 30 strains, recombination between the strain locus and the resistance loci did not have a significant impact on evolution of resistance. Each row represents a different strain mutation rate. Top row: Higher strain mutation rate, equal to genomic mutation rate (2.5 × 10-5) (A) Tubiq with recombination. (B) Tubiq with no recombination. Second row: default strain mutation rate (1 × 10-5) (C) Tubiq with recombination (repeated from Fig 3B in the main text). (D) Tubiq with no recombination. Third row: Reduced strain mutation rate 5 × 10-6 (E) Tubiq with recombination. (F) Tubiq with no recombination. Bottom row: Reduced strain mutation rate 1 × 10-6 (G) Tubiq with recombination. (H) Tubiq with no recombination. Tubiq.

https://doi.org/10.1371/journal.pcbi.1008577.s012

(TIF)

S10 Fig. Ranked frequency of gametocyte strains within a blood meal.

Representative blood meals were drawn from hosts in an equilibrium population and the frequencies of the strains it contained were ranked. With 300 vectors (left), most blood meals contained only a single strain, with minimal contribution from other strains. Diversity was higher with 1200 vectors (right), but the majority of the gametocytes within a single blood meal were still from a single strain.

https://doi.org/10.1371/journal.pcbi.1008577.s013

(TIF)

S11 Fig. Population-wide mean strain-specific immunity with a lower strain mutation rate.

Populations were initiated with maximum strain diversity and then allowed to evolve to equilibrium over 2000 timesteps using a strain mutation rate of 1 × 10-6, an order of magnitude lower than the default. As in Fig 2 in the main text, immunity to each strain was averaged over all hosts in a single representative population. Strain diversity was lower with lower strain mutation rate and was similar between 300 vectors (left) and 1200 vectors (right). This demonstrates that lower strain mutation rates do not generate enough antigenic diversity to cause broad population immunity or induce immune selection.

https://doi.org/10.1371/journal.pcbi.1008577.s014

(TIF)

S12 Fig. Evolutionary trajectory of costless resistance.

Following introduction of treatment at time 0, the prevalence of within-host frequencies of resistance was monitored over time. The top (purple) line indicates the prevalence of all treatment-resistant infections, defined as infections consisting of at least 10% resistant parasites. Moving down from the top, the prevalence of infections consisting of between 10% and 50% resistant parasites is shown in blue. Below that, the teal section shows the prevalence of infections in which the frequency of resistant parasites is between 50% and 90%. Finally, the bottom section in green shows the prevalence of infections in which resistance is at near fixation within the host, with a frequency exceeding 99.9%. Top row (A and B): 300 vectors. Bottom row (C and D): 1200 vectors. Left column (A and C): one strain. Right column (B and D): 30 strains. Patterns are qualitatively similar to costly resistance, in which mixed (i.e., sensitive and resistant) infections are rare within one strain conditions. Both 30 strain conditions show mixed infections at equilibrium, but, just as in costly resistance, mixed infections are common pre-ubiquity with 1200 vectors, indicating the greatest role for immune competition there. Data shows one representative simulation for each condition.

https://doi.org/10.1371/journal.pcbi.1008577.s015

(TIF)

S13 Fig. Frequency of resistance within each strain.

The frequency of resistance in each strain was monitored from treatment introduction to equilibrium. Left: 300 vectors. Most strains are fully resistant at equilibrium. Right: 1200 vectors. At equilibrium, partially resistant strains are common, suggesting influence from within-host competition. Representative populations are the same as those in Fig 6 of the main text. Gray indicates the strain is not present. The blue line indicates Tfail and the red line indicates Tubiq.

https://doi.org/10.1371/journal.pcbi.1008577.s016

(TIF)

S14 Fig. Evolution of antigenic diversity without immune selection.

Every individual in a population of naive hosts was inoculated with parasites from a single strain. Strain mutation was permitted at the default rate, but cross reactivity between strains was 100%, so there was no selection for antigenic novelty. Strain mutants were produced during peak infection, but they remained at very low frequencies (note scale of the legend) and were quickly lost due to genetic drift. This demonstrates that strain diversity is almost entirely the result of selection.

https://doi.org/10.1371/journal.pcbi.1008577.s017

(TIF)

S15 Fig. Population merozoite density.

In order to determine if differences in mutation supply affected the origin of resistance, the density of all merozoites in ten replicate host populations was monitored for the first ten time steps after treatment was introduced. Despite the differences in population size, Tfail was similar between conditions (shown in Fig 3A in the main text), indicating that mutation supply is not a limiting factor.

https://doi.org/10.1371/journal.pcbi.1008577.s018

(TIF)

Acknowledgments

We thank Madeline Peters and Megan Greischar for their insightful comments during the manuscript preparation.

Citation: Whitlock AOB, Juliano JJ, Mideo N (2021) Immune selection suppresses the emergence of drug resistance in malaria parasites but facilitates its spread. PLoS Comput Biol 17(7): e1008577. https://doi.org/10.1371/journal.pcbi.1008577

References

1. World Malaria Report. World Health Organization. 2019.

2. Paget-Mcnicol S, Saul A. Mutation rates in the dihydrofolate reductase gene of Plasmodium falciparum. Parasitology. 2001;122(5):497–505. pmid:11393822

3. White NJ, Pongtavornpinyo W. The de novo selection of drug-resistant malaria parasites. Proc R Soc B Biol Sci. 2003;270(1514):545–554.

4. Roper C, Pearce R, Nair S, Sharp B, Nosten F, Anderson T. Intercontinental spread of pyrimethamine-resistant malaria. Science. 2004;305(5687):1124. pmid:15326348

5. Cortese JF, Caraballo A, Contreras CE, Plowe CV. Origin and dissemination of Plasmodium falciparum drug-resistance mutations in South America. J Infect Dis. 2002;186(7):999–1006. pmid:12232841

6. Nair S, Williams JT, Brockman A, Paiphun L, Mayxay M, Newton PN, et al. A selective sweep driven by pyrimethamine treatment in Southeast Asian malaria parasites. Mol Biol Evol. 2003;20(9):1526–1536. pmid:12832643

7. Roper C, Pearce R, Bredenkamp B, Gumede J, Drakeley C, Mosha F, et al. Antifolate antimalarial resistance in southeast Africa: A population-based analysis. Lancet. 2003;361(9364):1174–1181. pmid:12686039

8. Wootton JC, Feng X, Ferdig MT, Cooper RA, Mu J, Baruch DI, et al. Genetic diversity and chloroquine selective sweeps in Plasmodium falciparum. Nature. 2002;418(6895):320–323. pmid:12124623

9. Imwong M, Suwannasin K, Kunasol C, Sutawong K, Mayxay M, Rekol H, et al. The spread of artemisinin-resistant Plasmodium falciparum in the Greater Mekong subregion: a molecular epidemiology observational study. Lancet Infect Dis. 2017;17(5):491–497. pmid:28161569

10. Chen N, Kyle DE, Pasay C, Fowler EV, Baker J, Peters JM, et al. pfcrt Allelic types with two novel amino acid mutations in chloroquine-resistant Plasmodium falciparum isolates from the Philippines. Antimicrob Agents Chemother. 2003;47(11):3500–3505. pmid:14576108

11. Trape JF. The public health impact of chloroquine resistance in Africa. In: Am. J. Trop. Med. Hyg. vol. 64. American Society of Tropical Medicine and Hygiene; 2001. p. 12–17. pmid:11425173

12. Hamilton WL, Amato R, van der Pluijm RW, Jacob CG, Quang HH, Thuy-Nhien NT, et al. Evolution and expansion of multidrug-resistant malaria in southeast Asia: a genomic epidemiology study. Lancet Infect Dis. 2019;19(9):943–951. pmid:31345709

13. van der Pluijm RW, Imwong M, Chau NH, Hoa NT, Thuy-Nhien NT, Thanh NV, et al. Determinants of dihydroartemisinin-piperaquine treatment failure in Plasmodium falciparum malaria in Cambodia, Thailand, and Vietnam: a prospective clinical, pharmacological, and genetic study. Lancet Infect Dis. 2019;19(9):952–961. pmid:31345710

14. Babiker HA, Creasey AM, Fenton B, Bayoumi RAL, Arnot DE, Walliker D. Genetic diversity of Plasmodium falciparum in a village in eastern Sudan. 1. Diversity of enzymes, 2D-PAGE proteins and antigens. Trans R Soc Trop Med Hyg. 1991;85(5):572–577. pmid:1780978

15. Paul REL, Hackford I, Brockman A, Muller-Graf C, Price R, Luxemburger C, et al. Transmission intensity and Plasmodium falciparum diversity on the northwestern border of Thailand. Am J Trop Med Hyg. 1998;58(2):195–203. pmid:9502604

16. Tanner M, Beck HP, Felger I, Smith T. The epidemiology of multiple Plasmodium falciparum infections. Trans R Soc Trop Med Hyg. 1999;93(SUPPL. 1):1–2. pmid:10450418

17. A-Elbasit IE, ElGhazali G, A-Elgadir TME, Hamad AA, Babiker HA, Elbashir MI, et al. Allelic polymorphism of MSP2 gene in severe P. falciparum malaria in an area of low and seasonal transmission. Parasitol Res. 2007;102(1):29–34. pmid:17768641

18. Babiker HA, Lines J, Hill WG, Walliker D. Population structure of Plasmodium falciparum in villages with different malaria endemicity in east Africa. Am J Trop Med Hyg. 1997;56(2):141–147. pmid:9080871

19. Arnot D. Clone multiplicity of Plasmodium falciparum infections in individuals exposed to variable levels of disease transmission. Trans R Soc Trop Med Hyg. 1998;92(6):580–585. pmid:10326095

20. Juliano JJ, Porter K, Mwapasa V, Sem R, Rogers WO, Ariey F, et al. Exposing malaria in-host diversity and estimating population diversity by capture-recapture using massively parallel pyrosequencing. Proc Natl Acad Sci U S A. 2010;107(46):20138–20143. pmid:21041629

21. Mideo N, Bailey JA, Hathaway NJ, Ngasala B, Saunders DL, Lon C, et al. A deep sequencing tool for partitioning clearance rates following antimalarial treatment in polyclonal infections. Evol Med Public Heal. 2016;2016(1):21–36. pmid:26817485

22. Wargo AR, Huijben S, De Roode JC, Shepherd J, Read AF. Competitive release and facilitation of drug-resistant parasites after therapeutic chemotherapy in a rodent malaria model. Proc Natl Acad Sci U S A. 2007;104(50):19914–19919. pmid:18056635

23. de Roode JC, Helinski MEH, Anwar MA, Read AF. Dynamics of multiple infection and within-host competition in genetically diverse malaria infections. Am Nat. 2005;166(5):531–542. pmid:16224719

24. Bushman M, Morton L, Duah N, Quashie N, Abuaku B, Koram KA, et al. Within-host competition and drug resistance in the human malaria parasite Plasmodium falciparum. Proc R Soc B Biol Sci. 2016;283(1826):20153038. pmid:26984625

25. Doolan DL, Dobaño C, Baird JK. Acquired immunity to malaria. Clin Microbiol Rev. 2009;22(1):13–36. pmid:19136431

26. Hay SI, Rogers DJ, Toomer JF, Snow RW. Annual Plasmodium falciparum entomological inoculation rates (EIR) across Africa: Literature survey, internet access and review. Trans R Soc Trop Med Hyg. 2000;94(2):113–127. pmid:10897348

27. Ataide R, Ashley EA, Powell R, Chan JA, Malloy MJ, O’Flaherty K, et al. Host immunity to Plasmodium falciparum and the assessment of emerging artemisinin resistance in a multinational cohort. Proc Natl Acad Sci U S A. 2017;114(13):3515–3520. pmid:28289193

28. Babiker HA, Hastings IM, Swedberg G. Impaired fitness of drug-resistant malaria parasites: Evidence and implication on drug-deployment policies. Expert Rev Anti Infect Ther. 2009;7(5):581–593. pmid:19485798

29. Huijben S, Nelson WA, Wargo AR, Sim DG, Drew DR, Read AF. Chemotherapy, within-host ecology and the fitness of drug-resistant malaria parasites. Evolution. 2010;64(10):2952–68. pmid:20584075

30. McCollum AM, Schneider KA, Griffing SM, Zhou Z, Kariuki S, Ter-Kuile F, et al. Differences in selective pressure on dhps and dhfr drug resistant mutations in western Kenya. Malar J. 2012;11:77. pmid:22439637

31. Nair S, Li X, Arya GA, McDew-White M, Ferrari M, Nosten F, et al. Fitness costs and the rapid spread of kelch13-C580Y substitutions conferring artemisinin resistance. Antimicrob Agents Chemother. 2018;62(9). pmid:29914963

32. Walliker D, Hunt P, Babiker H. Fitness of drug-resistant malaria parasites. Acta Trop. 2005;94(3 SPEC. ISS.):251–259. pmid:15845348

33. Nair S, Miller B, Barends M, Jaidee A, Patel J, Mayxay M, et al. Adaptive Copy Number Evolution in Malaria Parasites. PLOS Genet. 2008;4(10):e1000243. pmid:18974876

34. Straimer J, Gnädig NF, Stokes BH, Ehrenberger M, Crane AA, Fidock DA. Plasmodium falciparum K13 mutations differentially impact ozonide susceptibility and parasite fitness in vitro. MBio. 2017;8(2).

35. Mackinnon MJ, Marsh K. The selection landscape of malaria parasites. Science. 2010;328(5980):866–71. pmid:20466925

36. Bushman M, Antia R, Udhayakumar V, de Roode JC. Within-host competition can delay evolution of drug resistance in malaria. PLOS Biol. 2018;16(8):e2005712. pmid:30130363

37. Huijben S, Chan BHK, Nelson WA, Read AF. The impact of within-host ecology on the fitness of a drug-resistant parasite. Evol Med Public Heal. 2018;2018(1):127. pmid:30087774

38. Scott N, Ataide R, Wilson DP, Hellard M, Price RN, Simpson JA, et al. Implications of population-level immunity for the emergence of artemisinin-resistant malaria: A mathematical model. Malar J. 2018;17(1):279. pmid:30071877

39. Early AM, Lievens M, MacInnis BL, Ockenhouse CF, Volkman SK, Adjei S, et al. Host-mediated selection impacts the diversity of Plasmodium falciparum antigens within infections. Nat Commun. 2018;9(1):1381. pmid:29643376

40. Klein EY, Smith DL, Boni MF, Laxminarayan R. Clinically immune hosts as a refuge for drug-sensitive malaria parasites. Malar J. 2008;7:67. pmid:18439283

41. Antao T, Hastings IM. Environmental, pharmacological and genetic influences on the spread of drug-resistant malaria. Proc Biol Sci. 2011;278(1712):1705–1712. pmid:21084349

42. Brock AR, Ross JV, Greenhalgh S, Durham DP, Galvani A, Parikh S, et al. Modelling the impact of antimalarial quality on the transmission of sulfadoxine-pyrimethamine resistance in Plasmodium falciparum. Infect Dis Model. 2017;2(2):161–187. pmid:29928735

43. Hastings IM. A model for the origins and spread of drug-resistant malaria. Parasitology. 1997;115(2):133–141. pmid:10190169

44. Schneider KA, Kim Y. An analytical model for genetic hitchhiking in the evolution of antimalarial drug resistance. Theor Popul Biol. 2010;78(2):93–108. pmid:20600206

45. Kim Y, Escalante AA, Schneider KA. A population genetic model for the initial spread of partially resistant malaria parasites under anti-malarial combination therapy and weak intrahost competition. PLOS One. 2014;9(7):e101601. pmid:25007207

46. Dang Nguyen T, Olliaro P, Dondorp AM, Baird K, Lam HM, Farrar J, et al. Optimum population-level use of artemisinin combination therapies: a modelling study. Lancet Glob Heal. 2015;3:e758–e766.

47. Legros M, Bonhoeffer S. A combined within-host and between-hosts modelling framework for the evolution of resistance to antimalarial drugs. J R Soc Interface. 2016;13(117):20160148. pmid:27075004

48. Chang HH, Childs LM, Buckee CO. Variation in infection length and superinfection enhance selection efficiency in the human malaria parasite. Sci Rep. 2016;6(1):26370. pmid:27193195

49. Ayala MJC, Villela DAM. Early transmission of sensitive strain slows down emergence of drug resistance in Plasmodium vivax. PLOS Comput Biol. 2020;16(6). pmid:32555701

50. Chang HH, Moss EL, Park DJ, Ndiaye D, Mboup S, Volkman SK, et al. Malaria life cycle intensifies both natural selection and random genetic drift. Proc Natl Acad Sci U S A. 2013;110(50):20129–34. pmid:24259712

51. Chang HH, Hartl DL. Recurrent bottlenecks in the malaria life cycle obscure signals of positive selection. Parasitology. 2015;142 Suppl(Suppl 1):S98–S107. pmid:24560397

52. Ochola LI, Tetteh KKA, Stewart LB, Riitho V, Marsh K, Conway DJ. Allele frequency-based and polymorphism-versus-divergence indices of balancing selection in a new filtered set of polymorphic genes in Plasmodium falciparum. Mol Biol Evol. 2010;27(10):2344–2351. pmid:20457586

53. Amambua-Ngwa A, Tetteh KKA, Manske M, Gomez-Escobar N, Stewart LB, Deerhake ME, et al. Population genomic scan for candidate signatures of balancing selection to guide antigen characterization in malaria parasites. PLOS Genet. 2012;8(11):e1002992. pmid:23133397

54. Mobegi VA, Duffy CW, Amambua-Ngwa A, Loua KM, Laman E, Nwakanma DC, et al. Genome-wide analysis of selection on the malaria parasite Plasmodium falciparum in West African populations of differing infection endemicity. Mol Biol Evol. 2014;31(6):1490–1499. pmid:24644299

55. Kassegne K, Komi Koukoura K, Shen HM, Chen SB, Fu HT, Chen YQ, et al. Genome-Wide Analysis of the Malaria Parasite Plasmodium falciparum Isolates From Togo Reveals Selective Signals in Immune Selection-Related Antigen Genes. Front Immunol. 2020;11. pmid:33193320

56. Jalloh A, Jalloh M, Matsuoka H. T-cell epitope polymorphisms of the Plasmodium falciparum circumsporozoite protein among field isolates from Sierra Leone: Age-dependent haplotype distribution? Malar J. 2009;8(1):120. pmid:19500348

57. Miotto O, Amato R, Ashley EA, Macinnis B, Almagro-Garcia J, Amaratunga C, et al. Genetic architecture of artemisinin-resistant Plasmodium falciparum. Nat Genet. 2015;47(3):226–234. pmid:25599401

58. Iles MM, Walters K, Cannings C. Recombination can evolve in large finite populations given selection on sufficient loci. Genetics. 2003;165(4):2249–2258. pmid:14704200

59. Aly ASI, Vaughan AM, Kappe SHI. Malaria parasite development in the mosquito and infection of the mammalian host. Annu Rev Microbiol. 2009;63:195–221. pmid:19575563

60. Li X, Kumar S, McDew-White M, Haile M, Cheeseman IH, Emrich S, et al. Genetic mapping of fitness determinants across the malaria parasite Plasmodium falciparum life cycle. PLOS Genet. 2019;15(10). pmid:31609965

61. Moser KA, Madebe RA, Aydemir O, Chiduo MG, Mandara CI, Rumisha SF, et al. Describing the current status of Plasmodium falciparum population structure and drug resistance within mainland Tanzania using molecular inversion probes. Mol Ecol. 2020;Accepted Author Manuscript(n/a). pmid:33107096

62. Uwimana A, Legrand E, Stokes BH, Ndikumana JLM, Warsame M, Umulisa N, et al. Emergence and clonal expansion of in vitro artemisinin-resistant Plasmodium falciparum kelch13 R561H mutant parasites in Rwanda. Nat Med. 2020;26(10):1602–1608. pmid:32747827

63. Balikagala B, Sakurai-Yatsushiro M, Tachibana SI, Ikeda M, Yamauchi M, Katuro OT, et al. Recovery and stable persistence of chloroquine sensitivity in Plasmodium falciparum parasites after its discontinued use in Northern Uganda. Malar J. 2020;19(1):76. pmid:32070358

64. Boni MF, Smith DL, Laxminarayan R. Benefits of using multiple first-line therapies against malaria. Proc Natl Acad Sci U S A. 2008;105(37):14216–14221. pmid:18780786

65. Casares S, Richie TL. Immune evasion by malaria parasites: a challenge for vaccine development. Curr Opin Immunol. 2009;21(3):321–330. pmid:19493666

66. Kengne-Ouafo JA, Sutherland CJ, Binka FN, Awandare GA, Urban BC, Dinko B. Immune responses to the sexual stages of Plasmodium falciparum parasites. Front Immunol. 2019;10:136. pmid:30804940

67. Efficacy and safety of RTS,S/AS01 malaria vaccine with or without a booster dose in infants and children in Africa: Final results of a phase 3, individually randomised, controlled trial. Lancet. 2015;386(9988):31–45. pmid:25913272

68. Olotu A, Fegan G, Wambua J, Nyangweso G, Leach A, Lievens M, et al. Seven-Year Efficacy of RTS,S/AS01 Malaria Vaccine among Young African Children. N Engl J Med. 2016;374(26):2519–2529. pmid:27355532

69. Killeen GF, McKenzie FE, Foy BD, Schieffelin C, Billingsley PF, Beier JC. A simplified model for predicting malaria entomologic inoculation rates based on entomologic and parasitologic parameters relevant to control. Am J Trop Med Hyg. 2000;62(5):535–544. pmid:11289661

70. Anderson RA, Koella JC, Hurd H. The effect of Plasmodium yoelii nigeriensis infection on the feeding persistence of Anopheles stephensi Liston throughout the sporogonic cycle. Proc R Soc B Biol Sci. 1999;266(1430):1729–1733. pmid:10518321

71. Koella JC, Rieu L, Paul REL. Stage-specific manipulation of a mosquito’s host-seeking behavior by the malaria parasite Plasmodium gallinaceum. Behav Ecol. 2002;13(6):816–820.

72. Cator LJ, George J, Blanford S, Murdock CC, Baker TC, Read AF, et al. ‘Manipulation’ without the parasite: Altered feeding behaviour of mosquitoes is not dependent on infection with malaria parasites. Proc R Soc B Biol Sci. 2013;280 (1763).

73. Jeffery GM, Young MD, Burgess RW, Eyles DE. Early activity in sporozoite-induced Plasmodium falciparum infections. Ann Trop Med Parasitol. 1959;53(1):51–58. pmid:13650489

74. Sama W, Dietz K, Smith T. Distribution of survival times of deliberate Plasmodium falciparum infections in tertiary syphilis patients. Trans R Soc Trop Med Hyg. 2006;100(9):811–816. pmid:16451806

75. Bretscher MT, Maire N, Chitnis N, Felger I, Owusu-Agyei S, Smith T. The distribution of Plasmodium falciparum infection durations. Epidemics. 2011;3(2):109–118. pmid:21624782

76. Ashley EA, White NJ. The duration of Plasmodium falciparum infections. Malar J. 2014;13(1):500. pmid:25515943

77. Childs LM, Buckee CO. Dissecting the determinants of malaria chronicity: why within-host models struggle to reproduce infection dynamics. J R Soc Interface. 2015;12(104):20141379–20141379. pmid:25673299

78. World Health Organization. Severe malaria. Trop Med Int Heal. 2014;19 Suppl 1:7–131.

79. Jakeman GN, Saul A, Collins WE. Anaemia of acute malaria infections in non-immune patients primarily results from destruction of uninfected erythrocytes. Parasitology. 1999;119:127–133. pmid:10466119

80. Metcalf CJE, Long GH, Mideo N, Forester JD, Bjørnstad ON, Graham AL. Revealing mechanisms underlying variation in malaria virulence: effective propagation and host control of uninfected red blood cell supply. J R Soc Interface. 2012;9(76):2804. pmid:22718989

81. Fonseca LL, Alezi HS, Moreno A, Barnwell JW, Galinski MR, Voit EO. Quantifying the removal of red blood cells in Macaca mulatta during a Plasmodium coatneyi infection. Malar J. 2016;15(1):410. pmid:27520455

82. Sherman IW, Eda S, Winograd E. Erythrocyte aging and malaria. Cell Mol Biol (Noisy-le-grand). 2004;50(2):159–169. pmid:15095786

83. Pollitt LC, Mideo N, Drew DR, Schneider P, Colegrave N, Reece SE. Competition and the evolution of reproductive restraint in malaria parasites. Am Nat. 2011;177(3):358–367. pmid:21460544

84. Greischar MA, Mideo N, Read AF, Bjørnstad ON. Quantifying transmission investment in malaria parasites. PLOS Comput Biol. 2016;12(2):e1004718. pmid:26890485

85. Wargo AR, de Roode JC, Huijben S, Drew DR, Read AF. Transmission stage investment of malaria parasites in response to in-host competition. Proc R Soc B Biol Sci. 2007;274(1625):2629–2638. pmid:17711832

86. Miller LH, Good MF, Milon G. Malaria pathogenesis. Science. 1994;264(5167):1878–1883. pmid:8009217

87. Kyes S, Horrocks P, Newbold C. Antigenic variation at the infected red cell surface in malaria. Annu Rev Microbiol. 2001;55:673–707. pmid:11544371

88. Klein EY, Graham AL, Llinás M, Levin S. Cross-reactive immune responses as primary drivers of malaria chronicity. Infect Immun. 2014;82(1):140–51. pmid:24126530

89. Metcalf CJE, Graham AL, Huijben S, Barclay VC, Long GH, Grenfell BT, et al. Partitioning regulatory mechanisms of within-host malaria dynamics using the effective propagation number. Science. 2011;333(6045):984–8. pmid:21852493

90. Sinden RE, Dawes EJ, Alavi Y, Waldock J, Finney O, Mendoza J, et al. Progression of Plasmodium berghei through Anopheles stephensi is density-dependent. PLOS Pathog. 2007;3(12):e195. pmid:18166078

91. Barnes KI, White NJ. Population biology and antimalarial resistance: The transmission of antimalarial drug resistance in Plasmodium falciparum. Acta Trop. 2005;94(3 SPEC. ISS.):230–240. pmid:15878154

92. Vaughan JA, Noden BH, Beier JC. Population dynamics of Plasmodium falciparum sporogony in laboratory-infected Anopheles gambiae. J Parasitol. 1992;78(4):716. pmid:1635032

93. Mmbando BP, Lusingu JP, Vestergaard LS, Lemnge MM, Theander TG, Scheike TH. Parasite threshold associated with clinical malaria in areas of different transmission intensities in north eastern Tanzania. BMC Med Res Methodol. 2009;9(1):75. pmid:19909523

94. Dollat M, Talla C, Sokhna C, Sarr FD, Trape JF, Richard V. Measuring malaria morbidity in an area of seasonal transmission: Pyrogenic parasitemia thresholds based on a 20-year follow-up study. PLOS One. 2019;14(6). pmid:31246965

95. White NJ. Assessment of the pharmacodynamic properties of antimalarial drugs in vivo. Antimicrob Agents Chemother. 1997;41(7):1413–1422. pmid:9210658

96. Ihaka R, Gentleman R. R: a language for data analysis and graphics. J Comp Graph Stat. 1996;5(3):299–314.

97. Vargha A, Delaney HD. A Critique and Improvement of the CL Common Language Effect Size Statistics of McGraw and Wong. J Educ Behav Stat. 2000;25(2):101–132.