Epidemiological models for the spread of anti-malarial resistance

Background The spread of drug resistance is making malaria control increasingly difficult. Mathematical models for the transmission dynamics of drug sensitive and resistant strains can be a useful tool to help to understand the factors that influence the spread of drug resistance, and they can therefore help in the design of rational strategies for the control of drug resistance. Methods We present an epidemiological framework to investigate the spread of anti-malarial resistance. Several mathematical models, based on the familiar Macdonald-Ross model of malaria transmission, enable us to examine the processes and parameters that are critical in determining the spread of resistance. Results In our simplest model, resistance does not spread if the fraction of infected individuals treated is less than a threshold value; if drug treatment exceeds this threshold, resistance will eventually become fixed in the population. The threshold value is determined only by the rates of infection and the infectious periods of resistant and sensitive parasites in untreated and treated hosts, whereas the intensity of transmission has no influence on the threshold value. In more complex models, where hosts can be infected by multiple parasite strains or where treatment varies spatially, resistance is generally not fixed, but rather some level of sensitivity is often maintained in the population. Conclusions The models developed in this paper are a first step in understanding the epidemiology of anti-malarial resistance and evaluating strategies to reduce the spread of resistance. However, specific recommendations for the management of resistance need to wait until we have more data on the critical parameters underlying the spread of resistance: drug use, spatial variability of treatment and parasite migration among areas, and perhaps most importantly, cost of resistance.


Background
The treatment of malaria infections and the control of the spread of malaria are made increasingly difficult by drug resistant strains of the parasite. Resistance to chloroquine, an effective and safe anti-malarial that formed the first line of treatment, emerged more than 30 years ago [1], and since then malaria parasites have developed resist-ance against most of the widely used anti-malarials, including sulfadoxine-pyrimethamine [2,3], mefloquine [4] and quinine [5]. Indeed, it is now generally believed that the widespread use of a new drug will inevitably be followed by the appearance of resistance in the parasite population. There is, therefore, increasing emphasis in malaria control on rational management of drug use in an attempt, if not to prevent, then at least to delay the spread of resistance in the parasite population. One approach is to use combination therapy, where the simultaneous use of several drugs is hoped to delay the emergence of resistance [6,7]. Another approach could be to use drugs in a mosaic fashion (temporally or spatially), so that resistance to individual drugs is counter-selected where or when the drug is not used, as is discussed for the management of antibiotic resistance [8]. However, predicting the effect of different strategies for drug use on the spread of resistance to the drug is not straightforward because of the complex interactions that determine the transmission of the parasite between human and mosquito populations. Mathematical models have in the past provided a valuable framework for analysing the transmission dynamics of malaria [9,10]. These models have been widely used to consider the effect of different strategies such as vector control [11] and the use of future vaccines on the transmission dynamics of malaria [12,13].
In this paper, we extend these models to investigate the spread of drug resistance in malaria. We develop models of malaria transmission that allow us to examine the simplest case: the spread of resistance to a single drug, encoded at a single genetic locus. We modify standard epidemiological models of malaria transmission to incorporate drug sensitive and drug resistant parasite populations, and use this model to investigate the conditions under which drug resistance will spread in the population of hosts, as well as to predict the rate of spread of resistance. We then examine how superinfection or heterogeneity in treatment will modify the earlier results.
Our epidemiological models provide a counterweight to the recent models based on population genetic considerations [14][15][16]. One of the purposes of the paper is thus to discuss the similarities and differences in the predictions of the genetic and epidemiological models. But our main aims are to find which parameters would be necessary to predict the spread of resistance and to emphasize our lack of knowledge about any of these parameters.

