Sustained oscillations in an evolutionary epidemiological model of influenza A drift

Pierre Magal, Shigui Ruan

Abstract

Understanding the seasonal/periodic reoccurrence of influenza will be very helpful in designing successful vaccine programmes and introducing public-health interventions. However, the reasons for seasonal/periodic influenza epidemics are still not clear, even though various explanations have been proposed. In this paper, we study an age-structured type evolutionary epidemiological model of influenza A drift, in which the susceptible class is continually replenished because the pathogen changes genetically and immunologically from one epidemic to the next, causing previously immune hosts to become susceptible. Applying our recently established centre manifold theory for semi-linear equations with non-dense domain, we show that Hopf bifurcation occurs in the model. This demonstrates that the age-structured type evolutionary epidemiological model of influenza A drift has an intrinsic tendency to oscillate owing to the evolutionary and/or immunological changes of the influenza viruses.

1. Introduction

Influenza viruses are negative-stranded RNA viruses that are divided into three types: A, B and C, based on antigenic differences in nucleoprotein and matrix protein (Earn et al. 2002; Schweiger et al. 2002; WHO 2009). Influenza A and B are currently associated with human diseases, while influenza C has only subclinical importance (Schweiger et al. 2002). Influenza A viruses are further classified into subtypes based on two different spike-like protein components, haemagglutinin (HA) and neuraminidase (NA), called antigens, on the surface of the virus (WHO 2009). There are 16 different HA subtypes and nine different NA subtypes. H1N1 (the subtype responsible for the 1918 Spanish pandemic) and H3N2 (the strain that caused the 1968 Hong Kong pandemic) are two influenza A subtypes that are important for human (WHO 2009), which, together with influenza B, are included in each year’s influenza vaccine (CDC 2009).

Genes of influenza viruses mutate with high frequency (Buonagurio et al. 1986; Fitch et al. 1997), which plays an important role in causing reccurent influenza epidemics (Pease 1987). Influenza type A viruses undergo two kinds of changes. Antigenic drift is a major process that accumulates point mutations at the antibody sites in the HA protein, leading to the emergence of immunologically distinct strains and enabling the viruses to evade recognition by hosts’ antibodies (Schweiger et al. 2002; Treanor 2004; Shih et al. 2007). The immunologically distinct influenza virus strains can reinfect hosts that are immune to the progenitor influenza strain and reinvade the communities that recently suffered an epidemic of the progenitor strain (Pease 1987). Antigenic shift is an abrupt major change resulting in new HA and/or new HA and NA proteins in influenza viruses that infect humans. When shift happens, most people have little or no protection against the new influenza A subtype virus. While influenza viruses are changing by antigenic drift every time, antigenic shift happens only occasionally and causes pandemics (CDC 2009; WHO 2009).

During antigenic drift, point mutations, including substitutions, deletions and insertions produce genetic variation in influenza viruses. These genetic changes encode amino acid changes on the surface proteins that permit the virus to escape neutralization using the antibody generated for previous strains (Cox & Subbarao 2000). It is reported that, for sites involved in antigen determination, amino acid substitutions are more frequent than synonymous substitutions (Ina & Gojobori 1994; Earn et al. 2002). Also, an approximately equal number of amino acid substitutions occur in the influenza each year between the immunizing and the challenge virus strains (Pease 1987; Smith et al. 2004; Koelle et al. 2006).

To study the evolutionary epidemiological mechanism of influenza A drift, Pease (1987) developed an age-structured type model in which the susceptible class is continually replenished, as the pathogen changes genetically and immunologically from one epidemic to the next, causing previously immune hosts to become susceptible. Conditions were derived to show how the equilibrium number of infected hosts, the interepidemic period and the probability that a host becomes reinfected depend on the rate of amino acid substitution in the pathogen, the effect of these substitutions on host immunity, the host size and the recovery rate. However, applied to influenza A, the model only predicted damped oscillations, while the observed oscillations in the incidence of influenza are undamped over a long term (Pease 1987; Dushoff et al. 2004; Viboud et al. 2006).

