Hopf bifurcation of an epidemic model with a nonlinear birth in population and vertical transmission

Applied Mathematics and Computation 230 (2014) 164–173

Contents lists available at ScienceDirect

Applied Mathematics and Computation journal homepage: www.elsevier.com/locate/amc

Hopf bifurcation of an epidemic model with a nonlinear birth in population and vertical transmission q Zhang Yinying, Jia Jianwen ⇑ School of Mathematics and Computer Science, Shanxi Normal University, Shanxi, Linfen 041004, PR China

a r t i c l e

i n f o

Keywords: Nonlinear birth Vertical transmission Hopf bifurcation Normal theory Center manifold theorem

a b s t r a c t In this paper, an epidemic model involving a nonlinear birth in population and vertical transmission was studied. When R0 < 1, the disease-free equilibrium was stable, while if R0 > 1, the disease-free equilibrium was unstable. We researched the existence of Hopf bifurcation and obtained the stability and direction of the Hopf bifurcation by using the normal theory and the center manifold theorem. Numerical simulations were carried out to illustrate the main theoretical results and a brief discussion was given to conclude this work. Ó 2013 Elsevier Inc. All rights reserved.

1. Introduction Mathematical epidemiology, i.e. the construction and analysis of mathematical models are one of the major areas of biology to describe the spread and control of infectious diseases. Since Kermack and Mckendrick constructed a system of ODE [1] to study epidemiology in 1927, the concept of ‘‘Compartment model’’ have been used until now. Most of the research literatures described the spread of a non-lethal disease in a large population by dividing the total population into three classes: the susceptible (S), the infectious (I), and the recovered (i.e. with a permanent or temporary acquired immunity) (R). We usually call these compartmental models SIR models or SIRS models, and each letter in SIR or SIRS denotes a compartment. And there must be an individual belong to one compartment. In real life, some diseases may be passed from one individual to another via vertical transmission. That is to say, vertical transmission of diseases refer to diseases are infected to the offspring by their infected parentage. In recent years, a few studies of vertical transmission have been conducted to describe the effects of various and demographical factors [2–5]. For example, Busenberg and Cooke [5] discussed a variety of diseases which contained both of horizontally and vertically transmitting, a comprehensive and formulation survey. They also provided the mathematical analysis of compartmental models including vertical transmission. Some examples of such diseases are AIDS, Rubella, Hepatitis, etc. 2. The model Classical epidemic models assume that the total population size is constant, and that concentrate on describing the spread of disease based on the population. In recent years more and more models pay attention to a variable population size, and then disease causing death for a longer time scale should be taken into account to reduce reproduction. For example, according to the paper [6], a nonlinear birth term BðNÞ is considered and we can also find that the form BðNÞN is important

q

This work is supported by Natural Science Foundation of Shanxi province (2013011002-2).

⇑ Corresponding author.

E-mail address: [email protected] (J. Jianwen). 0096-3003/$ - see front matter Ó 2013 Elsevier Inc. All rights reserved. http://dx.doi.org/10.1016/j.amc.2013.12.084

165

Z. Yinying, J. Jianwen / Applied Mathematics and Computation 230 (2014) 164–173

in determining the qualitative dynamics. In the absence of disease, the paper [7] assumes that the total population size N changes according to a population growth equation

dNðtÞ ¼ BðNÞN  dN: dt Here d > 0 is the death rate constant, and BðNÞN is a birth rate function with BðNÞ satisfying with following basic assumptions for N 2 ð0; 1Þ: (A1) BðNÞ > 0; (A2) BðNÞ is continuously differentiable with B0 ðNÞ < 0; (A3) Bð0þ Þ > d > Bð1Þ. Based on the above description, we choose a BðNÞ ¼ NA þ B. It is clear that BðNÞ meets (A1) and (A2). Because of (A3), we can choose B; l ¼ minfl1 ; l2 ; l3 g satisfying B < l. Under these assumptions, the following epidemic model is considered as one model with nonlinear birth in population and nonlinear incidence.

8 dSðtÞ ¼ A  ðl1  BÞSðtÞ þ BqIðtÞ þ BRðtÞ  bSðtÞIðtÞ þ cel3 s Iðt  sÞ; > 1þaIðtÞ > < dt dIðtÞ bSðtÞIðtÞ ¼ 1þaIðtÞ þ BpIðtÞ  ðc þ l2 ÞIðtÞ; dt > > : dRðtÞ ¼ cIðtÞ  l3 RðtÞ  cel3 s Iðt  sÞ; dt

ð2:1Þ

where A; B; li ði ¼ 1; 2; 3Þ; q; p; b; a; c are nonnegative. SðtÞ; IðtÞ; RðtÞ denote the number of susceptible, infected and recover population stage at time t, respectively. A is constant immigrants, B is the birth rate and li ði ¼ 1; 2; 3Þ are natural death rate. q is the probability that a child who is born from infectious mother is susceptible; p is the probability that a child who is born from infectious mother is infected, then p þ q ¼ 1. b is contact rate between the susceptible and the infection. c is the recovery rate. s is the time delay. The initial conditions for system (2.1) are

