Stability and Hopf Bifurcation in a Computer Virus Model with Multistate Antivirus

Abstract

By considering that people may immunize their computers with countermeasures in susceptible state, exposed state and using anti-virus software may take a period of time, a computer virus model with time delay based on an SEIR model is proposed. We regard time delay as bifurcating parameter to study the dynamical behaviors which include local asymptotical stability and local Hopf bifurcation. By analyzing the associated characteristic equation, Hopf bifurcation occurs when time delay passes through a sequence of critical value. The linerized model and stability of the bifurcating periodic solutions are also derived by applying the normal form theory and the center manifold theorem. Finally, an illustrative example is also given to support the theoretical results.

1. Introduction

As globalization and development of communication networks have made computers more and more present in our daily life, the threat of computer viruses also becomes an increasingly important issue of concern. In 2003, a virus, called worm king, rapidly spread and attacked the global world, which results the network of the internet to be seriously congested and server to be paralyzed [1]. In 2010, the report of pestilence about computer virus in China revealed that more than 90% computers in China are infected computer virus.

Computer viruses are small programs developed to damage the computer systems erasing data, stealing information. Their action throughout a network can be studied by using classical epidemiological models for disease propagation [2–6]. In [7–9], based on SIR classical epidemic model, Mark had proposed the dynamical models for the computer virus propagation, which provided estimations for temporal evolutions of infected nodes depending on network parameters [10–12]. In [13], Richard and Mark propose a modified propagation model named SEIR (susceptible-exposed-infected-recover) model to simulate virus propagation. In [14], on this basis of the SIR model, Yao et al. proposed a SIDQV model with time delay which add a quarantine state to clean the virus. However, both above models assume the viruses are cleaned in the infective state. In fact, in addition to clean viruses in state I, people may immunize their computers with countermeasures in state S and state E in the real world. Moreover there may be a time lag when the node uses antivirus software to clean the virus.

In this paper, in order to overcome the above-mentioned limitation, we present a new computer virus model with time delay which is depending on the SEIR model [15]; time delay can be considered the period of the node uses antivirus software to clean the virus. This model provides an opportunity for us to study the behaviors of virus propagation in the presence of antivirus countermeasures, which are very important and desirable for understanding of the virus spread patterns, as well as for management and control of the spread. The remainder of this paper is organized as follows. In Section 2, the stability of trivial solutions and the existence of Hopf bifurcation are discussed. In Section 3, a formula for determining the direction of Hopf bifurcation and the stability of bifurcating periodic solutions will be given by using the normal form and center manifold theorem introduced by Hassard et al. in [16]. In Section 4, numerical simulations aimed at justifying the theoretical analysis will be reported.

2. Mathematical Model Formulation

Our model is based on the traditional SEIR model [7–9, 15, 17]. The SEIR model has four states: susceptible, exposed (infected but not yet infectious), infectious, and recovered. Our assumptions on the dynamical model are as follows.(1)In the real world, in addition cleaning viruses in state I, people may immunize their computers with countermeasures in state S and state E (after virus being cleaned), which may result in new state transition paths in comparison with SIR model:S-R: using countermeasure of real-time immunization,E-R: using real-time immunization after virus codes cleaning.(2)In state S, when people install the antivirus software on their computer, we assume that their computer can be immunized at a unit time.(3)In state E, since the computer is infected by the virus, the antivirus software may use a period to search the document and clean the viruses.(4)Denote the period of time of killing viruses when users find that their computers are infected by viruses.(5)While the computer is installed the antivirus software, it will not be quarantine or replacement. On the basis of the above hypotheses (1)–(5), the dynamical model can be formulated by the following equations:𝑑𝑆(𝑡)(𝑡)𝑆(𝑡)−𝜌𝑑𝑡=𝑢𝑁−𝛽𝐼_𝑆𝑅𝑆(𝑡),(𝑡)+𝜇𝑑𝐸(𝑡)𝑆(𝑡)−𝐸(𝑡)𝑑𝑡=𝛽𝐼𝛼+𝜇−𝜌_𝐸𝑅𝐸(),(𝑡)𝑡−𝜏𝑑𝐼(𝑡)−𝐼(𝑡),(𝑡)𝑑𝑡=𝛼𝐸𝛾+𝜇𝑑𝑅𝑑𝑡=𝜌_𝑆𝑅𝑆(𝑡)+𝜌_𝐸𝑅𝐸(𝑡−𝜏)+𝛾𝐼(𝑡−𝜏)−𝜇𝑅(𝑡), (2.1) where 𝜌_𝑆𝑅 describes the impact of implementing real-time immunization, 𝜌_𝐸𝑅 describes the impact of cleaning the virus and immunizing the nodes, and 𝜇 describes the impact of quarantine or replacement. 𝛼 is the transition rate from 𝐸 to 𝐼, and 𝛾 is the recovery rate from 𝐼 to 𝑅. 𝜏 is the time delay that the node usees antivirus software to clean the virus. 𝛽 is the transition rate from 𝑆 to 𝐸.

3. Local Stability of the Equilibrium and Existence of Hopf Bifurcation

