# Dynamic Modeling and Analysis of  the Email Virus Propagation

Authors: Yihong Li; Jinxiao Pan; Zhen Jin
Journal: Discrete Dynamics in Nature and Society
 (2012)

---

## Abstract

A novel deterministic SEIS model for the transmission of email viruses in growing communication networks is formulated. Interestingly, the model is different from classical SEIS models not only in the form, but also in the dynamical features. We study the equilibria and their stability and analyse the bifurcation dynamics of the model. In particular, we find that the virus-free equilibrium is locally asymptotically stable for any parameter values, which may attribute to the absence of the basic reproduction number. It is shown that the model undergoes a saddle-node bifurcation and admits the bistable phenomenon. Moreover, on the basis of the Lyapunov function, the domains of attraction of equilibria are estimated by solving an LMI optimization problem. Based on the above theoretical results, some effective strategies are also provided to control the propagation of the email viruses. Additionally, our results are confirmed by numerical simulations.


---

## Body

## 1. Introduction

Currently, email has become one of the most basic applications in the Internet with the development of networked computer. As more and more people rely on email communications for business and everyday life, email viruses constitute one of the Internet major security threats. When an email virus infects a machine, it sends an infected email to the addresses in the computer's email address book. This self-broadcast mechanism allows for the virus's rapid reproduction and spread. Due to this facility, and the advantages of sending viruses through email [1], hackers mostly tend to choose the email as the vector of spreading their viruses [2]. Indeed, according to the Virus Bulletin [3] and National Computer Virus Emergency Response Center [4], the email viruses still account for a large share of the virus prevalence today.In order to eradicate viruses, as well as to control and limit the impact of an outbreak, a lot of efforts have been devoted to develop various mathematical models to simulate the real case of viruses’ spreading [5–9]. Unlike scanning worms such as Code Red or Slammer, email viruses spread over a logical network defined by email address relationship, making traditional epidemic models invalid for modelling the propagation of email virus.The spread of email virus is a complex process. An adequate modeling of this process requires both a correct description of the logical network along which viruses spread [10–12] and a quantitative formulation of various behavioral mechanisms of the users who motivate the spread of the virus [1, 2].Due to the specific characteristics of email viruses propagation, some researchers focus attention on the topological structures of email networks [13], while some researchers deal with the development of new techniques for the detection and elimination of viruses [14]. Immunization strategies in email networks are presented by experiments and simulations in [15] and [16]. In 2003, Zou et al. [1] presented an email viruses’ model that accounts for the behaviors of email users, such as email checking frequency and the probability of opening an email attachment. But they rely primarily on simulation rather than mathematical analysis. In [12], the simulations for a stochastic susceptible-hidden-infectious-recovered model and mathematical analysis for homogeneous and heterogeneous SIR model are provided. When considering immunization in the heterogeneous SIR model, the saddle-node bifurcation occurs. Here we must emphasize the hidden state which is a key characteristic that the email virus different from other viruses. So based on the spreading mechanism of email viruses, we proposed a SEIS model in growing homogeneous networks by stochastic approach. Relying primarily on mathematical analysis, we give some effective control strategies.The rest of the paper is organized as follows. In Section2, a discrete-time model of the email viruses’ spreading is proposed whose detailed process is also given. Then we deal with the model to derive the continuous-time equations for the dynamics of email viruses spreading. In Section 3, we present a qualitative analysis of the model. We show that the virus-free equilibrium is locally asymptotically stable for any parameter values and the system admits a saddle-node bifurcation. Section 4 estimates the domain of attraction of the equilibrium on the basis of the determined Lyapunov function. The subsequent sections are the numerical examples and the conclusion and discussion.


## 2. Mechanism Analysis and Model Establishment

In this paper, we consider email viruses that are only transferred through users' email address books. Thus email address relationship forms a logical network for email viruses. Email address book of a user usually contains the user's friends' or business partners' email addresses. Thus if user A has user B's address, user B probably also has user A's address in his own address book. Due to this, we model the email network as an undirected graph, although this may not always be true. Without loss of generality, we can safely assume that every distinct email address is associated with one host computer which is associated with a single user. Then the logical network can be seen as host computer network, in which the vertices are the computers and the links are the address relationships of the corresponding users. Here we consider the network as a homogeneous network that assumes each computer has same neighbors which are supposed ask; that is, an email user has k addresses in his own address book and his address is in k users' address books.According to the already studied propagation characteristics, the spreading mechanism of email viruses can be described as follows: a healthy computer, which links with infectious one, would receive virus emails. But it may not be infectious immediately. The email viruses are hiding in the email box. At this moment, the computer is infected but not yet infectious. Only if the user opens the virus attachments, the computer becomes infectious. Then if the computer has been found infectious, the user would take antivirus measures to clean the viruses. After that, the computer recovers from infectious state and returns to healthy state again.Obviously, the viruses’ spreading is just the process of computers switching between different states essentially. We call computers that have not received virus emails and could be infected easily after receiving virus emails as the susceptible computers. The computer that has received virus emails but the virus attachments have not been opened is called the exposed computer. And further, a computer that has been infected by the viruses and can transmit them to susceptible computers is called infectious computer. The numbers of susceptible, exposed, and infectious computers at timet are denoted as S(t), E(t), and I(t), respectively. Then the total number of computers at time t is N(t)=S(t)+E(t)+I(t). Here we assume the total number of computers is varying dynamically as time t; that is, the new computers can be connected to the network and some old computers can be removed from the network (including the obsolescence computers). Suppose A is the rate at which new computers are connected and d is the rate at which the computers are removed per unit time. The flow diagram of the state transition is depicted in Figure 1.Figure 1
The flowchart of email virus spreading.By the above analysis, we have already known the mechanism of the virus propagation. In the following, we will elaborate the probabilities of transforming fromS to E, E to I, and I to S, by assuming the number of infected contacts as a random variable that follows the binomial distribution.


