# Stability and Hopf Bifurcation in a Delayed SEIRS Worm Model in Computer Network

Authors: Zizhen Zhang; Huizhong Yang
Journal: Mathematical Problems in Engineering
 (2013)

---

## Abstract

A delayed SEIRS epidemic model with vertical transmission in computer network is considered. Sufficient conditions for local stability of the positive equilibrium and existence of local Hopf bifurcation are obtained by analyzing distribution of the roots of the associated characteristic equation. Furthermore, the direction of the local Hopf bifurcation and the stability of the bifurcating periodic solutions are determined by using the normal form theory and center manifold theorem. Finally, a numerical example is presented to verify the theoretical analysis.


---

## Body

## 1. Introduction

Since the pioneering work of Kephart and White [1, 2], classical epidemiological models have been used to understand and predict virus propagations in computer network by many authors [3–10]. In [3], Thommes and Coates proposed a modified version of the SEI model to predict the virus propagation in a network. In [5], Yuan and Chen proposed an e-SEIR model and studied the behavior of virus propagation with the presence of antivirus countermeasures. Yuan and Chen [5] supposed that the recovered computers have a permanent immunization period and can no longer be infected. Considering that there is no permanent recovery from the virus till a node is attached to the computer network, Mishra and Pandey [6] proposed the following epidemic model for the transmission of worms with vertical transmission:
(1)dS(t)dt=b-βI(t)S(t)-pbE(t)-qbI(t)-dS(t)+ζR(t),dE(t)dt=βI(t)S(t)+pbE(t)+qbI(t)-εE(t)-dE(t),dI(t)dt=εE(t)-γI(t)-dI(t)-ηI(t),dR(t)dt=γI(t)-dR(t)-ζR(t),
where S(t), E(t), I(t), and R(t) denote the numbers of nodes at time t in susceptible, exposed, infectious and recovered states, respectively. b is the recruitment rate of susceptible nodes to the computer network. d is the crashing rate of nodes due to the reason other than the attack of worms. p and q are the fractions of the new nodes from the exposed and the infectious classes, respectively, that are introduced into the exposed class. β is the transmission rate. ε, γ, and ζ are the state transition rates. η is the crashing rate of nodes due to the attack of worms. Mishra and Pandey [6] discussed the characteristic of the basic reproduction number and investigated the global stability of system (1) by constructing a new Liyapunov function.It is well known that time delays can play a complicated role on the dynamics of a system. They can cause the stable equilibrium of a system to become unstable and make a system bifurcate periodic solutions. Dynamical systems with delay have been studied by many scholars [11–20]. Motivated by the works above and considering that the antivirus software may use a period to clean the worms in one node and that the recovered nodes may have a temporary immunity period due to the antivirus software, we consider the following delayed system in this paper:
(2)dS(t)dt=b-βI(t)S(t)-pbE(t)-qbI(t)-dS(t)+ζR(t-τ2),dE(t)dt=βI(t)S(t)+pbE(t)+qbI(t)-εE(t)-dE(t),dI(t)dt=εE(t)-γI(t-τ1)-dI(t)-ηI(t),dR(t)dt=γI(t-τ1)-dR(t)-ζR(t-τ2),
where τ1 is the period to clean the worms in one node and τ2 is the temporary immunity period. For the convenience of analysis, throughout this paper, we assume that the period to clean the worms in one node and the temporary immunity period are the same.Letτ1=τ2=τ, and then system (2) becomes the following form:
(3)dS(t)dt=b-βI(t)S(t)-pbE(t)-qbI(t)-dS(t)+ζR(t-τ),dE(t)dt=βI(t)S(t)+pbE(t)+qbI(t)-εE(t)-dE(t),dI(t)dt=εE(t)-γI(t-τ)-dI(t)-ηI(t),dR(t)dt=γI(t-τ)-dR(t)-ζR(t-τ).The main purpose of this paper is to investigate the effects of the time delay on the dynamics of system (3). We study the stability of the positive equilibrium of system (3) and find the critical value of the time delay where a Hopf bifurcation occurs. We also study the properties of the Hopf bifurcation such as direction and stability.This paper is organized as follows. In Section2, we investigate local stability of the positive equilibrium and obtain the sufficient conditions for the existence of local Hopf bifurcation. In Section 3, we determine the direction and the stability of the Hopf bifurcation by using the normal form theory and center manifold theorem. In order to testify the theoretical analysis, a numerical example is presented in Section 4.


