Within-host competition can delay evolution of drug resistance in malaria

In the malaria parasite P. falciparum, drug resistance generally evolves first in low-transmission settings, such as Southeast Asia and South America. Resistance takes noticeably longer to appear in the high-transmission settings of sub-Saharan Africa, although it may spread rapidly thereafter. Here, we test the hypothesis that competitive suppression of drug-resistant parasites by drug-sensitive parasites may inhibit evolution of resistance in high-transmission settings, where mixed-strain infections are common. We employ a cross-scale model, which simulates within-host (infection) dynamics and between-host (transmission) dynamics of sensitive and resistant parasites for a population of humans and mosquitoes. Using this model, we examine the effects of transmission intensity, selection pressure, fitness costs of resistance, and cross-reactivity between strains on the establishment and spread of resistant parasites. We find that resistant parasites, introduced into the population at a low frequency, are more likely to go extinct in high-transmission settings, where drug-sensitive competitors and high levels of acquired immunity reduce the absolute fitness of the resistant parasites. Under strong selection from antimalarial drug use, however, resistance spreads faster in high-transmission settings than low-transmission ones. These contrasting results highlight the distinction between establishment and spread of resistance and suggest that the former but not the latter may be inhibited in high-transmission settings. Our results suggest that within-host competition is a key factor shaping the evolution of drug resistance in P. falciparum.

In the following equations and explanations, we use subscripts and to denote type-specific variables and parameters; thus ( , ) = (1,2) or (2,1). The differential equation for uninfected red blood cells is as follows: where is the rate of production of new RBCs, is the death rate of uninfected RBCs, and is the rate of infection of RBCs by free merozoites (parameter values can be found in Table 1).
Infected RBCs of type are described by the following equation: where is the background death rate of infected RBCs; in the absence of antimalarial drug treatment, = 0 (making the death rate ). If the host is being treated with antimalarial drugs, then = where represents the efficacy of drug treatment against type . is the per capita rate of gametocyte formation. and are the rates of killing by adaptive and innate immune responses, respectively.
is the proportion of (the adaptive immune response to type ) that is effective against type . The relationship of to the antigenic overlap between different strains is discussed later on. RBC  The equation for free merozoites of type is: where is the burst size (number of merozoites released by a single infected red blood cell), is the fitness cost of type (implemented as a reduction in burst size [1]), and is the death rate of free merozoites. , , and are as described above.
The equation for gametocytes of type is below: where is the death rate of mature gametocytes and and are as described above. Due to the scarcity of gametocytes in the human host, we assume that adaptive immune responses to gametocytes are negligible. Therefore, "natural" death and killing by innate immunity are the only mechanisms by which gametocytes are eliminated in the model. Innate immunity is described by the following equation: where is considered the fraction of a fixed pool of innate immune effectors that are currently "activated." is the activation rate of these effectors, and is the inactivation rate.
The dynamics of adaptive immunity to type are described by the following equations:

= +
The parameter is the maximum growth rate of the adaptive immune response and is the density of merozoites ( + ) at which the growth rate is 2 ⁄ (Figure 1). is the proportion of fixed (non-variant) antigens or epitopes that are shared between strains of types and . The contribution of type merozoites to stimulation of is proportional to this overlap.
Although both merozoites and infected RBCs will stimulate adaptive immune responses, the above equation is written such that only merozoites drive growth of adaptive immunity. This simplification is justified because infected RBCs and free merozoites maintain a relatively fixed ratio in the host, such that + ≈ ; thus, this ratio can simply be incorporated into the parameters and instead.  The equations for the adaptive immune responses each includes two decay terms, but the application of these terms depends on which type(s) are present in the host. The variable is definied such that = 1 if type is present, and = 0 otherwise. In addition, the variable ensures that does not decline below the baseline value : = 0 if ≤ , and = 1 otherwise.
The first decay term, with coefficient (1 − max( , )), is simply the slow, exponential decline of adaptive immunity in the absence of continued stimulation (the halflife being measured in years). If type is present, this term simplifies to zero. If both types are absent, the term simplifies to − . If type is absent but type is present, and as long Exposure to antigenic variants as > , the applicable decay term is − (1 − ). The parameter is the proportion of antigens that are shared betweens strains of type and type ; therefore, only the nonoverlapping proportion (1 − ) decays when type is absent but type is present. The second decay term, with coefficient max( 1 , 2 ), is applied whenever the host is infected with either type. When the host is infected (with either type, or both types), and as long as > , the applicable decay term is − (1 − ( + )). As described in more detail below, increases with time, the fraction ( + ) approaches 1, and the decay rate approaches zero. The relationship between and the decay rate is depicted in Figure 2. The function of this second decay term is to approximate the process of immune evasion through antigenic variant switching. Variant switching is thought to be stochastic in nature, although the degree of randomness is not known. For what follows, we assume that switching is at least approximately random (not heavily biased toward particular switching patterns), and that the sequential appearance of individual variants is driven by selection from adaptive immunity [2].
P. falciparum has a large, but finite, pool of variant antigens to switch through; for example, the size of the var gene repertoire is generally around 60 variants. If variant switching is approximately random, the time it takes to "find" a variant that is not recognized by the adaptive immune response is primarily a function of how many variants are already recognized. Early in the infection, almost any variant will not be recognized, so "escape" through switching should happen rapidly. However, when most variants have been seen by the immune system, it will take many more random switches to find one that has not been seen before. Assuming random switching, the number of switches to find a novel variant follows a geometric distribution, with mean 1− where is the proportion of variants that have not been seen yet ( Figure 3). Figure 3. The mean number of antigenic variants that must be 'tried' before finding a variant the immune system has not seen, as a function of the number of variants that have already been seen (out of 60 total variants).
Rather than explicitly model the dynamics of variants and variant-specific immune responses, we use this hypothesized relationship between the number of variants already seen and the time required to "find" a novel variant to implicitly model the process of antigenic variation. The switch to a novel variant impairs the ability of the adaptive immune system to recognize and kill parasites; this loss of effectiveness is mathematically indistinguishable from a loss of immune effectors, and can thus be represented by a decay term in the equation for adaptive immunity.
As described above, novel variants should be found rapidly at the start of an infection, but much more slowly as the pool of variants is exhausted. Therefore, the rate of decay of adaptive immunity should be high initially and decrease as the infection progresses. The variable exists to track the "progress" of an infection -i.e. how much of the variant repertoire has been "seen" by the adaptive immune system. We assume that only one variant is expressed at any given time, and therefore increases linearly with time. However, different strains can have variants in common, and any shared variant expressed by one has been "used up" for all. Therefore, type contributes to the increase of over time at a rate that is proportional to the overlap in the variant repertoires of strains of types and (the parameter ).

Parasite diversity and acquired immunity
For the purposes of this model, we assume that the parasite population is comprised of a virtually infinite pool of strains, such that every exposure is considered to be a new strain. Strains are classified phenotypically into drug-sensitive and drug-resistant 'types' but there is assumed to be no underlying population structure. Any two strains (whether of the same type or different types) are assumed to have a fixed amount of overlap in the proteins/antigens that are visible to the adaptive immune system; the amount of overlap determines the extent of cross-reactivity been strains. At the population level, greater cross-reactivity reduces the number of exposures required to reach a given 'degree' of acquired immunity. At the within-host level, greater cross-reactivity can increase the severity of immune-mediated 'apparent competition' in which the immune response generated by one strain nevertheless has a negative effect on both strains.
The overlap between strains is governed by two parameters. is the proportion of fixed (non-variant) antigens that are shared between any two strains, while is the proportion of variant antigens (such as PfEMP1) that are shared. There are two reasons that overlap of fixed antigens and overlap of variant antigens are considered separately. The first is simply that overlap in variant repertoires can be quite low (sometimes approaching zero). The second is that fixed and variant antigens have different effects on the dynamics of immunity. When fixed antigens are shared, it has the effect of boosting the immune response, whereas when variant antigens are shared, it hastens the exhaustion of each strain's variant repertoire. The logic is as follows: suppose two strains in the same host share a particular antigenic variant. When one of the strains expresses this variant, the adaptive immune system mounts a response against it. However, when the other strain switches to expressing this variant, the specific immunity acquired from previous exposure will not contribute much to control of parasite growth; instead, it will simply exert selection for other variants that are not yet recognized. Thus, any variant expressed by either strain has been 'used up' for both, which decreases the time until both strains run out of novel variants.
As mentioned above, we assume that any two strains share an equal proportion ( ) of their non-variant antigens; the proportion shared by strains is −1 . When a new strain infects a host that has previously encountered strains, the host's immune system will recognize a proportion (1 − (1 − ) ) of the new strain's antigens. Thus, in the model, every time a new strain of type is introduced to a host, the adaptive immune response is multiplied by the proportion of the new strain's antigens that are recognizable based on past exposures: × (1 − (1 − ) ) where is the number of past encounters with type .
Something similar is done for the variable , which tracks how much of the current strain's antigenic variant repertoire the immune system has seen. When a new strain of type is introduced, is multiplied by the proportion of the new strain's variants that have been seen before: × (1 − (1 − ) + ) where and are the number of past exposures to type and type , respectively.
Finally, the number of previous exposures affects the degree to which each type is affected by the acquired immune response to the other type. The rate of killing of type by acquired immunity to type is proportional to where = 1 − (1 − ) ( is as defined above).

