key: cord-0005335-8pqf7dnm
authors: Alexander, Murray E.; Moghadas, Seyed M.; Röst, Gergely; Wu, Jianhong
title: A Delay Differential Model for Pandemic Influenza with Antiviral Treatment
date: 2007-08-16
journal: Bull Math Biol
DOI: 10.1007/s11538-007-9257-2
sha: 03c9c4d06fe6dab0aeb8c53321bfd755bb42c4a8
doc_id: 5335
cord_uid: 8pqf7dnm

The use of antiviral drugs has been recognized as the primary public health strategy for mitigating the severity of a new influenza pandemic strain. However, the success of this strategy requires the prompt onset of therapy within 48 hours of the appearance of clinical symptoms. This requirement may be captured by a compartmental model that monitors the density of infected individuals in terms of the time elapsed since the onset of symptoms. We show that such a model can be expressed by a system of delay differential equations with both discrete and distributed delays. The model is analyzed to derive the criterion for disease control based on two critical factors: (i) the profile of treatment rate; and (ii) the level of treatment as a function of time lag in commencing therapy. Numerical results are also obtained to illustrate the feasible region of disease control. Our findings show that due to uncertainty in the attack rate of a pandemic strain, initiating therapy immediately upon diagnosis can significantly increase the likelihood of disease control and substantially reduce the required community-level of treatment. This suggests that reliable diagnostic methods for influenza cases should be rapidly implemented within an antiviral treatment strategy.

Influenza pandemics have historically been associated with substantial morbidity, mortality and socioeconomic upheaval (Cox et al., 2003; Taubenberger and Morens, 2006) . The most devastating pandemic recorded in the history of human diseases occurred in 1918 and was responsible for approximately 50 million deaths among the countless infected (Taubenberger and Morens, 2006) . Such high fatality may result from exposure to a virus with major alterations in the viral genome (due to re-assortment or recombination), or to a novel influenza viral strain introduced in humans by zoonotic switch. This poses the great threat of an imminent influenza pandemic through the evolution of highly pathogenic avian influenza A virus to one that is readily transmitted between humans (Gani et al., 2005) .