Pease’s model has been modified and improved by several researchers. For example, Thieme & Yang (2002) reparameterized Pease’s model by using recovery age, which is the duration of time that has passed since the moment of recovering from the disease, as a structural variable. They investigated uniform disease persistence and global stability of the endemic equilibrium in their model with variable reinfection rate. Andreasen et al. (1996) improved Pease’s model by including mutation as a diffusion process in a phenotype space of variants. The existence of travelling waves in this generalized susceptible-infected-recovered (SIR) model of influenza A drift was established by Lin et al. (2003). Andreasen et al. (1997) developed models that describe the dynamics of multiple influenza A strains conferring partial cross-immunity, and showed that such models exhibit sustained oscillations. See also Castillo-Chavez et al. (1989), Lin et al. (1999), Andreasen (2003), Nuño et al. (2005) and the references cited therein for other multiple-strain influenza models with cross-immunity or immune-selection. Note that Castillo-Chavez et al. (1989) observed that sustained oscillations do not seem possible for a single-strain model, even with age-specific mortalities.

Inaba (1998, 2002) modified Pease’s model slightly by assuming that the effect of amino acid substitutions on host immunity (the transmission rate) is a monotone increasing function with a finite upper bound, proved the existence and uniqueness of solutions of the model and studied the stability of the endemic steady state by using the semigroup approach. Inaba (2002) further conjectured that, for realistic parameter values, the steady state could be destabilized and lead to a periodic solution via Hopf bifurcation.

The existence of non-trivial periodic solutions in age-structured models has been a very interesting and difficult problem. It is believed that such periodic solutions are induced by Hopf bifurcation; however, there is no general Hopf bifurcation theorem available for age-structured models (Iannelli 1994; Inaba 2002). One of the difficulties is that, rewriting age-structured models as a semilinear equation, the domain of the linear operator may not be dense. Recently, we (Magal & Ruan 2009a) developed the centre manifold theory for semilinear equations with non-dense domain.

In this paper, following the pioneering work of Kermack & McKendrick (1927, 1932, 1933), we first modify the model of Pease (1987) and Inaba (2002) by using the age of infection (instead of the number of amino acid substitutions) as a variable. We then apply the centre manifold theorem in Magal & Ruan (2009a) to study Hopf bifurcation in the age-structured type evolutionary epidemiological model of influenza A drift. The results indicate that sustained oscillations can occur in such models with a single strain (Castillo-Chavez et al. 1989) and without the seasonal force (Dushoff et al. 2004). Thus, the periodic recurrence of influenza may be a result of the evolutionary and/or immunological changes of the influenza viruses.

2. The model

To model the evolutionary epidemiological mechanism of influenza A drift, Pease (1987) made the following assumptions: (i) there is an approximately linear relation between the probability of reinfection and the number of amino acid substitutions between the immunizing and the challenge virus strains, (ii) only one virus strain circulates in the human population at any one time, and (iii) random drift causes amino acid substitutions to occur in the influenza virus.

Suppose that the total host population size N is a constant. Let I(t) be the number of infected individuals at time t. Let a≥0 be the time since the last infection, that is, the duration of time since an individual has been susceptible. Assume that the average number of amino acid substitutions is a continuous variable. More precisely, let k>0 be the average number of amino acid substitutions per unit of time (i.e. the mutation rate). Then, the number of substitutions after a period of time a in the susceptible class is given by ka. Let s(t,a) be the density of uninfected hosts (structured with respect to a), so thatEmbedded Image is the number of uninfected hosts that were last infected by a virus which differed by more than ka0 and less than ka1 amino acid substitutions from the virus strain prevailing at time t. Here, we consider a modified version of the model considered by Pease (1987) and Inaba (1998, 2002). The model to describe the evolutionary epidemiological model of influenza A drift takes the following form:Embedded Image 2.1 where ν>0 is the recovery rate of the infected hosts and Embedded Image describes how amino acid substitutions affect the probability of reinfection.

First, we make the following assumptions.

Assumption 2.1.

Assume that k>0, ν>0, Embedded Image and Embedded Image

From now on, we setEmbedded Image

In this model, we do not have a disease-free equilibrium. Indeed, we observe that, if I0=0, thenEmbedded Image is a solution of equation (2.1). So we do not obtain an equilibrium solution in the usual sense. Nevertheless, the total number of susceptible individuals remains constant with time, that is, Embedded Image It is clear thatEmbedded Image so d[S(t)+I(t)]/dt=0. It follows that I(t)=NS(t),∀t≥0, where N is the total number of individuals in the population. With this, we can rewrite system (2.1) as follows:Embedded Image 2.2

In Pease’s original article (1987), the function γ(ka) is a linear function. To consider such unbounded functions, we must change the state space by introducing the weighted L1 space in order to discuss the existence of solutions. Since this is not the goal of this work, we make a simplified assumption that γ(ka) is bounded. More precisely, we consider the following step function. DefineEmbedded Image as the threshold of sensitivity, where τ will be discussed later. We make the following assumption.

Assumption 2.2.

