Visitors   Views   Downloads
PEER-REVIEWED

Introduction

The evolution of drug-resistance in many infectious diseases has proven to be one of the most challenging problems of human health in modern medicine. While a number of evolutionary mechanisms are generic for the rise of resistance, there are several processes which are specific to the drugs and treatment regimens (Zur Wiesch et al., 2011). These processes are often characterized by the competition between the resistant and wild pathogen subtypes. If no resistant-type exists before the start of treatment, the likelihood of an emergent resistant mutant out-competing the wild-type and dominating the pathogen population during treatment depends greatly on its fitness advantage over the wild-type (Lipsitch, 2001; Zur Wiesch et al., 2011). For a sufficiently low replicative fitness, the resistant-type may still be out-competed by the wild-type, even under a strong selection pressure of drugs (Ferguson et al., 2003; Moghadas, 2008). However, if the difference between the intrinsic fitness of the two pathogen subtypes is sufficiently small, then the selection pressure of drugs can overturn the competitive dynamics in favour of the resistant-type (Lipsitch, 2001; Ferguson et al., 2003; Moghadas, 2008; Zur Wiesch et al., 2011).

Within the fitness landscape of the pathogen subtypes, a key factor determining the competitive dynamics between the resistant-type and wild-type is the probability of developing resistance, which may be affected by the time for start of the treatment. For example, in influenza infection and in the absence of resistant mutants, early treatment inhibits viral replication, which will therefore minimize the likelihood of evolving resistant mutants. However, when resistant mutants are present, early treatment can provide an opportunity for the outgrowth of such mutants under the pressure of drugs (Alexander et al., 2007; Handel, Longini Jr & Antia, 2007). In this context, a major barrier to this growth is the generation of within-host adaptive immune responses. It has been shown that drug resistance is less likely to develop if the immune responses are maintained above a certain threshold during treatment (Lloyd & Wodarz, 2006; Wodarz, 2001; Wodarz & Lloyd, 2004), suggesting that a delay in start of the treatment may be beneficial in preventing resistance emergence given the timelines for developing immune responses (Wodarz, 2001). Acquisition of resistance during treatment is thus affected not just by the fitness of resistant-type, but also by the probability of resistant mutants occurring, which may be viewed as a function of delay in the initiation of treatment.

Previous studies illustrate how treatment coverage and efficacy influence the population dynamics of drug-resistance spread, while considering the fitness of resistant-type but leaving out the possibility of time-dependent probability of resistance development (Alexander et al., 2007; Huijben et al., 2013; Lipsitch et al., 2007; Legros & Bonhoeffer, 2016; Moghadas, 2008; Moghadas et al., 2008; Xiao, Brauer & Moghadas, 2016). Here we consider this probability as a function of delay in start of the treatment and discuss both short-term (epidemic) and long-term (endemic) states of the resistant-type and with-type infections. For the endemic state, we characterize the behaviour of disease equilibria and obtain a condition to minimize the fraction of population infectious at equilibrium with a delay in start of the treatment within the infectious period. For the epidemic state, we show that in contrast with previous work on the short-term disease dynamics (Lipsitch et al., 2007; Moghadas, 2008; Moghadas et al., 2008; Xiao, Brauer & Moghadas, 2016), the minimum epidemic size may be achieved with different treatment coverages (i.e., the fraction of infected population treated) depending on the functional form of the probability of resistance development.

Methods

The model

We consider a population that is fully susceptible to the invading pathogen. The model includes two pathogen subtypes, namely, the wild-type and resistant-type. We assume that the resistant-type may prevail during the treatment of individuals infectious with the wild-type. The probability of developing resistance is assumed to be a function of delay in start of the treatment, and we assume that this delay is the same for everyone in the population. To describe the transmission dynamics between susceptible and infectious individuals, we use a mass-action incidence in a homogeneously mixing population, and assume that the resistant-type emerges with a relative transmission fitness (δ) compared to that of the wild-type, where δ = 1 describes the absence of fitness advantage or disadvantage for direct transmission of resistance. For many diseases (e.g., influenza), resistance generally emerges with a negative fitness advantage (δ < 1); however, fitness costs of resistance can be ameliorated with compensatory mutations (Handel, Regoes & Antia, 2006; Rimmelzwaan et al., 2005). There are also diseases (e.g., tuberculosis) for which resistance may initially emerge with a positive fitness advantage (δ > 1) (Borrell & Gagneux, 2009; Casali et al., 2014). In the model presented here, we consider a delay τ in the treatment initiation during the infectious period, and denote the delay-dependent probability of developing resistance by q(τ). We also assume that the treatment is ineffective against the resistant-type infection.

