Population dynamics of epidemic and endemic states of drug-resistance emergence in infectious diseases

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 I_w and I_r, respectively, the dynamics of disease propagation is governed by the time-dependent differential equations system: (1)${\displaystyle \begin{array}{c}{S}^{\prime }\left(t\right)=\mu -\beta S\left(t\right)\left[I_w\left(t\right)+\delta I_r\left(t\right)\right]-\mu S\left(t\right),\hfill \\ I_w^{\prime }\left(t\right)=\beta S\left(t\right)I_w\left(t\right)-\beta S\left(t-\tau \right)I_w\left(t-\tau \right)e^{-\left(\mu +\gamma \right)\tau }-\gamma I_w\left(t\right)-\mu I_w\left(t\right),\hfill \\ I_r^{\prime }\left(t\right)=\delta \beta S\left(t\right)I_r\left(t\right)+q\left(\tau \right)\beta S\left(t-\tau \right)I_w\left(t-\tau \right)e^{-\left(\mu +\gamma \right)\tau }-\gamma I_r\left(t\right)-\mu I_r\left(t\right),\hfill \\ T^{\prime }\left(t\right)=\left(1-q\left(\tau \right)\right)\beta S\left(t-\tau \right)I_w\left(t-\tau \right)e^{-\left(\mu +\gamma \right)\tau }+\gamma \left[I_w\left(t\right)+I_r\left(t\right)\right]-\mu T\left(t\right).\hfill \end{array}}$ 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)I_w(t) gives the cohort of individuals who enter the class I_w 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 I_w class τ units of time after entering the class is e^−(μ+γ)τ. Thus, the transition out of the I_w class by means of treatment at any time t can be described by the term βS(t − τ)I_w(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 I_r 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 I_r (Fig. 1). Lastly, the class T represents all individuals who are recovered from infection or effectively treated.

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 I_w 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.

Equilibria

To characterize the equilibria of the model, we define R₀ = β∕(μ + γ), and apply the next generation method to express the effective reproduction number (R_eff) as the spectral radius of the matrix ${\displaystyle \left(\begin{array}{cc}{R}_1\hfill & 0\hfill \\ \delta R_2\hfill & R_3\hfill \end{array}\right),}$ where ${\displaystyle R_1=R_0\left(1-e^{-\left(\mu +\gamma \right)\tau }\right),R_2=q\left(\tau \right)R_0{e}^{-\left(\mu +\gamma \right)\tau },R_3=\delta R_0.}$

With this notation, R_eff = max{R₁, R₃}. Each of these quantities provide important information on the spread of disease by each pathogen subtype. Specifically, R₁ represents the number of secondary cases of the wild-type infection generated by an infectious case of the wild-type before developing resistance; δR₂ represents the number of secondary cases of the resistant-type infection generated by an infectious case of the wild-type after developing resistance; and R₃ is the number of secondary cases of the resistant-type infection generated by an infectious case of the resistant-type. In this formulation, R₁ approaches R₀ for large τ, and R₀ 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 E₀ = (0, 0) exists unconditionally. The resistant-type equilibrium $E_r=\left(0,\widehat{I_r}\right)$, at which the wild-type infection is absent, exists only if R₃ > 1. The cotype equilibrium $E^{\ast }=\left(I_w^{\ast },I_r^{\ast }\right)$ exists whenever both conditions R₁ > 1 and R₁ > R₃ hold.

Delay-dependent analysis

We focus our analysis on the equilibria when R₀ > 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 R₃ ≥ R₀), E_r 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 (2)${\displaystyle \widehat{I_r}=\frac{\mu }{\mu +\gamma }\left(1-\frac{1}{R_3}\right).}$ If the resistant-type emerges with a negative fitness advantage (i.e., δ < 1 or equivalently R₃ < R₀), then there is a critical value τ₀ =  − log(1 − δ)∕(μ + γ) at which R₁(τ₀) = R₃ (Fig. 2). For τ < τ₀ (R₃ > R₁ > 1), similar to the case of δ ≥ 1, E_r remains the only infection equilibrium of the model.