Assume thatEmbedded Image where δ>0 and ρ≥0.

The function γ(ka) describes the rate at which individuals become reinfected per average number of amino acid substitution ka. The above construction thus provides a threshold ρ that is the time necessary to be reinfected after one infection. This is also equivalent to assuming that it is necessary to reach a threshold value τ for the average number of amino acid substitutions before reinfection. When the threshold value ρ=0, the function γ(ka) becomes a constant δ, and S(t) satisfies the following ordinary differential equation:Embedded Image The following is an easy consequence of the fact that S(t) satisfies the above ordinary differential equation.

Lemma 2.3.

Let assumptions (2.1 ) and (2.2 ) be satisfied and assume that τ=0. Then, there are two cases

  • — if ν/δN and I(0)=NS(0)>0, then Embedded Image and Embedded Image as Embedded Image and

  • — if ν/δ<N, then system (2.1) has a unique equilibrium solution Embedded Image such that Embedded Image, Embedded Image and Embedded Image Moreover, if I(0)=NS(0)>0, then Embedded Image as Embedded Image in Embedded Image Furthermore, Embedded Image is globally asymptotically stable for system (2.1) in Embedded Image

The above lemma shows that, when τ=0, if there exists an endemic equilibrium, then it is locally asymptotically stable. Nevertheless, since the system does not depend smoothly on τ, we cannot conclude that, for τ>0 small enough, the endemic equilibrium Embedded Image is also asymptotically stable.

From now on, we assume that τ>0 and make the following change of variable. By consideringEmbedded Image we obtainEmbedded Image and Embedded Image Thus, with the above changes of variables, we obtain the following system:Embedded Image 2.3 whereEmbedded Image Since the case τ=0 has been fully described in lemma 2.3, in the sequel we will only consider the case τ>0. Moreover, by making the above changes of variables, and by replacing the parameters δ and ν, respectively, by Embedded Image and Embedded Image in system (2.2), we can assume, without loss of generality, thatEmbedded Image

3. Preliminary results

The global existence and positivity of solutions follow by considering the systemEmbedded Image 3.1 where Embedded Image

(a) Global existence

In order to give a sense to the solutions of problem (3.1) we will use integrated semigroup theory. We refer to Arendt (1987), Neubrander (1988), Kellermann & Hieber (1989), Thieme (1990b, 1997), Arendt et al. (2001) and Magal & Ruan (2007, 2009a,b) for more detailed results on the subject.

Consider the Banach space Embedded Image endowed with the usual product norm Embedded Image. Consider a linear operator Embedded Image defined byEmbedded Image Then, Embedded Image Set Embedded Image and Embedded Image We also consider the nonlinear operator Embedded Image, defined byEmbedded Image By identifying s(t,.) with Embedded Image, we can rewrite the system as the following abstract Cauchy problem:Embedded Image 3.2 The global existence and positivity of solutions of equation (3.2) now follow from the results in Thieme (1990a), Magal (2001) and Magal & Ruan (2007, 2009a). See also the results in Inaba (1998, 2002).

Lemma 3.1.

For each ν>0 and δ>0, there exists a unique continuous semiflow {U(t)}t≥0 on X0+, such that, for each xX0+, the map Embedded Image is an integrated solution of the Cauchy problem (3.2), i.e. Embedded Image satisfies, for each t≥0, that Embedded Image andEmbedded Image

(b) Equilibrium solutions

Notice thatEmbedded Image We obtain the following lemma.

Lemma 3.2.

If S(0)<1, then S(t)∈[0,1) for t≥0 and Embedded Image In particular, ifν/δ≥1, then Embedded Image

If Embedded Image is an equilibrium solution of the system, we must haveEmbedded Image This is equivalent toEmbedded Image with Embedded Image So we must have Embedded Image and obtainEmbedded Image 3.3 By integrating Embedded Image over Embedded Image we obtain Embedded Image which implies that Embedded Image Thus,Embedded Image 3.4 Moreover, Embedded Image if and only if δ>ν. Now we can state the following lemma.

Lemma 3.3.

System (3.1) has an equilibrium if and only if δ>ν. Moreover, when the equilibrium exists, it is unique and given by equation (3.3 ).

4. Linearized equation

Since Embedded Image the map Embedded Image is differentiable in a neighbourhood of Embedded Image and the linearized equation of (3.2) at Embedded Image is given byEmbedded Image 4.1 We may rewrite this as the following PDE:Embedded Image 4.2

(a) Associated linear operator

