Given the danger of an unprecedented spread of the highly pathogenic avian influenza strain H5N1 in humans, and great challenges to the development of an effective influenza vaccine, antiviral drugs will probably play a pivotal role in combating a novel pandemic strain. A critical limitation to the use of these drugs is the evolution of highly transmissible drug-resistant viral mutants. Here, we develop a mathematical model to evaluate the potential impact of an antiviral treatment strategy on the emergence of drug resistance and containment of a pandemic. The results show that elimination of the wild-type strain depends crucially on both the early onset of treatment in indexed cases and population-level treatment. Given the probable delay of 0.5–1 day in seeking healthcare and therefore initiating therapy, the findings indicate that a single strategy of antiviral treatment will be unsuccessful at controlling the spread of disease if the reproduction number of the wild-type strain exceeds 1.4. We demonstrate the possible occurrence of a self-sustaining epidemic of resistant strain, in terms of its transmission fitness relative to the wild-type, and the reproduction number . Considering reproduction numbers estimated for the past three pandemics, the findings suggest that an uncontrollable pandemic is likely to occur if resistant viruses with relative transmission fitness above 0.4 emerge. While an antiviral strategy is crucial for containing a pandemic, its effectiveness depends critically on timely and strategic use of drugs.
The threat of an impending influenza pandemic, possibly through mutation of the present deadly avian strain H5N1, has galvanized global efforts to understand the potential benefits and limitations of mitigation strategies. Previous modelling studies have considered pharmaceutical and non-pharmaceutical interventions (Ferguson et al. 2003, 2005, 2006; Longini et al. 2004, 2005; Gani et al. 2005; Germann et al. 2006) and rationalized the use of antiviral drugs as the first-line defence against a new pandemic strain. The effects of these drugs are twofold: (i) they reduce the infectivity and duration of infectiousness by inhibiting virus replication and (ii) reduce susceptibility; these will in turn decelerate the spread of infection in the population to afford time for development of new vaccine candidates.
There are two groups of antiviral drugs available for treatment and prophylaxis of influenza: M2 inhibitors (amantadine and rimantadine) and neuraminidase inhibitors (oseltamivir and zanamivir). Despite the effectiveness of these drugs in reducing influenza-related morbidity and mortality, the emergence of drug resistance poses a critical limitation on their application. Incidence of viral resistance to M2 inhibitors has been associated with an increasing rate in seasonal influenza, possibly through widespread or indiscriminate use of the drugs (Bright et al. 2005). Neuraminidase inhibitors are less prone to selecting for resistant mutations (Moscona 2005; Regoes & Bonhoeffer 2006), and therefore offer a better option for pandemic preparedness. However, recent emergence of oseltamivir resistance has raised concerns about our strengths in facing an influenza pandemic (Kiso et al. 2004; de Jong et al. 2005; Moscona 2005; Regoes & Bonhoeffer 2006).
The strategy of antiviral therapy raises a number of public health concerns regarding the optimal use of drugs for treatment, prophylaxis, or combination thereof, in order to not only minimize the short-term impact of the virus on the population, but also account for the longer-term consequences of the evolutionary responses of the virus. This is particularly important for preventing pandemic waves of infection caused by the emergence of resistant viral mutants. These concerns can be addressed when appropriate models of evolutionary epidemiological aspects of the disease are employed. Despite several recent modelling efforts (Ferguson et al. 2003, 2005, 2006; Longini et al. 2004, 2005; Gani et al. 2005; Germann et al. 2006; Wu et al. 2006; Colizza et al. 2007), the interplay between these aspects and its consequences for containment of a pandemic are poorly exploited.
This study undertakes to evaluate the merit of the application of antiviral drugs from a different angle, through a modelling approach that provides a link between viral dynamics at the individual level and disease spread in the population. Central to our model is the inclusion of infectious compartments according to the stage progression of the disease. This allows us to monitor the density of infected individuals in terms of the time elapsed since the onset of clinical disease, explicitly as an independent structure variable. To capture the dynamics of the emergence of drug resistance within the limited window of opportunity for commencing antiviral therapy, we consider an evolutionary rate of emergence of resistant mutants that increases linearly with the outgrowth of viral replication. Our approach introduces a systematic way to account for the effect of delay in initiating treatment of indexed cases on the emergence of drug resistance and spread of the disease. We derive a criterion for the control of influenza infection and demonstrate the possible scenarios of disease outbreak, including the possibility of a self-sustaining epidemic of resistant viruses. Although this study does not address the influence of immunological/epidemiological characteristics of the individuals/population, these signatures remain crucial for gaining insights into understanding the complexity of transmission routes and identifying means by which the invasion of a pathogen in the population can be effectively contained. In the following, we detail the model structure, discuss our results and their epidemiological implications, and place them in the context of public health.
2. Model structure
The development of the model, based on clinical aspects of influenza infection, involves the processes of disease progression at both the individual and population levels. At the population level, we consider classes of susceptible, exposed (infected but not infectious), and both asymptomatic and symptomatic infected individuals. It is assumed that the transmission of infection occurs through contacts between susceptibles and infected individuals, where for simplicity, mass action incidence is used (Arino et al. 2006; Gardam et al. 2007). Denoting the class of susceptibles by S, and classes of exposed individuals with the wild-type and resistant viruses by E and Er, respectively, we have(2.1)where β is the baseline transmission rate of the wild-type strain; 1/μE represents the mean latent period (assumed to be the same for E and Er classes); and βQ(t)=β(Qs+Qr) is the force of infection, yet to be formulated, in which Qs and Qr comprise the infectious compartments of the wild-type and resistant strains, respectively. Exposed individuals in the E and Er classes cannot transmit the disease in the latent period, during which viral titres increase to detectable and transmissible levels (Baccam et al. 2006). Also, the inclusion of these classes more realistically represents the epidemiology of disease, due to the fact that treatment of infected individuals is not feasible until after the latent period has elapsed.
After the latent period, exposed individuals either develop clinical disease or undergo an asymptomatic phase without showing symptoms for the entire course of infection. Let A and Ar denote the classes of asymptomatic individuals corresponding to the wild-type and resistant strains, who can shed the virus during their infectious period 1/μA, respectively. Assuming a probability p∈(0, 1) for an exposed individual to develop clinical disease, we obtain(2.2)
To formulate the equations for symptomatic infectious classes, we need to consider the progression of disease at the individual level (figure 1). The clinical course of infection may be divided into three stages: (i) pre-symptomatic infection, during which transmission can occur before symptoms appear, (ii) primary stage of symptomatic infection after the onset of symptoms, which also represents the window of opportunity for initiating an effective course of antiviral therapy, and (iii) secondary stage of symptomatic infection following the window of opportunity. Furthermore, we assume that clinical cases who have not started antiviral therapy within this window will progress to the secondary stage without receiving any treatment.
Although antiviral drugs have been shown to reduce both the infectivity and the period of infectiousness (Nelson et al. 2004), their overall effectiveness may be limited by the emergence and spread of drug-resistant viral mutants (Moscona 2005; Regoes & Bonhoeffer 2006; Lipsitch et al. 2007). The likelihood of an emergent resistant mutant out-competing a wild-type virus and dominating the viral population depends greatly on its replicative fitness relative to that of the wild-type strain. In the absence of treatment, such dominance is improbable, due to a sufficient difference in the intrinsic fitness between resistant mutants and wild-types. However, treatment can result in a significant reduction in fitness of the wild-type, which may in turn lead to the resistant viral mutants prevailing (Handel et al. 2006). Without pre-existing resistant mutants, the inhibition of viral replication by the use of antiviral drugs will minimize the likelihood of evolving such mutants during a course of therapy. On the other hand, the probability of emergence of resistant mutants increases with the outgrowth of viral replication, and therefore the development of drug resistance becomes more probable with increasing delay in the onset of therapy.
To model the evolution of drug resistance, let r(a) represent the treatment rate at elapsed time a since an exposed individual becomes infectious, and iU(t, a) and iT(t, a) be the densities, with respect to a, of untreated and treated infected individuals, respectively, with the wild-type strain at current time t (figure 2). We define the rate ρ(a)=κθ(a) for the emergence of drug resistance within the window of opportunity (Stilianakis et al. 1998; Ferguson et al. 2003), where κ is the rate at which resistant mutants evolve during treatment and θ(a) is the likelihood of developing drug resistance with delay a in the onset of therapy. Let iU,r(t, a) and iT,r(t, a) denote the densities of untreated and treated individuals with resistant strain, respectively (figure 2). Then, with n representing the size of the window of opportunity for initiating therapy, we have (Webb 1985; Metz & Diekmann 1986), for a∈[0,n](2.3)(2.4)(2.5)(2.6)subject to the following boundary conditions
Let IU(t) and IT(t) denote the total number of untreated and treated individuals with the wild-type strain when a≥n, respectively. Thus, we obtain(2.11)(2.12)where 1/μU and 1/μT are the mean infectious periods of symptomatic infection (secondary stage) and dU and dT represent wild-type disease-induced mortality rates of untreated and treated individuals, respectively. The parameter αT denotes the rate at which infected individuals under treatment develop drug resistance during the secondary stage of symptomatic infection (Stilianakis et al. 1998; Ferguson et al. 2003). This rate is relatively high compared with ρ at the early stages following the onset of symptoms, as resistant mutants in viruses isolated from patients treated with oseltamivir were mostly detected 3 days after start of treatment (Kiso et al. 2004).
To express the equations for the total number of untreated and treated individuals with resistant viral strain for a≥n, denoted by IU,r(t) and IT,r(t), respectively, we need to solve equations (2.5) and (2.6). Substituting the solution of (2.5) into (2.6) and solving, we obtain(2.13)(2.14)for t≥a. These solutions at a=n can be used to derive(2.15)(2.16)where dU,r represents the disease-induced mortality rate of the resistant strain.
In order to gain a better understanding of the parameter v, we considered a simple scenario in which a constant treatment coverage q is assumed as a result of a pulse treatment rate. Thus,and therefore equation (2.12) reduces towhere E(t)=0 for t∈[−n,0]. Taking into account the second term in equation (2.14), the parameter v may be interpreted as the fraction of infected individuals undergoing therapy without developing drug resistance within the window of opportunity.
We are now in a position to define the force of infection βQ=β(Qs+Qr). Considering the infectious compartments with the wild-type strain, we have(2.17)where δA, δP and δU represent the relative transmission rates of the wild-type strain for asymptomatic, pre-symptomatic and the secondary stage of symptomatic infection without treatment, respectively. The parameter δT is the relative infectiousness of a treated clinical case with the wild-type strain. Assuming that the treatment has no effect in reducing the infectivity of drug-resistant cases (Stilianakis et al. 1998), and denoting the relative transmission fitness of the resistant strain by δr, the expression for Qr is given by(2.18)
Summarizing, the above derivations build our model in the form of a system of delay differential equations (see electronic supplementary material), where estimates of its parameter values from the published literature are given in table 1.
3. The control reproduction number
We derived the control reproduction number (Rc) in order to evaluate the potential impact of an antiviral strategy during a pandemic. In the absence of treatment, this quantity reduces to the basic reproduction number, defined as the number of new infectious cases generated by a single infected individual introduced into a wholly susceptible population during the course of infection (Diekmann & Heesterbeek 2000; van den Driessche & Watmough 2002).
To compute Rc, we assumed that the population is entirely susceptible to the emerging pandemic strain with no pre-existing immunity. Taking into account that treatment is feasible after the onset of symptoms, we first calculated the average number of drug-sensitive cases () emanating from the introduction of an infectious individual with the wild-type strain (see electronic supplementary material). In the absence of therapy (q≡1), reduces to the basic reproduction number of the wild-type strain given by(3.1)where S0 is the initial size of the susceptible population and other parameters are defined in table 1.
An individual infected with the wild-type virus who has started a course of therapy within the window of opportunity may develop drug resistance, and therefore generate a number of new resistant cases (; see electronic supplementary material). We also considered the scenario in which an individual exposed to a resistant viral strain is introduced into the population, and assumed that treatment does not reduce viral replication in individuals who shed resistant viruses (Stilianakis et al. 1998). The basic reproduction number of the resistant strain is given by(3.2)which determines whether an outbreak of resistant cases due to direct transmission occurs. We then calculated as the control reproduction number and showed that the disease outbreak of both wild-type and resistant strains can be prevented whenever Rc<1 (see electronic supplementary material).
In the absence of control measures, the final size relation (Arino et al. 2006; Gardam et al. 2007)(3.3)can be used to estimate based on the given clinical attack rate of the wild-type strain, defined as the fraction of the initial susceptible population S0 which develops clinical disease, and is given by p(1−S∞/S0), where S∞ is the size of susceptible population when the epidemic dies out. Then, equation (3.1) provides a range of values for the parameter βS0, and therefore the baseline transmission rate of the wild-type virus (see electronic supplementary material).
4. Rates of treatment and evolution of drug resistance
Influenza A resistant viruses can emerge during treatment, but they have only rarely been isolated from specimens collected as part of routine influenza surveillance, reporting less than 1% development of drug resistance from untreated infected cases (Ziegler et al. 1999). We therefore considered the scenario in which drug resistance may develop as a result of treatment.
To investigate the feasibility of containing a pandemic using antiviral agents, we prescribed a profile of treatment rate that decreases linearly towards the end of the window of opportunity. This corresponds to a strategy that gives priority to those who are clinically diagnosed early in this window by initiating a course of therapy immediately (figure 3a). Such priority is justified by two major epidemiological considerations: (i) early treatment reduces the level/duration of infectiousness, and therefore transmissibility of the virus, which affords a greater opportunity for containing the wild-type strain and (ii) it limits viral replication, and therefore inhibits formation of resistant mutants, which is crucial for preventing the emergence of drug resistance.
Owing to probable delay in seeking healthcare, we assumed that treatment commences with a time lag a0 after the onset of clinical symptoms and defined(4.1)where b is the slope of the treatment rate . Thus, those who have not started therapy within the window of opportunity will progress to the secondary stage of symptomatic infection without receiving treatment. Letting c denote the treatment level of infected cases so that q(n)=1−c, and using equation (2.7), we find b=rmax/(n−a0), where rmax is the maximum rate of treatment given by
The functional form of q(a; a0) with delay may then be expressed as
With the initiation of therapy, the population is threatened by the emergence of drug-resistant cases and transmission of resistant viruses. To capture the dynamics of emergence of resistant viral mutants, we considered a linearly increasing rate ρ(a) with delay in the onset of therapy within the window of opportunity (figure 3b), and defined bywhere ρmax is the maximum rate of emergence of drug resistance, corresponding to a high level of viral titre (Baccam et al. 2006). The linearly increasing rate ρ(a) provides a good approximation to the multi-cycle replication of neuraminidase inhibitor-resistant mutants (E116D, E116A, E116G, E116V, R149K, R291K and R292K) observed in clinical investigations (Jackson et al. 2005; Yen et al. 2005).
Experimental studies of influenza A infection provide a range of 1–18% incidence of neuraminidase resistance during a standard 5-day treatment with oseltamivir (Kiso et al. 2004; Yen et al. 2005). Taking into account the corresponding estimates of neuraminidase mutation rate in treated infected individuals, we set ρmax=κθ(n)=0.036 per day, and assumed the same rate of drug-resistant mutation during the secondary stage of symptomatic infection (Regoes & Bonhoeffer 2006).
5. Numerical experiments and results
In this section, we determine the feasible region for disease control under the prescribed profile of treatment, using parameter values given in table 1. We evaluated the impact of three important factors on the control reproduction number: (i) the level of treatment (c), (ii) the delay in initiating a course of therapy (a0), and (iii) the development of drug resistance during treatment.
Figure 4 illustrates boundaries of the domain , within which the wild-type strain is eliminated for different reproduction numbers. For these illustrations, we assumed that the resistant virus is fivefold less transmissible than that of the wild-type (Stilianakis et al. 1998), though neuraminidase inhibitor-resistant mutants may differ substantially in their fitness and transmissibility (Yen et al. 2005). While R292K and H274Y neuraminidase mutations are recognized to have compromised growth and transmissibility, the E119V mutant virus may exhibit a comparable fitness to that of the wild-type virus (Yen et al. 2005). With δr=0.2, the basic reproduction number of the resistant virus varies between 0.28 and 0.4, corresponding to the range 1.4–2 of . In this case, since (see electronic supplementary material), the occurrence of an outbreak can be prevented by reducing below unity. However, a substantial improvement in the level of treatment is required as the delay in initiating the course of treatment increases. This leads to a more restricted region of , with increasingly stringent requirements of treatment level as the transmissibility (and therefore ) of the wild-type strain increases. This is evident from figure 5 which illustrates the contour plots for , corresponding to , for the range of . Below each contour, for a particular value of , an outbreak of the wild-type strain occurs (), while above the contour wild-type infection can be contained ().
To demonstrate the effect of treatment on the time courses of wild-type infections, we simulated the model with a single infected case introduced into the susceptible population of size S0=100 000, and treatment of indexed cases beginning 1 day after the onset of clinical disease. For the particular values and , figure 6 clearly indicates a substantial reduction in the magnitude of disease outbreak due to treatment. Assuming that 60% of exposed individuals develop clinical disease (see electronic supplementary material), the clinical attack rate of the wild-type strain is reduced from 38% (without treatment) to 25% (with 60% treatment coverage) when . However, the mitigation impact of treatment is less pronounced in reducing the clinical attack rate (from 48 to 39%) when . We also evaluated the effect of transmission fitness of the resistant strain on the number of clinical infections generated through direct transmission of resistant viruses. With and δr=0.2 and 0.6, numerical results show that IU,r does not exceed one at any time during the entire course of the epidemic (figure 7). While this holds true also for and δr=0.2, a small outbreak of resistant infections (with a total number of 93 asymptomatic and 140 clinical infections) occurs when δr increases to 0.6, for which .
Figure 8 illustrates the threshold curve , above which an outbreak of resistant cases due to direct transmission occurs (shaded area). As is evident, the emergence of resistant viruses for which is more probable with higher reproduction numbers of the wild-type strain. This has important implications for pandemic plans involving the application of antiviral drugs, in order to not only prevent the development of drug resistance, but also block transmission of resistant viruses using appropriate control measures. If a pandemic were to occur with the reproduction number within the estimated ranges of the 1918, 1957 and 1968 pandemics (Mills et al. 2004; Gani et al. 2005; Viboud et al. 2006), the results indicate that the emergence of resistant mutants with a relative transmission fitness above 0.4 (for which ) could potentially lead to an uncontrollable pandemic (figure 8). While resistant mutants may emerge with lower transmission fitness during treatment (), the possible compensatory mutations may improve their relative fitness sufficiently to cause (Handel et al. 2006).
Antiviral therapy has been rationalized as the primary public health mitigation strategy in response to a newly emergent pandemic virus. Although this frontline defence is crucial in the absence of a proven vaccine (Ferguson et al. 2003, 2005, 2006; Longini et al. 2004, 2005; Gani et al. 2005; Germann et al. 2006), it can potentially lead to the emergence of drug-resistant viral strains by suppressing replication of the wild-type strain. Drug resistance may in turn limit the effective application of antiviral agents through fixation of mutants in a treated patient and transmission of these mutant viruses in the population as a whole. We incorporated these two mechanisms in a delay differential epidemic model to address the impact of three important factors influencing an antiviral strategy: (i) the delay in the onset of therapy, (ii) the population level of treatment required for the control of wild-type strain, and (iii) the emergence and spread of drug resistance at the population level.
Our results show that control of the wild-type strain is significantly affected by delay in initiating a course of therapy, regardless of the evolutionary characteristics of the virus. For sufficiently low transmission fitness of resistant viruses (), our simulations, over a range of values of the reproduction number of the wild-type strain (Longini et al. 2004, 2005; Gani et al. 2005), reveal that the containment of a pandemic is more likely to be achieved, and with the less stringent requirements on the population level of treatment, the shorter the delay in onset of therapy (figure 4). This underscores the importance of early detection of clinical infections through rapid implementation of reliable diagnostic tests for influenza cases. To emphasize this point, we performed a comparative analysis of different antiviral strategies in which the treatment rate is either constant or linearly increasing towards the end of the window of opportunity (see electronic supplementary material). Such analysis clearly demonstrates that the mitigation impact of the treatment rate discussed above is considerably higher than other treatment rates, offering a greater region for possible containment of pandemic in terms of timing and level of treatment (fig. 4 in electronic supplementary material). Given the probable delay of 0.5–1 day in seeking healthcare and therefore commencing treatment, the findings indicate that a single strategy of antiviral treatment would be unsuccessful at controlling the spread of disease if the transmissibility of the wild-type strain results in above 1.4 (figure 5).
While the early application of antiviral drugs is crucial for reducing the pandemic burden, extreme caution is required to prevent the emergence of drug-resistant viral mutants. The evolution of host–pathogen systems occurs on a short time-scale, due in part to the short generation times and rapid adaptation of viral strains under strong immunological pressure. Resistant strains with sufficiently small reproductive ratio are soon out-competed, and effective treatment may therefore result in disease elimination. However, it is possible that a resistant strain will undergo further mutations that compensate for the replication fitness comparable to that of the wild-type strain. Under these evolutionary changes, even a low level of treatment could lead to an outbreak of resistant cases (Handel et al. 2006), through a sufficient increase in transmission fitness of resistant viruses for which the reproduction number exceeds its threshold (). In our model, this is reflected in the relative transmission fitness of the resistant strain (δr), which determines whether it goes extinct or causes an uncontrollable pandemic.
The model discussed here incorporates the use of antiviral drugs as a single strategy which involves treatment of only clinically diagnosed cases. In the context of drug resistance, our study may be extended to the application of antiviral prophylaxis, which reduces the susceptibility to the wild-type strain (Regoes & Bonhoeffer 2006; Lipsitch et al. 2007), and therefore also the competitive interference between drug-sensitive and drug-resistant viruses. For resistant mutants with sufficiently high viral fitness, antiviral prophylaxis may result in a dramatic increase in the number of resistant cases, thereby establishing a self-sustaining epidemic of viral resistance. In previous work (Gardam et al. 2007), we have shown that in the absence of resistant mutants, targeted prophylaxis of healthcare workers within an antiviral strategy can significantly reduce morbidity and mortality of the general population. While maintaining a healthcare work force in place is critical, this strategy may enhance the spread of drug-resistant viruses through intense contacts with patients and colleagues in the healthcare setting and contacts outside of the workplace.
The modelling efforts in this paper are based on the homogenous mixing assumption, which may result in high sensitivity of the epidemic size to the particular conditions at the early stages of an outbreak, during which stochastic effects play a dominant role in disease propagation. However, simulating large-scale outbreaks may require the development of models based on detailed mobility patterns of interconnected heterogeneous subpopulations (Regoes & Bonhoeffer 2006). For the results reported here, we assumed that the use of antiviral drugs would be as effective as their application in seasonal influenza. We also used estimates for the parameters associated with neuraminidase inhibitor resistance in treated patients. Although the model has not addressed the choice of drugs, recent studies indicate a significantly lower incidence of resistance following treatment with neuraminidase inhibitors than with M2 inhibitors (Kiso et al. 2004; Yen et al. 2005). We should, however, emphasize that the incidence of resistance, and therefore the number of emergent cases (), depends greatly on several factors including the duration of treatment, the rate of treatment, the delay in onset of therapy and the rate at which de novo resistant mutations occur. Finally, our model assumed that the supply of antiviral drugs is secured for the entire course of a pandemic. In the case of scarce resources, the treatment rate would need to be modified to take into account depletion of antiviral drug supply, which would in turn affect the population level of treatment, and therefore the outcome of an antiviral strategy. This merits further investigation to determine an optimal treatment rate that minimizes the emergence of drug resistance while maximizing the benefit of drugs at the population level. The effectiveness and optimal use of drugs will then be significantly affected by limited supply, limited production capacity and the surge in demand for treatment. These constitute important public health concerns and should be integrated in devising any antiviral strategy.
This research was partially supported by the Natural Sciences and Engineering Research Council of Canada (NSERC), Canada Research Chairs Program (CRCP), Mathematics of Information Technology and Complex Systems (MITACS), Emergency Management Unit of the Ontario MOHLTC and the Hungarian Scientific Research Fund, OTKA, Grant no. T049516. The authors would like to thank Theresa Tam and Beni Sahai for their critical comments on the preliminary version of this paper that was presented during the ‘Mathematical Modelling Workshop’ for transmission and control of drug-resistant strains in HIV and influenza, organized by the Public Health Agency of Canada in Ottawa, on 15–16 January 2007. The authors also would like to thank reviewers for their comments that have greatly improved the paper.