Epidemiological framework
Our models are based on previously published extensions of the classical Macdonald-Ross model, which take into account acquired immunity by humans and allow the duration of immunity to depend on exposure (i.e. that immunity is boosted by additional infections) [17,18]. Though the assumptions about the development of immunity are undoubtedly overly simple, the model can predict several important epidemiological patterns, e.g. that the crude prevalence of malaria in a population may decrease at very high levels of transmission [19,20]. Furthermore, variations of the model (where immunity is life-long or where the duration of immunity is constant, rather than being boosted by new infections) gave similar predictions, so that our results are robust with respect to details of the dynamics of immunity.
The human population is split into a proportion x of susceptibles, a proportion y of infecteds and a proportion z of immunes (Fig. 1). The changes in the proportions of these three categories are described by three differential equations: Susceptibles become infected at the inoculation rate h (which is derived below). Infected individuals recover at a rate r to become immune, and immunes loose their immunity and become susceptible again at a rate ρ. Deaths occur at a rate δ (i.e. the life-expectancy is 1/δ) and are not affected by disease status. Deaths are balanced by births into the susceptible class, so that the population size remains constant.
If the duration of immunity is ω years in the absence of new infections and immunity is boosted by re-exposure,

Figure 1
Schematic illustration of the basic malaria model. As discussed more fully in the text, the diagram shows the categories of host and vector individuals considered in the model, the transitions between the categories and the transmission processes. reference [20] shows that the rate at which immunity is lost can be set to Consider next the infection rate of mosquitoes. The mosquito population is split into three categories: a proportion v of infected mosquitoes that have not yet developed infectious sporozoites, a proportion w of infectious individuals with sporozoites in the salivary glands, and a pro- Because the mosquito dynamics operate on a much faster time-scale than the human dynamics described by equation (1), the mosquito population can be considered to be at equilibrium with respect to changes in the human population. This allows its dynamics to be collapsed into the inoculation rate h: where m denotes the number of mosquitoes per human.
As is customary for epidemiological models, transmission is summarized by the basic reproductive number, which gives the number of secondary cases following introduction of a single infected individual into a susceptible population.
As an indicator of the intensity of transmission, R o will tend to be increased by natural selection, and parasites with higher R o will invade populations of parasites with lower reproductive numbers.

Spread of resistance
The basic model can easily be extended to describe the spread of resistant parasites. We assume that resistance is conferred by an allele at one locus, so that we can consider the spread of resistance as competition among two strains of parasites, an anti-malarial sensitive and a resistant one. We further assume that a fraction f of the infected individuals is treated. With these assumptions, the compartment model of Figure 1 is changed to allow for three sub-compartments of infected individuals: hosts that are infected with the sensitive strain but that are not treated, hosts infected with the sensitive strain and treated, and hosts infected with the resistant strain, which is not affected by treatment (Fig. 2). With the subscripts S and R denoting sensitive and resistant strains, respectively, and the subscripts U and T denoting untreated and treated infections, equations (1) are modified to where is the change of transmission due to a change from sensitivity to resistance; it denotes the transmission cost of resistance. Note that, although we have modelled the transmission cost in terms of infectivity from the human to the mosquito, the cost can also be expressed in terms of other transmission parameters. Mosquitoes harbouring resistant parasites, for example, may bite less frequently or experience higher mortality than those infected with sensitive parasites.
We further assume that the differences in recovery rate are due to treatment and to a cost of resistance within the human host. Therefore We can estimate the relative transmission of the sensitive and resistant parasites with their basic reproductive numbers. Letting equation (5) denote the basic reproductive number of the untreated, sensitive parasites and assuming that the average life-span of humans is considerably greater than the duration of infection (1/δ >> 1/r SU ), the basic reproductive number of the sensitive parasites (in the absence of resistant ones) becomes where R o equals the basic reproductive number defined in equation (5), i.e. for sensitive parasites in the absence of treatment. The basic reproductive number of the sensitive infections declines linearly with increasing treatment. In the absence of resistance treatment can eradicate the disease, if it drives the basic reproductive number to less than unity, i.e. if the fraction of infections treated is larger than a critical value f e given by: If we introduce resistant parasites, the outcome depends on the relative magnitude of the basic reproductive numbers of the sensitive and resistant infections. The latter is given by where Γ = γ t r R /r SU describes the reduction of the intensity of transmission due to a switch from sensitivity to resistance and can therefore be defined as the total cost of resistance. The long-term outcome will depend on whether R R is larger than R S , which in turn depends on whether the proportion of individuals receiving treatment exceeds a threshold f c defined as If the proportion of treated individuals is less than the threshold f c , R S >R R and all parasites in the population remain sensitive. If the proportion of treated individuals is higher than f c , all parasites will eventually become resistant. Note that this proportion depends only on two parameters: the cost of resistance and the ratio of the periods of infection in treated and untreated individuals. In particular, the threshold proportion is independent of the Modification of the illustration shown in Fig. 1 to incorporate resistant and sensitive parasites. The infected humans are split into 3 categories, representing untreated-sensitive, treated-sensitive, and resistant infections. The mosquitoes are split into two cohorts, representing sensitive and resistant infections. intensity of transmission (measured as R o , the basic reproductive number of the sensitive strain in the absence of treatment). Similarly, the rate at which resistance spreads, once the threshold proportion of treatment is exceeded, depends largely on the frequency of treatment and the cost of treatment, but is almost independent of the intensity of transmission (Fig. 3).
Thus, in this model, there is competitive exclusion between resistant and sensitive parasites. Whether sensitivity or resistance predominates in the population depends on three parameters: (i) the proportion of infected individuals that receive treatment, (ii) the effectiveness of treatment defined in terms of period of infection as ε = 1 -r SU / r ST , and (iii) the cost of resistance.

