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.

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 [59]. 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 [1012] 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 Section 2, 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 as 𝑘; that is, an email user has 𝑘 addresses in his own address book and his address is in 𝑘 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 time 𝑡 are denoted as 𝑆(𝑡), 𝐸(𝑡), and 𝐼(𝑡), respectively. Then the total number of computers at time 𝑡 is 𝑁(𝑡)=𝑆(𝑡)+𝐸(𝑡)+𝐼(𝑡). Here we assume the total number of computers is varying dynamically as time 𝑡; 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 𝐴 is the rate at which new computers are connected and 𝑑 is the rate at which the computers are removed per unit time. The flow diagram of the state transition is depicted in Figure 1.

By the above analysis, we have already known the mechanism of the virus propagation. In the following, we will elaborate the probabilities of transforming from 𝑆 to 𝐸, 𝐸 to 𝐼, and 𝐼 to 𝑆, 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 computer 𝑈 is in the susceptible state at time 𝑡. By the usual homogeneous mixing approximation, Θ(𝑡)=𝐼(𝑡)/𝑁(𝑡) denotes the rate that a neighbor of the computer 𝑈 is in the infectious state at time 𝑡. Suppose that 𝑘 neighbors of computer 𝑈 are independent with each other, we have 𝑘 independent events occurring with the same probability Θ(𝑡). Then the number of the infected neighbors of computer 𝑈 is a random variable that follows the binomial distribution 𝐵[𝑘,Θ(𝑡)]. And the probability of the number of the infected neighbors of computer 𝑈 equal to 𝑘𝑖 at time 𝑡 can be written as follows: 𝑃𝑋=𝑘𝑖=𝑘𝑘𝑖Θ(𝑡)𝑘𝑖(1Θ(𝑡))𝑘𝑘𝑖.(2.1) 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 𝑘 neighbors, the susceptible computers would have chance to turn into exposed. Let 𝑝𝑆𝑆 denote the probability that computer 𝑈 still stays in the susceptible state in the time interval [𝑡,𝑡+Δ𝑡], then 1𝑝𝑆𝑆 is the probability that it enters into the exposed state. 𝑘𝑖 is the number of the infected neighbors of computer 𝑈, and if 𝑝 denotes the probability that the computer 𝑈 received virus emails from each infected neighbor per unit time, then we have 𝑝𝑆𝑆𝑘𝑖,𝑡=(1Δ𝑡𝑝)𝑘𝑖.(2.2)

Because the infected number 𝑘𝑖 of the computer's 𝑘 neighbors can change from 0 to 𝑘 and the corresponding 𝑝𝑆𝑆(𝑘𝑖,𝑡) changes from 1 to (1Δ𝑡𝑝)𝑘, using (2.1) and (2.2), we obtain the transition probability 𝑝𝑆𝑆(𝑡) averaged over all possible values of 𝑘𝑖 as 𝑝𝑆𝑆𝑝(𝑡)=𝐸𝑆𝑆𝑘𝑖=,𝑡𝑘𝑘𝑖=0(1Δ𝑡𝑝)𝑘𝑖𝑘𝑘𝑖Θ(𝑡)𝑘𝑖(1Θ(𝑡))𝑘𝑘𝑖=(1Δ𝑡𝑝Θ(𝑡))𝑘.(2.3) Then in the interval [𝑡,𝑡+Δ𝑡], the probability of a susceptible computer transforming to the exposed state is 1(1Δ𝑡𝑝Θ(𝑡))𝑘.

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. Let 𝑝𝑜 be the probability that the email user opens a virus attachment (0<𝑝𝑜<1); thus 1𝑝𝑜 is the probability that the email user does not open a virus attachment. Let the number of the virus emails that the computer 𝑈 received during the time interval [𝑡,𝑡+Δ𝑡] be 𝑛𝑒; then the probability that all 𝑛𝑒 virus attachments would not be opened is (1𝑝𝑜)𝑛𝑒. Therefore, the probability of the computer infected is 1(1𝑝𝑜)𝑛𝑒.

In addition, from (2.1), we can obtain that the average infected number of the computer's neighbors is 𝑘𝑖=𝑘𝑘𝑖=0𝑘𝑖𝑃𝑋=𝑘𝑖=𝑘𝑘𝑖=0𝑘𝑖𝑘𝑘𝑖Θ(𝑡)𝑘𝑖(1Θ(𝑡))𝑘𝑘𝑖=𝑘Θ(𝑡).(2.4) Suppose the computer 𝑈 received 𝑎 virus emails from each infected neighbor per unit time, then in the time interval [𝑡,𝑡+Δ𝑡], the computer 𝑈 received 𝑎𝑘Θ(𝑡)Δ𝑡 virus emails on average; that is, 𝑛𝑒=𝑎𝑘Θ(𝑡)Δ𝑡.

