# Bifurcation Analysis of a Delayed Worm Propagation Model with Saturated Incidence

Authors: Zizhen Zhang; Yougang Wang; Luca Guerrini
Journal: Advances in Mathematical Physics
 (2018)

---

## Abstract

This paper is concerned with a delayed SVEIR worm propagation model with saturated incidence. The main objective is to investigate the effect of the time delay on the model. Sufficient conditions for local stability of the positive equilibrium and existence of a Hopf bifurcation are obtained by choosing the time delay as the bifurcation parameter. Particularly, explicit formulas determining direction of the Hopf bifurcation and stability of the bifurcating periodic solutions are derived by using the normal form theory and the center manifold theorem. Numerical simulations for a set of parameter values are carried out to illustrate the analytical results.


---

## Body

## 1. Introduction

Worms, as one kind of malicious codes, have become one of the main threats to the security of networks. Since the first Morris worm in 1998, new worms have come into networks frequently, including Slammer worm [1], Commwarrior worm [2], Cabir worm [3], and Chameleon worm [4]. Each of them can cause enormous financial losses and social panic [5–7]. Therefore, it is significant to explore effective methods to counter against worms. To this end, we need to accurately understand the dynamic behaviors of worm propagation in networks. Considering that the process of worm propagation in networks is similar to that of biological virus propagation in the population, mathematical models have been important tools used to analyze the propagation and control of worms based on the theory of Kermack and McKendrick [8].In [9], Kim et al. proposed the SIS (Susceptible-Infectious-Susceptible) model in order to analyze the dynamical behaviors of worm propagation on Internet. However, the SIS model neglects the effect of the antivirus software. Thus, the SIR (Susceptible-Infectious-Recovered) model is proposed [9]. Although SIR model considered the immunity of the nodes in which the worms have been cleaned, however, it assumes that the recovered hosts have permanent immunity. This is not consistent with the reality in networks, because they may be infected by some new emerging worms again. To overcome this drawback of the SIR model, Wang et al. investigated the SIRS (Susceptible-Infectious-Recovered-Susceptible) mode for analyzing the dynamics of worm propagation in networks [10–12]. It should be pointed out that both the SIR mode and the SIRS model assume that the susceptible nodes become infectious instantaneously. As we know, worms usually have a latent period. Based on this consideration, the SEIR (Susceptible-Exposed-Infectious-Recovered) model [13, 14] and the SEIRS (Susceptible-Exposed-Infectious-Recovered-Susceptible) model [11, 15] are proposed to describe the dynamics of worm propagation in networks. Considering influence of the quarantine strategy and the vaccination strategy on the propagation of worms, some worm models with quarantine strategy [16–19] and vaccination strategy [20–25] are formulated and analyzed.It should be pointed out that all the models above use the bilinear incidence rateβSI. As stated in [26], the dynamics of a model system heavily depends on the choice of the incidence rate. Gan et al. have considered the different incidence rate functions βSI/f(I) in their work [27, 28]. It was found that the saturated incidence rate βSI/(1+ηI) is more general than the bilinear incidence rate βSI. Based on this, Wang et al. [29] proposed the following model with partial immunization to defend against worms:(1)dStdt=1-pA-βStIt1+ηIt-μSt+γVt,dVtdt=pA-σβVtIt-μ+γVt,dEtdt=βStIt1+ηIt+σβVtIt-μ+ϖEt,dItdt=ϖEt-μ+α+δIt,dRtdt=δIt-μRt,where S(t), V(t), E(t), I(t), and R(t) present numbers of the susceptible, vaccinated, exposed, infectious, recovered hosts at time t, respectively. The meanings of more parameters are described and shown in “Parameters of the Model and Their Meanings” section. Wang et al. [29] investigated the stability of system (1).One of the significant features of computer viruses is their latent characteristics [30, 31]. In addition, time delays of one type or another could cause the numbers of hosts in system (1) to fluctuate. And worm propagation models with time delay have been investigated by some scholars [14, 17, 19]. Based on above discussions, in this paper, we extend system (1) by incorporating the time delay due to the latent period of the worms in the exposed hosts into system (1) and obtain the following delayed worm propagation model:(2)dStdt=1-pA-βStIt1+ηIt-μSt+γVt,dVtdt=pA-σβVtIt-μ+γVt,dEtdt=βStIt1+ηIt+σβVtIt-μEt-ϖEt-τ,dItdt=ϖEt-τ-μ+α+δIt,dRtdt=δIt-μRt,where τ is the latent period of the worms in the exposed nodes.The remainder of this paper is organized as follows. Local stability of the positive equilibrium and existence of a Hopf bifurcation at the positive equilibrium are analyzed in the next section. Properties of the Hopf bifurcation such as direction and stability are investigated in Section3. Numerical simulations are carried out in Section 4 to support the obtained theoretical results. Finally, conclusions are given in Section 5 to end our work.