We consider a linear operator Embedded Image, defined byEmbedded Image and a bounded linear operator Embedded Image, defined byEmbedded Image Then, Embedded Image and the Cauchy problem (4.1) can be rewritten asEmbedded Image By applying the results in Thieme (1990a) or Magal (2001) and Magal & Ruan (2007), we obtain the following result.

Lemma 4.1.

The resolvent set of B contains Embedded Image More precisely, for each Embedded ImageEmbedded Image Moreover, B is a Hille–Yosida operator andEmbedded Image with Embedded Image and Embedded Image

Recall that B0, the part of B in Embedded Image is defined byEmbedded Image andEmbedded Image So, we obtainEmbedded Image where Embedded Image with Embedded Image Since B is a Hille–Yosida operator, it follows that B0 is the infinitesimal generator of a strongly continuous semigroup of bounded linear operators {TB0(t)}t≥0 on X0. More precisely, we have the following result.

Lemma 4.2.

The linear operator B0 is the infinitesimal generator of a strongly continuous semigroup of bounded linear operators {TB0(t)}t≥0 on X0, and for each t≥0, the linear operator TB0(t) is defined byEmbedded Image where Embedded Image is a strongly continuous semigroup of bounded linear operators on Embedded Image andEmbedded Image Moreover, Embedded Image

Definition 4.3.

Let Embedded Image be the infinitesimal generator of a linear C0-semigroup {TL(t)}t≥0 on a Banach space X. Define the growth boundEmbedded Image of L byEmbedded Image The essential growth boundEmbedded Image of L is defined byEmbedded Image where ∥TL(t)∥ess is the essential norm of TL(t) defined byEmbedded Image Here, BX(0,1)={xX:∥xX≤1}, and for each bounded set Embedded ImageEmbedded Image is the Kuratovsky measure of non-compactness.

Since ∥TB(t)∥ess≤∥TB(t)∥, we deduce that ω0,ess(B)≤ω0(B). By using lemma 4.2, we have Embedded Image Since C is a bounded linear operator, it follows that Embedded Image is a Hille–Yosida operator. Thus, (B+C)0, the part of (B+C) on X0, is the infinitesimal generator of a strongly continuous semigroup of bounded linear operators {T(B+C)0(t)}t≥0 on X0. Moreover, since C is compact, by using the perturbation results in Thieme (1997) and Ducrot et al. (2008), it follows thatEmbedded Image 4.3

Next, we recall some results of spectral theory. In the following result, the existence of the projector was first proved by Webb (1985, 1987) and the fact that there is a finite number of points of the spectrum was proved by Engel & Nagel (2000).

Theorem 4.4.

Let Embedded Image be the infinitesimal generator of a linear C0-semigroup {TL(t)} on a Banach space X. Then,Embedded Image Assume in addition that ω0,ess(L)<ω0(L). Then, for each γ∈(ω0,ess(L),ω0(L)], Embedded Image is non-empty, finite and contains only poles of the resolvent of L. Moreover, there exists a finite rank bounded linear projector Embedded Image satisfying the following properties: (i ) Π(λL)−1=(λL)−1Π,∀λρ(L), (ii ) σ(LΠ(X))={λσ(L):Re(λ)≥γ}, and (iii ) σ(L(IΠ)(X))=σ(L)∖σ(LΠ(X)).

In theorem 4.4, the projector Π is the projection on the direct sum of the generalized eigenspaces of L associated to all points λσ(L) with Re(λ)≥γ. As a consequence of theorem 4.4, we have the following corollary.

Corollary 4.5.

Let Embedded Image be the infinitesimal generator of a linear C0-semigroup {TL(t)} on a Banach space X, and assume that ω0,ess(L)<ω0(L). Then,Embedded Image and each Embedded Image is a pole of the resolvent of L. That is, Embedded Image is isolated in σ(L), and there exists an integer k0≥1 (the order of the pole), such that Laurent’s expansion of the resolvent takes the following form:Embedded Image where Embedded Image are bounded linear operators on X, and the above series converges in the norm of operators whenever |λλ0| is small enough.

(b) Characteristic equation

By using lemma 2.1 in Magal & Ruan (2009a), we know that σ((B+C)0)=σ((B+C)). By using theorem 4.4, corollary 4.5 and equation (4.3), we deduce that, for each ε>0, Embedded Image is either empty or finite. Moreover, we have the following result.

Lemma 4.6.

For each Embedded Image with Embedded ImageEmbedded Image whereEmbedded Image Furthermore, if Embedded Image with Embedded Image and Δ(λ)≠0, thenEmbedded Image is equivalent toEmbedded Image whereEmbedded Image 4.4