ðw1 ðhÞ; w2 ðhÞ; w2 ðhÞÞ 2 C þ ¼ Cð½s; 0; R3þ Þ;

wi ð0Þ > 0; i ¼ 1; 2; 3;

ð2:2Þ

where

R3þ ¼ fðx1 ; x2 ; x3 Þ 2 R3 : xi P 0; i ¼ 1; 2; 3g:

Theorem 2.1. For any solution SðtÞ; IðtÞ; RðtÞ of system (2.1) with initial conditions (2.2), SðtÞ < M, IðtÞ < M; RðtÞ < M for all large A t, where M ¼ lB . h This paper is organized as following: in the next section, we obtain the basic reproduction number by the next generation method and the existence of equilibriums. We verify when R0 < 1, the disease-free equilibrium was stable, while if R0 > 1, the disease-free equilibrium was unstable. Then we focus on the local stability of the endemic equilibrium and the existence of the Hopf bifurcation. In Section 4, we obtain the stability and direction of the Hopf bifurcation by using the normal theory and the center manifold theorem. Numerical simulations are carried out in Section 5 to illustrate the main theoretical results and a brief discussion is given in last part to conclude this work. 3. The existence and stability of equilibria 3.1. The existence of equilibria and the stability of the disease-free equilibrium   It is easy to obtain the disease-free equilibrium E0 ¼ ðS0 ; 0; 0Þ ¼ l AB ; 0; 0 . Then, we define the basic reproduction num1 ber R0 of our model by directly using the next generation method presented in Diekmann et al. [8] and van den Driessche and Watmough [9]. Then

R0 ¼

bA : ðl1  BÞðl2 þ c  BpÞ

Theorem 3.1. When R0 < 1, there exists the unique disease-free equilibrium E0 . When R0 > 1, there also exist a endemic equilibrium

E ¼ ðS ðsÞ; I ðsÞ; R ðsÞÞ,

where

½bAðl1 BÞðl2 þcBpÞl3 al3 ðl1 BÞðl2 þcBpÞþbcBð1þel3 s Þþbl3 ðl2 BÞþbcl3 ð1el3 s Þ.

S ðsÞ ¼

ðl2 þcBpÞð1þaI Þ , b



R ðsÞ ¼ cbI3 ð1  el3 s Þ

and

I ðsÞ ¼

166

Z. Yinying, J. Jianwen / Applied Mathematics and Computation 230 (2014) 164–173

When

½bAðl1 BÞðl2 þcBpÞl3 s ¼ 0; I ð0Þ ¼ al3 ðl1 BÞð l2 þcBpÞþ2bcBþbl3 ðl2 BÞ.

Then we find the characteristic equation of any equilibrium E ¼ ðS; I; RÞ as follows:

   k þ l1  B þ 1þbIaI   detðkI  JðEÞÞ ¼   1þbIaI    0

bS ð1þaIÞ

2

 Bq  ceðl3 þkÞs

k þ l2 þ c  Bp 

bS ð1þaIÞ

2

ceðl3 þkÞs  c

  B     0 ¼0   k þ l3 

ð3:1Þ

Theorem 3.2. When R0 < 1, the disease-free equilibrium E0 is locally asymptotically stable for all unstable for all s P 0. h

s P 0; When R0 > 1; E0 is

3.2. The locally stability of the endemic equilibrium and the Hopf bifurcation In this section we will analyze the local stability of the positive equilibrium E . It is not easy to find rigorously local stability condition of E . In the following, using stability switch criteria, we try to analyze local stability of E . The criteria can indicate that the stability of a given steady state is simply determined by the graph which is expressed as some functions of s and thus can be easily depicted by Matlab. We will show the stability switch criteria for E . The characteristic of the system (2.1) near the infected equilibrium E is given by

Pðk; sÞ þ Q ðk; sÞeks ¼ 0;

ð3:1Þ

where

Pðk; sÞ ¼ k3 þ A1 ðsÞk2 þ A2 ðsÞk þ A3 ðsÞ;

ð3:2Þ

Q ðk; sÞ ¼ A4 ðsÞk þ A5 ðsÞ and

bI ðsÞ bS ðsÞ ;   1 þ aI ðsÞ ð1 þ aI ðsÞÞ2 !   bS ðsÞ bI ðsÞ bI ðsÞ þ l3 l1  B þ þ ðl2 þ c  BÞ ; A2 ðsÞ ¼ ðl1 þ l3  BÞ l2 þ c  Bp   2  1 þ aI ðsÞ 1 þ aI ðsÞ ð1 þ aI ðsÞÞ ! bS ðsÞ bI ðsÞ cbI ðsÞ þ A3 ðsÞ ¼ l3 ðl1  BÞ l2 þ c  Bp  l ð l  BÞ l  BÞ ; þ ð  3 2 3 2 1 þ aI ðsÞ 1 þ aI ðsÞ ð1 þ aI ðsÞÞ

A1 ðsÞ ¼ l1 þ l2 þ l3 þ c  Bp  B þ