### 2.1. The Transition from Susceptible State to Exposed State

Consider now a computerU is in the susceptible state at time t. By the usual homogeneous mixing approximation, Θ(t)=I(t)/N(t) denotes the rate that a neighbor of the computer U is in the infectious state at time t. Suppose that k neighbors of computer U are independent with each other, we have k independent events occurring with the same probability Θ(t). Then the number of the infected neighbors of computer U is a random variable that follows the binomial distribution B[k,Θ(t)]. And the probability of the number of the infected neighbors of computer U equal to ki at time t can be written as follows:
(2.1)P(X=ki)=(kki)Θ(t)ki(1-Θ(t))k-ki.
When a computer is infected, it will send a virus email to all email addresses in the user's email address book. So as long as there are infected computers in the k neighbors, the susceptible computers would have chance to turn into exposed. Let pSS denote the probability that computer U still stays in the susceptible state in the time interval [t,t+Δt], then 1-pSS is the probability that it enters into the exposed state. ki is the number of the infected neighbors of computer U, and if p denotes the probability that the computer U received virus emails from each infected neighbor per unit time, then we have
(2.2)pSS(ki,t)=(1-Δtp)ki.Because the infected numberki of the computer's k neighbors can change from 0 to k and the corresponding pSS(ki,t) changes from 1 to (1-Δtp)k, using (2.1) and (2.2), we obtain the transition probability p¯SS(t) averaged over all possible values of ki as
(2.3)p¯SS(t)=E(pSS(ki,t))=∑ki=0k‍(1-Δtp)ki(kki)Θ(t)ki(1-Θ(t))k-ki=(1-ΔtpΘ(t))k.
Then in the interval [t,t+Δt], the probability of a susceptible computer transforming to the exposed state is 1-(1-ΔtpΘ(t))k.


### 2.2. The Transition from Exposed State to Infectious State

An exposed computer being infected or not depends on whether the user opens the virus emails. As long as the email user opens a virus attachment (one virus email only has one virus attachment) arbitrarily, the computer will be infected, so only if the user does not open all virus emails that he has received, the computer will stay in the exposed state. Letpo be the probability that the email user opens a virus attachment (0<po<1); thus 1-po is the probability that the email user does not open a virus attachment. Let the number of the virus emails that the computer U received during the time interval [t,t+Δt] be ne; then the probability that all ne virus attachments would not be opened is (1-po)ne. Therefore, the probability of the computer infected is 1-(1-po)ne.In addition, from (2.1), we can obtain that the average infected number of the computer's neighbors is
(2.4)〈ki〉=∑ki=0k‍kiP(X=ki)=∑ki=0k‍ki(kki)Θ(t)ki(1-Θ(t))k-ki=kΘ(t).
Suppose the computer U received a virus emails from each infected neighbor per unit time, then in the time interval [t,t+Δt], the computer U received akΘ(t)Δt virus emails on average; that is, ne=akΘ(t)Δt.Substitutingne into the above formula, the probability of an exposed computer becomes infected is 1-(1-po)kaΘ(t)Δt in the time interval Δt.


### 2.3. The Transition from Infectious State to Susceptible State

An infected computer recovering from the infection and reentering into the susceptible state is dependent on the computer users. If the computers have been found virus-infected, the users would take antivirus measures; as a result the viruses could be removed and the computers will enter into susceptible state. But if the users do not do this, the computers will stay in the infected state. So, we suppose the computers return to the susceptible state at the probabilityδ per unit time.Based on the above analyses and Figure1, we can indeed write the spread mechanism of email viruses as the following deterministic model:
(2.5)S(t+Δt)=S(t)+ΔtA-[1-(1-ΔtpI(t)N(t))k]S(t)+ΔtδI(t)-ΔtdS(t),E(t+Δt)=E(t)+[1-(1-ΔtpI(t)N(t))k]S(t)-[1-(1-po)ka(I(t)/N(t))Δt]E(t)-ΔtdE(t),I(t+Δt)=I(t)+[1-(1-po)ka(I(t)/N(t))Δt]E(t)-ΔtδI(t)-ΔtdI(t).
In the limit Δt→0, the above model can be simplified by keeping only the two first terms in the development of the binomial. This leads to the following continuous-time model:
(2.6)dS(t)dt=A-αS(t)I(t)N(t)+δI(t)-dS(t),dE(t)dt=αS(t)I(t)N(t)-βE(t)I(t)N(t)-dE(t),dI(t)dt=βE(t)I(t)N(t)-dI(t)-δI(t),
where α=kp, β=kaln(1/(1-po)).Summing up the three equations in (2.6), we obtain
(2.7)dN(t)dt=A-dN(t),
when t→+∞, N(t)→A/d≜N0. It is easy to see that system (2.6) can be shown to be mathematically well posed in the positive invariant region 𝔻={(S,E,I)∣0≤S+E+I≤N0} and solution in 𝔻 exists for all positive time.Hence, (2.6) has the following limit system:
(2.8)dE(t)dt=αN0[N0-E(t)-I(t)]I(t)-βN0E(t)I(t)-dE(t),dI(t)dt=βN0E(t)I(t)-(d+δ)I(t).
Rescaling (2.8) by x1=(α/N0(d+δ))E, x2=(α/N0(d+δ))I, τ=(d+δ)t and still writing τ as t, we obtain
(2.9)dx1dt=(B-x1-x2)x2-mx1x2-nx1,dx2dt=mx1x2-x2,
where B=α/(d+δ),m=β/α,n=d/(d+δ), and the corresponding positive invariant set is 𝕏={(x1,x2)∣0≤x1+x2≤B}.