Proof.

For Embedded Image with Embedded Image since λIB is invertible, it follows that λI−(B+C) is invertible if and only if IC(λIB)−1 is invertible andEmbedded Image Note thatEmbedded Image 4.5 whereEmbedded Image So we must haveEmbedded Image By using the second equation of system (4.5), we obtainEmbedded Image Thus, we haveEmbedded Image 4.6 whereEmbedded Image and Embedded Image is given by equation (4.3). Moreover,Embedded Image 4.7 where Embedded Image Combining the first equation of system (4.5) with equations (4.6) and (4.7), we obtain Embedded Image and Embedded Image So,Embedded Image andEmbedded Image The characteristic function is given byEmbedded Image It follows thatEmbedded Image 4.8 By using the second equation of (4.5), Embedded Image together with equations (4.7) and (4.8), we obtainEmbedded Image Hence, Embedded Image Finally, we deduce thatEmbedded Image and the explicit formula for the resolvent of B+C follows by using the explicit formula for the resolvent of B.

Now assume that Δ(λ)=0. It remains to show that λI−(C+B) is not invertible, so it is sufficient to find Embedded ImageD(A), such thatEmbedded Image which is, in turn, equivalent toEmbedded Image If we fix Embedded Image then we obtainEmbedded Image Since Δ(λ)=0, we deduce that Embedded Image The result follows. ■

Lemma 4.7.

Δ(0)=1+ν>0. So 0∉σ(B+C).

Proof.

We have Δ(0)=1+νJ1J2, withEmbedded Image andEmbedded Image But, Embedded Image so we haveEmbedded Image which implies that Δ(0)=1+ν>0. This completes the proof. ■

Lemma 4.8.

For each Embedded Image with Embedded Image and λ≠0, we haveEmbedded Image

Proof.

We have Δ(λ)=1+νI1I2, withEmbedded Image andEmbedded Image The result thus follows. ■

The following lemma demonstrates that the purely imaginary eigenvalues of B+C are simple.

Lemma 4.9.

Assume that Embedded Image for some Embedded Image with ω>0. Then, Embedded Image is a simple eigenvalue of B+C.

Proof.

Assume that Embedded Image with ω>0 is a solution of Embedded Image Assume first thatEmbedded Image 4.9 Then, by using the explicit formula for the resolvent of B+C in lemma 4.6, we deduce that the following limit exists:Embedded Image Moreover, we haveEmbedded Image So, Embedded Image is a pole of order 1 of the resolvent of B+C, and Embedded Image the projector on the generalized eigenspace of B+C associated to Embedded Image has rank 1. It follows that Embedded Image is a simple eigenvalue of B+C.

It remains to show equation (4.8). We have thatEmbedded Image is equivalent toEmbedded Image 4.10 Notice thatEmbedded Image Therefore,Embedded Image and by using equation (4.10) we obtainEmbedded Image which is equivalent to Embedded Image. Since Embedded Image with ω>0, we obtain Embedded Image, which is impossible since the imaginary part is non-null. ■

5. Stability

In this section, we investigate the stability of the endemic equilibrium Embedded Image. We observe that the linear operator Embedded Image can be decomposed as the sum of a linear operator Embedded Image, defined byEmbedded Image and a bounded linear operator Embedded Image, defined byEmbedded Image That is, we have Embedded Image. Observe that Embedded Image. In order to have the existence of the endemic equilibrium (i.e. Embedded Image or Embedded Image), we must impose δ>ν. Fix δ>0. Then, we have Embedded Image as Embedded Image. Hence, Embedded Image as Embedded Image. Moreover, Embedded Image, the part of Embedded Image on X0, is the infinitesimal generator of the strongly continuous semigroup of bounded linear operators Embedded Image on X0, defined byEmbedded Image whereEmbedded Image Since Embedded Image it follows that Embedded Image By using the standard perturbation argument (see Thieme 1990a or Magal & Ruan 2009b), we obtain the following stability result.

Lemma 5.1.

Let δ>0 be fixed. Then, there exists Embedded Image, such that, for each Embedded Image the equilibrium Embedded Image is locally asymptotically stable.

By using a similar argument we also have the following lemma.

Lemma 5.2.

There exists ε>0, such that, if δ+νε, then the equilibrium Embedded Image is locally asymptotically stable.

Remark 5.3.