We may see that the first three equations in (2.1) are independent of the fourth equation, and therefore, the fourth equation can be omitted without loss of generality. Hence, model (2.1) can be rewritten as𝑑𝑆(𝑡)(𝑡)𝑆(𝑡)−𝜌𝑑𝑡=𝑢𝑁−𝛽𝐼_𝑆𝑅𝑆(𝑡),(𝑡)+𝜇𝑑𝐸(𝑡)𝑆(𝑡)−𝐸(𝑡)𝑑𝑡=𝛽𝐼𝛼+𝜇−𝜌_𝐸𝑅𝐸(),(𝑡)𝑡−𝜏𝑑𝐼(𝑡)−𝐼(𝑡).𝑑𝑡=𝛼𝐸𝛾+𝜇 (3.1) For the convenience of description, we define the basic reproduction number of the infection as𝑅₀=𝜇𝑁𝛽𝛼𝜌_𝑆𝑅+𝜇𝛼+𝜌_𝐸𝑅.+𝜇𝛾+𝜇 (3.2) Clearly, we have the following results with respect to the stable state of system (3.1). Here, the proof is omitted (see [17] for the details).

Theorem 3.1. If 𝑅₀<1, system (3.1) has only the disease-free equilibrium 𝐸₀=(𝜇𝑁/(𝜌_𝑆𝑅+𝜇),0,0) and is globally asymptotically stable. If 𝑅₀>1, 𝐸₀ becomes unstable and there exists a unique positive equilibrium 𝐸_𝑣𝑒, where 𝐸_𝑣𝑒=(𝜇𝑁/(𝜌_𝑆𝑅+𝜇)𝑅₀,𝜇𝑁(𝑅₀−1)/𝑅₀(𝛼+𝜇+𝜌_𝐸𝑅),𝛼𝐸^∗/(𝛾+𝜇)). Furthermore, for any 𝜏>0, 𝐸₀ is asymptotically stable if 𝑅₀<1 and unstable if 𝑅₀>1.

To investigate the qualitative properties of the positive equilibrium 𝐸^∗ with 𝜏>0, it is necessary to make the following assumption:

(H1) 𝑅₀>1.

Under hypothesis (H1), the Jacobian matrix of the system (3.1) about 𝐸_𝑣𝑒 is given by𝐽𝐸_𝑣𝑒=−𝑎₁0−𝑎₂𝑎₃−𝑎₄−𝑎₇𝑒^−𝜆𝜏𝑎₂0𝑎₅−𝑎₆, (3.3) where 𝑎₁=𝛽𝐼+𝜌_𝑆𝑅+𝜇,𝑎₂=𝛽𝑆,𝑎₃=𝛽𝐼,𝑎₄=𝛼+𝜇,𝑎₅=𝛼,𝑎₆=𝛾+𝜇,𝑎₇=𝜌_𝐸𝑅.

We can obtain the following characteristic equation:𝜆³+𝑏₁𝜆²+𝑏₂𝜆+𝑏₃+𝑒^−𝜆𝜏𝑏₄𝜆²+𝑏₅𝜆+𝑏₆=0, (3.4) where𝑏₁=𝑎₁+𝑎₄+𝑎₆,𝑏₂=𝑎₁𝑎₆+𝑎₄𝑎₁+𝑎₆−𝑎₂𝑎₅,𝑏₃=𝑎₁𝑎₄𝑎₆−𝑎₁𝑎₂𝑎₅+𝑎₂𝑎₃𝑎₅,𝑏₄=𝑎₇,𝑏₅=𝑎₇𝑎₁+𝑎₆,𝑏₆=𝑎₁𝑎₆𝑎₇. (3.5) If 𝑖𝜔(𝜔>0) is a root of (3.4), then−𝑖𝜔³−𝑏₁𝜔²+𝑏₂𝑖𝜔+𝑏₃+𝑒^{−𝑖𝜔𝜏}−𝜔²𝑏₄+𝑏₅𝑖𝜔+𝑏₆=0. (3.6) Separating the real and imaginary parts of (3.6), we have𝑏₅𝜔sin𝜔𝜏+𝑏₆−𝑏₄𝜔²cos𝜔𝜏=𝑏₁𝜔²−𝑏₃,𝑏₅𝜔cos𝜔𝜏−𝑏₆−𝑏₄𝜔²sin𝜔𝜏=𝜔³−𝑏₂𝜔. (3.7) Adding up the squares of (3.7) yields 𝜔⁶+𝑏²₁−2𝑏₂−𝑏²₄𝜔⁴+𝑏²₂−2𝑏₁𝑏₃+2𝑏₄𝑏₆−𝑏²₅𝜔²+𝑏²₃−𝑏²₆=0. (3.8) Letting 𝑧=𝜔²,𝑐₁=𝑏²₁−2𝑏₂−𝑏²₄,𝑐₂=𝑏²₂−2𝑏₁𝑏₃+2𝑏₄𝑏₆−𝑏²₅,𝑐₃=𝑏²₃−𝑏²₆, then (3.8) becomes𝑧³+𝑐₁𝑧²+𝑐₂𝑧+𝑐₃=0. (3.9) Letting 𝑧^∗=(1/3)(−𝑐₁+𝑐²₁−3𝑐₂),ℎ(𝑧^∗)=(𝑧^∗)³+𝑐₁(𝑧^∗)²+𝑐₂𝑧^∗+𝑐₃, then we have the following results (see [18–22] for details) about the distributions of the positive roots of (3.9).