Denoting the class of susceptible individuals by S, and the classes of individuals infectious with the wild-type and resistant-type by Iw and Ir, respectively, the dynamics of disease propagation is governed by the time-dependent differential equations system: S t = μ β S t I w t + δ I r t μ S t , I w t = β S t I w t β S t τ I w t τ e μ + γ τ γ I w t μ I w t , I r t = δ β S t I r t + q τ β S t τ I w t τ e μ + γ τ γ I r t μ I r t , T t = 1 q τ β S t τ I w t τ e μ + γ τ + γ I w t + I r t μ T t . where 1∕μ is the average lifetime of the host population, and 1∕γ is the average infectious period (assumed to be the same for both wild-type and resistant-type infections). Here, we omit the disease-induced death rate and assume that the disease has no effect on the average lifetime of the infectious individuals. A more general case where the average lifetime of the infectious individuals differs from that of the healthy individuals is considered in the Supplemental Information.

In this model, we assume that each infectious individual will receive treatment precisely τ units of time after becoming infected. The term βS(t)Iw(t) gives the cohort of individuals who enter the class Iw at any time t. Since the times until recovery and death are exponentially distributed with means 1∕γ and 1∕μ, respectively, the probability for an individual to remain in the Iw class τ units of time after entering the class is e−(μ+γ)τ. Thus, the transition out of the Iw class by means of treatment at any time t can be described by the term βS(t − τ)Iw(t − τ)e−(μ+γ)τ, accounting for the cohort of individuals who were infected by the wild-type pathogen at time t − τ and therefore receive treatment at time t (i.e., τ units of time after becoming infected). Treated individuals will either develop resistance with the probability q(τ) (and move to the Ir class), or become effectively treated (and move to the T class). Since the treatment is ineffective against resistance, for simplicity we included treated and untreated individuals infectious with the resistant-type in the same class Ir (Fig. 1). Lastly, the class T represents all individuals who are recovered from infection or effectively treated.

Schematic diagram of the transitions between model compartments.

Figure 1: Schematic diagram of the transitions between model compartments.

A careful bookkeeping of transitions between the compartments allows us to describe the dynamics by a system of differential equations, presented in (1). As described above, the system dynamics at any time t is not only determined by the present number of individuals in each compartment, but it also requires the infection rate of the susceptible individuals in the past time t − τ to account for those who receive treatment at the present time t. Therefore, we need to keep track of the history of the S and Iw compartments, which can be done both analytically and computationally (Smith, 2010). The computational codes for the system dynamics were developed in Matlab and are provided as additional Supplemental Information. For the analysis presented here, since re-infection is precluded, we omit the class (T) of individuals who are recovered from infection or effectively treated. This model indicates that with the introduction of a resistant-type infection (either initially or during treatment), a self-sustaining epidemic of resistant infections can occur through direct transmission. A more detailed description of the model is provided in the Supplemental Information.

Parameter space for the existence of model equilibria in terms of the reproduction numbers and delay in start of the treatment.

Figure 2: Parameter space for the existence of model equilibria in terms of the reproduction numbers and delay in start of the treatment.

Equilibria

To characterize the equilibria of the model, we define R0 = β∕(μ + γ), and apply the next generation method to express the effective reproduction number (Reff) as the spectral radius of the matrix R 1 0 δ R 2 R 3 , where R 1 = R 0 1 e μ + γ τ , R 2 = q τ R 0 e μ + γ τ , R 3 = δ R 0 .