## 2.1. The Transition from Susceptible State to Exposed State

Consider now a computerU is in the susceptible state at time t. By the usual homogeneous mixing approximation, Θ(t)=I(t)/N(t) denotes the rate that a neighbor of the computer U is in the infectious state at time t. Suppose that k neighbors of computer U are independent with each other, we have k independent events occurring with the same probability Θ(t). Then the number of the infected neighbors of computer U is a random variable that follows the binomial distribution B[k,Θ(t)]. And the probability of the number of the infected neighbors of computer U equal to ki at time t can be written as follows:
(2.1)P(X=ki)=(kki)Θ(t)ki(1-Θ(t))k-ki.
When a computer is infected, it will send a virus email to all email addresses in the user's email address book. So as long as there are infected computers in the k neighbors, the susceptible computers would have chance to turn into exposed. Let pSS denote the probability that computer U still stays in the susceptible state in the time interval [t,t+Δt], then 1-pSS is the probability that it enters into the exposed state. ki is the number of the infected neighbors of computer U, and if p denotes the probability that the computer U received virus emails from each infected neighbor per unit time, then we have
(2.2)pSS(ki,t)=(1-Δtp)ki.Because the infected numberki of the computer's k neighbors can change from 0 to k and the corresponding pSS(ki,t) changes from 1 to (1-Δtp)k, using (2.1) and (2.2), we obtain the transition probability p¯SS(t) averaged over all possible values of ki as
(2.3)p¯SS(t)=E(pSS(ki,t))=∑ki=0k‍(1-Δtp)ki(kki)Θ(t)ki(1-Θ(t))k-ki=(1-ΔtpΘ(t))k.
Then in the interval [t,t+Δt], the probability of a susceptible computer transforming to the exposed state is 1-(1-ΔtpΘ(t))k.


## 2.2. The Transition from Exposed State to Infectious State

An exposed computer being infected or not depends on whether the user opens the virus emails. As long as the email user opens a virus attachment (one virus email only has one virus attachment) arbitrarily, the computer will be infected, so only if the user does not open all virus emails that he has received, the computer will stay in the exposed state. Letpo be the probability that the email user opens a virus attachment (0<po<1); thus 1-po is the probability that the email user does not open a virus attachment. Let the number of the virus emails that the computer U received during the time interval [t,t+Δt] be ne; then the probability that all ne virus attachments would not be opened is (1-po)ne. Therefore, the probability of the computer infected is 1-(1-po)ne.In addition, from (2.1), we can obtain that the average infected number of the computer's neighbors is
(2.4)〈ki〉=∑ki=0k‍kiP(X=ki)=∑ki=0k‍ki(kki)Θ(t)ki(1-Θ(t))k-ki=kΘ(t).
Suppose the computer U received a virus emails from each infected neighbor per unit time, then in the time interval [t,t+Δt], the computer U received akΘ(t)Δt virus emails on average; that is, ne=akΘ(t)Δt.Substitutingne into the above formula, the probability of an exposed computer becomes infected is 1-(1-po)kaΘ(t)Δt in the time interval Δt.


## 2.3. The Transition from Infectious State to Susceptible State