Lemma 3.2 (see [18–22]). (i) If 𝑐₃<0, then (3.9) has at least one positive root.
(ii) If 𝑐₃≥0 and 𝑐²₁−3𝑐₂≤0, then (3.9) has no positive root.
(iii) If 𝑐₃≥0 and 𝑐²₁−3𝑐₂>0, then (3.9) has positive roots if and only if 𝑧^∗>0 and ℎ(𝑧^∗)≤0.

Suppose (3.9) has positive roots; without loss of generality, we assume that it has three positive roots defined by 𝜔_𝑘=√𝑧_𝑘, 𝑘=1,2,3. By (3.7), we havecos𝜔_𝑘𝜏=𝑏₁𝜔²_𝑘−𝑏₃𝑏₆−𝑏₄𝜔²_𝑘+𝑏₅𝜔²_𝑘𝜔²_𝑘−𝑏₂𝑏²₅𝜔²_𝑘+𝑏₆−𝑏₄𝜔²_𝑘². (3.10) Thus, denoting𝜏^𝑗_𝑘=1𝜔_𝑘𝑏arccos¹𝜔²_𝑘−𝑏₃𝑏₆−𝑏₄𝜔²_𝑘+𝑏₅𝜔²_𝑘𝜔²_𝑘−𝑏₂𝑏²₅𝜔²_𝑘+𝑏₆−𝑏₄𝜔²_𝑘²+2𝑗𝜋𝜔_𝑘, (3.11) where 𝑘=1,2,3; 𝑗=0,1,…, then ±𝑖𝜔 is a pair of purely imaginary roots of (3.4) with 𝜏^𝑗_𝑘. Define 𝜏₀=𝜏⁰_𝑘0=min_𝑘=1,2,3𝜏⁰_𝑘0,𝜔₀=𝜔_𝑘0. (3.12) Note that when 𝜏=0, (3.4) becomes𝜆³+𝑏₁+𝑏₄𝜆²+𝑏₂+𝑏₅𝜆+𝑏₃+𝑏₆=0. (3.13) In addition, Routh-Hurwitz criterion [13] implies that, if the following condition holds, then all roots of (3.13) have negative real parts.

(H2) (𝑏₁+𝑏₄)>0,(𝑏₁+𝑏₄)(𝑏₂+𝑏₅)−(𝑏₃+𝑏₆)>0.

Till now, we can employ a result from Ruan and Wei [23] to analyze (3.4), which is, for the convenience of the reader, stated as follows.

Lemma 3.3 (see [23]). Consider the exponential polynomial 𝑃𝜆,𝑒^−𝜆𝜏,…,𝑒^{−𝜆𝜏𝑚}=𝜆^𝑛+𝑝₁⁽⁰⁾𝜆^𝑛−1+⋯+𝑝⁽⁰⁾_𝑛−1𝜆+𝑝_𝑛⁽⁰⁾+𝑝₁⁽¹⁾𝜆^𝑛−1+⋯+𝑝⁽¹⁾_𝑛−1𝜆+𝑝_𝑛⁽¹⁾𝑒^−𝜆𝜏+⋯+𝑝₁^(𝑚)𝜆^𝑛−1+⋯+𝑝^(𝑚)_𝑛−1𝜆+𝑝_𝑛^(𝑚)𝑒^{−𝜆𝜏𝑚}, (3.14) where 𝜏_𝑖≥0 (𝑖=1,2,…,𝑚) and 𝑝_𝑗^(𝑖)(𝑗=1,2,…,𝑚) are constants. As(𝜏₁,𝜏₂,…,𝜏_𝑚) vary, the sum of the order of the zeros of 𝑃(𝜆,𝑒^−𝜆𝜏,…,𝑒^{−𝜆𝜏𝑚}) on the open right half plane can change only if a zero appears on or crosses the imaginary axis.

Using Lemmas 3.2 and 3.3 we can easily obtain the following results on the distribution of roots of the transcendental (3.4).

Lemma 3.4. (2.1) If 𝑐₃>0 and 𝑐²₁−3𝑐₂≤0, then all roots with positive real parts of (3.4) have the same sum as those of the polynomial (3.13) for all 𝜏≥0.
(3.1) If either 𝑐₃<0 or 𝑐₃≥0 and 𝑐²₁−3𝑐₂>0, 𝑧^∗>0,ℎ(𝑧^∗)≤0, then all roots with positive real parts of (3.4) have the same sum as those of the polynomial (3.13) for 𝜏∈[0,𝜏₀).

Lemma 3.5. If 3𝑤⁴_𝑘+𝑐₁𝑤¹_𝑘+𝑐₂≠0, then the following transversality condition holds: sgnRe𝑑𝜆𝑑𝜏⁻¹≠0when𝜏=𝜏₀. (3.15)