## 2. Existence of Hopf Bifurcation

By direct computation, we know that if the condition(H1): (μ+ϖ)(μ+α+δ)(σβI∗+μ+γ)>ϖσβpA holds, then system (2) has a positive equilibrium P∗(S∗,V∗,E∗,I∗,R∗), where(3)S∗=1+ηI∗μ+ϖμ+α+δσβI∗+μ+γ-ϖσβpAβϖσβI∗+μ+γ,V∗=pAσβI∗+μ+γ,E∗=μ+α+δϖI∗,R∗=δμI∗.And I∗ is the positive root of the following equation:(4)P2x2+P1x+P0=0,where(5)P0=μμ+ϖμ+γμ+α+δ-Aϖβγ+pμδ+1-pμ,P1=μ+ϖμ+α+δσβμ+μ+γμη+β-Aϖσββ+pμη,P2=βσμ+ϖμ+α+δβ+μη.The Jacobi matrix of system (2) about P∗(S∗,V∗,E∗,I∗,R∗) is given by(6)JP∗=m11m120m1400m220m240m31m32m33+n33e-λτm34000n43e-λτm440000m54m55,where(7)m11=-βI∗1+ηI∗+μ,m12=γ,m14=-βS∗1+ηI∗2,m22=-μ+γ+σβI∗,m24=-σβV∗,m31=βI∗1+ηI∗,m32=σβI∗,m33=-μ,m34=βS∗1+ηI∗2+σβV∗,m44=-μ+α+δ,m54=δ,m55=-μ,n33=-ϖ,n43=ϖ.The characteristic equation of that matrix (6) is(8)λ5+M4λ4+M3λ3+M2λ2+M1λ+M0+N4λ4+N3λ3+N2λ2+N1λ+N0e-λτ=0,with(9)M0=-m11m22m33m44m55,M1=m11m22m33m44+m55m11m22m33+m44+m33m44m11+m22,M2=-m11m22m33+m44+m33m44m11+m22-m55m11m22+m33m44+m11+m22m33+m44,M3=m11m22+m33m44+m11+m22m33+m44+m55m11+m22+m33+m44,M4=-m11+m22+m33+m44+m55,N0=m11m22m55n43m34-m44+m55n43m14m22m31-m11m24m32-m12m24m31,N1=n43m11m22m44+m55+m44m55m11+m22+m12m24m31+m24m32n43m11+m55-m14m31n43m22+m55-m34n43m11m22+m11m55+m22m55,N2=n43m14m31-m24m32+m34n43m11+m22+m55-n43m11m22+m44m55+m11+m22m44+m55,N3=n43m11+m22+m44+m55-m34,N4=-n43.When τ=0, (8) becomes(10)λ5+M04λ4+M03λ3+M02λ2+M01λ+M00=0,where(11)M00=M0+N0,M01=M1+N1,M02=M2+N2,M03=M3+N3,M04=M4+N4.Thus,P∗(S∗,E∗,I∗,R∗,V∗) is locally asymptotically stable when τ=0 if the condition (H2) is satisfied and (H2) is defined as follows:(12)M00>0,M04>0,M03M04>M02,M02M01+M03M04>M01M042+M022,M02M03M00+M01+M02+2M00M01M04>M002+M01M022+M04M012M04+M00M032.For τ>0, let λ=iω(ω>0) be the root of (8). Then, we have(13)N1ω-N3ω3sin⁡τω+N4ω4-N2ω2+N0cos⁡τω=M2ω2-M4ω4-M0,N1ω-N3ω3cos⁡τω-N4ω4-N2ω2+N0sin⁡τω=M3ω3-ω5-M1ω.Thus, we can get the following equation:(14)ω10+h4ω8+h3ω6+h2ω4+h1ω2+h0=0,where(15)h0=M02-N02,h1=M12-2M0M2-N12+2N0N2,h2=M22+2M0M4-2M1M3+2N1N3-N22-2N0N4,h3=M32-2M2M4+2M1-N32+2N2N4,h4=M42-2M3-N42.Let v=ω2; then (14) becomes(16)v5+h4v4+h3v3+h2ω2+h1v+h0=0.Based on the discussion about the distribution of the roots of (16) in [32], we suppose that (H3): (16) has at least one positive root v0.If the condition(H3) holds, then (16) has a positive root ω0=v0 and (8) has a pair of purely imaginary roots ±iω0. For ω0, we have(17)τ0=1ω0×arccos⁡G1ωG2ω,with(18)G1ω=N3-M4N4ω8+M2N4+M4N2-M3N3-N1ω6+M3N1+M1N3-M0N4-M2N2-M4N0ω4+M0N2+M2N0-M1N1ω2+M0N0,G2ω=N42ω8+N32-2N2N4ω6+N22+2N0N4-2N-1N3ω4+N12-2N0N2ω2+N02.Differentiating on both sides of (8) with respect to τ, we can obtain(19)dλdτ-1=5λ4+4M4λ3+3M3λ2+2M2λ+M1eλτλN4λ4+N3λ3+N2λ2+N1λ+N0-τλ.Further, we have(20)Re⁡dλdττ=τ0-1=f′v0G2ω0,where f(v)=v5+h4v4+h3v3+h2ω2+h1v+h0.Obviously, if the condition(H4):f′(v0)≠0 is satisfied, then Re[dλ/dτ]τ=τ0-1≠0. Based on the discussion above and the Hopf bifurcation theorem in [33], we have the following results.Theorem 1.
For system (2), if the conditions (H1)–(H4) hold, then the positive equilibrium P∗(S∗,V∗,E∗,I∗,R∗) is locally asymptotically stable when τ∈[0,τ0); system (2) undergoes a Hopf bifurcation at the P∗(S∗,V∗,E∗,I∗,R∗) when τ=τ0 and a family of periodic solutions bifurcate from P∗(S∗,V∗,E∗,I∗,R∗).