An infected computer recovering from the infection and reentering into the susceptible state is dependent on the computer users. If the computers have been found virus-infected, the users would take antivirus measures; as a result the viruses could be removed and the computers will enter into susceptible state. But if the users do not do this, the computers will stay in the infected state. So, we suppose the computers return to the susceptible state at the probabilityδ per unit time.Based on the above analyses and Figure1, we can indeed write the spread mechanism of email viruses as the following deterministic model:
(2.5)S(t+Δt)=S(t)+ΔtA-[1-(1-ΔtpI(t)N(t))k]S(t)+ΔtδI(t)-ΔtdS(t),E(t+Δt)=E(t)+[1-(1-ΔtpI(t)N(t))k]S(t)-[1-(1-po)ka(I(t)/N(t))Δt]E(t)-ΔtdE(t),I(t+Δt)=I(t)+[1-(1-po)ka(I(t)/N(t))Δt]E(t)-ΔtδI(t)-ΔtdI(t).
In the limit Δt→0, the above model can be simplified by keeping only the two first terms in the development of the binomial. This leads to the following continuous-time model:
(2.6)dS(t)dt=A-αS(t)I(t)N(t)+δI(t)-dS(t),dE(t)dt=αS(t)I(t)N(t)-βE(t)I(t)N(t)-dE(t),dI(t)dt=βE(t)I(t)N(t)-dI(t)-δI(t),
where α=kp, β=kaln(1/(1-po)).Summing up the three equations in (2.6), we obtain
(2.7)dN(t)dt=A-dN(t),
when t→+∞, N(t)→A/d≜N0. It is easy to see that system (2.6) can be shown to be mathematically well posed in the positive invariant region 𝔻={(S,E,I)∣0≤S+E+I≤N0} and solution in 𝔻 exists for all positive time.Hence, (2.6) has the following limit system:
(2.8)dE(t)dt=αN0[N0-E(t)-I(t)]I(t)-βN0E(t)I(t)-dE(t),dI(t)dt=βN0E(t)I(t)-(d+δ)I(t).
Rescaling (2.8) by x1=(α/N0(d+δ))E, x2=(α/N0(d+δ))I, τ=(d+δ)t and still writing τ as t, we obtain
(2.9)dx1dt=(B-x1-x2)x2-mx1x2-nx1,dx2dt=mx1x2-x2,
where B=α/(d+δ),m=β/α,n=d/(d+δ), and the corresponding positive invariant set is 𝕏={(x1,x2)∣0≤x1+x2≤B}.


## 3. Dynamic Analysis of System (2.9)