Substituting 𝑛𝑒 into the above formula, the probability of an exposed computer becomes infected is 1(1𝑝𝑜)𝑘𝑎Θ(𝑡)Δ𝑡 in the time interval Δ𝑡.

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 Figure 1, we can indeed write the spread mechanism of email viruses as the following deterministic model: 𝑆(𝑡+Δ𝑡)=𝑆(𝑡)+Δ𝑡𝐴11Δ𝑡𝑝𝐼(𝑡)𝑁(𝑡)𝑘𝑆(𝑡)+Δ𝑡𝛿𝐼(𝑡)Δ𝑡𝑑𝑆(𝑡),𝐸(𝑡+Δ𝑡)=𝐸(𝑡)+11Δ𝑡𝑝𝐼(𝑡)𝑁(𝑡)𝑘𝑆(𝑡)11𝑝𝑜𝑘𝑎(𝐼(𝑡)/𝑁(𝑡))Δ𝑡𝐸(𝑡)Δ𝑡𝑑𝐸(𝑡),𝐼(𝑡+Δ𝑡)=𝐼(𝑡)+11𝑝𝑜𝑘𝑎(𝐼(𝑡)/𝑁(𝑡))Δ𝑡𝐸(𝑡)Δ𝑡𝛿𝐼(𝑡)Δ𝑡𝑑𝐼(𝑡).(2.5) In the limit Δ𝑡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) where 𝛼=𝑘𝑝, 𝛽=𝑘𝑎ln(1/(1𝑝𝑜)).

Summing up the three equations in (2.6), we obtain 𝑑𝑁(𝑡)𝑑𝑡=𝐴𝑑𝑁(𝑡),(2.7) when 𝑡+, 𝑁(𝑡)𝐴/𝑑𝑁0. It is easy to see that system (2.6) can be shown to be mathematically well posed in the positive invariant region 𝔻={(𝑆,𝐸,𝐼)0𝑆+𝐸+𝐼𝑁0} and solution in 𝔻 exists for all positive time.

Hence, (2.6) has the following limit system: 𝑑𝐸(𝑡)=𝛼𝑑𝑡𝑁0𝑁0𝛽𝐸(𝑡)𝐼(𝑡)𝐼(𝑡)𝑁0𝐸(𝑡)𝐼(𝑡)𝑑𝐸(𝑡),𝑑𝐼(𝑡)=𝛽𝑑𝑡𝑁0𝐸(𝑡)𝐼(𝑡)(𝑑+𝛿)𝐼(𝑡).(2.8) Rescaling (2.8) by 𝑥1=(𝛼/𝑁0(𝑑+𝛿))𝐸, 𝑥2=(𝛼/𝑁0(𝑑+𝛿))𝐼, 𝜏=(𝑑+𝛿)𝑡 and still writing 𝜏 as 𝑡, we obtain 𝑑𝑥1=𝑑𝑡𝐵𝑥1𝑥2𝑥2𝑚𝑥1𝑥2𝑛𝑥1,𝑑𝑥2𝑑𝑡=𝑚𝑥1𝑥2𝑥2,(2.9) where 𝐵=𝛼/(𝑑+𝛿),𝑚=𝛽/𝛼,𝑛=𝑑/(𝑑+𝛿), and the corresponding positive invariant set is 𝕏={(𝑥1,𝑥2)0𝑥1+𝑥2𝐵}.

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 𝐷𝑥1,𝑥2=1𝑥2,(3.1) where 𝑥1>0, 𝑥2>0.
Set 𝐵𝑥1𝑥2𝑥2𝑚𝑥1𝑥2𝑛𝑥1𝑥=𝑃1,𝑥2,𝑚𝑥1𝑥2𝑥2𝑥=𝑄1,𝑥2;(3.2) we have 𝜕(𝐷𝑃)𝜕𝑥1+𝜕(𝐷𝑄)𝜕𝑥2𝑛=1𝑚𝑥2<0.(3.3) 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 𝐸0=(0,0), and the Jacobian matrix is 𝐽𝐸0=𝑛𝐵01. Obviously, (0,0) is locally asymptotically stable for any parameter values.