## 3. Properties of the Hopf Bifurcation

Letu1(t)=S(t)-S∗, u2(t)=V(t)-V∗, u3(t)=E(t)-E∗, u4(t)=I(t)-I∗, and u5(t)=R(t)-R∗, and normalize the time delay with the scaling t→(t/τ). Let τ=τ0+ϱ(ϱ∈R); then ϱ=0 is the Hopf bifurcation value of system (2). System (2) can be transformed into the following form:(21)u˙t=Lϱut+Fϱ,ut,where u(t)=(u1,u2,u3,u4,u5)T∈C=C([-1,0],R5) and Lϱ: C→R5 and F: R×C→R5 are given, respectively, by(22)Lϱϕ=τ0+ϱM1ϕ0+M2ϕ-1,Fϱ,ϕ=m15ϕ10ϕ40+m16ϕ420+m17ϕ10ϕ420+m18ϕ430+⋯m25ϕ20ϕ40m35ϕ10ϕ40+m36ϕ420+m37ϕ10ϕ420+m38ϕ430+m39ϕ20ϕ40⋯00,with(23)M1=m11m120m1400m220m240m31m32m33m340000m440a5100m54m55,M2=000000000000n330000n430000000,m15=-β1+ηI∗2,m16=βηS∗1+ηI∗3,m17=βη1+ηI∗3,m18=-βη2S∗1+ηI∗4,m25=-σβ,m35=β1+ηI∗2,m36=-βηS∗1+ηI∗3,m37=-βη1+ηI∗3,m38=βη2S∗1+ηI∗4,m39=σβ.According to the Riesz representation theorem, there exists a5×5 matrix function η(θ,ϱ): θ∈[-1,0]→R5 such that(24)Lϱϕ=∫-10dηθ,ϱϕθ,forϕ∈C.In fact, choosing(25)ηθ,ϱ=τ0+ϱAmax⁡δθ+Bmax⁡δθ+1and δ(θ) is the Dirac delta function.Forϕ∈C([-1,0],R5), define(26)Aϱϕ=dϕθdθ,-1≤θ<0,∫-10dηθ,ϱϕθ,θ=0,Rϱϕ=0,-1≤θ<0,Fϱ,ϕ,θ=0.Then system (21) can be transformed into the following operator equation:(27)u˙t=Aϱut+Rϱut.For φ∈C1([0,1],(R5)∗), we further define the adjoint operator(28)A∗φ=-dφsds,0<s≤1,∫-10dηTs,0φ-s,s=0and the bilinear inner product as follows:(29)φs,ϕθ=φ¯0ϕ0-∫θ=-10∫ξ=0θφ¯ξ-θdηθϕξdξ,where η(θ)=η(θ,0).Based on the discussion above, we know that±iω0τ0 are eigenvalues of A(0). Thus, they are also eigenvalues of A∗. Let q(θ)=(1,q2,q3,q4,q5)Teiτ0ω0θ be the eigenvector of A(0) corresponding to +iω0τ0 and q∗(s)=(1/D)(1,q2∗,q3∗,q4∗,q5∗)Teiτ0ω0s be the eigenvectors of A∗ corresponding to -iτ0ω0. By direct computation, we can obtain(30)q2=m24iω0-m11m12m24+m14iω0-m22,q3=iω0-m44n43e-iτ0ω0q4,q4=iω0-m22m24q2,q5=m54iω0-m55q4,q2∗=-m12+m32q3iω0+m22,q3∗=-iω0+m11m31,q4∗=-iω0+m33+n33eiτ0ω0n43eiτ0ω0,q5∗=-m24q2∗+m34q3∗+iω0+m11q4∗m54,D¯=1+q2q¯2∗+q3q¯3∗+q4q¯4∗+q5q¯5∗+τ0e-iτ0ω0q3n33q¯3∗+n43q¯4∗.Then we have 〈q∗,q〉=1 and 〈q∗,q¯〉=0.Next, we can obtain the coefficients which can determine the properties of the Hopf bifurcation atτ0 by following the algorithms given in [33] and using the computation process similar to those in [34–36]:(31)g20=2τ0D¯m15q4+m16q42+m25q2¯∗q2q4+q¯3∗m35q4+m36q42+m39q2q4,g11=τ0D¯m15q4+q¯4+2m16q4q¯4+m25q¯2∗q2q¯4+q¯2q4+q¯3∗m35q4+q4¯+2m36q4q¯4+m39q2q¯4+q¯2q4,g02=2τ0D¯m15q¯4+m16q¯42+m25q2¯∗q¯2q¯4+q¯3∗m35q¯4+m36q¯42+m39q¯2q¯4,g21=2τ0D¯m1512W2040+W1140+12W2010+W1110q¯4+m16W2040q¯4+2W1140q4+m17q42+2q4q¯4+3m18q42q¯4+m25q¯2∗12W2020q¯4+W1120q4+12W2040q¯2+W1140q2+q¯3∗m3512W2040+W1140+12W2010+W1110q¯4+m36W2040q¯4+2W1140q4+m37q42+2q4q¯4+3m38q42q¯4+m3912W2020q¯4+W1120q4+12W2040q¯2+W1140q2,with(32)W20θ=ig20ρ0τ0ω0eiτ0ω0θ+ig¯02ρ¯03τ0ω0e-iτ0ω0θ+E1e2iτ0ω0θ,W11θ=-ig11ρ0τ0ω0eiτ0ω0θ+ig¯11ρ¯0τ0ω0e-iτ0ω0θ+E2,E1=2m11′-m120-m1400m22′0-m240-m31-m32m33′-m34000-n43e-2iτ0ω0m44′0000-m54m55′-1×E11E12E1300,E2=-m11m120m1400m220m240m31m32m33+n33m34000n43m440000m54m55-1×E21E22E2300,where(33)m11′=2iω0-m11,m22′=2iω0-m22,m33′=2iω0-m33-n33e-2iτ0ω0,m44′=2iω0-m44,m55′=2iω0-m55,E11=m15q4+m16q42,E12=m25q2q4,E13=m35q4+m36q42+m39q2q4,E21=m15q4+q¯4+2m16q4q¯4,E22=m25q2q¯4+q¯2q4,E23=m35q4+q4¯+2m36q4q¯4+m39q2q¯4+q¯2q4.Then, one can obtain(34)C10=i2τ0ω0g11g20-2g112-g0223+g212μ2=-Re⁡C10Re⁡λ′τ0,β2=2Re⁡C10,T2=-Im⁡C10+μ2Im⁡λ′τ0τ0ω0.In conclusion, we have the following results.Theorem 2.
For system (2), μ2 determines the direction of the Hopf bifurcation: if μ2>0(μ2<0), then the Hopf bifurcation is supercritical (subcritical); β2 determines the stability of the bifurcating periodic solution: the bifurcating periodic solutions are stable (unstable) if β2<0(β2>0); T2 determines the period of the bifurcating periodic solutions: the period increases (decreases) if T2>0(T2<0).