Proof. Differentiating (3.4) with respect to 𝜏 yields 3𝜆²+2𝑏₁𝜆+𝑏₂+2𝑏₄𝜆+𝑏₅−𝜏𝑏₄𝜆²+𝑏₅𝜆+𝑏₆𝑒^−𝜆𝜏𝑑𝜆𝑏𝑑𝜏=𝜆₄𝜆²+𝑏₅𝜆+𝑏₆𝑒^−𝜆𝜏. (3.16) For the sake of simplicity, denoting 𝜔₀ and 𝜏₀ by 𝜔,𝜏 respectively, then 𝑑𝜆𝑑𝜏⁻¹=3𝜆²+2𝑏₁𝜆+𝑏₂𝜆𝑏₄𝜆²+𝑏₅𝜆+𝑏₆𝑒^−𝜆𝜏+2𝑏₄𝜆+𝑏₅𝜆𝑏₄𝜆²+𝑏₅𝜆+𝑏₆−𝜏𝜆=−2𝜆³+𝑏₁𝜆²−𝑏₃𝜆²𝜆³+𝑏₁𝜆²+𝑏₂𝜆+𝑏₃+𝑏₄𝜆²−𝑏₆𝜆𝑏₄𝜆²+𝑏₅𝜆+𝑏₆−𝜏𝜆=−2𝑖𝜔³−𝑏₁𝜔²−𝑏₃𝜔²𝑏₃−𝑏₁𝜔²−𝑖𝜔³−𝑏₂𝜔+𝑏₄𝜔²+𝑏₆𝜔²𝑏₆−𝑏₃𝜔²+𝑏₄−𝜏𝑖𝜔.𝑖𝜔 (3.17) Then we get Re𝑑𝜆𝑑𝜏⁻¹1=−𝜔²𝑏³₃−2𝜔⁶−𝑏²₁−2𝑏₂𝜔⁴𝑏₃−𝑏₁𝜔²²+𝜔³−𝑏₂𝜔²+𝑏²₆−𝑏²₄𝜔⁴𝑏₆−𝑏₄𝜔²²+𝑏²₅𝜔²=2𝜔⁶+𝑐₁𝜔⁴−𝑐₃𝜔²𝑏₆−𝑏₄𝜔²²+𝑏²₅𝜔²=3𝜔⁴+𝑐₁𝜔²+𝑐₂𝑏₆−𝑏₄𝜔²²+𝑏²₅𝜔². (3.18) Then, if 3𝜔⁴+𝑐₁𝜔²+𝑐₂≠0, we have sgn{Re{(𝑑𝜆/𝑑𝜏)⁻¹}}≠0, we complete proof.

Thus from Lemmas 3.2, 3.3, 3.4, and 3.5, and we have the following.

Theorem 3.6. Suppose that (H1) and (H2) hold, then the following results hold.(1)The positive equilibrium of (3.1) is asymptotically stable, if 𝑐₃>0 and 𝑐²₁−3𝑐₂≤0;(2)if either 𝑐₃<0 or 𝑐₃≥0 and 𝑐²₁−3𝑐₂>0, 𝑧^∗>0,ℎ(𝑧^∗)≤0, system (3.1) is asymptotically stable for 𝜏∈[0,𝜏₀) and system (3.1) undergoes a Hopf bifurcation at the origin when 𝜏=𝜏₀.

4. Direction of the Hopf Bifurcation

In this section, we derive explicit formulae for computing the direction of the Hopf bifurcation and the stability of bifurcation periodic solution at critical values 𝜏₀ by using the normal form theory and center manifold reduction.

Letting 𝑥₁=𝑆−𝑆^∗,𝑥₂=𝐸−𝐸^∗,𝑥₃=𝐼−𝐼^∗,𝑥_𝑖(𝑡)=𝑥_𝑖(𝜏𝑡),𝜏=𝜏₀+𝜇, and dropping the bars for simplification of notation, system (3.1) is transformed into an FDE as𝑥̇(𝑡)=𝐿_𝜇𝑥_𝑡+𝑓𝜇,𝑥_𝑡, (4.1) with𝐿_𝜇𝜑=𝜏₀+𝜇𝐵₁𝜑(0)+𝐵₂𝜑(−1), (4.2) where𝐵₁=−𝑎₁0−𝑎₂𝑎₃−𝑎₄𝑎₂0𝑎₅−𝑎₆,𝐵₂=0000−𝑎₇𝑒^−𝜆𝜏0,𝑓=𝜏000𝜇,𝜑₀+𝜇−𝛽𝜑₁(0)𝜑₂(0)𝛽𝜑₁(0)𝜑₂(0)0. (4.3) Using the Riesz representation theorem, there exists a function 𝜂(𝜃,𝜇) of bounded variation for 𝜃∈[−1,0], such that𝐿_𝜇𝜑=⁰₋₁𝑑𝜂𝜃,𝜇𝜑(𝜃)𝜑∈𝐶. (4.4) In fact, we can choose𝜂𝜃,𝜇=𝜏₀+𝜇𝐵₁𝛿(𝜃)+𝐵₂𝛿(𝜃+1), (4.5) where 𝛿(𝜃) is Dirac delta function.