To find the positive equilibrium, set 𝐵𝑥1𝑥2𝑥2𝑚𝑥1𝑥2𝑛𝑥1=0,𝑚𝑥11=0,(3.4) which yields the equation of 𝑥2: 𝑚𝑥22+(1+𝑚𝐵𝑚)𝑥2+𝑛=0.(3.5) We can obtain that(1)there is no positive equilibrium if 𝐵<2𝑛/𝑚+1/𝑚+1;(2)there is one positive equilibrium if 𝐵=2𝑛/𝑚+1/𝑚+1;(3)there are two positive equilibria if 𝐵>2𝑛/𝑚+1/𝑚+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)). Suppose 𝐵=2𝑛/𝑚+1/𝑚+1; system (2.9) has one positive equilibrium 𝐸=𝑥1,𝑥2=1𝑚,121𝐵𝑚;1(3.6) the Jacobian matrix at 𝐸 is 𝐽𝐸=(1+𝑚)𝑥2𝑛0𝑚𝑥20.(3.7) Clearly, the eigenvalues of matrix 𝐽𝐸 are 𝜎1=0 and 𝜎2=(1+𝑚)𝑥2𝑛; 𝐸 is a nonhyperbolic point. Therefore, we investigate the dynamics near 𝐸 by the center manifold theorem [17, 18].
Firstly, shift 𝐸 to the origin via 𝑦1=𝑥1𝑥1 and 𝑦2=𝑥2𝑥2; system (2.9) can be transformed into 𝑑𝑦1=𝑑𝑡𝑛(1+𝑚)𝑥2𝑦1(1+𝑚)𝑦1𝑦2𝑦22,𝑑𝑦2𝑑𝑡=𝑚𝑥2𝑦1+𝑚𝑦1𝑦2.(3.8) Secondly, define the transformation 𝑦1𝑦2𝑧=𝐻1𝑧21,𝐻=01𝑚𝑥2𝜎2,(3.9) which transformed system (3.8) into the following standard form: 𝑑𝑧1𝑑𝑡=𝑓1𝑧1,𝑧2,𝑑𝑧2𝑑𝑡=𝜎2𝑧2+𝑓2𝑧1,𝑧2,(3.10) where 𝑓1=(1+𝑚)(𝑚𝑥2/𝜎2)(𝑧1𝑧2+(𝑚𝑥2/𝜎2)𝑧22)+(𝑚𝑥2/𝜎2)(𝑧1+(𝑚𝑥2/𝜎2)𝑧2)2+𝑚(𝑧1𝑧2+(𝑚𝑥2/𝜎2)𝑧22) and 𝑓2=(1+𝑚)(𝑧1𝑧2+(𝑚𝑥2/𝜎2)𝑧22)(𝑧1+(𝑚𝑥2/𝜎2)𝑧2)2.
By the existence theorem [19], there exists a center manifold for system (3.10), which can be expressed locally as follows: 𝑊𝑐𝑧(0)=1,𝑧22𝑧2𝑧=1,𝑧1,<𝜌,(0)=0,𝐷(0)=0𝜌>0(3.11) with 𝜌 sufficiently small and 𝐷 is the derivative of with respect to 𝑧1.
To compute the center manifold 𝑊𝑐(0), we suppose (𝑧1) has the form 𝑧1=1𝑧21+2𝑧31+3𝑧41+4𝑧51+.(3.12) By the local center manifold theorem, the center manifold (3.12) satisfies 𝐷𝑓1𝑧1,𝑧2𝜎2𝑧2𝑓2𝑧1,𝑧2=0,(3.13) where 𝐷(𝑧1)=21𝑧1+32𝑧21+43𝑧31+54𝑧41+.
Rewrite 𝑓1 and 𝑓2 as 𝑓1=𝑎1𝑧21+𝑎2𝑧1𝑧2+𝑎3𝑧22 and 𝑓2=𝑏1𝑧21+𝑏2𝑧1𝑧2+𝑏3𝑧22, respectively, where 𝑎1=𝑚𝑥2𝜎2,𝑎2=(1+𝑚)𝑚𝑥2𝜎2+2𝑚𝑥2𝜎22𝑎+𝑚,3=(1+𝑚)𝑚𝑥2𝜎22+𝑚𝑥2𝜎23+𝑚2𝑥2𝜎2,𝑏1=1,𝑏2=(1+𝑚)2𝑚𝑥2𝜎2,𝑏3=(1+𝑚)𝑚𝑥2𝜎2𝑚𝑥2𝜎22.(3.14) Substituting (3.10) into (3.13) and then equating coefficients on each power of 𝑧1 to zero yields 1=1𝜎2,2=2𝑎11𝑏21𝜎2,3=2𝑎121+3𝑎12𝑏22𝑏321𝜎2,4=5𝑎212+2𝑎331+4𝑎13𝑏232𝑏312𝜎2.(3.15) Then, we get the approximation for : 1=𝜎2𝑧21+2𝑎11𝑏21𝑧31+2𝑎121+3𝑎12𝑏22𝑏321𝑧41+5𝑎212+2𝑎331+4𝑎13𝑏232𝑏312𝑧51+.(3.16)
Substituting (3.16) to the first equation of system (3.10), we achieve 𝑑𝑧1𝑑𝑡=𝑓1𝑧1,𝑧2=𝑚𝑥2𝜎2𝑧21+(1+𝑚)𝑚𝑥2𝜎2+2𝑚𝑥2𝜎221+𝑚𝜎2𝑧31+.(3.17) So from (3.17) (see [20, p. 338–340]), we get the following theorem.