With this notation, Reff = max{R1R3}. Each of these quantities provide important information on the spread of disease by each pathogen subtype. Specifically, R1 represents the number of secondary cases of the wild-type infection generated by an infectious case of the wild-type before developing resistance; δR2 represents the number of secondary cases of the resistant-type infection generated by an infectious case of the wild-type after developing resistance; and R3 is the number of secondary cases of the resistant-type infection generated by an infectious case of the resistant-type. In this formulation, R1 approaches R0 for large τ, and R0 represents the number of secondary cases generated by a single case of the wild-type infection in the absence of treatment. As illustrated in Fig. 2, we may characterize the existence of the model equilibria using the thresholds of persistence. For simplicity, we describe the equilibria by their infection components of the pathogen subtypes. The infection-free equilibrium E0 = (0, 0) exists unconditionally. The resistant-type equilibrium E r = 0 , I r ̂ , at which the wild-type infection is absent, exists only if R3 > 1. The cotype equilibrium E = I w , I r exists whenever both conditions R1 > 1 and R1 > R3 hold.

Delay-dependent analysis

We focus our analysis on the equilibria when R0 > 1 to establish a relationship between the delay in start of the treatment and the fraction of population infectious at equilibria. As we will show, this fraction is affected by two key parameters: the probability of developing resistance, and the relative transmission fitness of the resistant-type compared with that of the wild-type.

When the resistant-type emerges with a positive fitness advantage (i.e., δ ≥ 1 or equivalently R3 ≥ R0), Er is the only infection equilibrium of the model, at which the fraction of population infectious is independent of the delay in start of the treatment, given by I r ̂ = μ μ + γ 1 1 R 3 . If the resistant-type emerges with a negative fitness advantage (i.e., δ < 1 or equivalently R3 < R0), then there is a critical value τ0 =  − log(1 − δ)∕(μ + γ) at which R1(τ0) = R3 (Fig. 2). For τ < τ0 (R3 > R1 > 1), similar to the case of δ ≥ 1, Er remains the only infection equilibrium of the model.

At τ = τ0, the cotype equilibrium E emerges and exists along with the resistant-type equilibrium Er for τ ≥ τ0. We define the fraction of population infectious at the cotype equilibrium as a function of τ, and let G τ = I w τ + I r τ represent this function. At E, the value of G(τ) depends on the probability of developing resistance during treatment in addition to the transmission fitness of the resistance-type. Here we obtain the condition to characterize the behaviour of G(τ) at τ0. Taking the implicit derivative of G(τ) at τ0 gives G τ = d d τ I w τ + I r τ = d d τ R 1 R 3 R 2 I r τ + 1 + R 1 R 3 R 2 d d τ I r τ = R 1 R 2 R 1 R 3 R 2 R 2 2 I r τ + 1 + R 1 R 3 R 2 d d τ I r τ .

Since R 1 τ 0 R 2 τ 0 = μ + γ q τ 0 , we have G τ 0 = μ + γ q τ 0 I r τ 0 + d d τ I r τ 0 . Taking the derivative of infection component with the resistant-type at τ0 (see Supplemental Information) gives G τ 0 = μ q τ 0 1 δ 1 q τ 0 1 + 1 R 3 . This implies that the fraction of population infectious at E decreases if q τ 0 < 1 1 R 3 and increases otherwise. As we show in the Supplemental Information, the condition (4) also holds true for the more general case with the inclusion of disease-induced death rate in the Iw and Ir classes. When (4) holds, the start of the treatment with delay τ > τ0 reduces the fraction of population infectious at the cotype equilibrium. However, for sufficiently large τ, G(τ) increases and saturates at G = lim τ G τ = μ μ + γ 1 1 R 0 . One can show that G(τ) < G for any τ (see Supplemental Information), indicating that when δ < 1, treatment reduces the population level of infection at the cotype equilibrium, regardless of the probability of developing resistance. We also note that G > I r ̂ . This suggests that, given (4), there is an intermediate interval for delay in start of the treatment after τ0, in which the magnitude of infection at the cotype equilibrium reduces below that of the resistant-type equilibrium (Fig. 3A).