Superinfection of hosts
In the previous model we assumed that hosts cannot be superinfected by a second parasite, which would lead to mixture of sensitive and resistant parasites in the same patient. In this section we relax this assumption and allow superinfection of hosts, assuming that sensitive and resistant parasites develop independently within their host. We further assume that the treatment has a short half-life, so that the decision to treat superinfections is independent of any earlier treatment. In other words, a proportion f of all superinfections will be treated, while a proportion 1 -f will not be treated, irrespective of whether the first infection was treated or not. (We also analyzed a model where the treatment has a long-life and where the superinfecting parasite therefore receives the same treatment as the first infection. In this case, treated individuals cannot be superinfected with a sensitive parasite. The conclusions are the same as the ones reached with the equations given here). With these assumptions equations (6) can be modified to where the subscript M denotes mixed infections. The inoculation rates of sensitive and resistant parasites are changed accordingly as where y S = y ST + y SU and y M = y MT + y MU .
With this model the parasite population switches from being predominantly sensitive to being predominantly resistant at a similar frequency as in the model without superinfection (defined by equation 13). The main difference due to superinfection is that competitive exclusion is not always achieved. At intermediate levels of treatment resistant and sensitive parasites can coexist in the same population (Fig. 4). The area where coexistence is possible increases as the cost of resistance decreases. Furthermore, if the cost is sufficiently low, resistance can invade a sensitive population even when there is no treatment, so that  (6) and (7) were integrated with only sensitive parasites until an equilibrium was reached. Then a resistant parasites were introduced at a frequency of 10 -4 and the equations were integrated until the frequency of resistance reached 50%. The parameters used to integrate the equations were: δ = 0.01 year -1 , µ = 36.5 year -1 (ca. 0.1 day -1 ), τ = 0.03 years (ca. 10 days), a = 100 year -1 (ca. 1/3 day -1 ), ω = 2 years, r SU = 1 year -1 , r ST = 10 year -1 (ca. 1 month -1 ) and r R = 5 year -1 . The mosquito density m was calculated from equation (5) and the given parameters. For these parameters the cost of resistance Γ is 5 and the threshold frequency f o is 0.89. The possibility of coexistence of two parasite strains comes as no surprise. In many models of epidemiology, superinfection allows the coexistence of many strains, and in particular the coexistence of strains with considerable differences of transmission success [21].