## 4. Numerical Simulation

In this section, some numerical simulations are carried out for qualitative analysis by using Matlab software package. By extracting some values from [29] and considering the conditions for the existence of the Hopf bifurcation, we choose a set of parameters as follows: A=100, p=0.5, α=0.0001, β=0.003, γ=0.0001, δ=0.4, η=1, σ=0.05ϖ=0.02, and μ=0.001. Then, we can get the following specific case of system (2):(35)dStdt=50-0.003StIt1+It-0.001St+0.0001Vt,dVtdt=50-1.5000e-004VtIt-0.0011Vt,dEtdt=0.003StIt1+It+1.5000e-004VtIt-0.001Et-0.02Et-τ,dItdt=0.02Et-τ-0.4011It,dRtdt=0.4It-0.0001Rt.By some computations, we can obtain the following equation with respect to I:(36)5.0539e-009I2-1.0117e-006I-4.7907e-006=0.It follows that system (35) has a unique positive equilibrium P∗(12723,1571.3,4107.5,204.8103,81924) and we can verify that P∗(12723,1571.3,4107.5,204.8103,81924) is locally asymptotically stable when τ=0. Further, we have ω0=0.5508 and τ0=69.6986. According to Theorem 1, it can be concluded that P∗(12723,1571.3,4107.5,204.8103,81924) is locally asymptotically stable when τ∈[0,τ0=69.6986). This property can be shown as in Figures 1 and 2. However, a Hopf bifurcation will occur and a family of periodic solutions bifurcate from P∗(12723,1571.3,4107.5,204.8103,81924) when the value of τ passes through the Hopf bifurcation value τ0, which can be illustrated by Figures 3 and 4.Figure 1
Dynamic behavior of system (35): projection on S-E-R with τ=65.85<τ0.Figure 2
Dynamic behavior of system (35): projection on V-E-I with τ=65.85<τ0.Figure 3
Dynamic behavior of system (35): projection on S-E-R with τ=76.65>τ0.Figure 4
Dynamic behavior of system (35): projection on V-E-I with τ=76.65>τ0.In addition, we obtainC1(0)=-4.3990+2.9057i, λ′(τ0)=0.7014-0.0212i by some complicate computations. Thus, we get μ2=7.8776>0, β2=-8.798<0, and T2=-0.0713<0 based on (34). It follows from Theorem 2 that the Hopf bifurcation is supercritical and the bifurcating periodic solutions are stable and decrease. Since the bifurcating periodic solutions are stable, then the five classes of hosts in system (35) may coexist in an oscillatory mode from the view of the biological point, which is not welcome in networks.


## 5. Conclusions

In this study, the dynamical behaviors of a delayed SVEIR worm propagation model with saturated incidence are discussed based on the work in literature [29]. The dynamical behaviors of the model are investigated from the point of view of local stability and Hopf bifurcation both analytically and numerically. The threshold of the time delay τ0 at which the model causes a Hopf bifurcation is obtained by using eigenvalue method. We found that characteristics of the propagation of worms in the model can be predicted and controlled when the value of delay is suitably small (τ∈[0,τ0)). However, propagation of the worms in the model will be out of control once the value of the time delay is above the threshold value τ0. Accordingly, we can know that the propagation of worms in the model can be controlled by postponing occurrence of the Hopf bifurcation. Moreover, the properties of the Hopf bifurcation are investigated by applying the normal form theory and the center manifold theorem. Numerical simulations are also presented in order to testify our obtained theoretical results.



---

*Source: 7619074-2018-05-15.xml*