Now we return to the original system (2.1) and assume that the parameters ν>0, δ>0 and N>0 of system (2.1) are fixed such that ν/δ<N. Then, we can regard τ>0, defined in assumption 2.2, as a parameter of system (2.1). By the changes of variables described in §2, it corresponds to some parameters Embedded Image and Embedded Image of system (3.1) with Embedded Image and Embedded Image for some c1>c2>0. By lemma 5.2, we deduce that, if τ>0 is small enough, then the endemic equilibrium of the original system (2.1) is stable (see the numerical simulations in §7).

We have Embedded Image. Let Embedded Image be defined byEmbedded Image 5.1 Then, the characteristic equation can be rewritten as 1+ν[h(ν)/λ(λ+h(ν))](1−eλ)=0. Since Re(λ)>−h(ν), the characteristic equation is equivalent toEmbedded Image 5.2

Lemma 5.4.

Assume that Embedded Image with Re(λ)≥0 satisfies the characteristic equation. Then, we have the estimate |λ|≤2ν+2h(ν).

Proof.

From equation (5.2), we have [λ+h(ν)]2=h(ν){ν(eλ−1)+[λ+h(ν)]}. So [(λ/h(ν))+1]2=(ν/h(ν))(eλ−1)+[(λ/h(ν))+1]. It follows that |(λ/h(ν))+ 1|2 ≤ (ν/h(ν))(eRe(λ) + 1) + |(λ/h(ν)) + 1|. Since Re(λ) ≥ 0, we obtain that Embedded Image, it follows that |λ+h(ν)|≤2ν+h(ν) and |λ|≤|λ+h(ν)|+|h(ν)|≤2ν+2h(ν). This proves the lemma. ■

Now, we are in the position to state and prove the main result in this section.

Theorem 5.5 (stability)

Let δ>0 be fixed. Then, there exists Embedded Image, such that, for each Embedded Image the equilibrium Embedded Image of system (3.1) is locally asymptotically stable.

Proof.

Assume that there exist two sequences {νn} with Embedded Image and Embedded Image such that Re(λn)≥0 and Δ(λn)=0. Then, by lemma 5.4, we have |λn|≤2νn+2h(νn), ∀n≥0. We can find a converging subsequence {λnp} with Embedded Image. By using the characteristic equation, we obtainEmbedded Image and Embedded Image. So λ*=0. Notice thatEmbedded Image so we obtain Embedded Image. Thus, for p≥0 large enough, we have Re(λnp/h(νnp))<−(1+δ/2), which implies, for each p≥0 large enough, thatEmbedded Image a contradiction. ■

6. Hopf bifurcation

Recall that h(ν)=(δν/1+ν), we have Embedded Image. So h is strictly decreasing and h(0)=δ and h(δ)=0. Thus, Embedded Image. The characteristic equation becomes equivalent to find Embedded Image with Re(λ)>−h(ν), such that one of the following is satisfied:Embedded Image

(a) Existence of purely imaginary solutions

Assume that λ=iω, with ω>0, is a solution of the characteristic equation. We must haveEmbedded Image 6.1 This implies that Embedded Image and Embedded Image. So we obtain Embedded Image. Since we look for ω>0, the above equality is possible only if (h(ν)/ν)<2. Thus,Embedded Image

As a consequence, we obtain the following lemma.

Lemma 6.1.

Assume that there exists λ*=iω* for some ω*>0 and Δ(λ*)=0. Then,Embedded Image with (h(ν)/ν)<2. Moreover, Δ(iω)≠0, Embedded Image.

By considering the real and the imaginary parts of equation (6.1) with ω>0, we obtain Embedded Image and Embedded Image. Since ω>0 and ν>0, we must have Embedded Image,Embedded Image 6.2 The definition (5.1) of h(ν) yieldsEmbedded Image 6.3 By using equations (6.2) and (6.3), we obtain a family of curves in the (ν,δ)-plane parameterized by Embedded Image (see figure 1).

Figure 1.

Bifurcation diagrams. (a) A couple of bifurcation curves given by equation (6.3) in the (ν,δ)-plane with small (ν,δ) values. (b) A family of bifurcation curves defined by equation (6.3) in the (ν,δ)-plane.

Proposition 6.2.

Assume that δ=cν for some constant c>1. Then, we can find a sequence {νn} with Embedded Image as Embedded Image such that, for each νn, there exists ωn>0 so that λn=iωn satisfies the characteristic equationEmbedded Image

Proof.

By using equations (6.2) and (6.3) for Embedded Image we haveEmbedded Image Since, by the assumption, we have δ=cν with c>1, it follows thatEmbedded Image A simple analysis shows that both curves intersect infinitely many times in Embedded Image. ■

Remark 6.3.