A4 ð s Þ ¼ 

cbI ðsÞ l3 s cbI ðsÞ l3 s e ; A5 ðsÞ ¼ ðl3  BÞ e :  1 þ aI ðsÞ 1 þ aI ðsÞ

Theorem 3.3. When R0 > 1; E is locally asymptotically stable for Proof. When

s ¼ 0.

s ¼ 0, the characteristic equation at E is of the form:

k3 þ m1 ð0Þk2 þ m2 ð0Þk þ m3 ð0Þ ¼ 0; where

m1 ð0Þ ¼ A1 ð0Þ; m2 ð0Þ ¼ A2 ð0Þ þ A4 ð0Þ ¼ ðl1 þ l3  BÞ m3 ð0Þ ¼ A3 ð0Þ þ A5 ð0Þ ¼ l3 ðl1  BÞ For the fact that

M1 ¼ m1 ð0Þ;

bS ð0Þ ð1þaI ð0ÞÞ2

l2 þ c  Bp 

l2 þ c  Bp 

!

bS ð0Þ ð1 þ aI ð0ÞÞ ! bS ð0Þ 

2

ð1 þ aI ð0ÞÞ

2

þ l3



l1  B þ

þ l3 ðl2  BÞ

 bI ð0Þ bI ð0Þ þ ðl2  BÞ ;  1 þ aI ð0Þ 1 þ aI ð0Þ

bI ð0Þ ; 1 þ aI ð0Þ

2 þcBp ¼ l1þ aI ð0Þ < l2 þ c  Bp, we can obtain mk ð0Þ > 0ðk ¼ 1; 2; 3Þ. Then

M2 ¼ m1 ð0Þm2 ð0Þ  m3 ð0Þ:

We obtain Mk > 0ðk ¼ 1; 2Þ by directly computing. From the well-known Routh–Hurwitz criterion, we can obtain that when R0 > 1, the endemic equilibrium E is locally asymptotically stable for all s ¼ 0. h

167

Z. Yinying, J. Jianwen / Applied Mathematics and Computation 230 (2014) 164–173

In the following, we investigate the existence of purely imaginary roots k ¼ ixðx > 0Þ to Eq. (3.1). For simplify, we will use the abbreviations S ; I ; R instead of S ðsÞ; I ðsÞ; R ðsÞ. Eq. (3.1) takes the form of a third-degree exponential polynomial in k, with all the coefficients of P and Q depending on s. Beretta and Kuang [10] established a geometrical criterion which gives the existence of purely imaginary root of a characteristic equation with delay dependent coefficients. In order to apply the criterion due to Beretta and Kuang [10], we need to verify the following properties for all s 2 ½0; smax Þ, where smax is the maximum value which E exists. (a) Pð0; sÞ þ Q ð0; sÞ – 0. (b) Pðix; sÞ þ Q ðix; sÞ – 0. ðk;sÞ (c) lim supfj QPðk; sÞ j : jkj ! 1; Rek P 0g < 1.

(d) Fðx; sÞ ¼ jPðix; sÞj2  jQ ðix; sÞj2 has a finite number of zeros. (e) Each positive root xðsÞ of Fðx; sÞ ¼ 0 is continuous and differentiable in

s whenever it exists.

Here, Pðk; sÞ and Q ðk; sÞ are defined as in (3.2). Let s 2 ½0; smax Þ. It is easy to see that Pðix; sÞ þ Q ðix; sÞ ¼ A3 ðsÞ þ A5 ðsÞ – 0. This implies that (a) is satisfied. And (b) is obviously true because

    Pðix; sÞ þ Q ðix; sÞ ¼ A1 ðsÞx2 þ A3 ðsÞ þ A5 ðsÞ þ ix x2 þ A2 ðsÞ þ A4 ðsÞ – 0: From (3.2) we know that

          Q ðk; sÞ A4 ðsÞk þ A5 ðsÞ A4 ðsÞ     ¼ lim  lim  ¼ lim  ¼0   3 2 2 jkj!1 Pðk; sÞ jkj!1k þ A1 ðsÞk þ A2 ðsÞk þ A3 ðsÞ jkj!13k þ 2A1 ðsÞk þ A2 ðsÞ Therefore (c) follows. Let F be defined as in (d). From 2

2

jPðix; sÞj2 ¼ ðx2 þ A2 ðsÞxÞ þ ðA1 ðsÞx2 þ A3 ðsÞÞ ¼ x6 þ ½A21 ðsÞ  2A2 ðsÞx4 þ ½A22 ðsÞ  2A1 ðsÞA3 ðsÞx2 þ A23 ðsÞ and

jQðix; sÞj2 ¼ A24 ðsÞx2 þ A25 ðsÞ; we have

Fðx; sÞ ¼ x6 þ pðsÞx4 þ qðsÞx2 þ rðsÞ; where pðsÞ ¼ A21 ðsÞ  2A2 ðsÞ; qðsÞ ¼ A22 ðsÞ  2A1 ðsÞA3 ðsÞ  A24 ðsÞ; rðsÞ ¼ A23 ðsÞ  A25 ðsÞ. Let x2 ¼ h, then we have that 3