Spatial variability in treatment rates
In this section, we relax a further assumption of the basic model: that drug treatment is applied uniformly to a population of hosts. In a natural situation, it is more likely that treatment levels will vary among host individuals and among host populations. We therefore ask, to what extent will gene flow via migration counteract the selective forces due to treatment in a local environment? To do so, we extend the model defined by equations (6) and (7) by assuming that the gene flow is due to migration of the mosquito vectors infected by sensitive or resistant parasites.
We will consider two cases. First, following work by Comins [22] we assume that there is a reservoir of untreated humans, where all parasites are sensitive. In the treated area, infection will be due partly to parasites surviving the treatment in this area and partly to parasites that have migrated from the untreated into the treated area. Second, we consider several human populations with differing treatment rates. Infected vectors migrate among neighbouring populations, so that they can migrate from the treated area to the untreated one.
In the first of these possibilities we assume that the epidemiology in the untreated area follows the basic model formulated in equations (1)-(4), and that the epidemiology in the treated area is determined by the competition between the resistant and the sensitive strains according to equations (6) and (7). The equations are modified, so that a proportion m of the infected mosquitoes migrates from the untreated area to the treated one. In this case the migration of sensitive parasites into the treated area prevents the selection of resistance, if the migration rate and the cost of resistance are high enough (Fig. 5). For lower migration rates, resistance in the treated area increases towards an intermediate level, which is often close to fixation. These results are similar to those obtained previously [22] for the evolution of insecticide resistance.
In the more general case, migration from the treated area into the untreated zone will not be negligible. We model this situation by incorporating two areas that are identical except for their levels of treatment; in one area all infections are treated, in the other none are. In both, the epidemiology is modelled according to equations (6) and (7), and a proportion m (which is independent of area or resistance) of the infected mosquitoes migrates from one area to the other. At sufficiently low migration rates the two areas behave independently, with competitive exclusion within each area; sensitive parasites predominate in the untreated area (Fig. 6a) and resistant ones predominate in the treated area (Fig. 6b). At intermediate migration rates, competitive exclusion between sensitive and resistant parasites may break down. The level of resistance maintained in either area depends on the combination of the migration rate and the cost of resistance. A migration rate of 0.5 corresponds to a well-mixed population where half of the infections are treated. If migration rates exceed this value, the population as a whole is either completely sensitive or completely resistant; competitive exclusion operates at the level of the meta-population.

Summary of results
Our models enable us to explore the consequences of changing the level of drug use on the frequency of resistant malaria parasites. In the simplest model, we find a

Figure 4
Spread of resistance in the model that allows superinfection of individual hosts by resistant and sensitive parasites. The graph shows the proportion of parasites that are resistant when the population has reached equilibrium as a function of the cost of resistance (defined by equation 12) and the proportion of infections treated. The mosquito density m was set to 1, the recovery rate of resistant infections r R was determined by the cost of resistance, and the other parameters were as given in Fig. 3. threshold level of drug use, below which resistance cannot spread and above which resistance spreads and eventually becomes fixed. The threshold level is determined by the extent to which treatment of sensitive parasites reduces their duration of infection and by the cost of resistance, which is defined by the transmissibility of resistant parasites and by their duration of infection relative to untreated infections with sensitive parasites. This competitive exclusion between sensitive and resistant infections arises, because the parasite with the higher basic reproductive number R o drives the number of susceptible hosts below the threshold density that would be necessary to maintain the parasite with the lower basic reproductive number [23]. This result is the counterpart to earlier studies on the spread of resistance in directly transmitted infections [8]. It reflects that resistance is detrimental to the transmission of the parasite in untreated hosts, but beneficial in treated individuals, so that increasing drug use beyond a critical threshold tips the balance from 0% to 100% prevalence of resistance.
Allowing infection of hosts with both types of parasite gives similar results. The main differences are that, under some combinations of parameters, the resistant and sensitive parasites can coexist and that, for a low cost of resistance, resistant parasites can coexist with sensitive ones in the absence of treatment. This reflects the general conclusion that superinfection can maintain several parasite strains with different basic reproductive numbers [21]. The parasite with low transmission success can infect hosts that are already infected by the more rapidly transmitted parasite. Therefore, evolution does not lead to competitive exclusion and does not maximize transmission, in contrast to models without superinfection.
Spatial variability in treatment rates can have a larger effect on the spread of resistance. High migration between areas with and without treatment can swamp the local effects of selection. Thus, when the cost of resistance is high, sensitive parasites can dominate in areas where every host receives treatment; when the cost of resistance is low, all parasites in areas where no treatment is available can become resistant. At intermediate costs of resistance and migration rates, coexistence of sensitive and resistant parasites is possible within each area. This coexistence arises as a balance between selection within the populations and migration between populations. Migrants (with the inappropriate type) can persist for some time before they are removed by selection (drug pressure or the cost of resistance) occurring within each population. Therefore, at some ratios of migration and selection, both types of parasites can coexist.
Thus, from the small number of models described here, it is clear that knowledge about several parameters is critical for our understanding of the spread of resistance: the drug use, the spatial variability of treatment and parasite migration among areas, and perhaps most importantly, the cost of resistance.