Theorem 3.2. If 𝐵=2𝑛/𝑚+1/𝑚+1, the system (2.9) has one positive equilibrium 𝐸 and 𝐸 is a saddle node.

Case 2 (two positive equilibria in system (2.9)). Suppose 𝐵>2𝑛/𝑚+1/𝑚+1, there are two positive equilibria 𝐸1=𝑥1,𝑥21=1𝑚,𝐵1/𝑚1(𝐵1/𝑚1)24𝑛/𝑚2,𝐸2=𝑥1,𝑥22=1𝑚,𝐵1/𝑚1+(𝐵1/𝑚1)24𝑛/𝑚2.(3.18) We first determine the stability of 𝐸1. The Jacobian matrix at 𝐸1 is 𝐽𝐸1=(1+𝑚)𝑥21𝑛1𝐵𝑚12𝑛4𝑚𝑚𝑥210.(3.19) Obviously, det(𝐽𝐸1)=𝑚𝑥21(𝐵1/𝑚1)24𝑛/𝑚<0, so 𝐸1 is a saddle.
Now we analyze the stability of the second positive equilibrium 𝐸2. The Jacobian matrix at 𝐸2 is 𝐽𝐸2=(1+𝑚)𝑥221𝑛𝐵𝑚12𝑥22𝑚𝑥220.(3.20) By a similar argument as above, we obtain that 𝐽det𝐸2=𝑚𝑥221𝐵𝑚12𝑥22=𝑚𝑥221𝐵𝑚12𝑛4𝑚𝐽>0,trace𝐸2=(1+𝑚)𝑥22𝑛<0,(3.21)Δ=trace2𝐽𝐸2𝐽4det𝐸2=(1+𝑚)2𝑥222+2𝑛(1+𝑚)𝑥22+𝑛2+4(𝐵𝑚𝑚1)𝑥22𝑥8𝑚222.(3.22) After some algebra calculations, we obtain Δ>0; 𝐸2 is a stable node. So we have the following theorem.

Theorem 3.3. If 𝐵>2𝑛/𝑚+1/𝑚+1, the system (2.9) has two positive equilibria 𝐸1 and 𝐸2; 𝐸1 is a saddle and 𝐸2 is a stable node.

Case 3 (saddle-node bifurcation). From Cases 1 and 2, we can see system (2.9) experiences a saddle-node bifurcation. In this part, we give the proof.
Let us consider 𝐵 as a bifurcation parameter, define 𝐵0=2𝑛/𝑚+1/𝑚+1 and 𝑥0=𝐸. Rewrite system (2.9) as ̇𝑥=𝑓(𝑥,𝐵)=𝑛𝑥1+𝐵𝑥2(1+𝑚)𝑥1𝑥2𝑥22𝑥2+𝑚𝑥1𝑥2.(3.23) Then we have 𝑥𝐷𝑓0,𝐵0=(1+𝑚)𝑥2𝑛0𝑚𝑥20,𝑓𝐵𝑥0,𝐵0=𝑥20.(3.24)𝐷𝑓(𝑥0,𝐵0) has a simple eigenvalue 𝜆=0 with eigenvector 𝑣=[0,1]𝑇, and that 𝐷𝑇𝑓(𝑥0,𝐵0) has an eigenvector 𝑤=[𝑚𝑥2/((1+𝑚)𝑥2+𝑛),1]𝑇 corresponding to 𝜆=0. Furthermore, the following conditions are satisfied: 𝑤𝑇𝑓𝐵𝑥0,𝐵0=𝑚𝑥22(1+𝑚)𝑥2+𝑛0,𝑤𝑇𝐷2𝑓𝑥0,𝐵0(𝑣,𝑣)=2𝑚𝑥2(1+𝑚)𝑥2+𝑛0.(3.25) 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 𝑥0=𝐸 as the parameter 𝐵 passes through the bifurcation value 𝐵=𝐵0.