2

FðhÞ ¼ h þ pðsÞh þ qðsÞh þ rðsÞ ¼ 0:

ð3:3Þ

Theorem 3.4. Suppose that (3.3) has no positive roots, then when R0 > 1; E is locally asymptotically stable for all s P 0. h Suppose that (3.3) has positive roots, then (d) is satisfied. Let ðx0 ; s0 Þ be a point of its domain of definition such that Fðx0 ; s0 Þ ¼ 0. We know the partial derivatives F x and F s exist and are continuous in a certain neighborhood of ðx0 ; s0 Þ, and F x ðx0 ; s0 Þ – 0. By Implicit Function Theorem, (e) is also satisfied. Now let k ¼ ixðx > 0Þ be a root of Eq. (3.1), and from which we have that

ix3  A1 ðsÞx2 þ iA2 ðsÞx þ A3 ðsÞ þ ½iA4 ðsÞx þ A5 ðsÞeixs ¼ 0: Hence, we have that

x3  A2 ðsÞx ¼ A4 ðsÞx cosðxsÞ  A5 ðsÞ sinðxsÞ; A1 x2  A3 ðsÞx ¼ A4 ðsÞx sinðxsÞ þ A5 ðsÞ cosðxsÞ:

ð3:4Þ

From (3.4) it follows that

sinðxsÞ ¼

cosðxsÞ ¼

½A1 ðsÞA4 ðsÞ  A5 ðsÞx3 þ ½A2 ðsÞA5 ðsÞ  A3 ðsÞA4 ðsÞx A24 ðsÞx2 þ A25 ðsÞ A4 ðsÞx4 þ ½A1 ðsÞA5 ðsÞ  A2 ðsÞA4 ðsÞx2  A3 ðsÞA5 ðsÞ A24 ðsÞx2 þ A25 ðsÞ

;

ð3:5aÞ

:

ð3:5bÞ

By the definitions of Pðk; sÞ; Q ðk; sÞ as in (3.2), and applying the property (a), (3.5a) and (3.5b) can be written as

168

Z. Yinying, J. Jianwen / Applied Mathematics and Computation 230 (2014) 164–173

sinðxsÞ ¼ Im

Pðix; sÞ Q ðix; sÞ

ð3:6aÞ

and

cosðxsÞ ¼ Re

Pðix; sÞ ; Q ðix; sÞ

ð3:6bÞ

which yields

jPðix; sÞj2 ¼ jQ ðix; sÞj2 : Assume that I 2 Rþ0 is the set where xðsÞ is a positive root of

Fðx; sÞ ¼ jPðix; sÞj2  jQ ðix; sÞj2 and for s: 2 I , xðsÞ is not defined. Then for all

s in I, xðsÞ satisfied

Fðx; sÞ ¼ 0:

ð3:7Þ

Assume that Eq. (3.3) has pffiffiffiffiffi ffi only one positive real root, we denote by hþ this positive real root. Thus, Eq. (3.7) has only one positive real root x ¼ hþ . And the critical values of s and xðsÞ are impossible to solve explicitly, so we shall use the procedure described in Beretta and Kuang [10]. According to this procedure, we define hðsÞ 2 ½0; 2pÞ such that sin hðsÞ and cos hðsÞ are given by the right-hand sides of (3.5a) and (3.5b), respectively, with hðsÞ given by (3.6) And the relation between the argument h and xs in (3.6) for s > 0 must be

xs ¼ h þ 2np; n ¼ 0; 1; 2; . . . Hence we can define the maps

sn ðsÞ :¼

hðsÞ þ 2np ; xðsÞ

ð3:8Þ

sn : I ! Rþ0 given by

sn > 0; n ¼ 0; 1; 2; . . . ;

where a positive root xðsÞ of (3.3) exists in I. Let us introduce the functions Sn ðsÞ : I ! R,

Sn ðsÞ ¼ s 

hðsÞ þ 2np ; xðsÞ

n ¼ 0; 1; 2; . . . ;

that are continuous and differentiable in s. Thus, we give the following theorem which is due to Beretta and Kuang [10]. Theorem 3.5. Assume that xðsÞ is a positive root of (3.1) defined for s 2 I, I # Rþ0 , and at some s 2 I; Sn ðs Þ ¼ 0 for some n 2 N 0 . Then a pair of simple conjugate pure imaginary roots k ¼ ix exists at s ¼ s which crosses the imaginary axis from left to right if dðs Þ > 0 and crosses the imaginary axis from right to left if dðs Þ < 0, where



 dSn ðsÞ : dðs Þ ¼ signfF 0x ðxs ; s Þgsign  ds s¼s Applying Theorems 3.3, 3.4 and 3.5 and the Hopf bifurcation theorem for functional differential equation [11], we can conjecture the existence of a Hopf bifurcation as stated in Theorem 3.6. Theorem 3.6 (A conjecture). For system (2.1), then there exists s 2 I, such that the equilibrium E is asymptotically stable for 0 6 s < s , and becomes unstable for s staying in some right neighborhood of s , with a Hopf bifurcation occurring when s ¼ s . 4. Direction and stability of the Hopf bifurcation In this section, we shall study the direction of the Hopf bifurcations and stability of bifurcating periodic solutions by applying the normal theory and the center manifold theorem introduced by Hassard et al. [12]. e i ðtÞ ¼ ui ðstÞði ¼ 1; 2; 3Þ, s ¼ m þ s and dropping the bars for simplification of Let u1 ¼ S  S ; u2 ¼ I  I ; u3 ¼ R  R ; u notations, system (2.1) becomes an functional differential equation in C ¼ Cð½1; 0; R3 Þ as