The objective of this section is to perform a qualitative analysis of system (2.9). We first consider the nonexistence of limit cycle. Then we study the equilibria and their stability. Finally, we prove the occurrence of saddle-node bifurcation.Theorem 3.1.
For system (2.9), there is no limit cycle.Proof.
Take the Dulac function(3.1)D(x1,x2)=1x2,
where x1>0, x2>0.
Set(3.2)(B-x1-x2)x2-mx1x2-nx1=P(x1,x2),mx1x2-x2=Q(x1,x2);
we have
(3.3)∂(DP)∂x1+∂(DQ)∂x2=-1-m-nx2<0.
Therefore, there is not a closed orbit inside the first quadrant.Next, we will study the equilibrium situation of system (2.9). The system always has the virus-free equilibrium E0=(0,0), and the Jacobian matrix is JE0=[-nB0-1]. Obviously, (0,0) is locally asymptotically stable for any parameter values.To find the positive equilibrium, set(3.4)(B-x1-x2)x2-mx1x2-nx1=0,mx1-1=0,
which yields the equation of x2:
(3.5)mx22+(1+m-Bm)x2+n=0.
We can obtain that(1)
there is no positive equilibrium ifB<2n/m+1/m+1;(2)
there is one positive equilibrium ifB=2n/m+1/m+1;(3)
there are two positive equilibria ifB>2n/m+1/m+1.In the following, we analyse the stability of the positive equilibria in different situations, respectively, and prove the existence of saddle-node bifurcation, to elaborate from three cases.Case 1 (only one positive equilibrium in system (2.9)).
SupposeB=2n/m+1/m+1; system (2.9) has one positive equilibrium
(3.6)E*=(x1*,x2*)=(1m,12(B-1m-1));
the Jacobian matrix at E* is
(3.7)JE*=[-(1+m)x2*-n0mx2*0].
Clearly, the eigenvalues of matrix JE* are σ1=0 and σ2=-(1+m)x2*-n; E* is a nonhyperbolic point. Therefore, we investigate the dynamics near E* by the center manifold theorem [17, 18].
Firstly, shiftE* to the origin via y1=x1-x1* and y2=x2-x2*; system (2.9) can be transformed into
(3.8)dy1dt=[-n-(1+m)x2*]y1-(1+m)y1y2-y22,dy2dt=mx2*y1+my1y2.
Secondly, define the transformation
(3.9)[y1y2]=H[z1z2],H=[011mx2*σ2],
which transformed system (3.8) into the following standard form:
(3.10)dz1dt=f1(z1,z2),dz2dt=σ2z2+f2(z1,z2),
where f1=(1+m)(mx2*/σ2)(z1z2+(mx2*/σ2)z22)+(mx2*/σ2)(z1+(mx2*/σ2)z2)2+m(z1z2+(mx2*/σ2)z22) and f2=-(1+m)(z1z2+(mx2*/σ2)z22)-(z1+(mx2*/σ2)z2)2.
By the existence theorem [19], there exists a center manifold for system (3.10), which can be expressed locally as follows:
(3.11)Wc(0)={(z1,z2)∈R2∣z2=h(z1),‖z1‖<ρ,h(0)=0,Dh(0)=0∣ρ>0},
with ρ sufficiently small and Dh is the derivative of h with respect to z1.
To compute the center manifoldWc(0), we suppose h(z1) has the form
(3.12)h(z1)=h1z12+h2z13+h3z14+h4z15+⋯.
By the local center manifold theorem, the center manifold (3.12) satisfies
(3.13)Dh⋅f1(z1,z2)-σ2z2-f2(z1,z2)=0,
where Dh(z1)=2h1z1+3h2z12+4h3z13+5h4z14+⋯.
Rewritef1 and f2 as f1=a1z12+a2z1z2+a3z22 and f2=b1z12+b2z1z2+b3z22, respectively, where
(3.14)a1=mx2*σ2,a2=(1+m)mx2*σ2+2(mx2*σ2)2+m,a3=(1+m)(mx2*σ2)2+(mx2*σ2)3+m2x2*σ2,b1=-1,b2=-(1+m)-2mx2*σ2,b3=-(1+m)mx2*σ2-(mx2*σ2)2.
Substituting (3.10) into (3.13) and then equating coefficients on each power of z1 to zero yields
(3.15)h1=1σ2,h2=2a1h1-b2h1σ2,h3=2a1h12+3a1h2-b2h2-b3h12σ2,h4=5a2h1h2+2a3h13+4a1h3-b2h3-2b3h1h2σ2.
Then, we get the approximation for h:
(3.16)h=1σ2[z12+(2a1h1-b2h1)z13+(2a1h12+3a1h2-b2h2-b3h12)z14+(5a2h1h2+2a3h13+4a1h3-b2h3-2b3h1h2)z15]+⋯.
Substituting (3.16) to the first equation of system (3.10), we achieve
(3.17)dz1dt=f1(z1,z2)=mx2*σ2z12+[(1+m)mx2*σ2+2(mx2*σ2)2+m]1σ2z13+⋯.
So from (3.17) (see [20, p. 338–340]), we get the following theorem.Theorem 3.2.
IfB=2n/m+1/m+1, the system (2.9) has one positive equilibrium E* and E* is a saddle node.Case 2 (two positive equilibria in system (2.9)).
SupposeB>2n/m+1/m+1, there are two positive equilibria
(3.18)E1*=(x1*,x21*)=(1m,B-1/m-1-(B-1/m-1)2-4n/m2),E2*=(x1*,x22*)=(1m,B-1/m-1+(B-1/m-1)2-4n/m2).
We first determine the stability of E1*. The Jacobian matrix at E1* is
(3.19)JE1*=[-(1+m)x21*-n(B-1m-1)2-4nmmx21*0].
Obviously, det(JE1*)=-mx21*(B-1/m-1)2-4n/m<0, so E1* is a saddle.
Now we analyze the stability of the second positive equilibriumE2*. The Jacobian matrix at E2* is
(3.20)JE2*=[-(1+m)x22*-nB-1m-1-2x22*mx22*0].
By a similar argument as above, we obtain that
(3.21)det⁡(JE2*)=-mx22*(B-1m-1-2x22*)=mx22*(B-1m-1)2-4nm>0,trace(JE2*)=-(1+m)x22*-n<0,(3.22)Δ=trace2(JE2*)-4det⁡(JE2*)=(1+m)2(x22*)2+2n(1+m)x22*+n2+4(Bm-m-1)x22*-8m(x22*)2.
After some algebra calculations, we obtain Δ>0; E*2 is a stable node. So we have the following theorem.Theorem 3.3.
IfB>2n/m+1/m+1, the system (2.9) has two positive equilibria E1* and E2*; E1* is a saddle and E2* is a stable node.Case 3 (saddle-node bifurcation).
From Cases1 and 2, we can see system (2.9) experiences a saddle-node bifurcation. In this part, we give the proof.
Let us considerB as a bifurcation parameter, define B0=2n/m+1/m+1 and x0=E*. Rewrite system (2.9) as
(3.23)x˙=f(x,B)=[-nx1+Bx2-(1+m)x1x2-x22-x2+mx1x2].
Then we have
(3.24)Df(x0,B0)=[-(1+m)x2*-n0mx2*0],fB(x0,B0)=[x2*0].Df(x0,B0) has a simple eigenvalue λ=0 with eigenvector v=[0,1]T, and that DTf(x0,B0) has an eigenvector w=[mx2*/((1+m)x2*+n),1]T corresponding to λ=0. Furthermore, the following conditions are satisfied:
(3.25)wTfB(x0,B0)=mx2*2(1+m)x2*+n≠0,wTD2f(x0,B0)(v,v)=-2mx2*(1+m)x2*+n≠0.
According to the theorem (Sotomayor; see [21, p.148]), we get the following result.Theorem 3.4.
System (2.9) experiences a saddle-node bifurcation at the equilibrium x0=E* as the parameter B passes through the bifurcation value B=B0.


## 4. Estimation of the Domain of Attraction