Any data?
Unfortunately, data on any of these parameters are very limited, and almost completely restricted to chloroquine. Even data on the use of chloroquine, which may seem to be the easiest data to get, are poor. Though we may be able to estimate the total drug use in an area (e.g. 191 tons of chloroquine in Africa in 1988 [24]), though we may also obtain the consumption of the drug in an area (e.g. in endemic areas, large proportions of children have chloroquine residues in their blood at any given time [25]) and estimate the number of patients with clinical malaria receiving the drug, the more critical questions are: What proportion of infected people (as opposed to sick people) receive treatment, and is the dose of the drug sufficient to clear the parasites?  That resistance, at least against chloroquine, is costly in terms of transmission is suggested by the observation that resistance has dropped significantly in several areas when choloroquine was no longer used. In contrast, however, some data suggest that resistance is associated with a benefit, rather than a cost, in terms of transmission from humans to mosquitoes [26][27][28]. This is supported by in vitro studies of parasite development showing that resistant parasites grow more slowly than sensitive parasites, but develop their transmission stages earlier and at a higher rate [29]. However, additional factors must balance and outweigh this apparent benefit of resistance, otherwise resistance would have spread before the wide-spread use of drugs. The genetic basis for the spread is given. Even before the use of a drug, resistance is present at low levels in many populations, in particular for quinine and proguanil [30]. This probably applies for chloroquine as well, although the data are less convincing [31]). Possible factors leading to a cost of resistance might be related to the mortality of infected mosquitoes, since the parasite's transmission rate is more sensitive to mosquito mortality than to the effectiveness of transmission from humans to mosquitoes [9]. Thus, the more effective transmission from humans to mosquitoes could be more than compensated by the higher mortality of heavily infected mosquitoes [32,33]. Overall, there is some qualitative evidence for a cost of resistance, although quantitative measures of such a cost are absent.