In the next, for 𝜑∈[−1,0], we define𝐴𝜇𝜑=𝑑𝜑[),𝑑𝜃,𝜃∈−1,0₀₋₁𝜑(𝜃)𝑅𝜇[],𝑓𝑑𝜂𝜃,𝜇,𝜃=0,𝜑=0,𝜃∈−1,0𝜇,𝜑,𝜃=0. (4.6)(4.7) Then system (4.2) can be rewritten as𝑥̇(𝑡)=𝐴𝜇𝑥_𝑡+𝑅𝜇𝑥_𝑡, (4.8) where 𝑥_𝑡(𝜃)=𝑥(𝑡+𝜃).

The adjoint operator 𝐴^∗ of 𝐴 is defined by𝐴^∗𝜇−(𝑠)𝜓=𝑑𝜓(],𝑑𝜃,𝑠∈0,1₀₋₁𝑑𝜂^𝑇(𝑡,0)𝜓(−𝑡),𝑠=0, (4.9) where 𝜂^𝑇 is the transpose of the matrix 𝜂.

For 𝜑∈𝐶¹[−1,0] and 𝜓∈𝐶¹[0,1], we define𝜓,𝜑=𝜓(0)⋅𝜑(0)−⁰_𝜃=−1^𝜃_𝜉=0𝜓(𝜉−𝜃)𝑑𝜂(𝜃)𝜑(𝜉)𝑑𝜉, (4.10) where 𝜂(𝜃)=𝜂(𝜃,0). We know that ±𝑖𝜏₀𝜔₀ is an eigenvalue of 𝐴(0), so ±𝑖𝜏₀𝜔₀ is also an eigenvalue of 𝐴^∗(0). We can get𝑞(𝜃)=1𝑞₁𝑞₂𝑒^{𝑖𝜏0𝜔0𝜃},−1<𝜃≤0. (4.11) From the above discussion, it is easy to know that𝐴𝑞(0)=𝑖𝜏₀𝜔₀𝑞(0). (4.12) Hence we obtain𝑞₁=𝑖𝜔₀𝑞₂𝑎₅,𝑞₂=−𝑖𝜔₀+𝑎₁𝑎₂. (4.13) Suppose that the eigenvector 𝑞^∗ of 𝐴^∗ is𝑞^∗(𝑠)=1𝑞₁^∗𝑞₂^∗𝑒^{𝑖𝜏0𝜔0𝑠}, (4.14) Then the following relationship is obtained:𝐴^∗𝑞(0)=−𝑖𝜏₀𝜔₀𝑞^∗(0). (4.15) Hence we obtain𝑞^∗₁=𝑎₁−𝑖𝜔₀𝑎₃,𝑞^∗₂=𝑎₄+𝑎₇𝑒^{𝑖𝜔0𝜏0}𝑎₅𝑞^∗₁. (4.16) Let𝑞^∗,𝑞=1. (4.17) One can obtain𝑞^∗,𝑞=𝑞^∗(0)⋅𝑞(0)−⁰_𝜃=−1^𝜃_𝜉=0𝑞^∗𝑇()(𝜃)𝜑(𝜉)=1𝜉−𝜃𝑑𝜂𝑑𝜉𝜌1+𝑞₁𝑞₁^∗+𝑞₂𝑞₂^∗−⁰_𝜃=−1^𝜃_𝜉=0𝜏₀1𝜌1𝑞₁^∗𝑞₂^∗−𝑎₁0−𝑎₂𝑎₃−𝑎₄𝑎₂0𝑎₅−𝑎₆𝛿(𝜃)+0000−𝑎₇0𝛿()1𝑞000𝜃+1₁𝑞₂𝑒^{𝑖𝜏0𝜔0𝜃}=1𝑑𝜉𝑑𝜃𝜌1+𝑞₁𝑞₁^∗+𝑞₂𝑞₂^∗−1𝜌𝜏₀𝑒^{−𝑖𝜔0𝜏0}𝑎₇𝑞₁𝑞₁^∗. (4.18)

Hence we obtain𝜌=1+𝑞₁𝑞₁^∗+𝑞₂𝑞₂^∗−𝜏₀𝑒^{−𝑖𝜔0𝜏0}𝑎₇𝑞₁𝑞₁^∗. (4.19)

In the remainder of this section, by using the same notations as in Hassard et al. [16], we first compute the coordinates for describing the center manifold 𝐶₀ at 𝜇=0. Leting 𝑥_𝑡 be the solution of (4.1) with 𝜇=0, we define𝑧(𝑡)=𝑞^∗,𝑥_𝑡,𝑊(𝑡,𝜃)=𝑥_𝑡−2Re𝑧(𝑡)𝑞(𝜃). (4.20) On the center manifold 𝐶₀ we have𝑊(𝑡,𝜃)=𝑊(𝑧,𝑧,𝑡), (4.21) where𝑊(𝑧,𝑧,𝑡)=𝑊₂₀(𝜃)𝑧²2+𝑊₁₁(𝜃)𝑧𝑧2+𝑊₀₂(𝜃)𝑧²2+⋯. (4.22)

