Population-Wide Emergence of Antiviral Resistance during Pandemic Influenza

Background The emergence of neuraminidase inhibitor resistance has raised concerns about the prudent use of antiviral drugs in response to the next influenza pandemic. While resistant strains may initially emerge with compromised viral fitness, mutations that largely compensate for this impaired fitness can arise. Understanding the extent to which these mutations affect the spread of disease in the population can have important implications for developing pandemic plans. Methodology/Principal Findings By employing a deterministic mathematical model, we investigate possible scenarios for the emergence of population-wide resistance in the presence of antiviral drugs. The results show that if the treatment level (the fraction of clinical infections which receives treatment) is maintained constant during the course of the outbreak, there is an optimal level that minimizes the final size of the pandemic. However, aggressive treatment above the optimal level can substantially promote the spread of highly transmissible resistant mutants and increase the total number of infections. We demonstrate that resistant outbreaks can occur more readily when the spread of disease is further delayed by applying other curtailing measures, even if treatment levels are kept modest. However, by changing treatment levels over the course of the pandemic, it is possible to reduce the final size of the pandemic below the minimum achieved at the optimal constant level. This reduction can occur with low treatment levels during the early stages of the pandemic, followed by a sharp increase in drug-use before the virus becomes widely spread. Conclusions/Significance Our findings suggest that an adaptive antiviral strategy with conservative initial treatment levels, followed by a timely increase in the scale of drug-use, can minimize the final size of a pandemic while preventing large outbreaks of resistant infections.