8 h     i   bu1 u2 þabI u1 u2 abS u22 du1 ðtÞ l3 ðmþs Þ > > ¼ ðs þ mÞ  l1  B þ 1þbIaI u1 ðtÞ þ Bq  ð1þbSaI Þ2 u2 ðtÞ þ Bu3 ðtÞ þ ð1þ þ c e u 2 ðt  1Þ ;  2  > dt aI Þ ð1þaI þau2 Þ > < h     i abS u22 abI u1 u2 bu1 u2 du2 ðtÞ bI bS  ¼ ð s þ m Þ þ Bp  l  c ðtÞ þ u ;  u1 ðtÞ þ 2  2  2  2 1þ a I dt > ð1þ a I Þ ð1þ a I Þ ð1þ a I þ a u Þ 2 > > > du3 ðtÞ    : ¼ ðs þ mÞ cu2 ðtÞ  l3 u3 ðtÞ  cel3 ðmþs Þ u2 ðt  1Þ : dt Expanding the nonlinear part by Taylor expansion, we obtain

169

Z. Yinying, J. Jianwen / Applied Mathematics and Computation 230 (2014) 164–173

8 h  i    abS u22 ðtÞ abu1 ðtÞu22 ðtÞ du1 ðtÞ bu1 u2 bI bS   l ð mþ s  Þ > > u2 ðt  1Þ ; > dt ¼ ðs þ mÞ  l1  B þ 1þaI u1 ðtÞ þ Bq  ð1þaI Þ2 u2 ðtÞ þ Bu3 ðtÞ  ð1þaI Þ2 þ ð1þaI Þ3 þ ð1þaI Þ3 þ    þ ce 3 > < h  i    abS u22 ðtÞ abu1 ðtÞu22 ðtÞ du2 ðtÞ 1 ðtÞu2 ðtÞ ¼ ðs þ mÞ 1þbIaI u1 ðtÞ þ ð1þbSaI Þ2 þ Bp  l2  c u2 ðtÞ þ buð1þ þ    ;  2   3   3 dt > aI Þ ð1þaI Þ ð1þaI Þ > > > : du3 ðtÞ   l 3 ð mþ s  Þ ¼ ðs þ mÞ½cu2 ðtÞ  l3 u3 ðtÞ  ce u2 ðt  1Þ: dt

We denote the above system by

_ uðtÞ ¼ Lm ðut Þ þ f ðm; ut Þ;

ð4:1Þ T

where uðtÞ ¼ ðu1 ðtÞ; u2 ðtÞ; u3 ðtÞÞ 2 R3 , and Lm : C ! R3 , f : R  C ! R3 are given, respectively, by

0 B Lm ð/Þ ¼ ðs þ mÞB @



ðl1  B þ 1þbIaI Þ bI 1þaI

0

10 1 0 10 1  / ð0Þ /1 ð1Þ 0 cel3 ðs þmÞ 0 CB 1 C B C B C   bS 0 0 A@ /2 ð1Þ A þ Bp  ðl2 þ cÞ 0 C A@ /2 ð0Þ A þ ðs þ mÞ@ 0 ð1þaI Þ2  /3 ð0Þ /3 ð1Þ 0 cel3 ðs þmÞ 0 c l3 

Bq  ð1þbSaI Þ2

B

ð4:2Þ and

80 9 1 0 abS /22 ð0Þ 1 0 ab/1 ð0Þ/22 ð0Þ 1 1 ð0Þ/2 ð0Þ > >  b/ð1þ > >  2  3  3 > > ð1þaI Þ aI Þ > @ ð1þaI Þ3 A @ ð1þaI Þ > > ð1þaI Þ3 A > > : ; 0 0 0

ð4:3Þ

By the Riesz representation theorem, there exists a function gðh; mÞ of bounded variation for h 2 ½1; 0, such that

Lm ð/Þ ¼

Z

0

dgðh; mÞ/ðhÞ;

for / 2 C:

ð4:4Þ

1

In fact, we can choose

0     l1  B þ 1þbIaI B bI gðh; mÞ ¼ ðs þ mÞB B 1þaI @ 0



Bq  ð1þbSaI Þ2 bS ð1þaI Þ2

þ Bp  ðl2 þ cÞ

l3

c

1

0  0 cel3 ðs þmÞ C C B  dðhÞ  ðs þ mÞ@ 0 0 0 C A  0 cel3 ðs þmÞ B