At τ = τ₀, the cotype equilibrium E^∗ emerges and exists along with the resistant-type equilibrium E_r for τ ≥ τ₀. We define the fraction of population infectious at the cotype equilibrium as a function of τ, and let $G\left(\tau \right)=I_w^{\ast }\left(\tau \right)+I_r^{\ast }\left(\tau \right)$ 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 τ₀. Taking the implicit derivative of G(τ) at τ₀ gives (3)${\displaystyle G^{\prime }\left(\tau \right)=\frac{\mathrm{d}}{\mathrm{d}\tau }\left(I_w^{\ast }\left(\tau \right)+I_r^{\ast }\left(\tau \right)\right)=\frac{\mathrm{d}}{\mathrm{d}\tau }\left(\frac{R_1-R_3}{R_2}\right)I_r^{\ast }\left(\tau \right)+\left(1+\frac{R_1-R_3}{R_2}\right)\frac{\mathrm{d}}{\mathrm{d}\tau }{I}_r^{\ast }\left(\tau \right)=\left(\frac{R_1^{\prime }{R}_2-\left(R_1-R_3\right)R_2^{\prime }}{R_2^2}\right)I_r^{\ast }\left(\tau \right)+\left(1+\frac{R_1-R_3}{R_2}\right)\frac{\mathrm{d}}{\mathrm{d}\tau }{I}_r^{\ast }\left(\tau \right).}$

Since $R_1^{\prime }\left({\tau }_0\right)∕R_2\left({\tau }_0\right)=\left(\mu +\gamma \right)∕q\left({\tau }_0\right)$, we have ${\displaystyle G^{\prime }\left({\tau }_0\right)=\frac{\left(\mu +\gamma \right)}{q\left({\tau }_0\right)}{I}_r^{\ast }\left({\tau }_0\right)+\frac{\mathrm{d}}{\mathrm{d}\tau }{I}_r^{\ast }\left({\tau }_0\right).}$ Taking the derivative of infection component with the resistant-type at τ₀ (see Supplemental Information) gives ${\displaystyle G^{\prime }\left({\tau }_0\right)=\frac{\mu }{q\left({\tau }_0\right)}\left(\frac{1}{\delta }-1\right)\left(q\left({\tau }_0\right)-1+\frac{1}{R_3}\right).}$ This implies that the fraction of population infectious at E^∗ decreases if (4)${\displaystyle q\left({\tau }_0\right)<1-\frac{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 I_w and I_r classes. When (4) holds, the start of the treatment with delay τ > τ₀ reduces the fraction of population infectious at the cotype equilibrium. However, for sufficiently large τ, G(τ) increases and saturates at ${\displaystyle G_{\infty }=\underset{\tau \to \infty }{lim}G\left(\tau \right)=\frac{\mu }{\mu +\gamma }\left(1-\frac{1}{R_0}\right).}$ 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_{\infty }>\widehat{I_r}$. This suggests that, given (4), there is an intermediate interval for delay in start of the treatment after τ₀, in which the magnitude of infection at the cotype equilibrium reduces below that of the resistant-type equilibrium (Fig. 3A).

When q(τ₀) > 1 − 1∕R₃, the magnitude of infection at the cotype equilibrium increases for τ > τ₀ and sufficiently close to τ₀. The continual increase in the population level of infection for τ > τ₀ 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 τ₀. 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) ${\displaystyle F=\frac{\gamma \left(1-q\left(\tau \right)e^{-\gamma \tau }\right)}{1-e^{-\gamma \tau }}{\int }_0^{\infty }{I}_w\left(t\right)\mathrm{d}t+\gamma {\int }_0^{\infty }{I}_r\left(t\right)\mathrm{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.