Considering the catastrophic results of previous pandemics, the development of effective control strategies is undoubtedly a high priority of global public health efforts. Previous studies have rationalized the use of antiviral agents for both treatment and prophylaxis as a primary control measure during early stages of a pandemic (Ferguson et al., 2005 (Ferguson et al., , 2006 Gani et al., 2005; Germann et al., 2006; Longini et al., 2004 Longini et al., , 2005 . It has been shown that effective antiviral treatment can substantially reduce viral shedding and the infectiousness period of a clinical case. However, to achieve such effectiveness, therapy should be initiated as early as possible within the 48 hours of appearance of clinical symptoms, referred to as the window of opportunity for effective treatment (Ferguson et al., 2006; Linde, 2001) . In most cases, there has already been a time delay for clinical diagnosis, and therefore, also the onset of therapy, which could greatly jeopardize the impact of the treatment strategy. This point has been discussed in previous stochastic models employed to simulate potential outbreaks of the present deadly H5N1 strain in rural Southeast Asia and assess the effectiveness of antiviral prophylaxis (Ferguson et al., 2005; Longini et al., 2005) . Based on social networks and close contacts in household clusters, schools, workplaces, and other social settings, they predicted that a pandemic could be contained at the origin through a combination of targeted blanket prophylaxis and social distancing. A key assumption is that the virus remains less transmissible than pandemic viruses of the last century, so that the reproduction number stays below 2.4. The predictions of these models, however, depend strongly on the specific location of an initial outbreak, patterns of exposure to infection in localities, and how quickly infected cases are diagnosed and their contacts offered prophylaxis. In view of considerable uncertainties in social settings and parameters characterizing the nature of an impending epidemic, it is suggested that compartmental models with relatively few parameters may still be a viable approach to the general understanding of determinants in the course of an outbreak and the impact of interventions (Arino et al., 2006) .

In this paper, we develop a theoretical framework that uses the time elapsed since the onset of symptoms explicitly as an independent structure variable, and thus accurately accounts for the effect of time delay in starting therapy within the limited window of opportunity for antiviral treatment. The resulting epidemic model, formulated as a system of delay differential equations, incorporates a parameter for the level of treatment which is determined by the functional form of the treatment rate varying with age of infection. The importance of such variation and early application of antiviral drugs are highlighted in a recent experimental study of the kinetics of influenza A virus infection in humans (Baccam et al., 2006) . Thus, control of an epidemic is affected by the delay in initiating a course of treatment, as well as the prescribed rate at which diagnosed clinical cases are given priority and treated. The model is analyzed to obtain a criterion for containing the epidemic in the presence of the treatment strategy. Using different profiles of treatment rate, numerical results are obtained to display the feasible region of disease control as a function of the level of treatment and the time delay in the onset of therapy. The results are placed in the context of public health, in terms of epidemiological implications and development of effective treatment strategies.

To represent the clinical features of influenza infection, we divide the whole population into classes of susceptible, exposed (infected but not infectious), asymptomatic and symptomatic infected individuals. We assume that the spread of infection occurs through contacts between susceptible and infected individuals, where mass action incidence is used for sake of simplicity. This simplification can still provide a good approximation to a more realistic situation in which other incidence functions are used (Arino et al., 2006) . Since the period of a pandemic is expected to be relatively short, we ignore the effect of birth and natural death rates on the transmission dynamics of infection. Denoting susceptible and exposed classes by S and E, we have for t ≥ 0

where β is the baseline transmission rate, 1/μ E represents the incubation period, and Q(t) is the force of infection yet to be formulated. The probability for an exposed individual to develop symptoms is assumed to be a constant p ∈ (0, 1). Let A denote the class of asymptomatic individuals who shed virus during their infectious period (1/μ A ). Then

Since treatment of infected individuals is feasible only after the symptoms appear, we may divide the infectious period of the symptomatic infection into two stages. The primary stage, which represents the window of opportunity for initiating therapy, terminates at time t = n. Individuals who have not received antiviral drugs during this window will progress to the secondary stage and receive no treatment for the entire course of symptomatic infection.

The equations governing the disease dynamics within and beyond the primary interval are different. Let a represent the time elapsed since the onset of symptoms, and i U (t, a) and i T (t, a) be the densities, with respect to age a, of untreated and treated individuals, respectively, at current time t . Assuming that the treatment rate at time a after the onset of symptoms is r(a), we have for a ∈ [0, n] (Metz and Diekmann, 1986; Webb, 1985) , that

subject to the following boundary conditions

and initial conditions

Let I U (t) and I T (t) denote, respectively, the total number of untreated and treated individuals when a ≥ n. Thus, we obtain for t ≥ 0

where 1/μ U and 1/μ T are the mean infectious periods, and d U and d T represent diseaseinduced mortality rates of untreated and treated individuals, respectively. Solving (4) subject to (6), we obtain

Substituting this into (5) and using (6) leads to

where

This function describes the proportion of infected individuals which remains untreated until the age of clinical disease a. Therefore, Eqs. (8) and (9) can be written as (t ≥ n)

It should be mentioned that, for 0 ≤ t ≤ a ≤ n, we have from (4), (5) and (7)

and the equations for I U (t) and I T (t) with t ∈ [0, n] are given in (13), (14).

We are now in a position to formulate the force of infection Q(t). The model developed here does not explicitly incorporate a pre-symptomatic stage that corresponds to a short period of time during which infected individuals can transmit the disease with a lower rate of infectiousness, without showing clinical symptoms. We may include the pre-symptomatic stage as part of the symptomatic infection by considering a lag time (τ ) for the appearance of the symptoms. It is, therefore, assumed that antiviral treatment is not administered during the pre-symptomatic stage, and hence r(a) = 0 for 0 ≤ a ≤ τ . Thus, q(a) ≡ 1 for 0 ≤ a ≤ τ , and the force of infection may be expressed as

where δ A , δ P and δ U are the reduction factors in transmission rate during asymptomatic, pre-symptomatic, and symptomatic infection when a ≥ n, respectively, and δ T is the reduction factor in infectiousness (viral shedding) due to the treatment. A schematic diagram of the model is given in Fig. 1 , where R and D denote classes of recovered and dead individuals, respectively. Summarizing the above details, the model equations with both discrete and distributed delays may be written as (for t ≥ n)

where " " represents the derivative of the variables with respect to the time, q n := q(n) with q(a) defined in (12), and

We conclude this section with a final remark. In what follows, we will concentrate on the system (18-22) with Q(t) defined in (23). It should be noted, however, that the above system describes the transmission dynamics of disease only for t ≥ n. For t ∈ [0, n], the system consists of (18-20) together with (8), (9) and (17), for which the value of i U (t, a) and i T (t, a) (for t ≤ a) are given in terms of the initial values i U 0 (a) and i T 0 (a), as discussed earlier.

In this section, we analyze the model for its qualitative behavior and compute the control reproduction number as a threshold quantity determining the feasibility of disease control in the presence of antiviral treatment.

Since the model describes the dynamics of an epidemic, we consider only biologically feasible initial values, namely, S(0), A(0), I U (0), I T (0) ≥ 0 and E(s) ≥ 0 for all s ∈ [−n, 0], where the initial value of E must be specified on the whole interval [−n, 0] and we assume this initial function is integrable. It is easy to show that S(t), E(t), A(t), I U (t), I T (t) ≥ 0 for all t > 0. Note that the model has a set of diseasefree equilibria (DFE) given by E 0 = (S(0), 0, 0, 0, 0) with arbitrary S(0). In the presence of infection, S(t) is strictly decreasing and hence the model exhibits no endemic equilibrium.

A key parameter in epidemic models is the basic reproduction number (R 0 ), 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 and Heesterbeek, 2000) . A related quantity, called the control reproduction number (R c ), can be defined when certain control measures (such as the use of antiviral drugs for treatment) are implemented, and used to determine whether the epidemic can be contained (Gumel et al., 2004) .

In what follows, we fix the initial size S(0) of the susceptible population. To compute R c , we assume that an infected individual is introduced into the E-class such that E(0) = 1, and let E(t) = 0 for t ∈ [−n, 0], and A(0) = I U (0) = I T (0) = 0. After a short incubation period (1/μ E ), the infected individual can transmit the disease through either asymptomatic or symptomatic infection. Considering the duration and transmission rates associated with asymptomatic, untreated symptomatic and treated symptomatic stages (see Fig. 1 ), the total number of secondary cases generated in the A, I U and I T classes is given by

The number of secondary cases generated by the infected individual during the window of opportunity is given by

Adding (24) and (25) gives the control reproduction number as

In the absence of antiviral treatment (q(a) ≡ 1), R c reduces to the basic reproduction number

Using (18), we obtain S(t) ≤ S(0) for all t ≥ 0, and therefore (19) becomes

Thus, the standard comparison argument can be applied to show that the disease components (E, A, I U , I T ) of a solution of the model (18-22) are bounded by the corresponding components of the following linear system

Substituting the ansatz we λt , where w = (E 0 , A 0 , I 0 U , I 0 T ), and assuming that E(t) = e λt is an eigenfunction with E 0 = 1, it follows that

Eliminating e λt from (33-35) gives

Substituting in (32) and using (17), a simple calculation yields that λ is a root of the characteristic function

Let λ = x + iy be a complex root with x ≥ 0. Thus, |e −λa | ≤ 1 for any a ≥ 0, and we have

Therefore, if R c < 1, then all complex roots of h(λ) have negative real parts. Let

It is now easy to show that, for real values of λ > M, dh(λ)/dλ < 0, and hence h is a monotone decreasing function on (M, +∞). Furthermore, h(M) = +∞ and h(∞) = −∞. Noting that h(0) = μ E (R c − 1), it follows that h(λ) has a unique negative real root if and only if R c < 1. This, together with the above discussion for complex roots, implies that all the roots of the characteristic function h(λ) have negative real parts, and consequently the zero solution of the system (28-31) is asymptotically stable.

Since the linear system (28-31) is a cooperative and irreducible system of delay differential equations (Smith, 1995) , it follows that the zero of the characteristic function h(λ) with the largest real part, denoted by λ 0 , is a real number (Smith, 1995 , Theorems 3.1, 5.1), and the corresponding eigenvectors are (e λ 0 s , A * , I * U , I * T ) for s ∈ [−n, 0], with positive components (1, A * , I * U , I * T ). Thus, if initially E(s) ≤ E * e λ 0 s , A(0) ≤ A * E * , I U (0) ≤ I * U E * , and I T (0) ≤ I * T E * for some constant E * > 0, then a standard comparison argument shows that (E(t), A(t), I U (t), I T (t)) ≤ (1, A * , I * U , I * T )E * e λ 0 t for all t ≥ 0. Therefore, the disease components are bounded by exponential functions with negative exponent λ 0 for t ≥ 0, and hence the occurrence of an outbreak is prevented if R c < 1.

In this section, we establish a relationship between R c and the final size of the susceptible population when the epidemic dies out. This can be used to estimate the basic reproduction number, and therefore the baseline transmission rate, using the clinical attack rate of an influenza pandemic strain. We take advantage of this estimation for the purpose of numerical simulations.

For simplicity of the analysis, we use the following notation for any continuous function F (t):

and (36), (20), (21) and (22) from 0 to ∞ gives, respectively,

and hence we havê

Using these relations, and integrating (18), it follows that

For the particular case where E(t) = 0 for −n ≤ t < 0, and A(0) = I U (0) = I T (0) = 0, it can be seen thatẼ(a) = 0 for all 0 ≤ a ≤ n. Thus, if E(0) is small compared with S(0), then from (37) we obtain the final size relation

In this section, we simulate the model by calculating the values of some critical parameters, such as the basic reproduction number and the baseline transmission rate. These values are obtained using the final size relation given by (38), and an assumed clinical attack rate within the estimated range of previous pandemic strains. The clinical attack rate is defined as the fraction of the initial susceptible population which develops clinical symptoms. With our notation, this fraction is given by p(1 − S(∞)/S (0)). The purpose of our simulations is to evaluate the impact of three important parameters on the control reproduction number: (i) the level of treatment; (ii) the delay in the time of initiating the course of treatment; and (iii) the profile of treatment rate, as defined by r(u), subject to certain constraints.

In the absence of treatment, Eq. (38) can be used to estimate the basic reproduction number. Using a 30% clinical attack rate (Gani et al., 2005; Longini et al., 2004 Longini et al., , 2005 , and assuming that 60% of exposed individuals develop clinical symptoms (Ferguson et al., 2005 (Ferguson et al., , 2006 Longini et al., 2004 Longini et al., , 2005 , we obtain S(∞)/S(0) = 0.5. This gives an estimate of R 0 = 1.3863, using the final size relation. Now, the baseline transmission rate (β) can be calculated from the expression (27), in which other parameters are obtained from the published literature (Table 1) . We estimated a range of β corresponding to the range of 20%-35% clinical attack rate, and assumed that the treatment is administered within 2 days of clinical infection developing, with a delay of at least 6 hours for the pre-symptomatic stage (τ = 0.25 day). This is consistent with the assumptions made in a recent modeling study for a 0.25-day delay between the onset of clinical symptoms and diagnosis of the infected case (Ferguson et al., 2005) . Table 1 Description of the model parameters with their values obtained from the published literature (Arino et al., 2006; Ferguson et al., 2005 Ferguson et al., , 2006 Jefferson et al., 2006; Longini et al., 2004 Longini et al., , 2005 We prescribe the profile of treatment rate as a function that increases linearly with time during the window of opportunity for treatment. Assuming that treatment commences with the onset of the clinical symptoms (optimal scenario), we define

where b is the slope of the treatment function r * . To determine the value of b, we let c denote the treatment level of clinical cases under the profile of r * . Thus, q n = 1 − c, which gives

where r max is the maximum rate of treatment. This gives b = r max /(n − τ ). Having defined r * , we considered a similar rate r a 0 (u) when treatment begins with a delay a 0 during the time interval [τ, n] (see Fig. 2(a) ). This profile of r a 0 (u) corresponds to a strategy that increases the treatment rate towards the end of the window of opportunity for initiating a course of antiviral drugs. It effectively gives priority to individuals who are clinically diagnosed close to the end of this window. Thus, the functional form of q(a) with delay may be expressed as Figure 3 illustrates the feasible region of disease control (R c < 1) for different clinical attack rates, when τ ≤ a 0 ≤ n and 0.5 ≤ c ≤ 1. These simulations clearly show that higher levels of treatment are required as the delay in starting the course of treatment increases. This leads to a restricted region in which R c < 1, with increasingly stringent requirements of treatment level as the clinical attack rate increases.

To illustrate the comparative outcome of a different treatment profile, we considered a functional form of r * (u) (Fig. 2b ) defined as

which corresponds to a strategy that decreases the treatment rate across the window of opportunity, thus effectively giving priority to those clinically diagnosed early in this window. A simple calculation yields

Using different clinical attack rates, we ran the simulations to illustrate the feasible region of disease control (R c < 1). Figure 4 shows that profile of r a 0 (u) can play a ma- Fig. 4 The feasible region of disease control (shaded area), as a function of treatment level and delay in initiating antiviral treatment, corresponding to the profile of r a 0 (u) in Fig. 2(b) and attack rates between 20% and 35%. Parameter values are: 1/μ E = 1.25 days, 1/μ U = 2.85 days, 1/μ T = 1.35 days, 1/μ A = 4.1 days, d U = 0.002 day −1 , d T = 0.001 day −1 , δ P = 0.286, δ A = 0.071, δ U = 0.143, δ T = 0.4, p = 0.6, τ = 0.25, and n = 2. jor role in determining the outcome of a treatment strategy. As is evident, disease can be controlled with a lower level of treatment and even longer delay than was the case for simulations in Fig. 3 . Therefore, the strategy that gives priority to early treatment has a better chance of controlling the disease, and confers greater overall benefit to the healthcare system by substantially reducing the required level of treatment. For example, with no delay and an attack rate of 30%, a minimum of 85% treatment level is required to ensure R c < 1 in the former strategy (Fig. 3) , while this level is reduced to 62% in the latter strategy (Fig. 4) . Correspondingly, the disease cannot be controlled, even with 100% coverage, if the delay in initiating the course of treatment is longer than 0.8 day in the former, and 1.1 days in the latter strategy.

Strategies for mitigating the impact of an influenza pandemic have been evaluated in several studies (Ferguson et al., 2005 (Ferguson et al., , 2006 Germann et al., 2006; Longini et al., 2004 Longini et al., , 2005 World Health Organization Writing Group, 2006) , with significant public health implications for global containment. These strategies include antiviral treatment, vaccination, and non-pharmaceutical measures (social distancing, quarantine/isolation, school/border closure). Although a combination of these control measures would dramatically increase the likelihood of containing a pandemic, there are several limitations associated with their application (Ferguson et al., 2005 (Ferguson et al., , 2006 Germann et al., 2006) . Recent evaluations identify the use of antiviral treatment as the most effective strategy early on in a pandemic response, with the stringent requirement of starting therapy within a day of clinical symptoms developing (Ferguson et al., 2006) .

In this study, we developed a mathematical model to evaluate the effect of delay in initiating a course of antiviral treatment on containing a pandemic. To illustrate the feasible region of disease control, where R c is reduced to below one, we examined two different strategies corresponding to treatment rates that increase (decrease) toward the end of the window of opportunity for starting treatment. Simulation results, over an estimated range of attack rates (Gani et al., 2005; Longini et al., 2004 Longini et al., , 2005 , reveal that control of the disease is more likely to be achieved, and with less stringent requirements on the community-level of treatment, as the delay in onset of therapy diminishes. Therefore, a strategy that favors therapy immediately following the onset of clinical symptoms is more effective and confers greater population-wide benefit. This highlights the importance of early diagnosis of clinical infections through rapid implementation of reliable evaluation tests for influenza cases.

Our model assumes the availability of an unlimited supply of antiviral drugs for the entire course of a pandemic. In the case of a limited supply, the treatment rate would need to be modified to take into account its depletion, which would in turn affect the total number of treated clinical cases. Such modification introduces the drug supply as a new parameter in the model, which may affect the control reproduction number, and consequently the results of simulations reported here. However, whether the supply is limited or not, extreme caution is required for preventing the evolution of drug-resistant viral strains through the widespread or indiscriminate use of antiviral drugs (Ferguson et al., 2003; Stilianakis et al., 1998) . This is crucial for reducing the likelihood of multiple pandemic waves caused by the emergence of resistant viral mutants. These considerations, though important in practice, fall outside the scope of the present study but merit further investigation.

Simple models for containment of a pandemic

Kinetics of influenza A virus infection in humans

Influenza pandemic planning

Mathematical Epidemiology of Infectious Diseases

A population-dynamic model for evaluating the potential spread of drug-resistant influenza virus infections during community-based use of antivirals

Strategies for containing an emerging influenza pandemic in Southeast Asia

Strategies for mitigating an influenza pandemic

Potential impact of antiviral drug use during influenza pandemic

Mitigation strategies for pandemic influenza in the United States

Modeling strategies for controlling SARS outbreaks

Antivirals for influenza in healthy adults: systematic review

The importance of specific virus diagnosis and monitoring for antiviral treatment

Containing pandemic influenza with antiviral agents

Containing pandemic influenza at the source

The Dynamics of Physiologically Structured populations

Monotone Dynamical Systems, an Introduction to the Theory of Competitive and Cooperative Systems

Emergence of drug resistance during an influenza epidemic: insights from a mathematical model

influenza: the mother of all pandemics

Theory of Nonlinear Age-Dependent Population Dynamics. Dekker, New York. World Health Organization Writing Group