From the previous study, we know, whenB>2n/m+1/m+1, E2* is a stable node, and at the same time, the virus-free equilibrium is also locally asymptotically stable, that is, the bistable state. In this case, the eventual behavior of the system is sensitive to the initial positions. So in this section, we study the domains of attraction about the two equilibria using the method in [22]. For this task, we make some preparations.Definition 4.1 (domain of attraction [23]).
The domain of attraction of the origin is given by(4.1)S={x0∈Rn∣lim⁡t→∞x(t,x0)⟶0}.
Consider the following system:(4.2)x˙=f(x),
where f:ℝn→ℝn is an analytical function with the following properties:(1)
f(0)=0; that is,  x=0 is an equilibrium point of (4.2);(2)
every eigenvalue of∂f/∂x(0), the Jacobian matrix of f at the origin, has negative real part; that is, x=0 is an asymptotically stable equilibrium point.
LetDA be the domain of attraction of the zero solution of (4.2). The following result provides a tool to determine DA.Lemma 4.2 (see [24]).
D
A coincides with the natural domain of analyticity of the unique V of the problem
(4.3)〈∇V,f〉=-‖x‖2,V(0)=0.
The function V is positive on DA and limx→x0V(x)=∞ for any x0∈∂DA, where ∂DA is the boundary of DA.Based on Lemma4.2, the problem of determining DA can be reduced to finding the natural domain of analyticity of the solution V of (4.3). The function V is called the optimal Lyapunov function for (4.2). Generally, it is not easy to construct the optimal Lyapunov function V and determine its domain of analyticity. But, in the diagonalizable case, we can determine the coefficients of the expansion of W=V∘H at 0, where H reduces ∂f/∂x(0) to the diagonal form. By the Cauchy-Hadamard-type theorem [25], we can obtain the domain of convergence D0 of the series W, and DA0=H-1(D0) is a part of the domain of attraction.For system (4.2), the following lemma holds.Lemma 4.3 (see [26]).
For each isomorphismH:ℂn→ℂn and g=H-1∘f∘H, the problem(4.4)〈∇W,g〉=-‖Hz‖2,W(0)=0
has a unique analytical solution; namely, W=V∘H, where V is the optimal Lyapunov function for (4.2).LetH:ℂn→ℂn be one isomorphism which reduces ∂f/∂x(0) to the diagonal form H-1∂f/∂x(0)H=diag(λ1,λ2,…,λn). Denote g=H-1∘f∘H and W=V∘H, where V is the optimal Lyapunov function for (4.2).For givenH, we consider the expansion of W at origin:
(4.5)W(z1,z2,…,zn)=∑m=2∞‍∑|j|=m‍Bj1j2⋯jnz1j1z2j2⋯znjn,
and the expansions of the scalar components gi of g at origin:
(4.6)gi(z1,z2,…,zn)=λizi+∑m=2∞‍∑|j|=m‍bj1j2⋯jniz1j1z2j2⋯znjn,
then the coefficients Bj1j2⋯jnof the development (4.5) are given by
(4.7)Bj1j2⋯jn={-12λi0∑i=1n‍sii02if|j|=ji0=2,-2λp+λq∑i=1n‍sipsiqif|j|=2,jp=jq=1,-1∑i=1njiλiif|j|≥3,×∑p=2|j|-1‍∑|k|=p,ki≤ji‍∑i=1n‍[(ji-ki+1)×bk1k2⋯kniBj1-k1⋯ji-ki+1⋯jn-kn],
where λi0 is the i0th eigenvalue of ∂f/∂x(0), i0=1,2,…,n, sip(siq) represents the entry of H which lies in the ith row and the pth(qth) column.LetVp, p≥2 be the Taylor polynomials of degree p associated to V at the origin. For Vp, we have the following results [27].Lemma 4.4.
For anyp≥2, there exists rp>0 such that for any x∈B¯(rp)∖{0} one has
(4.8)Vp(x)>0,〈∇Vp(x),f(x)〉<0,
where B(rp)={x∈ℝn∣∥x∥<rp} and B¯(rp) is its closure.Hence, for any positive integerp≥2, Vp is a Lyapunov function for system (4.2) in B(rp).Now, we construct the Lyapunov function at the virus-free equilibrium. Define the following transformation:(4.9)[z1z2]=H-1[x1x2],H=[1B1-n0-1]:=[s11s12s21s22].Then system (2.9) is transformed into the following form:
(4.10)dzdt=g(z),
where
(4.11)z=[z1z2],g(z)=[-nz1+(1+m-B1-n)z1z2+(B(1+m)1-n-1-(B1-n)2m)z22-z2+mz1z2+Bm1-nz22].
For the isomorphism H:ℂ2→ℂ2 and the function g:ℝ2→ℝ2, it follows from Lemma 4.3 that there exists a unique analytical function W=V∘H such that 〈∇W,g〉=-∥Hz∥2 and W(0)=0, where V is the optimal Lyapunov function for the reduced system (4.10).Let the expansion ofW at (0,0) be
(4.12)W(z1,z2)=∑m=2∞‍∑|j|=m‍Bj1j2z1j1z2j2,
and denote the components gi of g as
(4.13)gi(z1,z2)=λizi+∑m=2∞‍∑|j|=m‍bj1j2iz1j1z2j2,i=1,2.
The coefficients of expansion for W are given by
(4.14)Bj1j2={-12λi0∑i=12‍sii02if|j|=ji0=2,-2λp+λq∑i=12‍sipsiqif|j|=2,jp=jq=1,-1∑i=12jiλi∑p=2|j|-1‍∑|k|=p,ki≤ji‍Φif|j|≥3,
where Φ=(j1-k1+1)×bk1k21Bj1-k1+1j2-k2+(j2-k2+1)×bk1k22Bj1-k1j2-k2+1.Letλ1=-n, λ2=-1; we obtain
(4.15)B20=-12λ1(s112+s212)=12n,B11=-2λ1+λ2(s11s12+s21s22)=2B1-n2,B02=-12λ2(s122+s222)=12(B2(1-n)2+12).
Let W2(z1,z2) be the sum of the square terms of the expansion (4.12). Clearly,
(4.16)W2(z1,z2)=B20z12+B11z1z2+B02z22.
Correspondingly, by the reversibility of H, we obtain an analytical function
(4.17)V2(x1,x2)=(W2∘H-1)(x1,x2)=B20(x1+B1-nx2)2+B11(x1+B1-nx2)(-x2)+B02(-x2)2=k1x12+k2x1x2+k3x22,
where
(4.18)k1=12n>0,k2=Bn(1+n)>0,k3=B22n(1+n)+12>0.
By using Lemma 4.4, it is easy to see that there exists a domain of the origin and in which V2 is a positive definite function and further, it provides a Lyapunov function to guarantee that the trivial equilibrium E0 of system (2.9) is asymptotically stable. Obviously, this domain belongs to DA.Next step, we will use the Lyapunov functionV2 to estimate the domain of attraction for the virus-free equilibrium. For this, we first introduce the following results.Lemma 4.5 (see [28]).
LetV(x) be a Lyapunov function for system (4.2) in the domain
(4.19)Ωc={x∈Rn∣V(x)≤c,c>0}.
Assume, moreover, that Ωc is bounded and contains the origin. If V˙ is negative definite in Ωc, then the origin is asymptotically stable and every solution in Ωc tends to the origin as t→∞.To determinate the whole domain of attractionDA of a general nonlinear system is an extremely difficult problem. However, by Lemma 4.5 it is possible to compute subsets of the domain of attraction which are defined by (4.19). Our objective is to find the maximum value c* of c [29]. This c* is defined by the following optimization problem.
(4.20)c*=min⁡V(x),x∈{Rn∣dV(x)dt=0,x≠0}.Geometrically, this means that we seek the global minimum of the Lyapunov functionV(x) along the hypersurface {x∈ℝn∣dV(x)/dt=0,x≠0}.Correspondingly, for system (2.9) and the determined Lyapunov function V2, the optimization problem (4.20) should be reformulated as
(4.21)c*=min⁡V2(x),x∈{R2∣dV2(x)dt=0,x≠0},
where V2 is defined by (4.17).Clearly, (4.21) is a linear matrix inequality (LMI) relaxations of global optimization problem with multivariable real-valued polynomial criterion and constraints (see [30, 31] and the references therein for details).ThenΩc* is the estimation of domain of attraction for E0.LetΛ1=(-[(1+m)x22*+n]+Δ)/2 and Λ2=(-[(1+m)x22*+n]-Δ)/2, where Δ=(1+m)2(x22*)2+2n(1+m)x22*+n2+4(Bm-m-1)x22*-8m(x22*)2>0(defined by (3.22)), analogously, for the positive equilibrium E2*, after tedious computation, we can obtain the Lyapunov function:
(4.22)V2(x1,x2)=k1′(x1-x1*)2+k2′(x1-x1*)(x2-x22*)+k3′(x2-x22*)2,
where
(4.23)k1′=-12Δ(Λ1-Λ2)2Λ1+Λ2(mx22*)2+Λ1Λ2Λ1Λ2,k2′=1Δ(Λ1-Λ2)2Λ1Λ2,k3′=-12Δ(Λ1-Λ2)2Λ1+Λ2(mx22*)2((Λ1+Λ2)2+Λ1Λ2)+Λ12Λ22(mx22*)2Λ1Λ2.
From the discussion in Section 3, we can obtain
(4.24)Λ1+Λ2=-[(1+m)x22*+n]<0,Λ1Λ2=mx22*(B-1m-1)2-4nm>0;
hence k1′, k2′, and k3′ are all positive.Correspondingly, we can obtain the estimation of domain of attraction forE2*.


## 5. Numerical Examples

In this section, we will perform numerical examples to verify the results that obtained in previous sections, hence, provide an intuitive impression about saddle-node bifurcation, bistability, and the domain of attraction.Firstly, we present the saddle-node bifurcation. Here we emphasize that the virus-free equilibrium is locally asymptotically stable for any parameter values, and it can be observed in Figures2, 3, and 4 simultaneously (the green asterisk in Figures 2–4). In Figure 2, k=3, p=0.65, δ=0.6, d=0.05, po=0.35, and a=1; we have B=3, m=0.66274294783454, n=0.07692307692308, which satisfy the condition B<2n/m+1/m+1. In this case, system (2.9) has only virus-free equilibrium, and it is asymptotically stable.Figure 2
The phase diagram of system (2.9) when B<2n/m+1/m+1.Figure 3
The phase diagram of system (2.9) when B=2n/m+1/m+1(saddle node).Figure 4
The phase diagram of system (2.9) when B>2n/m+1/m+1(bistability).Settingk=4, p=0.4005, and keeping other parameter values unchanged, we have B=2.46455056685463, m=1.07561277426331, n=0.07692307692308, then B=2n/m+1/m+1. According to Theorem 3.2, system (2.9) has a positive equilibrium, and it is a saddle node. If p=0.52, we have B=3.2, m=0.82842868479318, n=0.07692307692308, then B>2n/m+1/m+1. According to Theorem 3.3, there have been two positive equilibria; one of them is a saddle and the other is a node.Above all, system (2.9) experiences a saddle-node bifurcation which can be seen obviously from Figures 2–4. After crossing the critical value, the saddle node (the red asterisk in Figure 3) becomes one saddle and one node (the red pentagrams in Figure 4). Consequently, the bistable case occurs.Secondly, we show the domain of attraction for the equilibria. From Figure4, we can see that system (2.9) has two locally asymptotically stable equilibria under the condition B>2n/m+1/m+1. Letting the parameters be the same as in Figure 4, we give the estimated results of domains of attraction (Figures 5 and 6). From (4.17), we get the Lyapunov function corresponding to E0 as follows:
(5.1)V2E0=6.5x12+38.62857143x1x2+62.30571430x22,
and its total derivative along the trajectory of system (2.9) has the form
(5.2)dV2E0dt=-x12-x22+8.231443707x12x2+19.60209380x1x22-38.62857143x23.
Then solving (4.21) by using the MATLAB tool, we obtain c1*=2.8989975698 for E0. The corresponding region is V2(x)E0=c1* which is shown in Figure 5.Figure 5
The white region is the estimation of domain of attraction forE0.Figure 6
The white region is the estimation of domain of attraction forE2*.In analogy with the virus-free equilibriumE0, we obtain the Lyapunov function of E2* as follows:
(5.3)V2E2*=0.5698434628(x1-x1*)2+1.327812551(x1-x1*)(x2-x22*)+2.081459239(x2-x22*)2,dV2E2*dt=-0.9616845362(x1-x1*)2-0.088571026(x1-x1*)(x2-x22*)-1.040809141(x2-x22*)2-0.983838261(x1-x1*)2(x2-x22*)-0.118816403(x1-x1*)(x2-x22*)2-1.327812551(x2-x22*)3,c2*=0.9898659789. The corresponding region is V2E2*=c2* which is shown in Figure 6.


## 6. Conclusion and Discussion

In this paper, we have established a model to study email virus propagation and analyzed its dynamical behaviors. It is worthwhile to note that we useβEI/N instead of βE to depict the transition from E to I, which is different from most SEIS infectious disease models. For example, Fan et al. [32] proposed an SEIS model as
(6.1)dS(t)dt=A-λS(t)I(t)+δI(t)-dS(t),dE(t)dt=λS(t)I(t)-ɛE(t)-dE(t),dI(t)dt=ɛE(t)-dI(t)-δI(t)-θI(t).It is shown that the global dynamics is completely determined by the basic reproduction number R0. If R0≤1, the disease-free equilibrium is globally stable. If R0>1, a unique endemic equilibrium is globally stable in the interior of the feasible region.The reason of usingβEI/N is that we have considered the number of virus emails that the user received in our model, which is proportional to the relative number of infected computers, that is, I/N. Due to this difference in the form, our model shows many different and complex dynamics. In our model, however, the basic reproduction number R0 does not appear; hence the virus-free equilibrium is locally asymptotically stable for any parameter values. And it is proved that the model experiences a saddle-node bifurcation as parameter varies through the bifurcation value, which also occurs in vector-borne disease model [33]. Furthermore, the bistability is obtained. In short, the spreading mechanism and dynamical features of email virus are indeed different from most biological virus. This can offer a valuable perspective to the future dynamics modeling. Especially, the amount of virus is considerable in the transition from E to I when modeling the infectious disease. By doing so, the models would be beneficial for better understanding the infectious disease propagation.The main aim of the model is to investigate how email viruses transmit and to identify effective strategies for their prevention and control. So based on the arguments above and the simulation results, it is straightforward to provide some control strategies.From the theoretical analyses and Figures2–4, we can take
(6.2)R=B-2nm-1m-1=kpd+δ-2dpaln(1/(1-po))(d+δ)-paln(1/(1-po))-1
as a threshold parameter. If R<0, the virus will die out; if R>0, virus may die out or outbreak.In the first case, to makeR below than zero, the most direct way is to change the parameter values which correspond to the related measures. Obviously, R increases with k, a, and po. The relationships between R and p, δ are shown in Figure 7; R increases with p, while decreases with δ. According to the sensitivity analysis, some measures can be taken to control virus propagation for instance, enhancing the people's understanding of email viruses to minimize the probability of opening a virus email (po) and appealing people to install the antivirus software and update in time to increase the ability of finding and removing the viruses (δ).Figure 7
The combined influence of parameters onR; R is plotted as a function of p and δ. Here, k=6, a=1, po=0.5, d=0.05.In the second case, from the definition of the domain of attraction, we know that starting from any initial points in the regionΩc* of an equilibrium, the solution of a dynamical system will converge to the equilibrium. As is well known, the solution tending to the virus-free equilibrium indicates that the virus will be extinct without taking any measures, whereas tending to the positive equilibrium indicates the outbreak of the virus; that is, a large number of computers will be infected. As a result, in order to control the propagation of the virus, we should take measures to restrict the initial value of the computers within the domain of attraction for virus-free equilibrium (the white region in Figure 5), such as improving the computers' detection ability and enhancing the people's awareness of email virus defense.It should be noted that the epidemiological model described here for email viruses could also be used to study other disease similar to email viruses from an epidemiological point of view. Moreover, in our model, we only consider some fundamental factors. In the future work; more factors will be incorporated into new and improved models for prediction and containment of email viruses.



---

*Source: 472072-2012-08-23.xml*