In fact, 𝑧 and 𝑧 are local coordinate for 𝐶₀ in the direction of 𝑞 and 𝑞^∗. Note that, if 𝑥_𝑡 is, we will deal with real solutions only. Since 𝜇=0𝑧̇(𝑡)=𝑞^∗,𝑥̇_𝑡=𝑞^∗,𝐴𝜇𝑥_𝑡+𝑅𝜇𝑥_𝑡=𝑞^∗,𝐴𝑥_𝑡+𝑞^∗,𝑅𝑥_𝑡=𝑖𝜏₀𝑤₀𝑧+𝑞^∗(0)⋅𝑓0,𝑊(𝑡,0)+2Re𝑧(𝑡)𝑞(0). (4.23) Rewrite (4.23) as𝑧̇(𝑡)=𝑖𝜏₀𝜔₀𝑧+𝑔(𝑧,𝑧), (4.24) where𝑔(𝑧,𝑧)=𝑔₂₀𝑧2²+𝑔₁₁𝑧𝑧+𝑔₀₂𝑧²2+𝑔₂₁𝑧²𝑧2+⋯. (4.25) From (4.1) and (4.24), we have𝑊̇=𝑥̇_𝑡−𝑧𝑞̇−𝑧̇𝑞̇=𝐴𝑊−2Re𝑞^∗(0)⋅𝑓(𝑧,𝑧)𝑞(𝜃),𝜃∈[−2𝜏,0),𝐴𝑊−2Re𝑞^∗(0)⋅𝑓(𝑧,𝑧)𝑞(𝜃)+𝑓₀(𝑧,𝑧),𝜃=0. (4.26) Let𝑊̇=𝐴𝑊+𝐻(𝑧,𝑧,𝜃), (4.27) where𝐻(𝑧,𝑧,𝜃)=𝐻₂₀(𝜃)𝑧²2+𝐻₁₁(𝜃)𝑧𝑧2+𝐻₀₂(𝜃)𝑧²2+⋯. (4.28) Expanding the above series and comparing the corresponding coefficients, we obtain𝐴−2𝑖𝑤₀𝑊₂₀(𝜃)=−𝐻₂₀(𝜃),𝐴𝑊₁₁(𝜃)=−𝐻₁₁(𝜃),𝐴+2𝑖𝑤₀𝑊₀₂(𝜃)=−𝐻₀₂(𝜃). (4.29) Since 𝑥_𝑡=𝑥(𝑡+𝜃)=𝑊(𝑧,𝑧,𝜃)+𝑧𝑞+𝑧⋅𝑞, we have𝑥_𝑡=𝑊⁽¹⁾(𝑧,)𝑊𝑧,𝜃₍₂₎(𝑧,)𝑊𝑧,𝜃₍₃₎(𝑧,)1𝑞𝑧,𝜃+𝑧₁𝑞₂𝑒^{𝑖𝜔0𝜃}+𝑧1𝑞₁𝑞₂𝑒^{−𝑖𝜔0𝜃}. (4.30) Thus, we can obtain𝜑₁(0)=𝑧+𝑧+𝑊⁽¹⁾₂₀(0)𝑧²2+𝑊⁽¹⁾₁₁(0)𝑧𝑧+𝑊⁽¹⁾₀₂(0)𝑧²2,𝜑₂(0)=𝑧𝑞₁+𝑧𝑞^∗₁+𝑊⁽¹⁾₂₀(0)𝑧²2+𝑊⁽¹⁾₁₁(0)𝑧𝑧+𝑊⁽¹⁾₀₂(0)𝑧²2. (4.31) So𝜑₁(0)𝜑₂(0)=𝑞₁𝑧²+𝑞₁𝑧²+𝑞₁+𝑞₂𝑧𝑧+𝑊⁽²⁾₁₁+12𝑊⁽²⁾₂₀+𝑊⁽¹⁾₁₁𝑞₁+12𝑊⁽¹⁾₂₀𝑞𝑧²𝑧. (4.32) It follows from (4.24) and (4.25) that𝑓=𝐾𝜑,𝜇¹¹𝑧²+𝐾₁₂𝑧𝑧+𝐾₁₃𝑧²+𝐾₁₄𝑧²𝑧𝐾₂₁𝑧²+𝐾₂₂𝑧𝑧+𝐾₂₃𝑧²+𝐾₂₄𝑧²𝑧0, (4.33) where𝐾₁₁=−𝛽𝑞₁,𝐾₁₂=−𝛽𝑞₁,𝐾₁₃=−𝛽𝑞₁+𝑞₁,𝐾₁₄=−𝛽𝑊⁽²⁾₁₁+12𝑊⁽²⁾₂₀+𝑊⁽¹⁾₁₁𝑞₁+12𝑊⁽¹⁾₂₀𝑞,𝐾₂₁=𝛽𝑞₁,𝐾₂₂=𝛽𝑞₁,𝐾₂₃=𝛽𝑞₁+𝑞₁,𝐾₂₄=𝛽𝑊⁽²⁾₁₁+12𝑊⁽²⁾₂₀+𝑊⁽¹⁾₁₁𝑞₁+12𝑊⁽¹⁾₂₀𝑞. (4.34) Since 𝑞^∗(0)=(1/𝜌)(1,𝑞^∗₁,𝑞^∗₂)^𝑇, we have𝑔(𝑧,𝑧)=1𝜌1,𝑞^∗₁,𝑞^∗₂𝐾₁₁𝑧²+𝐾₁₂𝑧𝑧+𝐾₁₃𝑧²+𝐾₁₄𝑧²𝑧𝐾₂₁𝑧²+𝐾₂₂𝑧𝑧+𝐾₂₃𝑧²+𝐾₂₄𝑧²𝑧0. (4.35) Comparing the coefficients of the above equation with those in (4.27), we have𝑔₂₀=1𝜌𝐾₁₁+𝐾₂₁𝑞^∗₁,𝑔₁₁=1𝜌𝐾₁₂+𝐾₂₂𝑞^∗₁,𝑔₀₂=1𝜌𝐾₁₃+𝐾₂₃𝑞^∗₁,𝑔₂₁=1𝜌𝐾₁₄+𝐾₂₄𝑞^∗₁. (4.36) In what follows, we focus on the computation of 𝑊₂₀(𝜃) and 𝑊₁₁(𝜃). For the expression of 𝑔₂₁, we have𝐻(𝑧,𝑧,𝜃)=−2Re𝑞^∗(0)⋅𝑓(𝑧,𝑧)𝑞(𝜃𝑔)=−₂₀𝑧2²+𝑔₁₁𝑧𝑧+𝑔₀₂𝑧²2+⋯𝑞(𝜃)−𝑔₂₀𝑧2²+𝑔₁₁𝑧𝑧+𝑔₀₂𝑧²2+⋯𝑞(𝜃). (4.37) Comparing the coefficients of the above equation, we can obtain that𝐻₂₀(𝜃)=−𝑔₂₀𝑞(𝜃)−𝑔₀₂𝑞(𝜃)[),𝐻,𝜃∈−1,0₁₁(𝜃)=−𝑔₁₁𝑞(𝜃)−𝑔₁₁𝑞(𝜃),𝜃∈[−1,0). (4.38)(4.39) Substituting (4.39) into (4.27) and (4.38) into (4.27), respectively, we get𝑊̇₂₀(𝜃)=2𝑖𝜏₀𝜔₀𝑊₂₀(𝜃)+𝑔₂₀𝑞(𝜃)+𝑔₂₀𝑞(𝜃),𝑊̇₁₁(𝜃)=+𝑔₁₁𝑞(𝜃)+𝑔₁₁𝑞(𝜃). (4.40) So𝑊₂₀(𝜃)=𝑖𝑔₂₀𝜏₀𝜔₀𝑞(0)𝑒^{𝑖𝜏0𝜔0𝜃}−𝑔₀₂3𝑖𝜏₀𝜔₀𝑞(0)𝑒^{−𝑖𝜏0𝜔0𝜃}+𝐸₁𝑒^{2𝑖𝜏0𝜔0𝜃},𝑊₁₁(𝜃)=𝑔₁₁𝑖𝜏₀𝜔₀𝑞(0)𝑒^{𝑖𝜏0𝜔0𝜃}−𝑔₁₁𝑖𝜏₀𝜔₀𝑞(0)𝑒^{−𝑖𝜏0𝜔0𝜃}+𝐸₂. (4.41) In the sequel, we will determine 𝐸₁ and 𝐸₂. Form the definition of 𝐴 in (4.8), we have𝑑𝜂(⁰₋₁𝜃)𝑊₂₀(𝜃)=2𝑖𝜏₀𝜔₀𝑊₂₀(0)−𝐻₂₀(0),₀₋₁𝑑𝜂(𝜃)𝑊₁₁(𝜃)=−𝐻₁₁(0). (4.42)(4.43) From (4.6) and (4.38)-(4.39), we have𝐻₂₀(𝜃)=−𝑔₂₀𝑞(𝜃)−𝑔₀₂𝑞(𝜃)+𝐾₁₁,𝐾₂₁,0^𝑇,𝐻₁₁(𝜃)=−𝑔₁₁𝑞(𝜃)−𝑔₁₁𝑞(𝜃)+𝐾₁₂,𝐾₂₂,0^𝑇. (4.44)(4.45) Substituting (4.41) and (4.44) into (4.42) and noticing that𝑖𝜔₀𝐼−⁰₋₁𝑒^{𝑖𝜔0𝜃}𝑑𝜂(𝜃)𝑞(0)=0,−𝑖𝜔₀𝐼−⁰₋₁𝑒^{−𝑖𝜔0𝜃}𝑑𝜂(𝜃)𝑞(0)=0, (4.46) we can obtain2𝑖𝜔₀𝐼−⁰₋₁𝑒^{2𝑖𝜏0𝜔0𝜃}𝑑𝜂(𝜃)𝐸₁=𝐾₁₁𝐾₂₁0^𝑇, (4.47) which leads to2𝑖𝜔₀+𝑎₁0𝑎₂−𝑎₃2𝑖𝜔₀+𝑎₄+𝑎₇𝑒^{−2𝑖𝜔0𝜏0}−𝑎₂0−𝑎₅2𝑖𝜔₀+𝑎₆𝐸₁⁽¹⁾𝐸₁⁽²⁾𝐸₁⁽³⁾=𝐾₁₁𝐾₂₁0,𝑎₁0𝑎₂−𝑎₃𝑎₄+𝑎₇−𝑎₂0−𝑎₅𝑎₆𝐸₂⁽¹⁾𝐸₂⁽²⁾𝐸₂⁽³⁾=𝐾₁₂𝐾₂₂0. (4.48) It follows that𝐸₁⁽¹⁾=𝐾₁₁−𝑎₂𝐸₁⁽³⁾2𝑖𝜔₀+𝑎₁,𝐸₁⁽²⁾=2𝑖𝜔₀+𝑎₆𝑎₅𝐸₁⁽³⁾,𝐸₁⁽³⁾=𝐾₂₂−𝐾₁₁/2𝑖𝜔₀+𝑎₁𝑎₂𝑎₃/2𝑖𝜔₀+𝑎₁+2𝑖𝜔₀+𝑎₄+𝑎₇𝑒^{−2𝑖𝜔0𝜏0}2𝑖𝜔₀+𝑎₆/𝑎₅−2𝑖𝜔₀−𝑎₆,𝐸₂⁽¹⁾=𝐾₁₂−𝑎₂𝐸₂⁽³⁾𝑎₁,𝐸₂⁽²⁾=𝑎₆𝑎₅𝐸₂⁽³⁾,𝐸₂⁽³⁾=𝑎₁𝑎₅𝐾₂₂−𝑎₅𝐾₁₂𝑎₂𝑎₃𝑎₅+𝑎₄+𝑎₇𝑎₁𝑎₆−𝑎₁𝑎₅𝑎₆. (4.49) Based on the above analysis, we can see each 𝑔_𝑖𝑗 in (4.37) is determined by parameters and delays in (3.1). Thus, we can compute the following quantities:𝜇₂=−Re𝐶₁(0)Re𝜆′𝜏₀,𝑇₂=−Im𝐶₁(0)+𝜇₂Im𝜆^′(0)𝜔₀,𝛽₂=−2Re𝐶₁(0). (4.50)

