# 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 that 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: 2.1 where ν>0 is the recovery rate of the infected hosts and describes how amino acid substitutions affect the probability of reinfection.

First, we make the following assumptions.

## Assumption 2.1.

Assume that k>0, ν>0, and

From now on, we set

In this model, we do not have a disease-free equilibrium. Indeed, we observe that, if I0=0, then 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, It is clear that 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: 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. Define as the threshold of sensitivity, where τ will be discussed later. We make the following assumption.

## Assumption 2.2.

Assume that 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: 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 and as and

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

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 is also asymptotically stable.

From now on, we assume that τ>0 and make the following change of variable. By considering we obtain and Thus, with the above changes of variables, we obtain the following system: 2.3 where 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 and in system (2.2), we can assume, without loss of generality, that

## 3. Preliminary results

The global existence and positivity of solutions follow by considering the system 3.1 where

### (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 endowed with the usual product norm . Consider a linear operator defined by Then, Set and We also consider the nonlinear operator , defined by By identifying s(t,.) with , we can rewrite the system as the following abstract Cauchy problem: 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 is an integrated solution of the Cauchy problem (3.2), i.e. satisfies, for each t≥0, that and

### (b) Equilibrium solutions

Notice that We obtain the following lemma.

### Lemma 3.2.

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

If is an equilibrium solution of the system, we must have This is equivalent to with So we must have and obtain 3.3 By integrating over we obtain which implies that Thus, 3.4 Moreover, 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 the map is differentiable in a neighbourhood of and the linearized equation of (3.2) at is given by 4.1 We may rewrite this as the following PDE: 4.2

### (a) Associated linear operator

We consider a linear operator , defined by and a bounded linear operator , defined by Then, and the Cauchy problem (4.1) can be rewritten as 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 More precisely, for each Moreover, B is a Hille–Yosida operator and with and

Recall that B0, the part of B in is defined by and So, we obtain where with 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 by where is a strongly continuous semigroup of bounded linear operators on and Moreover,

### Definition 4.3.

Let be the infinitesimal generator of a linear C0-semigroup {TL(t)}t≥0 on a Banach space X. Define the growth bound of L by The essential growth bound of L is defined by where ∥TL(t)∥ess is the essential norm of TL(t) defined by Here, BX(0,1)={xX:∥xX≤1}, and for each bounded set 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 Since C is a bounded linear operator, it follows that 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 that 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 be the infinitesimal generator of a linear C0-semigroup {TL(t)} on a Banach space X. Then, Assume in addition that ω0,ess(L)<ω0(L). Then, for each γ∈(ω0,ess(L),ω0(L)], is non-empty, finite and contains only poles of the resolvent of L. Moreover, there exists a finite rank bounded linear projector 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 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, and each is a pole of the resolvent of L. That is, 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: where 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, is either empty or finite. Moreover, we have the following result.

### Lemma 4.6.

For each with where Furthermore, if with and Δ(λ)≠0, then is equivalent to where 4.4

### Proof.

For with since λIB is invertible, it follows that λI−(B+C) is invertible if and only if IC(λIB)−1 is invertible and Note that 4.5 where So we must have By using the second equation of system (4.5), we obtain Thus, we have 4.6 where and is given by equation (4.3). Moreover, 4.7 where Combining the first equation of system (4.5) with equations (4.6) and (4.7), we obtain and So, and The characteristic function is given by It follows that 4.8 By using the second equation of (4.5), together with equations (4.7) and (4.8), we obtain Hence, Finally, we deduce that 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 D(A), such that which is, in turn, equivalent to If we fix then we obtain Since Δ(λ)=0, we deduce that The result follows. ■

### Lemma 4.7.

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

### Proof.

We have Δ(0)=1+νJ1J2, with and But, so we have which implies that Δ(0)=1+ν>0. This completes the proof. ■

### Lemma 4.8.

For each with and λ≠0, we have

### Proof.

We have Δ(λ)=1+νI1I2, with and The result thus follows. ■

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

### Lemma 4.9.

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

### Proof.

Assume that with ω>0 is a solution of Assume first that 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: Moreover, we have So, is a pole of order 1 of the resolvent of B+C, and the projector on the generalized eigenspace of B+C associated to has rank 1. It follows that is a simple eigenvalue of B+C.

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

## 5. Stability

In this section, we investigate the stability of the endemic equilibrium . We observe that the linear operator can be decomposed as the sum of a linear operator , defined by and a bounded linear operator , defined by That is, we have . Observe that . In order to have the existence of the endemic equilibrium (i.e. or ), we must impose δ>ν. Fix δ>0. Then, we have as . Hence, as . Moreover, , the part of on X0, is the infinitesimal generator of the strongly continuous semigroup of bounded linear operators on X0, defined by where Since it follows that 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 , such that, for each the equilibrium 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 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 and of system (3.1) with and 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 . Let be defined by 5.1 Then, the characteristic equation can be rewritten as 1+ν[h(ν)/λ(λ+h(ν))](1−eλ)=0. Since Re(λ)>−h(ν), the characteristic equation is equivalent to 5.2

## Lemma 5.4.

Assume that 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 , 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 , such that, for each the equilibrium of system (3.1) is locally asymptotically stable.

## Proof.

Assume that there exist two sequences {νn} with and 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 . By using the characteristic equation, we obtain and . So λ*=0. Notice that so we obtain . Thus, for p≥0 large enough, we have Re(λnp/h(νnp))<−(1+δ/2), which implies, for each p≥0 large enough, that a contradiction. ■

## 6. Hopf bifurcation

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

### (a) Existence of purely imaginary solutions

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

As a consequence, we obtain the following lemma.

### Lemma 6.1.

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

By considering the real and the imaginary parts of equation (6.1) with ω>0, we obtain and . Since ω>0 and ν>0, we must have , 6.2 The definition (5.1) of h(ν) yields 6.3 By using equations (6.2) and (6.3), we obtain a family of curves in the (ν,δ)-plane parameterized by (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 as such that, for each νn, there exists ωn>0 so that λn=iωn satisfies the characteristic equation

### Proof.

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

### Remark 6.3.

By remark 5.3 and proposition 6.2, we deduce that there is an infinite sequence {τn} with such that the linearized equation of system (3.1) around the endemic equilibrium 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 parameters 6.4 for some c≥1. Taking ν=ω/c, equations (6.2) and (6.3) become and for . We have the following lemma (figure 2).

### Lemma 6.4.

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

Similarly, for the case 6.5 we have , and the following lemma.

### Lemma 6.5.

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

Figure 2.

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

### (b) Transversality condition

For 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), where Hence, Δ1(ν,λ)=λ[λ+νκ(c)]+ν2κ(c)(1−eλ). Moreover, we knew that if λn=ωni and νn=cωn, then Δ1(νn,λn)=0, and and . So . Since ∂Δ1(ν,λ)/∂λ=2λ+νκ(c)+ν2κ(c)eλ, it follows that which is impossible since both the real and the imaginary parts are strictly positive.

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

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 , such that

For the symmetric case (6.5), we have and where Similarly, we have the following lemma.

### Lemma 6.7.

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

### (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, is a Hopf bifurcation point for the system 6.6 around the branch of equilibrium points . Moreover, the period of the bifurcating periodic orbits is close to .

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, is a Hopf bifurcation point for the system (6.6 ) around the branch of equilibrium points . Moreover, the period of the bifurcating periodic orbits is close to .

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

• Accepted November 4, 2009.

View Abstract