By remark 5.3 and proposition 6.2, we deduce that there is an infinite sequence {τn} with Embedded Image such that the linearized equation of system (3.1) around the endemic equilibrium Embedded Image has some purely imaginary eigenvalues in its spectrum. As a consequence, we deduce that the linearized equation of system (2.1) for the same values τn also has some purely imaginary eigenvalues. Therefore, when τ increases, the original system (2.1) may exhibit an infinite number of Hopf bifurcations.

In order to find an infinite number of purely imaginary eigenvalues, we consider the curves in the (ν,δ)-plane of parametersEmbedded Image 6.4 for some c≥1. Taking ν=ω/c, equations (6.2) and (6.3) become Embedded Image and Embedded Image for Embedded Image. We have the following lemma (figure 2).

Lemma 6.4.

There is an infinite sequence {ωn} defined by Embedded Image, such that, for n≥0 ifEmbedded Image then Δ(iωn)=0.

Similarly, for the caseEmbedded Image 6.5 we have Embedded Image, Embedded Image and the following lemma.

Lemma 6.5.

There is a second sequence Embedded Image given by Embedded Image, such that, for n≥0 ifEmbedded Image then Embedded Image.

Figure 2.

In this figure, we plot the curve Embedded Image for c=1 in red, for c=5 in blue and the curve Embedded Image for c=5 in green.

(b) Transversality condition

For Embedded Image with Re(λ)>−h(ν), the characteristic equation Δ(λ)=0 is equivalent to Δ1(ν,λ)=λ[λ+h(ν)]+νh(ν)(1−eλ)=0. Recall the form of h(ν), by equation (6.4), we have h(ν)=νκ(c), whereEmbedded Image Hence, Δ1(ν,λ)=λ[λ+νκ(c)]+ν2κ(c)(1−eλ). Moreover, we knew that if λn=ωni and νn=cωn, then Δ1(νn,λn)=0, and Embedded Image and Embedded Image. So Embedded Image. Since ∂Δ1(ν,λ)/∂λ=2λ+νκ(c)+ν2κ(c)eλ, it follows thatEmbedded Image which is impossible since both the real and the imaginary parts are strictly positive.

Next, we have ∂Δ1(ν,λ)/∂ν=λκ(c)+2νκ(c)[1−eλ]. Hence,Embedded Image Now, since ∂Δ1(νn,λn)/∂λ≠0, applying the implicit function theorem, we deduce that there exists a C1 map Embedded Image, such that Embedded Image and Embedded Image, ∀ν∈[νnεn,νn+εn]. Thus, we haveEmbedded Image Notice thatEmbedded Image and it follows from the fact νn=cωn thatEmbedded Image But, Embedded Image and κ(c)>0, so we haveEmbedded Image Set Embedded Image. Then,Embedded Image Since the map Embedded Image is strictly increasing on [0,1] and is equal to 1 for y=1, it implies that Embedded Image.

The above analysis can be summarized as the following lemma.

Lemma 6.6.

Assume that c≥1, then, for each n≥0, there exists a C1 map Embedded Image, such thatEmbedded Image

For the symmetric case (6.5), we have Embedded Image and Embedded Image whereEmbedded Image Similarly, we have the following lemma.

Lemma 6.7.

Assume that c≥1, then, for each n≥0, there exists a C1 map Embedded Image, such thatEmbedded Image

(c) Hopf bifurcation

By combining results on the essential growth rate of the linearized equations (equation (4.3)), the simplicity of the imaginary eigenvalues (lemma 4.9), the existence of purely imaginary eigenvalues (lemma 6.4) and the transversality condition (lemma 6.6), we are in a position to apply the centre manifold theorem 4.21 and proposition 4.22 in Magal & Ruan (2009a). Applying the Hopf bifurcation theorem proved in Hassard et al. (1981) to the reduced system, we obtain the following theorem.

Theorem 6.8.

Consider the curves defined by (6.4) in the (ν,δ)-plane for some c≥1. Then, for each n≥0,Embedded Image is a Hopf bifurcation point for the systemEmbedded Image 6.6 around the branch of equilibrium points Embedded Image. Moreover, the period of the bifurcating periodic orbits is close to Embedded Image.

By using the same arguments as above (i.e. combining equation (4.3), lemmas 4.9, 6.5 and 6.7), we also obtain the following theorem.

Theorem 6.9.

Consider the curves defined by equation (6.5 ) in the (ν,δ)-plane for some c≥1. Then, for each n≥0,Embedded Image is a Hopf bifurcation point for the system (6.6 ) around the branch of equilibrium points Embedded Image. Moreover, the period of the bifurcating periodic orbits is close to Embedded Image.