4. Estimation of the Domain of Attraction

From the previous study, we know, when 𝐵>2𝑛/𝑚+1/𝑚+1, 𝐸2 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 𝑥𝕊=0𝑛lim𝑡𝑥𝑡,𝑥00.(4.1)
Consider the following system: ̇𝑥=𝑓(𝑥),(4.2) where 𝑓𝑛𝑛 is an analytical function with the following properties:(1)𝑓(0)=0; that is, 𝑥=0 is an equilibrium point of (4.2);(2)every eigenvalue of 𝜕𝑓/𝜕𝑥(0), the Jacobian matrix of 𝑓 at the origin, has negative real part; that is, 𝑥=0 is an asymptotically stable equilibrium point.
Let 𝐷𝐴 be the domain of attraction of the zero solution of (4.2). The following result provides a tool to determine 𝐷𝐴.

Lemma 4.2 (see [24]). 𝐷𝐴 coincides with the natural domain of analyticity of the unique 𝑉 of the problem 𝑉,𝑓=𝑥2,𝑉(0)=0.(4.3) The function 𝑉 is positive on 𝐷𝐴 and lim𝑥𝑥0𝑉(𝑥)= for any 𝑥0𝜕𝐷𝐴, where 𝜕𝐷𝐴 is the boundary of 𝐷𝐴.

Based on Lemma 4.2, the problem of determining 𝐷𝐴 can be reduced to finding the natural domain of analyticity of the solution 𝑉 of (4.3). The function 𝑉 is called the optimal Lyapunov function for (4.2). Generally, it is not easy to construct the optimal Lyapunov function 𝑉 and determine its domain of analyticity. But, in the diagonalizable case, we can determine the coefficients of the expansion of 𝑊=𝑉𝐻 at 0, where 𝐻 reduces 𝜕𝑓/𝜕𝑥(0) to the diagonal form. By the Cauchy-Hadamard-type theorem [25], we can obtain the domain of convergence 𝐷0 of the series 𝑊, and 𝐷𝐴0=𝐻1(𝐷0) is a part of the domain of attraction.

For system (4.2), the following lemma holds.

Lemma 4.3 (see [26]). For each isomorphism 𝐻𝑛𝑛 and 𝑔=𝐻1𝑓𝐻, the problem𝑊,𝑔=𝐻𝑧2,𝑊(0)=0(4.4) has a unique analytical solution; namely, 𝑊=𝑉𝐻, where 𝑉 is the optimal Lyapunov function for (4.2).

Let 𝐻𝑛𝑛 be one isomorphism which reduces 𝜕𝑓/𝜕𝑥(0) to the diagonal form 𝐻1𝜕𝑓/𝜕𝑥(0)𝐻=diag(𝜆1,𝜆2,,𝜆𝑛). Denote g=𝐻1𝑓𝐻 and 𝑊=𝑉𝐻, where 𝑉 is the optimal Lyapunov function for (4.2).

For given 𝐻, we consider the expansion of 𝑊 at origin: 𝑊𝑧1,𝑧2,,𝑧𝑛=𝑚=2||𝑗||=𝑚𝐵𝑗1𝑗2𝑗𝑛𝑧𝑗11𝑧𝑗22𝑧𝑗𝑛𝑛,(4.5) and the expansions of the scalar components 𝑔𝑖 of 𝑔 at origin: 𝑔𝑖𝑧1,𝑧2,,𝑧𝑛=𝜆𝑖𝑧𝑖+𝑚=2||𝑗||=𝑚𝑏𝑖𝑗1𝑗2𝑗𝑛𝑧𝑗11𝑧𝑗22𝑧𝑗𝑛𝑛,(4.6) then the coefficients 𝐵𝑗1𝑗2𝑗𝑛of the development (4.5) are given by 𝐵𝑗1𝑗2𝑗𝑛=12𝜆𝑖0𝑛𝑖=1𝑠2𝑖𝑖0||𝑗||if=𝑗𝑖02=2,𝜆𝑝+𝜆𝑞𝑛𝑖=1𝑠𝑖𝑝𝑠𝑖𝑞||𝑗||if=2,𝑗𝑝=𝑗𝑞1=1,𝑛𝑖=1𝑗𝑖𝜆𝑖||𝑗||×if3,|𝑗|1𝑝=2||𝑘||=𝑝,𝑘𝑖𝑗𝑖𝑛𝑖=1𝑗𝑖𝑘𝑖+1×𝑏𝑖𝑘1𝑘2𝑘𝑛𝐵𝑗1𝑘1𝑗𝑖𝑘𝑖+1𝑗𝑛𝑘𝑛,(4.7) where 𝜆𝑖0 is the 𝑖0th eigenvalue of 𝜕𝑓/𝜕𝑥(0), 𝑖0=1,2,,𝑛, 𝑠𝑖𝑝(𝑠𝑖𝑞) represents the entry of 𝐻 which lies in the 𝑖th row and the 𝑝th(𝑞𝑡) column.