Theorem 4.1. In (4.50), the following results hold.(1)The sign of 𝜇₂ determines the directions of the Hopf bifurcation: if 𝜇₂>0(𝜇₂<0) then the Hopf bifurcation is forward (backward) and the bifurcating periodic solutions exist for 𝜏>𝜏₀(𝜏<𝜏₀).(2)The sign of 𝛽₂ determines the stability of the bifurcating periodic solutions: the bifurcating periodic solutions are stable (unstable) if 𝛽₂<0(𝛽₂>0).(3)The sign of 𝑇₂ determines the period of the bifurcating periodic solutions: the period increases (decreases) if 𝑇₂>0(𝑇₂<0).

5. Numerical Examples

In this section, some numerical results of system (3.1) are presented to justify the Previous theorem above. As an example, considering the following parameters: 𝜇=0.01,𝑁=10000,𝛾=0.08,𝛼=0.1,𝛽=0.01,𝜌_𝑆𝑅=0.2,𝜌_𝐸𝑅=0.2, then 𝑅₀=1.706,𝑐₃=−3.6284𝑒−5, and 𝐸^∗=(279,133.6,148.4). According to the Lemma 3.2, (3.9) has one positive real root 𝜔=0.1194. Correspondingly, by (3.13), we obtain 𝜏₀=14.05. First, we choose 𝜏=13<𝜏₀, the corresponding wave form and phase plots are shown in Figure 1; it is easy to see from Figure 1 that system (3.1) is asymptotically stable. Finally, we choose 𝜏=14.15>𝜏₀ the corresponding wave form and phase plots are shown in Figure 2; it is easy to see that Figure 2 undergoes a Hopf bifurcation.

(a)

(b)

(c)

(d)

Figure 1 

𝜏=13<𝜏0. The positive equilibrium 𝐸∗ of system (3.1) is asymptotically stable.

(a)

(b)

(c)

(d)

Figure 2 

𝜏=14.15>𝜏0. The bifurcation periodic solution for system (3.1) is stable.

6. Conclusions

In this paper, considering that in addition to cleaning viruses in state I, people may immunize their computers with countermeasures in state 𝑆 and state 𝐸, and since using antivirus software will take a period of time, we have constructed a computer virus model with time delay depending on the SEIR model. The theoretical analyses for the computer virus models are given. Furthermore, we have proved that when time cross through the critical value, the system exist a Hopf bifurcation. Finally, simulation clarifies our results.

Acknowledgments

This work was supported in part by the National Natural Science Foundation of China under Grant 60973114 and Grant 61170249, in part by the Research Fund of Preferential Development Domain for the Doctoral Program of Ministry of Education of China under Grant 20110191130005, in part by the Natural Science Foundation project of CQCSTC under Grant 2009BA2024, in part by Changjiang Scholars, and in part by the State Key Laboratory of Power Transmission Equipment & System Security and New Technology, Chongqing University, under Grant 2007DA10512711206.