## 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 2009*a*) 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 (2009*a*) 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 *k**a*. 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 *k**a*_{0} and less than *k**a*_{1} 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 *I*_{0}=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*)]/*d**t*=0. It follows that *I*(*t*)=*N*−*S*(*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 *γ*(*k**a*) is a linear function. To consider such unbounded functions, we must change the state space by introducing the weighted *L*^{1} space in order to discuss the existence of solutions. Since this is not the goal of this work, we make a simplified assumption that *γ*(*k**a*) 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 *γ*(*k**a*) describes the rate at which individuals become reinfected per average number of amino acid substitution *k**a*. 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 *γ*(*k**a*) 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)=*N*−*S*(0)>0,*then**and**as**and*— if

*ν*/*δ*<*N*,*then system (2.1) has a unique equilibrium solution**such that*,*and**Moreover, if I*(0)=*N*−*S*(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 (1990*b*, 1997), Arendt *et al.* (2001) and Magal & Ruan (2007, 2009*a*,*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 (1990*a*), Magal (2001) and Magal & Ruan (2007, 2009*a*). 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 X_{0+}, such that, for each x∈X_{0+}, 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 (1990*a*) 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 *B*_{0}, the part of *B* in is defined by
and
So, we obtain
where with Since *B* is a Hille–Yosida operator, it follows that *B*_{0} is the infinitesimal generator of a strongly continuous semigroup of bounded linear operators {*T*_{B0}(*t*)}_{t≥0} on *X*_{0}. More precisely, we have the following result.

### Lemma 4.2.

*The linear operator B_{0} is the infinitesimal generator of a strongly continuous semigroup of bounded linear operators {T_{B0}(t)}_{t≥0} on X_{0}, and for each t≥0, the linear operator T_{B0}(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 *C*_{0}-semigroup {*T*_{L}(*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 ∥*T*_{L}(*t*)∥_{ess} is the essential norm of *T*_{L}(*t*) defined by
Here, *B*_{X}(0,1)={*x*∈*X*:∥*x*∥_{X}≤1}, and for each bounded set
is the Kuratovsky measure of non-compactness.

Since ∥*T*_{B}(*t*)∥_{ess}≤∥*T*_{B}(*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 *X*_{0}, is the infinitesimal generator of a strongly continuous semigroup of bounded linear operators {*T*_{(B+C)0}(*t*)}_{t≥0} on *X*_{0}. 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 C_{0}-semigroup {T_{L}(t)} on a Banach space X. Then*,

*Assume in addition that*(i )

*ω*_{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:*Π*(

*λ*−

*L*)

^{−1}=(

*λ*−

*L*)

^{−1}

*Π*,∀

*λ*∈

*ρ*(

*L*), (ii )

*σ*(

*L*

_{Π(X)})={

*λ*∈

*σ*(

*L*):

*R*

*e*(

*λ*)≥

*γ*},

*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 *R**e*(*λ*)≥*γ*. As a consequence of theorem 4.4, we have the following corollary.

### Corollary 4.5.

*Let be the infinitesimal generator of a linear C_{0}-semigroup {T_{L}(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*k*_{0}≥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 (2009*a*), 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 *λ**I*−*B* is invertible, it follows that *λ**I*−(*B*+*C*) is invertible if and only if *I*−*C*(*λ**I*−*B*)^{−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+*ν**J*_{1}−*J*_{2}, 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+*ν**I*_{1}−*I*_{2}, 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 *X*_{0}, is the infinitesimal generator of the strongly continuous semigroup of bounded linear operators on *X*_{0}, defined by
where
Since it follows that By using the standard perturbation argument (see Thieme 1990*a* or Magal & Ruan 2009*b*), 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 *c*_{1}>*c*_{2}>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 *R**e*(*λ*)>−*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

*ν*+2

*h*(

*ν*).

## 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*(*ν*))(e^{−Re(λ)} + 1) + |(*λ*/*h*(*ν*)) + 1|. Since *R**e*(*λ*) ≥ 0, we obtain that , it follows that |*λ*+*h*(*ν*)|≤2*ν*+*h*(*ν*) and |*λ*|≤|*λ*+*h*(*ν*)|+|*h*(*ν*)|≤2*ν*+2*h*(*ν*). 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 *R**e*(*λ*_{n})≥0 and Δ(*λ*_{n})=0. Then, by lemma 5.4, we have |*λ*_{n}|≤2*ν*_{n}+2*h*(*ν*_{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 *R**e*(*λ*_{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 *R**e*(*λ*)>−*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).

### 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 .

### (b) Transversality condition

For with *R**e*(*λ*)>−*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}=*ω*_{n}*i* 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 *C*^{1} 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 C^{1} 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 C^{1} 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 (2009*a*). 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.

## 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, *γ*(*k**a*) 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 3*a*). When *ρ*=100 days, the endemic equilibrium becomes unstable and there are periodic solutions via Hopf bifurcation (figure 3*b*). The oscillations are amplified when *ρ*=200 days (figure 4*a*) and irregular when *ρ*=365 days (figure 4*b*).

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 2009*a*) 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.

- © 2009 The Royal Society