Let 𝑉𝑝, 𝑝2 be the Taylor polynomials of degree 𝑝 associated to 𝑉 at the origin. For 𝑉𝑝, we have the following results [27].

Lemma 4.4. For any 𝑝2, there exists 𝑟𝑝>0 such that for any 𝑥𝐵(𝑟𝑝){0} one has 𝑉𝑝(𝑥)>0,𝑉𝑝(𝑥),𝑓(𝑥)<0,(4.8) where 𝐵(𝑟𝑝)={𝑥𝑛𝑥<𝑟𝑝} and 𝐵(𝑟𝑝) is its closure.

Hence, for any positive integer 𝑝2, 𝑉𝑝 is a Lyapunov function for system (4.2) in 𝐵(𝑟𝑝).

Now, we construct the Lyapunov function at the virus-free equilibrium. Define the following transformation: 𝑧1𝑧2=𝐻1𝑥1𝑥21𝐵,𝐻=𝑠1𝑛01=11𝑠12𝑠21𝑠22.(4.9)

Then system (2.9) is transformed into the following form: 𝑑𝑧𝑑𝑡=𝑔(𝑧),(4.10) where 𝑧𝑧=1𝑧2,𝑔(𝑧)=𝑛𝑧1+𝐵1+𝑚𝑧1𝑛1𝑧2+𝐵(1+𝑚)𝐵1𝑛11𝑛2𝑚𝑧22𝑧2+𝑚𝑧1𝑧2+𝐵𝑚𝑧1𝑛22.(4.11) For the isomorphism 𝐻22 and the function 𝑔22, it follows from Lemma 4.3 that there exists a unique analytical function 𝑊=𝑉𝐻 such that 𝑊,𝑔=𝐻𝑧2 and 𝑊(0)=0, where 𝑉 is the optimal Lyapunov function for the reduced system (4.10).

Let the expansion of 𝑊 at (0,0) be 𝑊𝑧1,𝑧2=𝑚=2||𝑗||=𝑚𝐵𝑗1𝑗2𝑧𝑗11𝑧𝑗22,(4.12) and denote the components 𝑔𝑖 of 𝑔 as 𝑔𝑖𝑧1,𝑧2=𝜆𝑖𝑧𝑖+𝑚=2||𝑗||=𝑚𝑏𝑖𝑗1𝑗2𝑧𝑗11𝑧𝑗22,𝑖=1,2.(4.13) The coefficients of expansion for 𝑊 are given by 𝐵𝑗1𝑗2=12𝜆𝑖02𝑖=1𝑠2𝑖𝑖0||𝑗||if=𝑗𝑖02=2,𝜆𝑝+𝜆𝑞2𝑖=1𝑠𝑖𝑝𝑠𝑖q||𝑗||if=2,𝑗𝑝=𝑗𝑞1=1,2𝑖=1𝑗𝑖𝜆𝑖|𝑗|1𝑝=2||𝑘||=𝑝,𝑘𝑖𝑗𝑖||𝑗||Φif3,(4.14) where Φ=(j1𝑘1+1)×𝑏1𝑘1𝑘2𝐵𝑗1𝑘1+1𝑗2𝑘2+(𝑗2𝑘2+1)×𝑏2𝑘1𝑘2𝐵𝑗1𝑘1𝑗2𝑘2+1.