The fraction of population infectious at the equilibria.

Figure 3: The fraction of population infectious at the equilibria.

Parameters are: R0 = 3, μ = 1∕70 year−1, and γ = 1∕5 day−1. (A) δ = 0.4, τ0 = 2.55 days, τ1 = 3.4 days, q(τ) = 2e∕(1 + e−2(τ−1)); with a = 1.5 (solid magenta curve) and a = 0.5 (dashed magenta curve). (B): q(τ) = 0.6 for τ ≤ τ0 − 0.06 and q(τ) = 0.6e−2(ττ0+0.06) for τ > τ0 − 0.06, with δ = 0.38 (solid magenta curve), δ = 0.42 (dashed magenta curve), and δ = 0.46 (dotted magenta curve). Solid and dashed horizontal lines correspond to the fraction of population infectious at the resistant-type equilibrium.

When q(τ0) > 1 − 1∕R3, the magnitude of infection at the cotype equilibrium increases for τ > τ0 and sufficiently close to τ0. The continual increase in the population level of infection for τ > τ0 depends on the functional form of the probability of resistance development. As shown in our simulations (Fig. 3B), it is possible to reduce the fraction of population infectious at the cotype equilibrium below that of the resistant-type for some intermediate interval of τ that remains at positive distance from τ0. Characterizing q(τ) and the conditions for the existence of such an intermediate interval is not an easy task, but we show this possibility in our simulations.

When the demographics are omitted (μ = 0), the model represents an epidemic scenario without exhibiting any endemic states. In this case, a similar question on the fraction of population infected during the epidemic may be formulated in terms of the epidemic final size (i.e., the total number of infections throughout the epidemic), given by (see Supplemental Information) F = γ 1 q τ e γ τ 1 e γ τ 0 I w t d t + γ 0 I r t d t . Previous work on the population dynamics of drug-resistance emergence without demographics focuses on the fraction of infected population treated to quantify the epidemic size (Lipsitch et al., 2007; Moghadas, 2008; Moghadas et al., 2008; Xiao, Brauer & Moghadas, 2016). Summarizing the outcomes, these models demonstrate that there is an optimal treatment coverage at which the epidemic final size is minimized. The uniqueness of such optimal coverage has not been proven (Xiao, Brauer & Moghadas, 2016). Our simulations for the model presented here show that, for a given transmission fitness of the resistant-type, different delays in start of the treatment may lead to the same epidemic final size depending on the probability of resistance development at the time of treatment initiation (Figs. 4 and 5). This effectively means that the minimum epidemic final size may be achieved with different population levels of treatment.

Final size of the epidemic and resistant infections as functions of delay in start of the treatment and the transmission fitness of the resistant-type.

Figure 4: Final size of the epidemic and resistant infections as functions of delay in start of the treatment and the transmission fitness of the resistant-type.

Parameters are: R0 = 3, γ = 1∕5 day−1, and the functional forms of q(τ) in (A, D), (B, E) and (C, F) correspond to the red, black, and blue curves illustrated in Fig. S2.
Illustration of the parametric curves (F(τ), e−γτF(τ)) for the epidemic final size and the fraction of infected population treated, with τ as an independent variable.

Figure 5: Illustration of the parametric curves (F(τ), eγτF(τ)) for the epidemic final size and the fraction of infected population treated, with τ as an independent variable.

Parameters are: R0 = 2.2, γ = 1∕5 day−1, δ = 0.48, and the functional forms of q in black and red curves are given by black and red curves shown in Fig. S2.

Simulation Results

To illustrate the theoretical results, based on the competitive dynamics between the wild-type and resistant-type, we simulated the cotype equilibrium by varying the delay in start of the treatment within the average infectious period of 1∕γ = 5 days (Table 1). We assumed a reproduction number R0 = 3 in the absence of treatment, and fixed the transmission fitness of resistant-type at δ = 0.4, giving 1 − 1∕R3 = 0.1667. Figure 3A shows the curves of equilibria for two different functional forms of q(τ) = 2e∕(1 + e−2(τ−1)), as illustrated in Fig. S1. The solid curve represents the emergence of the cotype equilibrium at τ0 = 2.55, where a = 1.5 and q(τ0) = 0.0418 satisfies the condition (4), and therefore the fraction of population infectious at the cotype equilibrium reduces below that of the resistant-type equilibrium in the interval (τ0τ1) (Fig. 3A). The dashed curve corresponds to the case where a = 0.5 and q(τ0) = 0.5348 disobeys the condition (4).

