Drug resistance is a long-standing economic, veterinary and human health concern in human and animal populations. Efficacy of prophylactic drug treatments targeting a particular pathogen is often short-lived, as drug-resistant pathogens evolve and reach high frequency in a treated population. Methods to combat drug resistance are usually costly, including use of multiple drugs that are applied jointly or sequentially, or development of novel classes of drugs. Alternatively, there is growing interest in exploiting untreated host populations, refugia, for the management of drug resistance. Refugia do not experience selection for resistance, and serve as a reservoir for native, drug-susceptible pathogens. The force of infection from refugia may dilute the frequency of resistant pathogens in the treated population, potentially at an acceptable cost in terms of overall disease burden. We examine this concept using a simple mathematical model that captures the core mechanisms of transmission and selection common to many host–pathogen systems. We identify the roles of selection and gene flow in determining the utility of refugia.

1. Introduction

The emergence of drug-resistant pathogens is the natural outcome of the steady application of drugs to hosts. The emergence rate is governed by the population size of the pathogen, its mutation and recombination rates, drug efficacy, the intrinsic cost of resistance as it relates to the drug's biochemical targets—how mutation affects other aspects of pathogen physiology (epistasis or pleiotropy), and the intrinsic benefit of resistance in treated hosts [1,2]. Examples of drug resistance are numerous and the mechanisms contributing to resistance are diverse [1,3,4]. Consequently, the cost of maintaining effective interventions via the development and use of multiple drugs is considerable. These costs have driven the exploration of alternative strategies for integrated pathogen management, such as the use of refugia (untreated populations) for reducing levels of resistant pathogens in treated animals [5]. Importantly, refugia can serve to dilute drug-resistant pathogens through the force of infection from drug-susceptible pathogens.

Refugia may be exploited for managing drug resistance by enabling cross-host transmission between treated and untreated animal populations [6] and benefitting from inevitable cross-group transmission in the case of zoonoses and age-structured intervention strategies for human infectious diseases [7]. Refugia may increase the lifespan of drug products and classes at an acceptable cost in terms of overall pathogen pressure, and have been seriously considered in a range of systems including haemonchosis [6] and heartworm [8]. In the former, parasites may be shared between treated livestock and both untreated livestock and wildlife owing to free-living parasite stages whose transmission may be managed by co-grazing. In the latter, untreated canids may share parasites with treated companion animals via generalist-feeding mosquito species that vector the parasite, in which case understanding heterogeneity in mixing owing to spatial variation in vector and host abundance as well as drug compliance, may help to predict hotspots for emergence of drug resistance. However, there is little theory to guide predictions on short- and long-term consequences of refugia for levels of prevalence and resistance in target, treated host groups. Because the degree of transmission between treated and untreated groups is either partly controllable or likely to vary geographically [9], we give special attention to how degree of mixing and drug coverage affect prevalence and frequency of resistance in treated hosts via their modulation of the force of infection. The success and recommendation of refugia in managing drug resistance [6,10] serve as a motivation to understand the limits of this form of disease management, especially any counter-productive effects. Here, we present an illustrative, mathematical model for exploring the development of drug resistance in a host–pathogen system that includes refugia.

2. Model and analysis

(a) Model structure

We investigate the consequences of drug-susceptible and drug-resistant pathogen fitness differences in treated and untreated hosts using a transmission model for which there is a specified degree of cross-group transmission. Parameter ϕ determines the extent to which transmission is between-group rather than within-group. For ϕ = 0, all transmission is within-group (untreated hosts contact and infect only other untreated hosts, and similarly for treated hosts). Conversely, when ϕ = 0.5, transmission between-groups occurs at the same rate as transmission within group. Hypothetically, ϕ can take values greater than 0.5 (indicating more between-group transmission than within-group) but we do not consider this scenario in detail. We denote uninfected untreated and treated hosts by S and S_T, respectively. Each may become infected with drug-susceptible strain 1 or drug-resistant strain 2. Untreated infected hosts are denoted by I_X, where X = 1 refers to infection by the drug-susceptible strain and X = 2 by the drug-resistant strain. Similarly, T_X denotes treated infected hosts. Differences between strains in the two host environments are assumed to affect the transmission rates; transmission of each strain from untreated individuals occurs at rates β₁ and β₂, whereas treated individuals transmit at rates β_T1 and β_T2. In this way, E = β₁ – β_T1 measures the efficacy of prophylactic drug treatment, C = β₁ − β₂ measures the cost of resistance, and B = β_T2 − β_T1 measures the benefit of resistance (where E, C and B > 0). Typically, Inline Formula such that drug-susceptible strain 1 may be considered a specialist (high fitness in untreated hosts, low fitness in treated hosts) and drug-resistant strain 2 may be considered a generalist (intermediate fitness in the two host environments). The overall host population is capped at size θ/μ (the ratio of the birth rate to the per capita death rate), a constant fraction p is maintained in the treated group and all individuals recover back to their susceptible state at rate γ.

2.1

(b) Transient dynamics

We may rewrite equation (2.1) to focus on the epidemiological dynamics (where we sum infected types by treatment status such that Inline Formula and ):

2.2a

where

and

are the average transmission rates from untreated and treated hosts, respectively. The parameters Inline Formula

and

refer to the frequencies of resistant strains in untreated and treated hosts, respectively. Hence, evolution affects epidemiological dynamics via the change in strain frequencies, which is given by

2.2b