The Model Structure
We assume that the population is entirely susceptible to the emerging pandemic strain with no preexisting immunity. Let S denote the class of susceptible individuals who may become infected with either wild-type or resistant strains. Denoting the classes of individuals exposed to wild-type viruses by E, resistant strains with low fitness by E r , and resistant strains with high fitness by E rH , we have S ′ (t) = −βQ(t)S(t), where β is the baseline transmission rate of the wild-type strain, 1/µ E represents the mean latent period (assumed to be the same for wild-type and resistant infections), and βQ(t) = β(Q w + Q r + Q rH ) is the force of infection, yet to be formulated. Let p represent the probability of developing clinical disease after the latent period. Then, for corresponding classes of individuals with asymptomatic infections (i.e. those who are infectious without showing clinical symptoms, and therefore are not treated), we obtain 1 Simulations and sensitivity analyses of the model were performed using a solver for delay integro-differential equations (Paul, 1997 (Ferguson et al., 2005(Ferguson et al., , 2006(Ferguson et al., , 2003Halloran et al., 2006;Jefferson et al., 2006;Longini et al., 2004Longini et al., , 2005Regoes & Bonhoeffer, 2006 where 1/µ A is the infectiousness period. To derive the equations for the symptomatic infection, we assume a window of opportunity of two days for start of treatment following the onset of clinical symptoms. Using rates of treatment and emergence of drug-resistance described in Alexander et al. (2007), the corresponding equations for untreated and treated symptomatic infections are given by where n is the size of the window of opportunity for the start of treatment; 1−q represents the populationlevel of treatment, µ U and µ T are the recovery rates of untreated and treated symptomatic infections (during secondary stage), respectively; d U , d U,r , and d U,rH are the corresponding disease-induced mortality rates for untreated symptomatic infections; V represents the fraction of treated individuals which develops drug-resistance with low fitness during the window of opportunity; α T is the rate for developing drug-resistance during the secondary stage of symptomatic infection; and γ U and γ T are the rates of conversion between resistant mutants of untreated and treated symptomatic infections, respectively. To formulate the force of infection, let i U (t, a) and i T (t, a) be the densities of untreated and treated wild-type infections after a time a has elapsed since an exposed individual becomes infectious. Consider-ing the infectious compartments of the wild-type strain, and detailed description presented in Alexander et al. (2007Alexander et al. ( , 2008, we have where τ is the period of pre-symptomatic infection; δ A , δ P , and δ U represent the relative infectiousness of the wild-type strain for asymptomatic, pre-symptomatic, and the secondary stage of symptomatic infection without treatment, respectively; and δ T is the relative infectiousness of a treated clinical case with the wild-type strain. The treatment is assumed to have no effect on individuals infected with resistant strains. With the corresponding notation for resistant strains, we have and where δ r and δ r H represent the relative infectiousness of resistant strains with low and high transmission fitness, respectively. It is assumed that treatment of wild-type infection reduces transmissibility of the virus by 60% (through reduction factor δ T ), but has no effect in transmission of resistant strains. Summarizing, the above represents the model as a system of delay differential equations, where estimates of its parameter values from the published literature are given in Table 1 (see also Table 1 in the main text).

Analysis of Reproduction Numbers
In our analysis, we fixed the initial size S 0 of the susceptible population to compute the control reproduction number R w c when an individual infected with the wild-type strain is introduced into the E-class. We assumed that E(0) = 1, and let E(t) = 0 for t ∈ [−n, 0), and A(0) = I U (0) = I T (0) = 0. Considering the duration and transmission rates associated with asymptomatic, untreated and treated symptomatic infections (see Figure 1 in the main text), the total number of secondary cases generated in the A, I U and I T classes is given by We also calculated the number of new cases generated during the primary stage of symptomatic infection (window of opportunity for effective treatment), which involves the history of the E-class. Noting that q = 1 during pre-symptomatic infection (without treatment), this number is given by and therefore the control reproduction number of the wild-type strain can be expressed as In the absence of antiviral treatment (q ≡ 1 and V ≡ 0), R w c reduces to the basic reproduction number Since treatment has no effect on individuals infected with resistant strains, similar calculations to the above lead to the reproduction numbers of resistant strains with low fitness (R r 0 ) and high fitness (R rH 0 ) as The next generation matrix has the form and therefore the criterion for the control of disease, defined in terms of the spectrum of this matrix (Diekmann & Heesterbeek, 2000), is given by R c = max{R w c , R r 0 , R rH 0 }.

Sensitivity Analysis
To investigate the effect of parameter changes on the results shown by simulations using baseline values, we performed a sensitivity analysis by considering a sampling approach that allows for the simultaneous variations of the basic reproduction number, the rates of de novo resistant mutations, and the rates of conversion between resistant strains. Using the Latin Hypercube Sampling (LHS) technique (McKay et al., 1979), we generated samples of size n=100 in which each parameter is treated as a random variable and assigned a probability function. In this technique, the parameters are uniformly distributed and sampled within their respective ranges. The reproduction number R 0 was uniformly sampled from the range [1.4, 2], which includes the estimated ranges of reproduction numbers for the 1918, 1957, and 1968pandemics (Gani et al., 2005Viboud et al., 2006). The rates of de novo resistant mutations (ρ max , α T ) were sampled from the range [0.018, 0.072] (Regoes & Bonhoeffer, 2006;Débarre et al., 2007), corresponding to 5.8% − 17% incidence of resistance, which lies within the estimated range of neuraminidase resistance reported in clinical studies (Kiso et al., 2004;Ward et al., 2005;Yen et al., 2005). The corresponding ranges for the conversion rates of resistant strains (γ U , γ T ) was computed using the constraint that the fraction of treated individuals hosting resistance, which undergoes compensatory mutations and subsequently generates resistant strains with high fitness, lies between 1/5000 and 1/500 (Lipsitch et al., 2007). The baseline value of 1/2000 was used for simulations. Furthermore, we assumed that compensatory mutations are less likely to occur in the absence of treatment, and considered γ U to be 10−fold smaller than γ T . The values of other parameters are given in Table 1.
For the sensitivity analyses, we introduced the parameter T a to represent the "minimum final size" of the pandemic within an adaptive treatment strategy, when the initial treatment level is increased to 90% at a time t * . For each set of parameter values in the sample, we then computed the ratio T a /T c as a function of δ rH (the relative transmission of the resistant strain with high fitness), where T c is the final size of the pandemic when treatment is maintained constant at the corresponding optimal level (below 90%) at all times during the outbreak. The results of sensitivity analyses are illustrated in Figures 1a,  1b, 1c, when initial treatment levels in the adaptive treatment strategy are assumed to be 0%, 25%, and 50%, respectively. These figures indicate that for low values of δ rH (below ∼ 0.8), the risk of a resistant epidemic developing is small, and both treatment strategies are comparable in their effectiveness. However, as δ rH increases (above ∼ 0.8), self-sustaining epidemics of resistant viruses can be established, and the benefit of an adaptive treatment strategy becomes more pronounced. Using the above sample, we also projected the corresponding ranges of time t * as a critical parameter in this strategy. The results, depicted in Figures 2a, 2b, 2c, suggest that aggressive treatment should be further delayed (following the onset of the outbreak) for higher initial treatment levels, should resistant strains with high transmission fitness (above ∼ 0.8) emerge. Such high treatment levels decelerate the spread of the wild-type virus in the population, and therefore extend the time required for a sufficient drop in the number of susceptible individuals to prevent resistant outbreaks.  Figure 2: Sensitivity analyses showing box plots for the variations in the optimal transition time t * corresponding to the minimum total number of infections, as a function of δ rH , with other parameters sampled from their respective ranges, as described in the text. The solid curve passes through the median values of data points for t * , and each box contains 50% of data points between the first and third quartiles of the sampling distribution. The remaining 50% of data points are represented by whiskers. Initial treatment levels in adaptive antiviral strategy before transition time t * are: (a) 0%; (b) 25%; and (c) 50%.

The Effect of Delay in Start of Treatment
Not only the population level of drug use, but also early onset of treatment of indexed cases within the window of opportunity can significantly influence the outcome of an antiviral strategy. To demonstrate this, we compared the final size of clinical infections in three scenarios with different delays in the start of treatment following the onset of clinical disease. Figure 3a (dotted curve) shows that early treatment with 0.5-day delay results in smaller number of clinical infections (and therefore the minimum epidemic size is feasible with a lower level of drug use) than when treatment is initiated with 1-day delay (solid curve). This is mostly due to a greater reduction in transmission of the wild-type infection (Figure 3b). A more rapid decline in the final size of wild-type infections occurs when compensated mutants become the driving force for disease progression with increasing level of treatment. However, a more dramatic increase in the number of resistant infections is observed for higher levels of drug use with less delay in start of treatment (Figure 3a, dotted curve). Initiating treatment with a longer delay of 1.5 days has little impact on suppressing wild-type infection, even with high levels of treatment (Figure 3b, dashed curve). However, in this case, resistance emergence is limited due to the wide spread of the wild-type virus, thereby rapidly depleting the pool of susceptible hosts.