Table 1:
Parameter values used for simulations extracted from the published literature (Alexander et al., 2007; Moghadas, 2008; Moghadas et al., 2008; Xiao, Brauer & Moghadas, 2016).
Parameter Baseline values (range) Comments
μ 1/70 per year assumed
R0 2.2, 3 (>1) the value of R0 used for simulations are within the estimated ranges for influenza epidemics and pandemics
δ variable [0–1] varied in simulations
γ 0.2 per day the value of γ corresponds to an average infectious period of five days within the estimated range for influenza infection
τ variable [0–1∕γ] varied in simulations for the length of infectious period
q(τ) variable the probability of developing resistance during treatment was determined from the functional form of q(τ)
DOI: 10.7717/peerj.2817/table-1

We further simulated the cotype equilibrium with the probability of developing resistance that has the functional form q τ = 0 . 6 if τ < τ 0 0 . 06 0 . 6 e 2 τ τ 0 + 0 . 06 if τ τ 0 0 . 06 with τ0 = 2.55. Figure 3B shows the behaviour of the cotype equilibrium for different values of δ, and represents the possibility of an intermediate interval of τ in which the magnitude of infection at the cotype equilibrium is reduced below that of the resistant-type equilibrium (i.e., dotted curve for the case of δ = 0.46).

In the absence of demographics (μ = 0), we simulated the model with parameter values described in Table 1 to determine the final size of the epidemic with different functional forms of q(τ). Figure 4 shows heatmaps for the epidemic final size and the final size of the resistant infections as functions of delay in start of the treatment and transmission fitness of resistance. Clearly, the patterns of epidemic final size are affected by q(τ), which also determines how early the competitive balance is shifted in favour of the resistant-type. The functional forms of q(τ) are shown in Fig. S2. We observed that for a given transmission fitness δ, it is possible to obtain the same epidemic final size with several delays in start of the treatment. For R0 = 2.2 and δ = 0.48, Fig. 5 illustrates the change in the final size of epidemic with the fraction of infected population that is treated. The particular forms of q(τ) (Fig. S2) indicate that if delay in start of treatment exceeds a certain amount of time, the same final size may be achieved with different treatment levels (Fig. 5, black and red curves). For the black curve, the lowest final size is achieved when the treatment starts at the onset of infectious period, which corresponds to the highest fraction of infected individuals being treated. However, using a different functional form of q(τ) (Fig. S2, red curve), we observed that the lowest final size occurs with τ = 2.75 days delay in start of the treatment (within the infectious period), which corresponds to a lower treatment level of the infected population. Importantly, the lowest epidemic final sizes in black and red curves are approximately the same (represented by black and red boxed), but occurs at different delays in start of the treatment and at different fractions of infected population treated.

Discussion

The emergence and spread of drug-resistance have been studied in a number of epidemic and endemic models of infectious diseases (Lipsitch et al., 2007; Moghadas, 2008; Handel, Longini Jr & Antia, 2007; Legros & Bonhoeffer, 2016; Mills, Cohen & Colijn, 2013). A key parameter in these models is the probability of resistance development at the host level, which has been considered to be independent of the time for start of the treatment. In this study, we developed a simple two-subtype disease transmission model to study the population dynamics of drug-resistance spread when the probability of developing resistance is a function of delay in the treatment initiation. We discussed the transient and equilibrium states of the resistant and wild pathogen subtypes in short-term (epidemic) and long-term (endemic) scenarios. In contrast to previous work (D’Agata et al., 2008), our results show that an early initiation of treatment may not necessarily be the optimal strategy for preventing the emergence and spread of drug-resistance (Fig. 5). Our results also show that, depending on the probability of resistance development at the host level and the relative transmission fitness of the resistant-type at the population level, the minimum infection state of the system at equilibrium for both pathogen subtypes may occur with considerable delay in start of the treatment during the infectious period (Fig. 3). This is not necessarily consistent with current treatment practices, in which the management of infection and severe outcomes in patients takes precedence over the possible evolution and spread of drug-resistance under the selection pressure of drugs.