0

1

C 0 Adðh þ 1Þ;

0 ð4:5Þ

3

where d denote the Dirac delta function. For / 2 Cð½1; 0; R Þ, define

( d/ðhÞ dh

AðmÞ/ ¼

R0

1

h 2 ½1; 0Þ;

;

dgðh; mÞ/ðhÞ; h ¼ 0

and

RðmÞð/Þ ¼



0;

h 2 ½1; 0Þ;

f ðm; /Þ; h ¼ 0:

Then, system (4.1) is equivalent to

_ uðtÞ ¼ AðmÞut þ RðmÞut ;

ð4:6Þ

where ut ¼ uðt þ hÞ for h 2 ½1; 0.  For w 2 Cð½0; 1; ðR3 Þ Þ, define

( 

A wðsÞ ¼

 dwðsÞ ; s 2 ð0; 1; ds R0 wðtÞdgðt; 0Þ; s ¼ 0: 1

and a bilinear inner product

hwðsÞ; /ðhÞi ¼ wð0Þ/ð0Þ 

Z

0

1

Z

h

wðr  hÞdgðhÞ/ðrÞdr;

ð4:7Þ

r¼0

where gðhÞ ¼ gðh; 0Þ. Then Að0Þ and A are adjoint operators. By the discussion in Section 3.2, we know that in s are eigenvalues of Að0Þ. Hence, they are also eigenvalues of A . We first need to compute the eigenvectors of Að0Þ and A corresponding to in s and in s , respectively.   Assume that qðhÞ ¼ ð1; q1 ; q2 ÞT ein s h is the eigenvector of Að0Þ corresponding to in s , then Að0ÞqðhÞ ¼ in s qðhÞ. Then from   the definition of Að0Þ and (4.2), (4.4) and (4.5), and for qð1Þ ¼ qð0Þein s , we have

170

Z. Yinying, J. Jianwen / Applied Mathematics and Computation 230 (2014) 164–173

0 B B B @



l2  B þ 1þbIaI

bS ð1þaI Þ2

bI 1þaI

B

þ Bp  ðl2 þ cÞ 

c  ceðl3 þin Þs

0

1

0 1 0 1 1 1 C CB C C B 0 C@ q1 ð0Þ A ¼ in @ q1 ð0Þ A: A q2 ð0Þ q2 ð0Þ l3

   Bq  ð1þbSaI Þ2 þ ceðl3 þin Þs

Then we obtain

8 < q1 ¼ ðin þl

bI ð1þaI Þ

2 þcBpÞð1þaI

:q ¼ 2

 2

Þ bS

;

:   bcI ð1þaI Þð1eðl3 þin Þs Þ ðin þl3 Þ½ðin þl2 þcBpÞð1þaI Þ2 bS 

  Similarly, we can calculate the eigenvector q ðsÞ ¼ Dð1; q1 ; q2 ÞT ein s of A corresponding to in s . Where

(

 aI ÞþbI q1 ¼ ðin þl1 BÞð1þ ; bI

q2 ¼ inBþl : 3

We normalize q and q by the condition hq ðsÞ; qðhÞi ¼ 1. Clearly, hq ðsÞ; qðhÞi ¼ 0. In order to assure hq ðsÞ; qðhÞi ¼ 1, we need to determine the value of D. By (4.7), we have

! hq ðsÞ; qðhÞi ¼ Dð1; q1 ; q 2 Þð1; q1 ; q2 ÞT  ¼ Df1 þ

q1 q1

þ

q2 q2



Z

0

1

Z

0

1

Z

h

r¼0

    Dð1; q1 ; q2 Þein s ðrhÞ dgðhÞð1; q1 ; q2 ÞT ein s r dr

  ð1; q1 ; q2 Þhein s h d

¼ Df1 þ q1 q1 þ q2 q2 þ s q1 cel3 ðs

 þmÞ

gðhÞð1; q1 ; q2 ÞT g

  ð1  q2 Þein s g:

Therefore, we can choose D as

D¼

1 :  1 þ q1 q1 þ q2 q2 þ s q1 cel3 ðs þmÞ ð1  q2 Þein s

Following the algorithms given in [12] and using similar computation process in [13], we can get the coefficients which will be used to determine the important quantities

g 11 g 02 g 21

2bs Dðq1  1Þ

2abs S D

ð1  q1 Þq21 ; 3 ð1 þ aI Þ 2abs S D ¼ Refq1 g þ ð1  q1 Þjq1 j2 ; 3  2 ð1 þ aI Þ ð1 þ aI Þ 2bs Dðq1  1Þ 2abs S D ¼ q1 þ ð1  q1 Þq21 ; 3  2 ð1 þ aI Þ ð1 þ aI Þ i bs Dðq1  1Þ h ð1Þ ð1Þ ð2Þ ð2Þ ¼ W 20 ð0Þq1 þ 2q1 W 11 ð0Þ þ W 20 ð0Þ þ 2W 11 ð0Þ  2 ð1 þ aI Þ 2abs S D 2abs D ð2Þ ð2Þ þ ð1  q1 Þ½q1 W 20 ð0Þ þ 2q1 W 11 ð0Þ þ ð1  q1 Þ½q21 þ 2jq1 j2 ; 3  3 ð1 þ aI Þ ð1 þ aI Þ

g 20 ¼

 2

ð1 þ aI Þ 2bs Dðq1  1Þ

q1 þ

ð4:8Þ

where

W 20 ðhÞ ¼



ig 20 ig 02 g 20 g 02   ig   in s h in s h e2in s h ,  20 qð0Þein s h qð0Þe þ qð0Þe þ W ð0Þ þ qð0Þ þ qð0Þ 20 n s 3n s in s 3in s ns þ

ig 02     qð0Þein s h þ E1 e2in s h 3n s

ð4:9Þ

and

W 11 ðhÞ ¼

ig 11   qð0Þein s h þ E2 : n s

ð4:10Þ

Moreover E1 ; E2 satisfy the following equations, respectively,

0 B B B @



2in þ l1  B þ 1þbIaI 

   Bq þ ð1þbSaI Þ2  ceðl3 þin Þs 

 1þbIaI

2in þ ðl2 þ cÞ  Bp  ð1þbSaI Þ2

0

ceðl3 þin Þs  c



B 0 2in þ l3

1 C 2bs q1 2abs S q21 C T ð1; 1; 0Þ þ ð1; 1; 0ÞT CE1 ¼ 2 3 A ð1 þ aI Þ ð1 þ aI Þ

171

Z. Yinying, J. Jianwen / Applied Mathematics and Computation 230 (2014) 164–173

and

0 B B B @







l1  B þ 1þbIaI Bq þ ð1þbSaI Þ2  ceðl3 þin Þs B

C 2bs Refq1 g 2abs S jq1 j2 C T ð1; 1; 0Þ þ ð1; 1; 0ÞT : 0 CE2 ¼ 2 3 A ð1 þ aI Þ ð1 þ aI Þ





 1þbIaI

l2 þ c  Bp  ð1þbSaI Þ2

0

ceðl3 þin Þs  c

1



l3

Because each g ij is expressed by the parameters and display in (4.8), we can compute the following quantities:

c1 ð0Þ ¼

i jg j g ðg g  2jg 11 j2  02 Þ þ 21 ; 2n s 20 11 3 2

m2 ¼ 

Refc1 ð0Þg Refdkðdss Þg 

; ð4:11Þ



Refc1 ð0Þg þ m2 Refdkðdss Þg T2 ¼  : n  s

b2 ¼ 2Refc1 ð0Þg;

By the result of Hassard et al. [12], we have the following: Theorem 4.1. In (4.11), the following results hold: (i) the sign of m2 determines the directions of the Hopf bifurcation: if m2 > 0ðm2 < 0Þ, then the Hopf bifurcation is supercritical ðsubcriticalÞ and the bifurcating periodic solutions exist for s > s ðs < s Þ; (ii) the sign of b2 determines the stability of the bifurcating periodic solutions: the bifurcating periodic solutions are stable ðunstableÞ if b2 < 0ðb2 > 0Þ; (iii) the sign of T 2 determines the period of the bifurcating periodic solutions: the period increases ðdecreasesÞ if T 2 > 0ðT 2 < 0Þ. 5. Numerical example In this section, we present some numerical results of system (2.1). In Section 3, we have obtained that when R0 < 1, the disease-free equilibrium was stable, while if R0 > 1, the disease-free equilibrium was unstable. We have researched the existence of Hopf bifurcation and have obtained the stability and direction of the Hopf bifurcation according to the normal theory and the center manifold theorem. Our theoretical results show that the time delay s must be responsible for the observed regular cycles of disease incidence. We consider the following two set of the parameters: A ¼ 0:2; l1 ¼ 0:1; l2 ¼ 0:2; l3 ¼ 0:1; b ¼ 4; a ¼ 1; c ¼ 10; B ¼ 0:0001; q ¼ 0:2; p ¼ 0:8; s ¼ 5 and Sð0Þ ¼ 3; Ið0Þ ¼ 1; Rð0Þ ¼ 1. By directly computing, we obtain R0 ¼ 0:7851 < 1. According to Theorem 3, we know the ‘‘infectionfree’’ periodic solution of system (2.1) is locally asymptotically stable for this case (see Fig. 1). A ¼ 0:2; l1 ¼ 0:1; l2 ¼ 0:2; l3 ¼ 0:1; b ¼ 4; a ¼ 1; c ¼ 5; B ¼ 0:0001; q ¼ 0:2; p ¼ 0:8; s ¼ 5 and Sð0Þ ¼ 3; Ið0Þ ¼ 1; Rð0Þ ¼ 1. By directly computing, we obtain R0 ¼ 1:54 > 1. For this case, the diseases will be permanent (see Figs. 2 and 3).

(1)

(2)

(a) Time series of S(t) of the system.

(c) Time series of R(t) of the system.

(b) Time series of I(t) of the system. 1.4

3

2.8 1.2

2.5

2.7 1

2.6

2

2.4

R(t)

0.8 t

S(t)

2.5

2.3

1

2.2

0.4

2.1

0.5

0.2

2 1.9

1.5

0.6

500

1000 t

1500

2000

0

500

1000 I(t)

1500

2000

0

500

1000

1500

2000

t

Fig. 1. Time series on the threshold values R0 ¼ 0:7851 < 1. The ‘‘infection-free’’ periodic solution of system (2.1) is locally asymptotically stable.

172

Z. Yinying, J. Jianwen / Applied Mathematics and Computation 230 (2014) 164–173

(a) Time series of S(t) of the system.

(b) Time series of I(t) of the system.

3 2.8

2.2

0.25

2

0.15

1.8

2.2

0.1

1.6

2

0.05

R(t)

0.2

2.4

t

2.6

0

1.8

1.4 1.2

−0.05

1.6

1

−0.1 1.4

0.6

−0.2 0

1

2

3

4

5

t

6

x 10

2

2.5

4

3

3.5

4

I(t)

0.4 4

0

1

2

3

4

5

(a) Phase diagrams of S(t), I(t), R(t).

x 10

(b) Phase diagrams of S(t), I(t). 0.1 0.09 0.08

0.7

0.07

0.6

0.06

I(t)

0.8

0.5

0.05 0.04

0.4

0.03

0.3 0.1

0.02 1.6 0.05

0.01

1.4 1.2 0

I(t)

1

0 S(t)

1

(c) Phase diagrams of S(t), R(t).

0.75

0.75

0.7

0.7

0.65

0.65

0.6

0.6

R(t)

0.8

0.55

0.5

0.45

0.45

0.4

0.4

0.35

0.35

1.1

1.2

1.3 S(t)

1.4

1.3 S(t)

1.4

1.5

1.6

0.55

0.5

1

1.2

(d) Phase diagrams of I(t), R(t).

0.8

0.3

1.1

1.5

1.6

0.3

0

0.02

Fig. 3. Phase diagram on the threshold values R0 ¼ 1:54 > 1.

0.04

0.06 I(t)

6 4

t

x 10

Fig. 2. Time series on the threshold values R0 ¼ 1:54 > 1. The diseases will be permanent.

R(t)

1

0.8

−0.15

1.2

R(t)

S(t)

(c) Time series of R(t) of the system. 2.4

0.3

0.08

0.1

Z. Yinying, J. Jianwen / Applied Mathematics and Computation 230 (2014) 164–173

173

References [1] W.O. Kermack, A.G. Mckendrick, A contribution to the mathematical theory of epidemics, Proc. R. Soc. London A 115 (1927) 700–721. [2] S. Busenberg, K.L. Cooke, M.A. Pozio, Analysis of a model of a vertically transmitted disease, J. Math. Biol. 17 (1983) 305–329. [3] S.N. Busenberg, K.L. Cooke, Models of vertical transmitted diseases with sequential-continuous dynamics, in: V. Lakshmicantham (Ed.), Nonlinear Phenomena in Mathematical Sciences, Academic Press, New York, 1982, pp. 179–187. [4] K.L. Cook, S.N. Busenberg, Vertically transmitted diseases, in: V. Lakshmicantham (Ed.), Nonlinear Phenomena in Mathematical Sciences, Academic Press, New York, 1982, pp. 189–197. [5] S.N. Busenberg, K.L. Cooke, Vertically transmitted diseases models and dynamics, Biomathematics, vol. 23, Springer-Verlag, New York, 1993. [6] T.L. Zhang, J.L. Liu, Z.D. Teng, Stability of Hopf bifurcation of a delayed SIRS epidemic model with stage structure, Nonlinear Anal.: Real World Appl. 11 (2010) 293–306. [7] K. Cooke, P. van den Driessche, X. Zou, Interaction of maturation delay and nonlinear birth in population and epidemic models, J. Math. Biol. 39 (1999) 332–352. [8] O. Diekmann, J.A.P. Heesterbeek, J.A.J. Metz, On the definition and the computation of the basic reproduction ratio R0 in models for infectious diseases in heterogeneous populations, J. Math. Biol. 28 (1990) 365–382. [9] P. van den Driessche, J. Watmough, Reproduction numbers and sub-threshold endemic equilibria for compartmental models of disease transmission, Math. Biosci. 180 (2002) 29–48. [10] E. Beretta, Y. Kuang, Geometric stability switch criteria in delay differential systems with delay dependent parameters, SIAM J. Math. Anal. 33 (2002) 1144–1165. [11] J.K. Hale, S.M. Verduyn Lunel, Introduction to Functional Differential Equations, Springer, 1993. [12] B.D. Hassard, N.D. Kazarinoff, Y.H. Wan, Theory and Applications of Hopf Bifurcation, CUP Archive, 1981. [13] J.J. Wei, S.G. Ruan, Stability and bifurcation in a neural network model with two delays, Phys. D: Nonlinear Phenom. 130 (1999) 255–272.