Let 𝜆1=𝑛, 𝜆2=1; we obtain 𝐵201=2𝜆1𝑠211+𝑠221=1,𝐵2𝑛112=𝜆1+𝜆2𝑠11𝑠12+𝑠21𝑠22=2𝐵1𝑛2,𝐵021=2𝜆2𝑠212+𝑠222=12𝐵2(1𝑛)2+12.(4.15) Let 𝑊2(𝑧1,𝑧2) be the sum of the square terms of the expansion (4.12). Clearly, 𝑊2𝑧1,𝑧2=𝐵20𝑧21+𝐵11𝑧1𝑧2+𝐵02𝑧22.(4.16) Correspondingly, by the reversibility of 𝐻, we obtain an analytical function 𝑉2𝑥1,𝑥2=𝑊2𝐻1𝑥1,𝑥2=𝐵20𝑥1+𝐵𝑥1𝑛22+𝐵11𝑥1+𝐵𝑥1𝑛2𝑥2+𝐵02𝑥22=𝑘1𝑥21+𝑘2𝑥1𝑥2+𝑘3𝑥22,(4.17) where 𝑘1=12𝑛>0,𝑘2=𝐵𝑛(1+𝑛)>0,𝑘3=𝐵2+12𝑛(1+𝑛)2>0.(4.18) By using Lemma 4.4, it is easy to see that there exists a domain of the origin and in which 𝑉2 is a positive definite function and further, it provides a Lyapunov function to guarantee that the trivial equilibrium 𝐸0 of system (2.9) is asymptotically stable. Obviously, this domain belongs to 𝐷𝐴.

Next step, we will use the Lyapunov function 𝑉2 to estimate the domain of attraction for the virus-free equilibrium. For this, we first introduce the following results.

Lemma 4.5 (see [28]). Let 𝑉(𝑥) be a Lyapunov function for system (4.2) in the domain Ω𝑐={𝑥𝑛𝑉(𝑥)𝑐,𝑐>0}.(4.19) Assume, moreover, that Ω𝑐 is bounded and contains the origin. If ̇𝑉 is negative definite in Ω𝑐, then the origin is asymptotically stable and every solution in Ω𝑐 tends to the origin as 𝑡.

To determinate the whole domain of attraction 𝐷𝐴 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 𝑐 of 𝑐 [29]. This 𝑐 is defined by the following optimization problem. 𝑐=min𝑉(𝑥),𝑥𝑛𝑑𝑉(𝑥).𝑑𝑡=0,𝑥0(4.20)

Geometrically, this means that we seek the global minimum of the Lyapunov function 𝑉(𝑥) along the hypersurface {𝑥𝑛𝑑𝑉(𝑥)/𝑑𝑡=0,𝑥0}.

Correspondingly, for system (2.9) and the determined Lyapunov function 𝑉2, the optimization problem (4.20) should be reformulated as 𝑐=min𝑉2(𝑥),𝑥2𝑑𝑉2(𝑥),𝑑𝑡=0,𝑥0(4.21) where 𝑉2 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 Ω𝑐 is the estimation of domain of attraction for 𝐸0.

Let Λ1=([(1+𝑚)𝑥22+𝑛]+Δ)/2 and Λ2=([(1+𝑚)𝑥22+𝑛]Δ)/2, where Δ=(1+𝑚)2(𝑥22)2+2𝑛(1+𝑚)𝑥22+𝑛2+4(𝐵𝑚𝑚1)𝑥228𝑚(𝑥22)2>0(defined by (3.22)), analogously, for the positive equilibrium 𝐸2, after tedious computation, we can obtain the Lyapunov function: 𝑉2𝑥1,𝑥2=𝑘1𝑥1𝑥12+𝑘2𝑥1𝑥1𝑥2𝑥22+𝑘3𝑥2𝑥222,(4.22) where 𝑘11=Λ2Δ1Λ22Λ1+Λ2𝑚𝑥222+Λ1Λ2Λ1Λ2,𝑘21=ΔΛ1Λ22Λ1Λ2,𝑘31=Λ2Δ1Λ22Λ1+Λ2𝑚𝑥222Λ1+Λ22+Λ1Λ2+Λ21Λ22𝑚𝑥222Λ1Λ2.(4.23) From the discussion in Section 3, we can obtain Λ1+Λ2=(1+𝑚)𝑥22+𝑛<0,Λ1Λ2=𝑚𝑥221𝐵𝑚12𝑛4𝑚>0;(4.24) hence 𝑘1, k2, and 𝑘3 are all positive.

Correspondingly, we can obtain the estimation of domain of attraction for 𝐸2.

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 Figures 2, 3, and 4 simultaneously (the green asterisk in Figures 24). In Figure 2, 𝑘=3, 𝑝=0.65, 𝛿=0.6, 𝑑=0.05, 𝑝𝑜=0.35, and 𝑎=1; we have 𝐵=3, 𝑚=0.66274294783454, 𝑛=0.07692307692308, which satisfy the condition 𝐵<2𝑛/𝑚+1/𝑚+1. In this case, system (2.9) has only virus-free equilibrium, and it is asymptotically stable.