## 2. Stability and Existence of Local Hopf Bifurcation

It is not difficult to verify that if(4)R0=bβε+bdqε+bdp(d+γ+η)d(d+ε)(d+γ+η)>1,(d+ε)(d+γ+η)>pb(d+γ+η)+qbε,
system (2) has a unique positive equilibrium D*(S*,E*,I*,R*), where
(5)S*=(d+ε)(d+γ+η)-pb(d+γ+η)-qbεβε,I*=(d+ζ)[bβε+d(pb(d+γ+η)222222222+qbε-(d+ε)(d+γ+η))]×(β(d+ε)(d+ζ)(d+γ+η)-βγεζ)-1,E*=d+γ+ηεI*,R*=γd+ζI*.R
0 is called the basic reproduction number.The Jacobian matrix of system (3) about the positive equilibrium D* is(6)J(D*)=(λ-a11-a12-a12-b14e-λτ-a21λ-a22-a2300-a32λ-a33-b33e-λτ000-b43e-λτλ-a44-b43e-λτ),where
(7)a11=-(d+βI*),a12=-pb,a13=-(qb+βS*),a21=βI*,a22=pb-d-ε,a23=qb+βS*,a32=ε,a33=-(d+η),a44=-d,b14=ζ,b33=-γ,b43=γ,b44=-ζ.
Thus, the characteristic equation of system (2) at the positive equilibrium D* is
(8)λ4+A3λ3+A2λ2+A1λ+A0+(B3λ3+B2λ2+B1λ+B0)e-λτ+(C2λ2+C1λ+C0)e-2λτ=0,
where
(9)A0=a11a44(a22a33-a23a32)+a21a44(a13a32-a12a33),A1=(a33+a44)(a12a21-a11a22)+a23a32(a11+a44)-a33a44(a11+a22)+a13a21a32,A2=(a11+a22)(a33+a44)+a11a22+a33a44-a12a21-a23a32,A3=-(a11+a22+a33+a44),B0=(a11a22-a12a21)(a33b44+a44b33)+(a13a21-a11a23)a32b44,B1=(b33+b44)(a12a21-a11a22)+a23a32b44-(a11+a22)(a33b44+a44b33),B2=(a11+a22)(b33+b44)+a33b44+a44b33,B3=-(b33+b44),C2=b33b44,C1=b33b44(a11+a22),C0=(a11a22-a12a21)b33b44-a21a32b14b43.Multiplyingeλτ on both sides of (8), it is easy to get
(10)B3λ3+B2λ2+B1λ+B0+(λ4+A3λ3+A2λ2+A1λ+A0)eλτ+(C2λ2+C1λ+C0)e-λτ=0.Whenτ=0, (10) reduces to
(11)λ4+A13λ3+A12λ2+A11λ+A10=0,
where
(12)A10=A0+B0+C0,A11=A1+B1+C1,A12=A2+B2+C2,A13=A3+B3.By the Routh-Hurwitz criterion, the sufficient conditions for all roots of (11) to have a negative real part given in the following form:
(13)D1=A13>0,D2=det(A131A11A12)>0,D3=det(A1310A11A12A130A10A11)>0,(14)D4=det(A13100A11A12A1310A10A11A12000A10)>0.Thus, if the condition(H1) holds, which means that (13) and (14) are satisfied, the positive equilibrium D* is locally asymptotically stable in the absence of delay.Forτ>0, let λ=iω(ω>0) be a root of (10). Then, we can get
(15)(ω4-(A2+C2)ω2+A0+C0)cosτω-((A1-C1)ω-A3ω3)sinτω=B2ω2-B0,(ω4-(A2-C2)ω2+A0-C0)sinτω+((A1+C1)ω-A3ω3)cosτω=B3ω3-B1ω.
Then, we can get
(16)sinτω=p7ω7+p5ω5+p3ω3+p1ωω8+q6ω6+q4ω4+q2ω2+q0,cosτω=p6ω6+p4ω4+p2ω2+p0ω8+q6ω6+q4ω4+q2ω2+q0,
where
(17)p0=B0(C0-A0),p1=(A1+C1)B0-(A0+C0)B1,p2=(A2-C2)B0+(A0-C0)B2+(C1-A1)B1,p3=(A0+C0)B3-(A1+C1)B2+(A2+C2)B1-A3B0,p4=(A1-C1)B3-(A2-C2)B2+A3B1-B0,p5=A3B2-B1-(A2+C2)B3,p6=B2-A3B3,p7=B3,q0=A02-C02,q2=A12-C12-2A0A2+2C0C2,q4=A22-C22+2A0-2A1A3,q6=A32-2A2.It is well known thatsin2τω+cos2τω=1. Thus, we have
(18)ω16+e7ω14+e6ω12+e5ω10+e4ω8+e3ω6+e2ω4+e1ω2+e0=0,
where
(19)e0=q02-p02,e1=2q0q2-2p0p2-p12,e2=q22+2q0q4-p22-2p0p4-2p1p3,e3=2q0q6+2q2q4-p32-2p0p6-2p1p5-2p2p4,e4=q42+2q0+2q2q6-p42-2p1p7-2p2p6-2p3p5,e5=2q2+2q4q6-2p3p7-2p4p6-p52,e6=q62-p62+2q4-2p5p7,e7=2q6-p72.Letω2=v, and then (18) becomes
(20)v8+e7v7+e6v6+e5v5+e4v4+e3v3+e2v2+e1v+e0=0.In order to give the main results in this paper, we made the following assumption.(
H
2
) Equation (20) has at least one positive real root.Suppose that the condition(H2) holds. Without loss of generality, we suppose that (20) has eight positive real roots, which are denoted as v1,v2,…,v8, respectively. Then, (18) has eight positive roots ωk=vk,k=1,2,…,8. For every fixed ωk, the corresponding critical value of time delay is
(21)τk(j)=1ωkarccosp6ω6+p4ωk4+p2ωk2+p0ω8+q6ωk6+q4ωk4+q2ωk2+q0+2jπωk,22222222222222222k=1,2,…,8,j=0,1,2,….Let(22)τ0=min{τk(0)},k∈{1,2,…,8},ω0=ωk|τ=τ0.Taking the derivative of with respect to (10), it is easy to obtain
(23)[dλdτ]-1=(3B3λ2+2B2λ+B1+(2C2λ+C1)e-λτ+(4λ3+3A3λ2+2A2λ+A1)eλτ)×((C2λ3+C1λ2+C0λ)e-λτ-(λ5+A3λ4+A2λ3+A1λ2+A0λ)eλτ))-1-τλ.
Then, we have
(24)Re[dλdτ]τ=τ0-1=PRQR+PIQIQR2+QI2,PR=(A1+C1-3A3ω02)cosτ0ω0+2(2ω03-(A2-C2)ω0)sinτ0ω0-3B3ω02+B1,PI=(A1-C1-3A3ω02)sinτ0ω0+2(ω03-(A2+C2)ω0)cosτ0ω0+2B2ω0,QR=(A1ω02-A3ω04)cosτ0ω0+(ω05-A2ω03+A0ω0)sinτ0ω0,QI=(A1ω02-A3ω04)sinτ0ω0-(ω05-A2ω03+A0ω0)cosτ0ω0.Obviously, if the condition(H3)PRQR+PIQI≠0 holds, then Re[dλ/dτ]τ=τ0-1≠0. Thus, according to the Hopf bifurcation theorem in [21], we have the following results.Theorem 1.
Suppose that the conditions (H1)–(H3) hold. The positive equilibrium D*(S*,E*,I*,R*) of system (3) is asymptotically stable for τ∈[0,τ0). System (3) undergoes a Hopf bifurcation at the positive equilibrium D*(S*,E*,I*,R*) when τ=τ0 and a family of periodic solutions bifurcate from the positive equilibrium D*(S*,E*,I*,R*) near τ=τ0.