The above equations are fully equivalent to (2.1). However, distinguishing epidemiology (2.2a) and evolution (2.2b) better elucidates the forces affecting the dynamics of the system. For instance, in each type of host, three forces are acting on evolutionary dynamics: (i) the first term refers to the selection acting on the philopatric pathogens, (ii) the second term refers to selection coming from allopatric pathogens introduced by cross-group transmission, and (iii) the final term refers to the homogenizing effect of gene flow between the two types of hosts. The last term is likely the least intuitive because it depends heavily on epidemiology (i.e. the ratio of T/I and the differentiation Inline Formula between pathogen populations). These equations can be used to describe the transient dynamics but they also help explain the effect of the model parameters on the equilibrium state.

(c) Equilibrium state

Both the equilibrium frequency of resistance and equilibrium prevalence are strongly influenced by the degree of cross-group transmission (mixing, ϕ) and the proportion of all hosts treated (drug coverage, p) (figure 1).

Figure 1. Equilibrium frequency of resistance and prevalence in treated (T) (*a,b*) and untreated (U) (*c,d*) hosts from the deterministic model system (equation (2.1)) initialized with and parameter values β₁ = 1.0, β₂ = 0.8, β_T1 = 0.5, β_T2 = 0.8, ϕ = [0, 0.5], p = [0, 1], θ = 10.0, γ₁ = 1.0, γ₂ = 1.0 and μ = 1.0. (Online version in colour.)

When ϕ = 0.5, drug coverage has only a weak effect on prevalence of infection in treated hosts. Indeed, when coverage is very low, there is not drug resistance, but infections are maintained in treated hosts by cross-group transmission. When coverage is very high, drug resistance is fixed but the resistant strain is not a very good transmitter. Prevalence is minimized for intermediate levels of drug coverage.

Interestingly, when drug coverage is relatively low, a higher level of spatial clustering of host and/or pathogen philopatry (i.e. ϕ < 0.5) decreases prevalence in treated hosts. Indeed, lower values of ϕ boost the incidence of infection in the ‘good quality’ hosts (i.e. the untreated ones), which feeds back on the evolutionary dynamics. There is still selection for drug resistance in treated hosts, but gene flow introduces many drug-susceptible strains into the ‘bad quality’ hosts (i.e. the treated hosts). When drug coverage is sufficiently low this can prevent the spread of drug resistance, which has been demonstrated in spatially explicit models under certain assumptions of the distribution of treated hosts and pathogen dispersal [11].

When drug coverage is relatively high, however, selection for drug resistance overwhelms the evolutionary dynamics. The frequency of resistance is very high in treated hosts and very low in untreated hosts (because there is no longer selection from immigrants). Prevalence is maximized in treated hosts because this set of parameter values (high ϕ and high p) is creating the perfect habitat for the drug-resistant strain (i.e. a homogeneously resistant host population).

3. Discussion

Keeping pace with pathogens through the development and use of novel antimicrobials is increasingly challenging [12,13]. As a consequence, attention is turning to prudent use of drugs in management strategies designed to slow the development of drug resistance. These include diagnosing and treating only infected (or suspected infected) individuals [6], using dose control to limit the rate of drug resistance emergence within-hosts [14], personalized medicine [15], drug combinations [16] and refugia [10].

In seeking to establish the challenges and opportunities of using refugia to mitigate drug resistance, we have found great benefit in decoupling the epidemiological and evolutionary dynamics. For instance, malaria-specific models have identified threshold treatment levels to prevent fixation of drug resistance [17], but the evolutionary dynamics were not formally uncoupled to explain the threshold in terms of selection and gene flow. Our parsimonious refugia model, capturing essential features likely underlying many multi-host pathogen systems, demonstrates that refugia may lead to complex effects on the prevalence and frequency of resistance. Models tailored to specific host–pathogen systems will require extra complexity such as density-dependent pathogen regulation [18] and targeting only infected (and suspected infected) individuals for drug treatment [6].

The two main parameters that influence epidemiology and evolution are the drug coverage (p) and degree of mixing between treated and untreated hosts groups (ϕ). Optimizing the values of these two parameters requires a specified objective. If the objective is to reduce the spread of drug resistance, then one can derive the value of ϕ that allows maximizing drug coverage without the spread of drug resistance. Some segregation between untreated and treated populations is always favourable towards this goal (but not extreme values). Alternatively, if the objective is to reduce prevalence in treated hosts then this may be achieved with relatively low coverage and extreme segregation. In other words, a small zone of intensively treated individuals surrounded by untreated populations is best, but maybe not very realistic.

Beyond drug resistance, the use of refugia in Bacillus thuringiensis (Bt) transgenic crops has been extensively discussed and implemented [19]. To date, theoretical assessment of the utility of crop refuges (proximal crops supporting insecticide-susceptible pests) has been approached via computer simulations. Our analytical approach could be extended in this context, but would require the inclusion of diploid pests to account for the fitness differences between homozygous and heterozygous forms (a challenge that similarly exists for studying macroparasites), and potentially discrete-time models. Regardless of the likely increase in model complexity [20], we contend that there is much to be gained from teasing apart the epidemiological and evolutionary dynamics, and the guiding principles yielded from our simple model provide a good foundation for moving forward to specific systems.

Data accessibility

R code producing reported results on prevalence and frequency of drug resistance is provided in the electronic supplementary material.

Authors' contributions

A.W.P., J.H., R.K. and S.G. participated in the design of the study and drafted the manuscript; A.W.P., J.H. and S.G. performed the mathematical analysis and numerical simulation; all authors gave final approval for publication and agree to be held accountable for the content therein.