Critique of the models
Obviously, like other mathematical (or non-mathematical) models, our equations are simple caricatures of the epidemiological processes being considered. The success of a model will depend on the inclusion of the relevant processes. There are several reasons for believing that the models we propose capture some aspect of reality. Compartment models of infectious diseases in general have had considerable success in epidemiology [34]. In particular, compartment models similar to the ones used here can describe several epidemiological patterns of malaria [11], including, for example, the observation that increased transmission can decrease the malaria prevalence   [20]. Similar models are considered helpful in evaluating the sensitivity of malaria transmission to different control measures [9] and in predicting the effectiveness of vaccine programmes [13,35]. Compartment models have also been used to examine the epidemiology of antimicrobial resistance in acute infections [8] as well as hospital acquired infections [36] and specific infections such as tuberculosis [37] and HIV [38]. The current paper is a first attempt at applying them to anti-malarial resistance.
There are a number of specific problems with modelling the epidemiology of resistance to malaria, which we briefly discuss below. One set of problems concerns the epidemiological processes of malaria. In this category are problems associated with the lack of our understanding of the nature of immunity to malaria at the individual and epidemiological levels. Thus, the equations describing the development of immunity reflect the true process only crudely. However, simulations with two extreme alternative models (life-long immunity after the first infection or immunity without boosting, i.e. that is lost at a constant rate) suggest that our results are insensitive to the details of the immunological process. Furthermore, we have assumed that treatment does not affect the development of immunity. While we have no data on this issue for malaria, it is a reasonable assumption that is based on observations in other systems. The T-cell response to bacteria, for example, is programmed early in an infection (before the peak of the inflammatory response) and independent of the duration and severity of the bacterial infection. In other words, the infection starts the program, but has little influence on the further development of the immune response [39].
A second set of problems is associated with the genetics of malaria parasites. In this category falls the possibility that malaria is structured in distinct strains [40,41]. This could slow the evolution of drug resistance if resistance is due to several genetic loci [15]. This model and other models that take into account details of population genetics [14][15][16], as have population genetic models of the resistance of other parasites [34], show the importance and uses of population genetics. However, just as our equations neglect population genetics, these models have neglected the epidemiological processes underlying the evolutionary changes. From our models it appears that the neglect of epidemiology could qualitatively bias the results; the differences between competitive exclusion and coexistence, for example, can only be reached by explicitly contrasting different epidemiological assumptions and settings. Perhaps the best approach would be to merge epidemiological and genetic models. This approach has been taken several times in the recent literature with models taking into account individual infections [42][43][44], but it has some major disadvantages. Most importantly, because of the complexity of the models some important epidemiological features must be omitted. These models have not (or only in basic terms) taken account of acquired immunity, which is an important aspect of malaria epidemiology. If merging transmission and genetics is too complex, it remains to be seen whether focusing on the genetic side or on the epidemiological side will lead to better predictions. Certainly, where the genetics are simple (one locus), the epidemiological approach is more promising.

Conclusions
The approach offered here can provide insights into several aspects of the epidemiology of anti-malarial resistance. It remains unclear why chloroquine resistance has emerged much later in sub-Saharan Africa than in South-East Asia [24]. This fact is often attributed to the high transmission in Africa, with the assumption that high transmission (and therefore strong acquired and innate immunity) leads to a slow spread of resistance [31]. This assumption has indeed been found in some models of the spread of resistance, because intense transmission is associated with less inbreeding, which in turn can lead to slower spread of resistance under some conditions [15,16]. However, this prediction is not universal, but depends on the details of the assumptions underlying transmission. Thus our model predicts either little or the opposite effect of transmission rate, while another recent model [45] shows that the within-host dynamics of malaria help to determine whether resistance is more likely to emerge in areas with high or with low transmission. But perhaps the explanation is much simpler: the use of chloroquine was fairly low until the mid-70s [24], so that resistance would be maintained at a low level due to the cost of resistance. As the chloroquine consumption increased to a peak of 191 tons in 1988 [24], it brought a rapid increase of resistance. Unfortunately, we do not have the quantitative data required to explore this possibility in more detail.
A further puzzle is that pyrimethamine resistance often spreads to areas not submitted to pyrimethamine pressure [46]. One unlikely explanation is that resistance is associated with a benefit, rather than a cost. Another explanation could be that resistance spreads from areas with high drug consumption, provided the migration rate is sufficiently high and the cost of resistance is low.
Similarly, there appears to be a low, but non-negligible level of resistance in most geographical areas even before exposure to drugs [31]. A possible explanation for this is a balance between mutation and selection [16], but it might also be explained by superinfection, which allows resist-ance to be maintained at low levels in areas in the absence of treatment, if resistance is not too costly.
Overall, an epidemiological approach could lead to valuable insights into the spread of resistance. What is needed is, on the one hand, investigation of the processes determining the spread and, on the other hand, more detailed knowledge of the critical parameters involved, e.g. cost of resistance. Only this will enable us to give advice about rational drug management that could slow the spread of resistance.

Authors' contributions
The authors contributed equally to the formulation of the models. JK undertook the calculations. Both authors edited, read and approved the final manuscript. RA acknowledges support from the NIH (R29-GM-54268 and AI 49334).