## 3. Properties of the Hopf Bifurcation

In the previous section, we have obtained the conditions for the Hopf bifurcation to occur whenτ=τ0. In this section, we will study properties of the Hopf bifurcation such as direction and stability by using the normal form theory and the center manifold theorem in [21].Letu1(t)=S(t)-S*, u2(t)=E(t)-E*, u3(t)=I(t)-I*, u4(t)=R(t)-R*, and τ=τ0+μ,μ∈R and normalize t→(t/τ). Then system (3) can be transformed into the following form:
(25)u˙(t)=Lμut+F(μ,ut),
where ut=(u1(t),u2(t),u3(t),u4(t))T∈C=C([-1,0],R4) and Lμ:C→R4 and F:R×C→R4 are given, respectively, by
(26)Lμϕ=(τ0+μ)(a11a12a130a21a22a2300a32a330000a44)ϕ(0)-(τ0+μ)(000b14000000b33000b43b44)ϕ(-1),F(μ,ϕ)=(τ0+μ)(-βϕ1(0)ϕ3(0)βϕ1(0)ϕ3(0)00).By the Riesz representation theorem, there exists a4×4 matrix function η(θ,μ):[-1,0]→R4 whose elements are of bounded variation such that
(27)Lμϕ=∫-10‍dη(θ,μ)ϕ(θ),ϕ∈C.In fact, we choose(28)η(θ,μ)=(τ0+μ)(a11a12a130a21a22a2300a32a330000a44)δ(θ)-(τ0+μ)(000b14000000b33000b43b44)δ(θ+1),
where δ is the Dirac delta function.Forϕ∈C([-1,0],R4), we define
(29)A(μ)ϕ={dϕ(θ)dθ,-1≤θ<0,∫-10‍dη(θ,μ)ϕ(θ),θ=0,R(μ)ϕ={0,-1≤θ<0,F(μ,ϕ),θ=0.
Then system (25) can be transformed into the following operator equation:
(30)u˙(t)=A(μ)ut+R(μ)ut.The adjoint operatorA* of A is defined by
(31)A*(φ)={-dφ(s)ds,0<s≤1∫-10dηT(s,0)φ(-s),s=0,
associated with a bilinear form:
(32)〈φ(s),ϕ(θ)〉=φ-(0)ϕ(0)-∫θ=-10‍∫ξ=0θ‍φ-(ξ-θ)dη(θ)ϕ(ξ)dξ,
where η(θ)=η(θ,0).Letρ(θ)=(1,ρ2,ρ3,ρ4)Teiω0τ0θ be the eigenvector of A corresponding to iω0τ0 and let ρ*(s)=D(1,ρ2*,ρ3*,ρ4*)eiω0τ0s be the eigenvector of A* corresponding to -iω0τ0. From the definition of A(0) and A*(0), we can get
(33)ρ2=-iω0-a33-b33e-iτ0ω0(iω0-a22)(iω0-a33-b33e-iτ0ω0)-a23a32,ρ3=a21a32(iω0-a22)(iω0-a33-b33e-iτ0ω0),ρ4=b43e-τ0ω0iω0-a44+b33e-τ0ω0ρ3,ρ2*=-a11iω0+a21,r3*=-(iω0(a12-a11)+a12a21-a11a22a32(iω0+a21),ρ4*=b14eiτ0ω0iω0+a44+b44eiτ0ω0.From (32), we have
(34)〈ρ*,ρ〉=D-[1+ρ2ρ-2*+ρ3ρ-3*+ρ4ρ-4*+τ0e-iτ0ω0×(ρ3(b33ρ-3*+b43ρ-4*)+ρ4(b14+b44ρ-4*))τ0e-iτ0ω0].Then, we choose(35)D-=[1+ρ2ρ-2*+ρ3ρ-3*+ρ4ρ-4*+τ0e-iτ0ω0×(ρ3(b33ρ-3*+b43ρ-4*)+ρ4(b14+b44ρ-4*))τ0e-iτ0ω0]-1
such that 〈q*,q〉=1, 〈q*,q-〉=0.Next, we can obtain the coefficients which will be used to determine the properties of the Hopf bifurcation by following the algorithms introduced in [21] and using a computation process as in [22]:
(36)g20=2βτ0D-ρ3(ρ-2*-1),g02=2βτ0D-ρ-3(ρ-2*-1),g11=βτ0D-(r3+r-3)(r-2*-1)),g21=2βτ0D-(ρ-2*-1)(W11(1)(0)ρ3+12W20(1)(0)ρ-3+W11(3)(0)+12W20(3)(0)),
with
(37)W20(θ)=ig20q(0)τ0ω0eiτ0ω0θ+ig-02q-(0)3τ0ω0e-iτ0ω0θ+E1e2iτ0ω0θ,W11(θ)=-ig11q(0)τ0ω0eiτ0ω0θ+ig-11q-(0)τ0ω0e-iτ0ω0θ+E2,
where E1 and E2 can be determined by the following equations, respectively:(38)E1=2(2iω0-a11-a12-a13-b14e-2iω0τ0-a212iω0-a22-a2300-a322iω0-a33-b33e-2iω0τ0000-b43e-2iω0τ02iω0-a44-b44e-2iω0τ0)-1(E1(1)E1(2)00),E2=-(a11a12a13b14a21a22a2300a32a33+b33000b43a44+b44)-1(E2(1)E2(2)00),with
(39)E1(1)=-βρ3,E1(2)=βρ3,E2(1)=-β(ρ3+ρ-3),E2(2)=β(ρ3+ρ-3).Therefore, we can calculate the following values:(40)C1(0)=i2τ0ω0(g11g20-2|g11|2-|g02|23)+g212,μ2=-Re{C1(0)}Re{λ′(τ0)},β2=2Re{C1(0)},T2=-Im{C1(0)}+μ2Im{λ′(τ0)}τ0ω0.
Based on the discussion above, we can obtain the following results.Theorem 2.
For system (3), if μ2>0(μ2<0), the Hopf bifurcation is supercritical (subcritical). If β2<0(β2>0), the bifurcating periodic solutions are stable (unstable). If T2>0(T2<0) the period of the bifurcating periodic solutions increases (decreases).


## 4. Numerical Simulation

In this section, we present a numerical example to verify the theoretical analysis in Sections2 and 3. Let b=0.7, d=0.06, p=0.02, q=0.01, β=0.1, ε=0.2, γ=0.05, η=0.01, and ζ=0.1. We can get a particular case of system (3):
(41)dS(t)dt=0.7-0.1I(t)S(t)-0.014E(t)-0.007I(t)-0.06S(t)+0.1R(t-τ),dE(t)dt=0.1I(t)S(t)+0.014E(t)+0.007I(t)-0.2E(t)-0.06E(t),dI(t)dt=0.2E(t)-0.05I(t-τ)-0.06I(t)-0.01I(t),dR(t)dt=0.05I(t-τ)-0.06R(t)-0.1R(t-τ).Then, we can get thatR0=7.5774>1, (ε+d)(d+γ+η)=0.0312>pb(d+γ+η)+qbε=0.0031, and system (41) has a unique positive equilibrium D*(1.4050,1.5422,4.9350, and 2.9610). Further, we have D1=1.0795>0, D2=0.3360>0, D3=0.0141>0, and D4=2.8200e-004>0. That is, condition (H1) holds. Thus, we can obtain ω0=0.1245, τ0=13.6109, and λ′(τ0)=0.0073+0.0299i. Therefore, from Theorem 1, the positive equilibrium D* is asymptotically stable when τ∈[0,13.6109) and unstable when τ>13.6109. As can be seen from Figures 1, 2, and 3, when τ=12.5∈[0,13.6109), the positive equilibrium D* is asymptotically stable. However, when τ=14.95>13.6109, the positive equilibrium D* will lose its stability, a Hopf bifurcation occurs, and a family of periodic solutions bifurcate from the positive equilibrium D*. This property can be illustrated by Figures 4, 5, and 6. In addition, from (40), we can get C1(0)=-0.00904+0.6833i, μ2=12.3836>0, β2=-0.1808<0, and T2=-0.6217<0. Thus, according to Theorem 2, we know that the Hopf bifurcation is supercritical. The bifurcating periodic solutions are stable and the period of the periodic solutions decreases.Figure 1
The track of the statesS, E, I, and R for τ=12.5000<13.6109=τ0.Figure 2
The phase plot of the statesS, E, and I for τ=12.5000<13.6109=τ0.Figure 3
The phase plot of the statesE, I, and R for τ=12.5000<13.6109=τ0.Figure 4
The track of the statesS, E, I, and R for τ=14.9500>13.6109=τ0.Figure 5
The phase plot of the statesS, E, and I for τ=14.9500>13.6109=τ0.Figure 6
The phase plot of the statesE, I, and R for τ=14.9500>13.6109=τ0.


## 5. Conclusions

In this paper, we incorporate the time delay due to the period the antivirus software has to use to clean the worms in one node and the temporary immunity period of the recovered nodes into the model considered in the literature [6] and get a delayed SEIRS epidemic model for the transmission of worms in computer network through vertical transmission. The effects of the time delay on the dynamics of the model are investigated. It is found that the time delay can play a complicated role on the model by analyzing the distribution of the roots of the associated characteristic equation. When the time delay is suitable small, the positive equilibrium is asymptotically stable. However, a local Hopf bifurcation occurs and a branch of periodic solutions bifurcates from the positive equilibrium when the delay passes through the critical value τ0. Furthermore, the properties of the Hopf bifurcation such as direction and stability are determined by using the normal form theory and center manifold theorem. In order to verify the theoretical analysis, a numerical example is also included.



---

*Source: 319174-2013-11-20.xml*