We would like to point out that we have studied the existence of Hopf bifurcation, while the stability of the bifurcating periodic solutions remains open.

Figure 3.

(a) When ρ=30 days, δ=0.0033976 and 1/v=5 in days, the trajectories converge to the positive equilibrium value. (b) There are sustained oscillations when ρ=100 days, δ=0.0016988 and 1/v=5 in days, via Hopf bifurcation.

Figure 4.

(a) The oscillations are amplified when ρ=200 days, δ=0.0033976 and 1/v=5 in days. (b) Irregular oscillations when ρ=365 days, δ=0.0033976 and 1/v=5 in days.

7. Numerical simulations

In this section, we present some numerical simulations for model (3.1). We consider the percentage of the population in a community. Recall that, by assumption 2.2, γ(ka) is equal to δ if aρ and 0 if 0<a<ρ, where ρ=τ/k. Here, we fix τ=30,ν=1/5 days−1, δ=0.0034, and let ρ vary from 0 to 365 days.

When ρ=30 days, there is an endemic equilibrium that is asymptotically stable (figure 3a). When ρ=100 days, the endemic equilibrium becomes unstable and there are periodic solutions via Hopf bifurcation (figure 3b). The oscillations are amplified when ρ=200 days (figure 4a) and irregular when ρ=365 days (figure 4b).

Recall that τ is the average number of amino acid substitutions before reinfection and k is the average number of amino acid substitutions per unit of time, that is, the mutation rate. Thus, the threshold ρ represents the time necessary to be reinfected after one infection. The simulations in figures 3 and 4 indicate that, for a given mutation rate, oscillations can start as early as 100 days after the previous infection and become irregular after about a year. It should be emphasized that ρ depends on the average number of amino acid substitutions after one infection and the mutation rate.

8. Discussion

Understanding the seasonal/periodic reoccurrence of influenza will be very helpful in designing successful vaccine programmes and introducing public-health interventions. However, the reasons for seasonal/periodic influenza epidemics are still not clear, even though various explanations have been proposed (Fine 1982; Dushoff et al. 2004; Levin et al. 2004).

Deterministic SIR models for some infectious diseases with permanent immunity have been shown to have an intrinsic tendency to oscillate periodically (Hethcote et al. 1981; Hethcote & Levin 1989; Ruan & Wang 2003; Dushoff et al. 2004; Tang et al. 2008, etc.). However, this type of models do not apply directly to influenza since immunity is not permanent for influenza. Dushoff et al. (2004) modelled antigenic drift by allowing individuals to lose their resistance to the circulating virus and re-enter the susceptible class after a few years, and used a periodic incidence function to describe seasonal changes. They showed that the large oscillations in incidence may be caused by undetectably small seasonal changes in the influenza transmission rate that are amplified by dynamical resonance. Castillo-Chavez et al. (1989), Andreasen et al. (1997), Lin et al. (1999), Andreasen (2003) and Nuño et al. (2005) developed multiple influenza A strain models with partial cross-immunity, and showed that such models exhibit sustained oscillations. Note that Castillo-Chavez et al. (1989) observed that sustained oscillations do not seem possible for a single-strain model, even with age-specific mortalities. For sustained oscillations to occur, they required at least an age-structured population and two or more co-circulating viral strains.

 Pease (1987) used an age-structured type model to describe the evolutionary epidemiological mechanism of influenza A drift, in which the susceptible class is continually replenished because the pathogen changes genetically and immunologically from one epidemic to the next, causing previously immune hosts to become susceptible. However, applied to influenza A, the model only predicted damped oscillations. Inaba (2002) modified the model of Pease (1987) and conjectured that, for realistic parameter values, the modified model may exhibit sustained oscillations due to Hopf bifurcation. In this paper, we first modified the model of Pease (1987) and Inaba (2002) by using the age of infection (instead of the number of amino acid substitutions) as a variable. Then, we applied our recently established centre manifold theory for semilinear equations with non-dense domain (Magal & Ruan 2009a) to show that Hopf bifurcation occurs in the modified model. This demonstrates that the age-structured type evolutionary epidemiological model of influenza A drift has an intrinsic tendency to oscillate owing to the evolutionary and/or immunological changes of the influenza viruses.

Acknowledgements

The authors would like to thank the two anonymous referees for helpful comments and suggestions. Research was partially supported by NSF grant DMS-0715772.

Footnotes

    • Received August 18, 2009.
    • Accepted November 4, 2009.

References

View Abstract