Human-mosquito contact and parasite transmission
Every day, each human host is assigned to be bitten by a number of mosquitoes that is drawn from a Poisson distribution with mean . Each mosquito bites only one host per day, and which mosquitoes bite on any given day is random (mosquitoes can bite on sequential days but do not necessarily do so).
The probability that a mosquito is infected upon feeding on a host is determined by a function described by Churcher et al. If a mosquito is determined to be infected, the number of gametocytes of type picked up is drawn from a Poisson distribution with mean * where = type gametocytes/ and is the volume of a mosquito blood meal in . (Draws of zero gametocytes are disallowed because the number of mosquitoes infected is pre-determined by the gametocyte density-infectivity function shown above.) Infection in the mosquito has a latent period of days. After the latent period ends (simulating the appearance of sporozoites in the salivary glands), the mosquito becomes infective. When an infective mosquito bites a host, there is a fixed probability that sporozoites are transmitted to the host; if this occurs, the mosquito introduces sporozoites to the host. If gametocytes of type were originally, then drawing from a binomial distribution with size and probability 1 ( 1 + 2 ) ⁄ determines the number of sporozoites of each type transmitted. If parasites from blood meals have reached the infective stage, assuming is the number of type gametocytes acquired from blood meal (1 ≤ ≤ ), then a draw of size is made from a multinomial distribution with probability ( 1 ) ( + ) for type from blood meal . The sporozoites from different blood meals are considered separately for the purposes of tracking the human host's exposure to strains of each type. If the host receives sporozoites of type that were derived from 3 different blood meals, then the 'strain exposure count' for type increases by 3, even though the parasites were introduced by the same mosquito. However, when a mosquito acquires gametocytes from a host, the gametocytes of each type are considered to constitute one strain, even if they were derived from multiple introductions. Infection in the human host also has a latent period of days, which simulates the liver stage of the infection. At the end of the latent period, merozoites are released for each sporozoite introduced days before and are added to the circulating merozoites tracked by the within-host model. At this point, the host's 'exposure count' for each type is updated to reflect the number of strains represented among the newly-emerged parasites.

Populations and turnover
The human population consists of hosts with ages uniformly distributed between zero and the human lifespan, . A host that reaches age is replaced by a host of age zero with a 'clean slate' -no current infection or latent infection, no history of infection, and no immunity. The mosquito population is similar (except for having a much higher rate of turnover); the population consists of mosquitoes that are evenly divided between ages zero to (the mosquito lifespan), and each day the oldest mosquitoes are removed and replaced with new mosquitoes of age zero.

Treatment
The simulated use of antimalarial drugs is flexible in a few ways. Treatment can be made conditional on the total parasite density ( 1 + 2 ) exceeding a threshold, which can simulate treating only symptomatic infections or only infections detectable by standard diagnostic methods; alternatively, treatment can be administered to any infected host (by setting the threshold parasite density to zero) or to any host regardless of infection status (by setting the threshold to -1). Antimalarial drug use can be started in the middle of a simulation, to simulate introduction of a drug into a population at equilibrium. Treatment can also be restricted to start only on certain days, which can simulate mass drug administration (MDA) or mass screening and treatment (MSAT) where antimalarial drugs are administered en masse at regular intervals. Not all of these options are used in the simulations presented, but they provide opportunities to further explore the fate of drug resistance in scenarios not considered here. Maximum growth rate of adaptive immunity 10 3 Shape parameter (adaptive immunity growth curve) 0.1 Decay rate of adaptive immunity due to antigenic escape 120 Shape parameter (adaptive immunity decay due to antigenic escape) 8 Shape parameter (adaptive immunity decay due to antigenic escape) 10 −3 Background decay rate of adaptive immunity 10 −3 Starting value of adaptive immunity to each type Ω 10 −4 Extinction threshold (infected RBC density) 14 Duration of treatment (days) 14 Mosquito lifespan (days) 10 Latent period in mosquito (days) 12 Latent period (liver stage) in humans 2 Blood meal volume (μL) 12 Number of sporozoites introduced by each mosquito bite 10 4 Number of merozoites produced per sporozoite (liver stage) 3000 Human lifespan (days) 400 Human population size 1. 2