The probability of developing resistance at the host level is affected by several factors, including the frequency with which the resistant mutants arise, compensatory mutations that restore impaired fitness of resistance, the efficacy of drugs in preventing pathogen growth, the concentration of drugs, and the level and strength of the host immune responses (White, 2004). The time for start of the treatment can have a substantial influence on the underlying processes by which these factors are characterized. At the population level, the interplay between these factors manifests itself in the rise and fall of resistance (Lipsitch, 2001), making the identification of optimal treatment regimens extremely difficult to achieve in order to simultaneously minimize the incidence of disease and limit resistance emergence and spread (Colijn & Cohen, 2015; Day & Read, 2016). As has been demonstrated in recent studies, it is therefore important to develop multi-scale models that integrate both within-host infection dynamics and between-host disease transmission (Legros & Bonhoeffer, 2016). Our theoretical framework and analyses here are not meant to address the question of optimal treatment strategies in either the host or population level, but rather to underscore the complexity of drug-resistance dynamics in both the epidemic and endemic disease states.

In the model presented here, we have made several simplifying assumptions. We assumed that those who are effectively treated (without developing resistance) do not contribute to disease transmission following the initiation of treatment. This effect can be adjusted in the model by tuning other parameters such as the infectious period. We assumed the same infectious period for both the resistance-type and wild-type infections. Our distinct theoretical analyses show that similar analytical results can still be achieved, and we therefore chose this simplification. We also assumed that all infectious individuals will receive treatment with delay τ after the onset of infectious period, if they have not recovered until time τ. Realistically, in several diseases (e.g., influenza), a sizeable fraction of infectious individuals may not receive treatment for a number of reasons such as mild illness or asymptomatic infection. While the coverage of treatment can be included in the model as an independent parameter, we note that the effect of reduced treatment coverage can also be adjusted through a longer infectious period. Our model is based on a single treatment regimen which is assumed to be ineffective against resistant-type infection. In some infectious diseases (e.g., influenza, HIV/AIDS, tuberculosis) there are a number of drugs that may be used sequentially or in combination to prevent or manage resistance in treated patients. Although we have not considered multi-drug treatment, the effect of drugs on resistant-type infection may be included in the probability of developing resistance in our model. While these assumption can be relaxed and the model be extended to consider the aforementioned factors, we believe this simple structure underscores the importance of the probability of resistance development with delay in start of the treatment, and clearly demonstrates the complex dynamics of drug-resistance spread in the population. Finally, our results echo the take-home message drawn from some recent studies that the management of drug-resistance requires theoretical frameworks that combine the two scales of micro and macro dynamics to not only minimize the short-term impact of disease on the population, but also to address the long-term epidemiological consequences of the pathogen evolutionary responses under the pressure of drug treatment.

Supplemental Information

Revised Supplementary Material

DOI: 10.7717/peerj.2817/supp-1

Probability of developing resistance

DOI: 10.7717/peerj.2817/supp-2

Probability of developing resistance with delay

DOI: 10.7717/peerj.2817/supp-3

Illustration of parametric curves

DOI: 10.7717/peerj.2817/supp-4

Matlab code

The Matlab code to find the numerical solutions of the full model for a given set of parameters

DOI: 10.7717/peerj.2817/supp-5

Matlab Code for Fig. 3A of the main text

DOI: 10.7717/peerj.2817/supp-6

Matlab Code for Fig. 3B of the main text

DOI: 10.7717/peerj.2817/supp-7

Matlab Code for Fig. 4 of the main text

DOI: 10.7717/peerj.2817/supp-8

Matlab Code for Fig. 5 of the main text

DOI: 10.7717/peerj.2817/supp-9