Setting 𝑘=4, 𝑝=0.4005, and keeping other parameter values unchanged, we have 𝐵=2.46455056685463, 𝑚=1.07561277426331, 𝑛=0.07692307692308, then 𝐵=2𝑛/𝑚+1/𝑚+1. According to Theorem 3.2, system (2.9) has a positive equilibrium, and it is a saddle node. If 𝑝=0.52, we have 𝐵=3.2, 𝑚=0.82842868479318, 𝑛=0.07692307692308, then 𝐵>2𝑛/𝑚+1/𝑚+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 24. 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 Figure 4, we can see that system (2.9) has two locally asymptotically stable equilibria under the condition 𝐵>2𝑛/𝑚+1/𝑚+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 𝐸0 as follows: 𝑉𝐸02=6.5𝑥21+38.62857143𝑥1𝑥2+62.30571430𝑥22,(5.1) and its total derivative along the trajectory of system (2.9) has the form 𝑑𝑉𝐸02𝑑𝑡=𝑥21𝑥22+8.231443707𝑥21𝑥2+19.60209380𝑥1𝑥2238.62857143𝑥32.(5.2) Then solving (4.21) by using the MATLAB tool, we obtain 𝑐1=2.8989975698 for 𝐸0. The corresponding region is 𝑉2(𝑥)𝐸0=𝑐1 which is shown in Figure 5.

In analogy with the virus-free equilibrium 𝐸0, we obtain the Lyapunov function of 𝐸2 as follows: 𝑉𝐸22𝑥=0.56984346281𝑥12𝑥+1.3278125511𝑥1𝑥2𝑥22𝑥+2.0814592392𝑥222,𝑑𝑉𝐸22𝑥𝑑𝑡=0.96168453621𝑥12𝑥0.0885710261𝑥1𝑥2𝑥22𝑥1.0408091412𝑥222𝑥0.9838382611𝑥12𝑥2𝑥22𝑥0.1188164031𝑥1𝑥2𝑥222𝑥1.3278125512𝑥223,(5.3)𝑐2=0.9898659789. The corresponding region is 𝑉𝐸22=𝑐2 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 𝛽𝐸𝐼/𝑁 instead of 𝛽𝐸 to depict the transition from 𝐸 to 𝐼, which is different from most 𝑆𝐸𝐼𝑆 infectious disease models. For example, Fan et al. [32] proposed an 𝑆𝐸𝐼𝑆 model as 𝑑𝑆(𝑡)𝑑𝑡=𝐴𝜆𝑆(𝑡)𝐼(𝑡)+𝛿𝐼(𝑡)𝑑𝑆(𝑡),𝑑𝐸(𝑡)𝑑𝑡=𝜆𝑆(𝑡)𝐼(𝑡)𝜀𝐸(𝑡)𝑑𝐸(𝑡),𝑑𝐼(𝑡)𝑑𝑡=𝜀𝐸(𝑡)𝑑𝐼(𝑡)𝛿𝐼(𝑡)𝜃𝐼(𝑡).(6.1)It is shown that the global dynamics is completely determined by the basic reproduction number 𝑅0. If 𝑅01, the disease-free equilibrium is globally stable. If 𝑅0>1, a unique endemic equilibrium is globally stable in the interior of the feasible region.

The reason of using 𝛽𝐸𝐼/𝑁 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, 𝐼/𝑁. Due to this difference in the form, our model shows many different and complex dynamics. In our model, however, the basic reproduction number 𝑅0 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 𝐸 to 𝐼 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 Figures 24, we can take 𝑅=𝐵2𝑛𝑚1𝑚1=𝑘𝑝𝑑+𝛿2𝑑𝑝aln1/1𝑝𝑜𝑝(𝑑+𝛿)aln1/1𝑝𝑜1(6.2) as a threshold parameter. If 𝑅<0, the virus will die out; if 𝑅>0, virus may die out or outbreak.

In the first case, to make 𝑅 below than zero, the most direct way is to change the parameter values which correspond to the related measures. Obviously, 𝑅 increases with 𝑘, 𝑎, and 𝑝𝑜. The relationships between 𝑅 and 𝑝, 𝛿 are shown in Figure 7; 𝑅 increases with 𝑝, 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 (𝑝𝑜) and appealing people to install the antivirus software and update in time to increase the ability of finding and removing the viruses (𝛿).

In the second case, from the definition of the domain of attraction, we know that starting from any initial points in the region Ω𝑐 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.

Acknowledgments

This work is supported by the National Natural Science Foundation of China (11171314, 70871072, 61171179, and 11147015), the Program for Basic Research of Shan'xi province (2010011007 and 2010011002-1), research project supported by the Shanxi Scholarship Council of China (2010-074), the Foundation of the National Key Laboratory for Electronic Measurement Technology (9104c12040051010).