# Stability and Bifurcation Analysis of a Modified Epidemic Model for Computer Viruses

Authors: Chuandong Li; Wenfeng Hu; Tingwen Huang
Journal: Mathematical Problems in Engineering
 (2014)

---

## Abstract

We extend the three-dimensional SIR model to four-dimensional case and then analyze its dynamical behavior including stability and bifurcation. It is shown that the new model makes a significant improvement to the epidemic model for computer viruses, which is more reasonable than the most existing SIR models. Furthermore, we investigate the stability of the possible equilibrium point and the existence of the Hopf bifurcation with respect to the delay. By analyzing the associated characteristic equation, it is found that Hopf bifurcation occurs when the delay passes through a sequence of critical values. An analytical condition for determining the direction, stability, and other properties of bifurcating periodic solutions is obtained by using the normal form theory and center manifold argument. The obtained results may provide a theoretical foundation to understand the spread of computer viruses and then to minimize virus risks.


---

## Body

## 1. Introduction

Computer viruses arose in the 1980s with the widespread use of Internet in a variety of fields, such as communication, Internet business, and commercial system [1, 2]. With the development of hardware and software technology, computer viruses started to be a major threat to the information security. Usually, computer viruses can cause severe damage to the individuals and the corporations by different ways, including acquiring confidential data from network users, attacking the whole system, and even causing fatal damage to the hardware [3]. So the behaviors of computer viruses have attracted much attention from the fundamental researchers to the network security professional.It is well-known that computer virus is a malicious mobile code including virus, worm, Trojan horses, and logic bomb [4, 5]. Though different computer viruses vary in many respects, they all have many similar characteristics including infectivity, invisibility, latent, destructibility, and unpredictability [6]. The word “latent” means that the viruses hide themselves in the computers and spread them in the Internet through a period of time. Thus, in the construction of the system, the latent period should not be ignored [7–9].Because of a lot of similarities between the computer viruses and the infectious diseases, many researchers choose the epidemic models to find out the rule in the computer viruses [10, 11] and much attention is now paid to the effect of topological structure of the network on the spread of the viruses. Among the epidemic models, the SIR, SEI, and SEIR epidemic models are some of the most famous ones. Inspired by these epidemic models, the models of the computer virus have been proposed in recent years [12, 13].However, some shortcomings of such models arose due to the inevitable difference between computer viruses and infectious diseases. Consequently, the results obtained from the infectious diseases models cannot be carried over to computer viruses completely. As a result, we need to make some modifications in order to model the computer viruses.From the discussion above, we will propose a modified epidemic model for computer viruses and make some dynamic analysis on it. Specifically, we will extend the three-dimensional SIR model to four-dimensional SIRA model. However, because of the increased dimension, the complexity of the proposed model increases highly. We will present several theoretical results for its stability property and bifurcation dynamics by the rigorous mathematical analysis.


## 2. Model Description and Preliminaries

The present model is a modification of the original compartmental model [14]. Here, we assume that each node is denoted as one computer and the total population T can be divided into the following four groups by the state of each node:(1)
S
(
t
) is the number of noninfected computers subjected to possible infection;
(2)
I
(
t
) is the number of infected computers;
(3)
R
(
t
) is the number of removed ones due to infection or not;
(4)
A
(
t
) is the antidotal population representing computers equipped with fully effective antivirus programs [15].In the present paper, we use the antivirus distribution strategy; namely, we convert the susceptible into antidotal, which is proportional to the productS
I and with a controlled parameter a
S
A. In addition, the infected computers can be fixed by using antivirus programs, by which the infected computers can be converted into antidotal ones with a rate proportional to A
I and a proportion factor given by a
I
A, or we let the infected ones become useless and be removed with a rate controlled by δ because of the antivirus cost. Usually, the removed computers can be restored and converted into susceptible with a proportional factor σ. It is noticed that all the compartments have the mortality rate not due to the viruses. Here we assume all of them are the same and are denoted by the proportion coefficient μ. We further suppose that the influx rate N represents the incorporation of new computers to the network.It should be pointed out that the rate of the conversion from susceptible into infected ones is called incidence rate. It has been suggested by several authors that the viruses’ transmission process may have a nonlinear incidence rate. This allows one to include behavioral changes and prevent unbounded contact rates (e.g., [16]). In many epidemic models, the bilinear incidence rate β
S
I and the standard incidence rate β
S
I
/
N are frequently used. The bilinear incidence rate is based on the law of mass action. This contact law is more appropriate for communicable diseases such as influenza but not for point-to-point computer viruses. In the paper, we introduce the following saturated incidence rate g
(
I
)
S into models, where g
(
I
) tends to a saturation level when I becomes large:
(1)
g
(
I
)
=
β
I
(
t
-
τ
)
1
+
α
I
(
t
-
τ
)
,
where β
I
(
t
-
τ
) measures the infection force of the viruses and 1
/
1
+
α
I
(
t
-
τ
) measures the inhibition effect from the behavioral change of the susceptible ones when their number increases or from the crowding effect of the infected ones. It is noticed that τ represents the latent period; it means a fixed time during which some viruses develop in a susceptible computer and it is only after that time the susceptible one is converted to an infected one.Considering all these facts above, we can propose a new model with an economical use of the antivirus programs:(2)
S
˙
=
N
-
a
S
A
S
A
-
β
S
(
t
-
τ
)
I
(
t
-
τ
)
1
+
α
I
(
t
-
τ
)
-
μ
S
+
σ
R
,
I
˙
=
β
S
(
t
-
τ
)
I
(
t
-
τ
)
1
+
α
I
(
t
-
τ
)
-
a
I
A
A
I
-
δ
I
-
μ
I
,
R
˙
=
δ
I
-
σ
R
-
μ
R
,
A
˙
=
a
S
A
S
A
+
a
I
A
A
I
-
μ
A
.
For simplicity, let a
S
A
=
a
1, a
I
A
=
a
2, and a
S
A
+
a
I
A
=
a
A.


## 3. Mathematical Analysis

### 3.1. Virus-Free Equilibrium Point

Under the condition of virus-free, namely, we assumeI
=
0, there is no need to equip the computers with the antivirus programs, which means A
=
0. Then, bringing the equilibrium point E
*
=
[
S
*
,
I
*
,
R
*
,
A
*
] into (2), we get
(3)
N
-
a
1
S
*
A
*
-
β
S
*
I
*
1
+
α
I
*
-
μ
S
*
+
σ
R
*
=
0
,
β
S
*
I
*
1
+
α
I
*
-
a
2
A
*
I
*
-
δ
I
*
-
μ
I
*
=
0
,
δ
I
*
-
σ
R
*
-
μ
R
*
=
0
,
a
1
S
*
A
*
+
a
2
A
*
I
*
-
μ
A
*
=
0
.Based on (3), the virus-free equilibrium point can be calculated:
(4)
E
1
=
(
S
1
,
I
1
,
R
1
,
A
1
)
=
(
N
μ
,
0,0
,
0
)
.The characteristic equation of (2) at E
1 is given by the following:(5)
det
⁡
(
λ
+
μ
N
μ
β
e
-
λ
τ
-
σ
a
1
N
μ
0
λ
-
N
μ
β
e
-
λ
τ
+
(
δ
+
μ
)
0
0
0
-
δ
λ
+
σ
+
μ
0
0
0
0
λ
-
N
μ
a
1
+
μ
)
=
0
,which equals
(6)
(
λ
+
μ
)
(
λ
-
N
μ
β
e
-
λ
τ
+
(
δ
+
μ
)
)
×
(
λ
+
σ
+
μ
)
(
λ
-
N
μ
a
1
+
μ
)
=
(
λ
+
μ
)
(
λ
+
σ
+
μ
)
(
λ
-
N
μ
a
1
+
μ
)
×
[
λ
+
(
δ
+
μ
)
(
1
-
K
0
e
-
λ
τ
)
]
=
0
,
where K
0
=
β
N
/
μ
(
δ
+
μ
).Clearly, (6) always has two negative eigenvalues λ
1
=
-
μ, λ
2
=
-
σ
-
μ and two indefinite eigenvalues λ
3
=
(
N
/
μ
)
a
1
-
μ and λ
4
=
-
(
δ
+
μ
)
(
1
-
K
0
e
-
λ
τ
). We let K
1
=
a
1
N
/
μ
2.It is clear to see that ifK
1
<
1 and K
0
<
1, the other two eigenvalues must be negative too. So, the following theorem can be acquired.Theorem 1.
For the system (2), if K
0
=
β
N
/
μ
(
δ
+
μ
)
<
1 and K
1
=
a
1
N
/
μ
2
<
1 are satisfied, then the virus-free equilibrium point E
1
=
(
N
/
μ
,
0,0
,
0
) is locally asymptotically stable. Besides, if K
1
>
1, the equilibrium point E
1 is unstable.Remark 2.
Theorem1 investigates the local stability of the virus-free equilibrium point by analyzing the eigenvalues of the corresponding characteristic equation. We can see that when viruses do not appear, all computers in the network are subjected to possible infection.


### 3.2. Endemic Equilibrium Points

Endemic equilibrium points are characterized by the existence of infected ones in the network; that is,I
≠
0. First, we consider the case when the network has no antidotal node; namely, A
=
0. Then, it is not difficult to solve (3) when A
=
0, I
≠
0 and the solution is
(7)
E
2
=
(
S
2
,
I
2
,
R
2
,
A
2
)
=
(
S
2
,
I
2
,
R
2
,
0
)
,
where
(8)
I
2
=
(
σ
+
μ
)
[
β
N
-
μ
(
δ
+
μ
)
]
(
σ
+
μ
)
(
δ
+
μ
)
(
β
+
μ
α
)
-
β
σ
δ
,
S
2
=
(
δ
+
μ
)
(
1
+
α
I
2
)
β
,
R
2
=
δ
I
2
σ
+
μ
.It indicates that whenK
0
<
1, I
2
<
0. Consequently, the condition K
0
<
1 can avoid the existence of the equilibrium point E
2.Another more important case isA
≠
0 when I
≠
0. Again, calculating (3) and the endemic point is given by
(9)
E
3
=
(
S
3
,
I
3
,
R
3
,
A
3
)
=
(
μ
-
a
2
I
3
a
1
,
I
3
,
δ
I
3
σ
+
μ
,
β
(
μ
-
a
2
I
3
)
a
1
a
2
(
1
+
α
I
3
)
-
δ
+
μ
a
2
)
.Bringing theE
3 point into the first equation of (3), this leads to
(10)
N
-
a
1
S
3
A
3
-
β
S
3
I
3
1
+
α
I
3
-
μ
S
3
+
σ
R
3
=
N
-
a
1
μ
-
a
2
I
3
a
1
(
β
(
μ
-
a
2
I
3
)
a
1
a
2
(
1
+
α
I
3
)
-
δ
+
μ
a
2
)
-
β
I
3
1
+
α
I
3
·
μ
-
a
2
I
3
a
1
-
μ
μ
-
a
2
I
3
a
1
+
σ
δ
I
3
σ
+
μ
=
0
.Then, we expand (10) and merge the similar terms; we let I
=
I
3 without confusion:
(11)
p
1
(
1
+
α
I
)
+
p
2
I
(
1
+
α
I
)
+
p
3
(
μ
-
a
2
I
)
+
p
4
(
μ
-
a
2
I
)
(
1
+
α
I
)
=
0
,
where
(12)
p
1
=
a
1
a
2
N
(
σ
+
μ
)
,
p
2
=
a
1
a
2
σ
δ
,
p
3
=
-
μ
β
(
σ
+
μ
)
,
p
4
=
(
σ
+
μ
)
(
μ
a
1
+
δ
a
1
-
μ
a
2
)
.Furthermore, a quadratic equation of the variableI is obtained:
(13)
b
1
I
2
+
b
2
I
+
b
3
=
0
,
where
(14)
b
1
=
p
2
α
-
p
4
α
a
2
,
b
2
=
p
1
α
+
p
2
-
p
3
a
2
+
p
4
α
μ
-
p
4
a
2
,
b
3
=
p
1
+
p
3
μ
+
p
4
μ
.From (13), we have the following result.Theorem 3.
Suppose thatS, I, R, A are all positive, and K
2
>
1, K
3
>
1. Then (1)
ifb
1
=
0 and b
3
/
b
2
<
0 hold, there exists unique positive solution (9) for (3), and I
3
=
-
b
3
/
b
2;
(2)
whenΔ
>
0 always holds, then if b
2
/
b
1
<
0 and b
3
/
b
1
>
0, there exist two positive solutions, where I
3
(
1
)
=
I
+ and I
3
(
2
)
=
I
-; if b
3
/
b
1
<
0, there is only one positive solution I
3
=
I
+ with b
1
>
0 or I
3
=
I
- with b
1
<
0; if b
3
=
0 and b
2
/
b
1
<
0, there is only one positive solution I
3
=
-
b
2
/
b
1;
(3)
ifΔ
=
0 and b
2
/
b
1
<
0, there is only one positive solution I
3
=
-
b
2
/
2
b
1,where K
2
=
(
β
(
μ
-
a
2
I
)
)
/
(
a
1
(
1
+
α
I
)
(
δ
+
μ
)
), K
3
=
μ
/
a
2
I, Δ
=
b
2
2
-
4
b
1
b
3, I
+
=
(
-
b
2
+
b
2
2
-
4
b
1
b
3
)
/
2
b
1, and I
-
=
(
-
b
2
-
b
2
2
-
4
b
1
b
3
)
/
2
b
1.Remark 4.
Theorem3 mainly focuses on the existence of the positive equilibrium point of the system. We can see that all conditions in the theorem are easily verified.
When the equilibrium pointE
3 exists, the characteristic equation of (2) at the point is(15)
det
⁡
(
λ
+
a
1
A
3
+
β
I
3
e
-
λ
τ
1
+
α
I
3
+
μ
β
S
3
e
-
λ
τ
(
1
+
α
I
3
)
2
-
σ
a
1
S
3
-
β
I
3
e
-
λ
τ
1
+
α
I
3
λ
-
β
S
3
e
-
λ
τ
(
1
+
α
I
3
)
2
+
a
2
A
3
+
δ
+
μ
0
a
2
I
3
0
-
δ
λ
+
σ
+
μ
0
-
a
1
A
3
-
a
2
A
3
0
λ
)
=
0
.We setc
1
=
β
I
3
/
(
1
+
α
I
3
), c
2
=
β
S
3
/
(
1
+
α
I
3
)
2, c
3
=
a
1
A
3
+
μ, and c
4
=
a
2
A
3
+
δ
+
μ. For notational simplicity, we use the (
S
,
I
,
R
,
A
) in place of E
3. Then, the following four-degree exponential polynomial equation is obtained:
(16)
λ
4
+
d
1
λ
3
+
d
2
λ
2
+
d
3
λ
+
d
4
+
e
-
λ
τ
(
d
5
λ
3
+
d
6
λ
2
+
d
7
λ
+
d
8
)
=
0
,
where
(17)
d
1
=
c
3
+
c
4
+
σ
+
μ
,
d
2
=
a
1
2
A
S
+
a
2
2
A
I
+
c
3
c
4
+
(
σ
+
μ
)
(
c
3
+
c
4
)
,
d
3
=
a
1
2
A
S
(
c
4
+
σ
+
μ
)
+
a
2
2
A
I
(
c
3
+
σ
+
μ
)
+
c
3
c
4
(
σ
+
μ
)
,
d
4
=
(
σ
+
μ
)
(
a
1
2
c
4
A
S
+
a
2
2
c
3
A
I
)
+
a
1
a
2
σ
δ
A
I
,
d
5
=
c
1
-
c
2
,
d
6
=
c
1
(
c
4
+
σ
+
μ
)
-
c
2
(
c
3
+
σ
+
μ
)
,
d
7
=
A
(
a
2
I
+
a
1
S
)
(
c
1
a
2
-
c
2
a
1
)
-
c
1
σ
δ
+
(
c
1
c
4
-
c
2
c
3
)
(
σ
+
μ
)
,
d
8
=
A
(
a
2
I
+
a
1
S
)
(
c
1
a
2
-
c
2
a
1
)
(
σ
+
μ
)
.Multiplyinge
λ
τ on both sides of (16), it is obvious to get
(18)
J
=
(
λ
4
+
d
1
λ
3
+
d
2
λ
2
+
d
3
λ
+
d
4
)
e
λ
τ
+
(
d
5
λ
3
+
d
6
λ
2
+
d
7
λ
+
d
8
)
=
0
.Letλ
=
i
ω
0, τ
=
τ
0, and substituting this into (18), for the sake of simplicity, denote ω
0 and τ
0 by ω, τ, respectively; then (18) becomes
(19)
(
ω
4
-
i
d
1
ω
3
-
d
2
ω
2
+
i
d
3
ω
+
d
4
)
(
cos
⁡
ω
τ
+
i
sin
ω
τ
)
+
(
-
i
d
5
ω
3
-
d
6
ω
2
+
i
d
7
ω
+
d
8
)
=
0
.Separating the real and imaginary parts, we have(20)
(
ω
4
-
d
2
ω
2
+
d
4
)
cos
⁡
ω
τ
+
(
d
1
ω
3
-
d
3
ω
)
sin
ω
τ
=
d
6
ω
2
-
d
8
,
(
-
d
1
ω
3
+
d
3
ω
)
cos
⁡
ω
τ
+
(
ω
4
-
d
2
ω
2
+
d
4
)
sin
ω
τ
=
d
5
ω
3
-
d
7
ω
.By simple calculation, the following equations are obtained:(21)
cos
⁡
ω
τ
=
e
5
ω
6
+
e
6
ω
4
+
e
7
ω
2
+
e
8
ω
8
+
e
1
ω
6
+
e
2
ω
4
+
e
3
ω
2
+
e
4
,
(22)
sin
ω
τ
=
e
9
ω
7
+
e
10
ω
5
+
e
11
ω
3
+
e
12
ω
ω
8
+
e
1
ω
6
+
e
2
ω
4
+
e
3
ω
2
+
e
4
,
where
(23)
e
1
=
d
1
2
-
2
d
2
,
e
2
=
-
2
d
1
d
3
+
d
2
2
+
2
d
4
,
e
3
=
-
2
d
2
d
4
+
d
3
2
,
e
4
=
d
4
2
,
e
5
=
d
6
-
d
1
d
5
,
e
6
=
d
1
d
7
-
d
2
d
6
+
d
3
d
5
-
d
8
,
e
7
=
d
2
d
8
-
d
3
d
7
+
d
4
d
6
,
e
8
=
-
d
4
d
8
,
e
9
=
d
5
,
e
10
=
d
1
d
6
-
d
2
d
5
-
d
7
,
e
11
=
-
d
1
d
8
+
d
2
d
7
-
d
3
d
6
+
d
4
d
5
,
e
12
=
d
3
d
8
-
d
4
d
7
.As is known to all thatsi
n
2
ω
τ
+
cos
⁡
2
ω
τ
=
1, we get
(24)
ω
16
+
f
7
ω
14
+
f
6
ω
12
+
f
5
ω
10
+
f
4
ω
8
+
f
3
ω
6
+
f
2
ω
4
+
f
1
ω
2
+
f
0
=
0
,
where
(25)
f
7
=
2
e
1
-
e
9
2
,
f
6
=
e
1
2
+
2
e
2
-
e
5
2
-
2
e
9
e
10
,
f
5
=
2
e
1
e
2
+
2
e
3
-
2
e
5
e
6
-
2
e
9
e
11
-
e
10
2
,
f
4
=
2
e
1
e
3
+
e
2
2
+
2
e
4
-
2
e
5
e
7
-
e
6
2
-
2
e
9
e
12
-
2
e
10
e
11
,
f
3
=
2
e
1
e
4
+
2
e
2
e
3
-
2
e
5
e
8
-
2
e
6
e
7
-
2
e
10
e
12
-
e
11
2
,
f
2
=
2
e
2
e
4
+
e
3
2
-
2
e
6
e
8
-
e
7
2
-
2
e
11
e
12
,
f
1
=
2
e
3
e
4
-
2
e
7
e
8
-
e
12
2
,
f
0
=
e
4
2
-
e
8
2
.
Denote z
=
ω
2; (24) can be rewritten as
(26)
z
8
+
f
7
z
7
+
f
6
z
6
+
f
5
z
5
+
f
4
z
4
+
f
3
z
3
+
f
2
z
2
+
f
1
z
+
f
0
=
0
.
We suppose that(H1)
(26) has at least one positive real root.Without loss of generality, we can assume the equation hasl (1
≤
l
≤
8) positive real roots, which are represented as z
i (1
≤
i
≤
l); then ω
i
=
z
i.By (21), we get
(27)
τ
k
j
=
1
ω
k
{
arccos
(
e
5
ω
k
6
+
e
6
ω
k
4
+
e
7
ω
k
2
+
e
8
ω
k
8
+
e
1
ω
k
6
+
e
2
ω
k
4
+
e
3
ω
k
2
+
e
4
)
+
2
j
π
}
,
i
i
i
i
i
i
i
i
i
i
i
i
i
i
i
i
i
i
i
i
i
i
i
i
i
i
i
i
i
i
i
i
i
i
i
i
i
i
i
i
i
i
i
i
i
i
i
i
i
i
i
i
i
i
i
i
i
i
i
i
i
i
i
i
i
i
i
j
=
0,1
,
…
.From the early discussions, we know that the±
i
ω
k are a pair of purely imaginary roots of (16) with τ
k
j. Define
(28)
τ
0
=
τ
k
0
0
=
min
⁡
k
∈
{
1
,
…
,
l
}
{
τ
k
0
}
,
ω
0
=
ω
k
0
.It is noted that whenτ
=
0, (16) becomes
(29)
λ
4
+
d
1
λ
3
+
d
2
λ
2
+
d
3
λ
+
d
4
+
(
d
5
λ
3
+
d
6
λ
2
+
d
7
λ
+
d
8
)
=
λ
4
+
(
d
1
+
d
5
)
λ
3
+
(
d
2
+
d
6
)
λ
2
+
(
d
3
+
d
7
)
λ
+
(
d
4
+
d
8
)
=
0
.By virtue of the well-known Routh-Hurwitz criteria, a set of necessary and sufficient conditions for all roots of (29) to have the negative real part is given in the following form:
(30)
D
1
=
d
1
+
d
5
>
0
,
(31)
D
2
=
|
d
1
+
d
5
d
3
+
d
7
1
d
2
+
d
6
|
=
(
d
1
+
d
5
)
(
d
2
+
d
6
)
-
(
d
3
+
d
7
)
>
0
,
(32)
D
3
=
|
d
1
+
d
5
d
3
+
d
7
0
1
d
2
+
d
6
d
4
+
d
8
0
d
1
+
d
5
d
3
+
d
7
|
=
(
d
1
+
d
5
)
[
(
d
2
+
d
6
)
(
d
3
+
d
7
)
-
(
d
1
+
d
5
)
(
d
4
+
d
8
)
]
-
(
d
3
+
d
7
)
2
>
0
,
(33)
D
4
=
|
d
1
+
d
5
d
3
+
d
7
0
0
1
d
2
+
d
6
d
4
+
d
8
0
0
d
1
+
d
5
d
3
+
d
7
0
0
1
d
2
+
d
6
d
4
+
d
8
|
=
(
d
4
+
d
8
)
D
3
>
0
.If (30)–(33) hold, (29) has four roots with negative real parts, and therefore when τ
=
0, system (2) is stable near the equilibrium point E
3.In order to give the main results, it is necessary to make the following assumption:(H2)
Re
(
d
λ
d
τ
)
|
τ
=
τ
0
≠
0.In order to calculate the derivative ofλ with respect to τ in (18), it is followed by
(34)
d
λ
d
τ
=
-
∂
J
/
∂
τ
∂
J
/
∂
λ
=
-
(
λ
(
λ
4
+
d
1
λ
3
+
d
2
λ
2
+
d
3
λ
+
d
4
)
e
λ
τ
)
×
(
(
4
λ
3
+
3
d
1
λ
2
+
2
d
2
λ
+
d
3
)
e
λ
τ
+
τ
e
λ
τ
×
(
λ
4
+
d
1
λ
3
+
d
2
λ
2
+
d
3
λ
+
d
4
)
+
(
3
d
5
λ
2
+
2
d
6
λ
+
d
7
)
(
4
λ
3
+
3
d
1
λ
2
+
2
d
2
λ
+
d
3
)
e
λ
τ
+
τ
e
λ
τ
)
-
1
;
thus
(35)
(
d
λ
d
τ
)
-
1
=
(
(
4
λ
3
+
3
d
1
λ
2
+
2
d
2
λ
+
d
3
)
e
λ
τ
+
(
3
d
5
λ
2
+
2
d
6
λ
+
d
7
)
)
×
(
-
λ
e
λ
τ
(
λ
4
+
d
1
λ
3
+
d
2
λ
2
+
d
3
λ
+
d
4
)
)
-
1
-
τ
λ
=
(
(
4
λ
3
+
3
d
1
λ
2
+
2
d
2
λ
+
d
3
)
e
λ
τ
+
(
3
d
5
λ
2
+
2
d
6
λ
+
d
7
)
)
×
(
d
5
λ
4
+
d
6
λ
3
+
d
7
λ
2
+
d
8
λ
)
-
1
-
τ
λ
.Whenτ
=
τ
0, ±
i
ω
0 are a pair of purely imaginary roots of (16). Substituting the τ
0, i
ω
0 into (35) and denoting τ
0 and ω
0 by τ, ω for simplicity, then
(36)
(
d
λ
d
τ
)
-
1
|
λ
=
i
ω
0
,
τ
=
τ
0
=
(
(
-
3
d
1
ω
2
+
d
3
)
cos
⁡
ω
τ
+
(
4
ω
3
-
2
d
2
ω
)
sin
ω
τ
-
3
d
5
ω
2
+
d
7
)
×
(
(
d
5
ω
4
-
d
7
ω
2
)
+
i
(
-
d
6
ω
3
+
d
8
ω
)
)
-
1
+
i
(
(
-
4
ω
3
+
2
d
2
ω
)
cos
⁡
ω
τ
+
(
-
3
d
1
ω
2
+
d
3
)
sin
ω
τ
+
2
d
6
ω
)
×
(
(
d
5
ω
4
-
d
7
ω
2
)
+
i
(
-
d
6
ω
3
+
d
8
ω
)
)
-
1
-
τ
i
ω
.We setQ
=
(
d
5
ω
4
-
d
7
ω
2
)
2
+
(
-
d
6
ω
3
+
d
8
ω
)
2,
(37)
Q
·
Re
(
d
λ
d
τ
)
-
1
|
τ
=
τ
0
=
[
(
-
3
d
1
ω
2
+
d
3
)
cos
⁡
ω
τ
+
(
4
ω
3
-
2
d
2
ω
)
sin
ω
τ
-
3
d
5
ω
2
+
d
7
]
×
(
d
5
ω
4
-
d
7
ω
2
)
+
[
(
-
4
ω
3
+
2
d
2
ω
)
cos
⁡
ω
τ
+
(
-
3
d
1
ω
2
+
d
3
)
sin
ω
τ
+
2
d
6
ω
]
×
(
-
d
6
ω
3
+
d
8
ω
)
,
sgn
⁡
{
Re
(
d
λ
d
τ
)
|
τ
=
τ
0
}
=
sgn
⁡
{
Q
·
Re
(
d
λ
d
τ
)
-
1
|
τ
=
τ
0
}
.On the basis of (H2), we can know that the roots of characteristic equation (16) cross the imaginary axis as τ continuously varies from a number less than τ
0 to one greater than τ
0 by Rouche’s theorem [17]. Therefore, the transversality condition holds and the conditions for Hopf bifurcation are then satisfied at τ
=
τ
0.Lemma 5 (see [18]).
Consider the exponential polynomial(38)
P
(
λ
,
e
-
λ
τ
1
,
…
,
e
-
λ
τ
m
)
=
λ
n
+
p
1
(
0
)
λ
n
-
1
+
⋯
+
p
n
-
1
(
0
)
λ
+
p
n
(
0
)
+
[
p
1
(
1
)
λ
n
-
1
+
⋯
+
p
n
-
1
(
1
)
λ
+
p
n
(
1
)
]
e
-
λ
τ
1
+
⋯
+
[
p
1
(
m
)
λ
n
-
1
+
⋯
+
p
n
-
1
(
m
)
λ
+
p
n
(
m
)
]
e
-
λ
τ
m
,
where τ
i
≥
0
(
i
=
1,2
,
…
,
m
) and p
j
(
i
)
(
i
=
0,1
,
…
,
m
;
j
=
1,2
,
…
,
n
) are constants. As (
τ
1
,
τ
2
,
…
,
τ
m
) vary, the sum of the order of the zeros of P
(
λ
,
e
-
λ
τ
1
,
…
,
e
-
λ
τ
m
) on the open right plane can change only if a zero appears on or crosses the imaginary axis.
From Lemma5, it is easy to get the theorem.Theorem 6.
Suppose that (H1) and (H2) hold; then (1)
for system (2), the equilibrium point E
3 is asymptotically stable for τ
∈
[
0
,
τ
0
);
(2)
system (2) undergoes a Hopf bifurcation at the E
3 when τ
=
τ
0. That is to say, it has a branch of periodic solutions bifurcating from the E
3 near τ
=
τ
0, where τ
0 can be calculated by (28).
Furthermore, we can investigate the stability and direction of bifurcating periodic solutions by analyzing higher order terms according to Hassard et al. [19].Remark 7.
The analyses above study the existence of the Hopf bifurcation and obtain the critical value of the parameterτ. We can apply the theorem into the system; thus, through changing some parameters, some bad performances of the model can be avoided.Remark 8.
When we setN
=
0, μ
=
0, α
=
0, and τ
=
0, the present model can be transformed to the model in [14]. However, [14] only calculated the disease-free/endemic equilibrium points and analyzed the stability. This paper expands the model and makes a detailed analysis about the bifurcations.Remark 9.
In this paper, we propose a new group called antidotal population, that is,A
(
t
), which is more reasonable when we study the computer viruses. It is well known that when dealing with Hopf bifurcation, the complexity of computation increases significantly as the dimension of system increases. Sometimes, it even cannot be calculated, especially when nonlinear terms β
S
(
t
-
τ
)
I
(
t
-
τ
)
/
(
1
+
α
I
(
t
-
τ
)
) exist in our system. Partly because of this, the majority of the literatures published use the traditional SIR three-dimensional model to study the epidemic model or virus [7, 10, 20], which has some limitations.


### 3.3. Stability and Direction of the Hopf Bifurcation

We have obtained the conditions under which a family of periodic solutions bifurcate from the positive equilibriumE
3 at the critical value of τ
0. In this subsection, the formulae for determining the direction of Hopf bifurcation and stability of bifurcating periodic solutions of system at τ
0 will be presented by employing the normal form theory and the center manifold reduction [19, 21–24].For convenience, letx
1
=
S
-
S
3, x
2
=
I
-
I
3, x
3
=
R
-
R
3, x
4
=
A
-
A
3, x
-
i
(
t
)
=
x
i
(
τ
t
), and τ
=
τ
0
+
u; we drop the bars for simplification of notations. In the light of multivariate Taylor expansion, system (2) can be transformed into an FDE in C
=
C
(
[
-
1,0
]
,
R
3
) as
(39)
x
˙
(
t
)
=
L
u
(
x
t
)
+
f
(
u
,
x
t
)
,
where x
(
t
)
=
(
x
1
(
t
)
,
x
2
(
t
)
,
x
3
(
t
)
,
x
4
(
t
)
)
T
∈
R
4, and L
u
:
C
→
R, f
:
R
×
C
→
R are given, respectively, by
(40)
L
u
(
ϕ
)
=
(
τ
0
+
u
)
(
-
c
3
0
σ
-
a
1
S
0
-
c
4
0
-
a
2
I
0
δ
-
σ
-
μ
0
a
1
A
a
2
A
0
0
)
(
ϕ
1
(
0
)
ϕ
2
(
0
)
ϕ
3
(
0
)
ϕ
4
(
0
)
)
+
(
τ
0
+
u
)
(
-
c
1
-
c
2
0
0
c
1
c
2
0
0
0
0
0
0
0
0
0
0
)
(
ϕ
1
(
-
1
)
ϕ
2
(
-
1
)
ϕ
3
(
-
1
)
ϕ
4
(
-
1
)
)
=
(
τ
0
+
u
)
B
1
ϕ
(
0
)
+
(
τ
0
+
u
)
B
2
ϕ
(
-
1
)
,
f
(
u
,
ϕ
)
=
(
τ
0
+
u
)
×
(
l
1
ϕ
2
2
(
-
1
)
-
l
2
ϕ
1
(
-
1
)
ϕ
2
(
-
1
)
-
a
1
ϕ
1
(
0
)
ϕ
4
(
0
)
+
h
.
o
.
t
-
l
1
ϕ
2
2
(
-
1
)
+
l
2
ϕ
1
(
-
1
)
ϕ
2
(
-
1
)
-
a
2
ϕ
2
(
0
)
ϕ
4
(
0
)
+
h
.
o
.
t
0
a
1
ϕ
1
(
0
)
ϕ
4
(
0
)
+
a
2
ϕ
2
(
0
)
ϕ
4
(
0
)
)
,
where l
1
=
α
c
2
/
(
1
+
α
I
), l
2
=
c
2
/
S, and L
u is a one-parameter family of bounded linear operators in C
[
-
1,0
]
→
R
4.By the Riesz representation theorem, there exists a functionη
(
θ
,
u
) of bounded variation for θ
∈
[
-
1,0
], such that
(41)
L
u
ϕ
=
∫
-
1
0
d
η
(
θ
,
u
)
ϕ
(
θ
)
for
ϕ
∈
C
1
(
[
-
1,0
]
,
R
4
)
.
In fact, we can choose
(42)
η
(
θ
,
u
)
=
(
τ
0
+
u
)
B
1
δ
(
θ
)
-
(
τ
0
+
u
)
B
2
δ
(
θ
+
1
)
,
where δ
(
θ
) is the Dirac delta function.Forϕ
∈
C
1
(
[
-
1,0
]
,
R
4
), define
(43)
A
(
u
)
ϕ
=
{
d
ϕ
(
θ
)
d
θ
,
-
1
≤
θ
<
0
,
∫
-
1
0
d
η
(
θ
,
u
)
ϕ
(
θ
)
,
θ
=
0
,
R
(
u
)
ϕ
=
{
0
,
θ
∈
[
-
1,0
)
,
f
(
u
,
ϕ
)
,
θ
=
0
.Then, system (39) is equivalent to
(44)
x
˙
t
=
A
(
u
)
x
t
+
R
(
u
)
x
t
,
where x
t
(
θ
)
=
x
(
t
+
θ
), for θ
∈
[
-
1,0
), and x
=
(
x
1
,
x
2
,
x
3
,
x
4
)
T.Forψ
∈
C
1
(
[
0,1
]
,
(
R
4
)
*
), define
(45)
A
*
(
u
)
ψ
=
{
-
d
ψ
(
s
)
d
s
,
0
<
s
≤
1
,
∫
-
1
0
d
η
T
(
t
,
u
)
ϕ
(
-
t
)
,
s
=
0
.In order to normalize the eigenvectors of operatorA and adjoint operator A
*, the following bilinear inner product is needed to introduce
(46)
〈
ψ
(
s
)
,
ϕ
(
θ
)
〉
=
ψ
-
(
0
)
ϕ
(
0
)
-
∫
-
1
0
∫
ξ
=
0
θ
ψ
-
(
ξ
-
θ
)
d
η
(
θ
)
ϕ
(
ξ
)
d
ξ
,
where η
(
θ
)
=
η
(
θ
,
0
).From the early discussions, we know that±
i
τ
0
ω
0 are eigenvalues of A
(
0
) and other eigenvalues have strictly negative real parts; thus, they are also eigenvalues of A
*. Define that
(47)
q
(
θ
)
=
(
1
,
q
2
,
q
3
,
q
4
)
T
e
i
τ
0
ω
0
θ
,
-
1
<
θ
≤
0
,
which is the eigenvector of A
(
0
) belonging to the eigenvalue i
τ
0
ω
0; namely, A
(
0
)
q
(
θ
)
=
i
τ
0
ω
0
q
(
θ
). Then, we can easily obtain
(48)
τ
0
(
i
ω
0
+
c
3
+
c
1
e
-
i
τ
0
ω
0
c
2
e
-
i
τ
0
ω
0
-
σ
a
1
S
-
c
1
e
-
i
τ
0
ω
0
i
ω
0
+
c
4
-
c
2
e
-
i
τ
0
ω
0
0
a
2
I
0
-
δ
i
ω
0
+
σ
+
μ
0
-
a
1
A
-
a
2
A
0
i
ω
0
)
×
q
(
0
)
=
(
0
0
0
0
)
.
Hence, we obtain
(49)
q
2
=
i
ω
0
c
1
e
-
i
τ
0
ω
0
-
a
1
a
2
I
A
i
ω
0
(
i
ω
0
+
c
4
-
c
2
e
-
i
τ
0
ω
0
)
+
a
2
2
I
A
,
q
3
=
δ
(
i
ω
0
c
1
e
-
i
τ
0
ω
0
-
a
1
a
2
I
A
)
(
i
ω
0
+
σ
+
μ
)
[
i
ω
0
(
i
ω
0
+
c
4
-
c
2
e
-
i
τ
0
ω
0
)
+
a
2
2
I
A
]
,
q
4
=
a
1
A
i
ω
0
+
a
2
A
(
i
ω
0
c
1
e
-
i
τ
0
ω
0
-
a
1
a
2
I
A
)
-
ω
0
2
(
i
ω
0
+
c
4
-
c
2
e
-
i
τ
0
ω
0
)
+
i
ω
0
a
2
2
I
A
.
Suppose that the eigenvector q
* of A
* belonging to the eigenvalue -
i
τ
0
ω
0 is
(50)
q
*
(
s
)
=
1
ρ
(
1
,
q
2
*
,
q
3
*
,
q
4
*
)
e
i
τ
0
ω
0
s
,
0
≤
s
<
1
.
Similar to the calculation of (49), we can get
(51)
q
2
*
=
i
ω
0
(
i
ω
0
-
c
3
-
c
1
e
i
τ
0
ω
0
)
+
a
1
2
A
S
-
i
ω
0
c
1
e
i
τ
0
ω
0
-
a
1
a
2
I
A
,
q
3
*
=
σ
-
i
ω
0
+
σ
+
μ
,
q
4
*
=
a
1
S
i
ω
0
+
a
2
I
i
ω
0
(
i
ω
0
-
c
3
-
c
1
e
i
τ
0
ω
0
)
+
a
1
2
a
2
I
A
S
ω
0
2
c
1
e
i
τ
0
ω
0
-
i
ω
0
a
1
a
2
I
A
.
Let 〈
q
*
,
q
〉
=
1; then
(52)
〈
q
*
(
s
)
,
q
(
θ
)
〉
=
1
ρ
-
(
1
+
q
2
q
2
*
+
q
3
q
3
*
+
q
4
q
4
*
)
-
∫
θ
=
-
1
0
∫
ξ
=
0
θ
1
ρ
-
(
1
,
q
-
2
*
,
q
-
3
*
,
q
-
4
*
)
e
-
i
τ
0
ω
0
(
ξ
-
θ
)
d
η
(
θ
)
×
(
1
q
2
q
3
q
4
)
e
i
τ
0
ω
0
ξ
d
ξ
=
1
ρ
-
(
1
+
q
2
q
2
*
+
q
3
q
3
*
+
q
4
q
4
*
)
-
∫
-
1
0
1
ρ
-
(
1
,
q
-
2
*
,
q
-
3
*
,
q
-
4
*
)
θ
e
i
τ
0
ω
0
θ
d
η
(
θ
)
(
1
q
2
q
3
q
4
)
=
1
ρ
-
{
τ
0
e
-
i
τ
0
ω
0
(
1
+
q
2
q
2
*
+
q
3
q
3
*
+
q
4
q
4
*
)
+
τ
0
e
-
i
τ
0
ω
0
[
(
q
-
2
*
-
1
)
(
c
1
+
α
c
2
)
]
}
,
which leads to
(53)
ρ
-
=
(
1
+
q
2
q
2
*
+
q
3
q
3
*
+
q
4
q
4
*
)
+
τ
0
e
-
i
τ
0
ω
0
[
(
q
-
2
*
-
1
)
(
c
1
+
α
c
2
)
]
.In the following, we apply the method in [19] to compute the coordinates describing the center manifold Ω
0 near u
=
0. Let x
t be the solution of (39) when u
=
0. We define
(54)
z
(
t
)
=
〈
q
*
,
x
t
〉
,
W
(
t
,
θ
)
=
x
t
(
θ
)
-
2
Re
[
z
(
t
)
q
(
θ
)
]
.On the center manifoldΩ
0, we have W
(
t
,
θ
)
=
W
(
z
(
t
)
,
z
-
(
t
)
,
θ
), where
(55)
W
(
z
(
t
)
,
z
-
(
t
)
,
θ
)
=
W
20
(
θ
)
z
2
2
+
W
11
(
θ
)
z
z
-
+
W
02
(
θ
)
z
-
2
2
+
W
30
(
θ
)
z
3
6
+
⋯
.In fact,z and z
- are local coordinates of center manifold Ω
0 in the direction of q and q
*, respectively. For solution x
t
∈
Ω
0 of (44), since u
=
0, we get
(56)
z
˙
(
t
)
=
〈
q
*
,
x
˙
t
〉
=
i
τ
0
ω
0
z
(
t
)
+
q
-
*
(
0
)
×
f
(
0
,
W
(
t
,
0
)
+
2
Re
[
z
(
t
)
q
(
0
)
]
)
=
i
τ
0
ω
0
z
(
t
)
+
q
-
*
(
0
)
f
0
(
z
,
z
-
)
≜
i
τ
0
ω
0
z
(
t
)
+
g
(
z
,
z
-
)
.
It is noted that
(57)
g
(
z
,
z
-
)
=
q
-
*
(
0
)
f
0
(
z
,
z
-
)
=
g
20
z
2
2
+
g
11
z
z
-
+
g
02
z
-
2
2
+
g
21
z
2
z
-
2
+
⋯
.By (54), we know that
(58)
x
t
(
θ
)
=
(
x
1
t
(
θ
)
,
x
2
t
(
θ
)
,
x
3
t
(
θ
)
,
x
4
t
(
θ
)
)
=
W
(
t
,
θ
)
+
z
q
(
θ
)
+
z
-
·
q
-
(
θ
)
.Considering (47) and (55), we have
(59)
x
1
t
(
0
)
=
z
+
z
-
+
W
20
(
1
)
(
0
)
z
2
2
+
W
11
(
1
)
(
0
)
z
z
-
+
W
02
(
1
)
(
0
)
z
-
2
2
+
o
(
|
(
z
,
z
-
)
|
3
)
,
x
2
t
(
0
)
=
q
2
z
+
q
-
2
z
-
+
W
20
(
2
)
(
0
)
z
2
2
+
W
11
(
2
)
(
0
)
z
z
-
+
W
02
(
2
)
(
0
)
z
-
2
2
+
o
(
|
(
z
,
z
-
)
|
3
)
,
x
4
t
(
0
)
=
q
4
z
+
q
-
4
z
-
+
W
20
(
4
)
(
0
)
z
2
2
+
W
11
(
4
)
(
0
)
z
z
-
+
W
02
(
4
)
(
0
)
z
-
2
2
+
o
(
|
(
z
,
z
-
)
|
3
)
,
x
1
t
(
-
1
)
=
z
e
-
i
τ
0
ω
0
+
z
-
e
i
τ
0
ω
0
+
W
20
(
1
)
(
-
1
)
z
2
2
+
W
11
(
1
)
(
-
1
)
z
z
-
+
W
02
(
1
)
(
-
1
)
z
-
2
2
+
o
(
|
(
z
,
z
-
)
|
3
)
,
x
2
t
(
-
1
)
=
z
q
2
e
-
i
τ
0
ω
0
+
z
-
·
q
-
2
e
i
τ
0
ω
0
+
W
20
(
2
)
(
-
1
)
z
2
2
+
W
11
(
2
)
(
-
1
)
z
z
-
+
W
02
(
2
)
(
-
1
)
z
-
2
2
+
o
(
|
(
z
,
z
-
)
|
3
)
.It follows that(60)
g
(
z
,
z
-
)
=
q
-
*
(
0
)
f
0
(
z
,
z
-
)
=
τ
0
ρ
-
(
1
,
q
-
2
*
,
q
-
3
*
,
q
-
4
*
)
×
(
l
1
x
2
t
2
(
-
1
)
-
l
2
x
1
t
(
-
1
)
x
2
t
(
-
1
)
-
a
1
x
1
t
(
0
)
x
4
t
(
0
)
+
h
.
o
.
t
-
l
1
x
2
t
2
(
-
1
)
+
l
2
x
1
t
(
-
1
)
x
2
t
(
-
1
)
-
a
2
x
2
t
(
0
)
x
4
t
(
0
)
+
h
.
o
.
t
0
a
1
x
1
t
(
0
)
x
4
t
(
0
)
+
a
2
x
2
t
(
0
)
x
4
t
(
0
)
)
=
τ
0
ρ
-
{
1
2
[
l
1
(
1
-
q
-
2
*
)
q
2
2
e
-
2
i
τ
0
ω
0
+
l
2
(
q
-
2
*
-
1
)
q
2
e
-
2
i
τ
0
ω
0
+
a
1
(
q
-
4
*
-
1
)
q
4
+
a
2
(
q
-
4
*
-
q
-
2
*
)
q
2
q
4
l
2
(
q
-
2
*
-
1
)
q
2
e
-
2
i
τ
0
ω
0
]
z
2
+
[
l
1
(
1
-
q
-
2
*
)
q
-
2
2
e
2
i
τ
0
ω
0
+
l
2
(
q
-
2
*
-
1
)
q
-
2
e
2
i
τ
0
ω
0
+
a
1
(
q
-
4
*
-
1
)
q
-
4
+
a
2
(
q
-
4
*
-
q
-
2
*
)
q
-
2
q
-
4
l
2
(
q
-
2
*
-
1
)
q
2
e
-
2
i
τ
0
ω
0
]
z
-
2
+
[
l
1
(
1
-
q
-
2
*
)
2
q
2
q
-
2
+
l
2
(
q
-
2
*
-
1
)
(
q
2
+
q
-
2
)
+
a
1
(
q
-
4
*
-
1
)
(
q
4
+
q
-
4
)
+
a
2
(
q
-
4
*
-
q
-
2
*
)
(
q
2
q
-
4
+
q
-
2
q
4
)
]
z
z
-
+
[
1
2
q
-
2
W
20
(
4
)
(
0
)
l
1
(
1
-
q
-
2
*
)
×
(
2
q
2
W
11
(
2
)
(
-
1
)
e
-
i
τ
0
ω
0
+
q
-
2
W
20
(
2
)
(
-
1
)
e
i
τ
0
ω
0
)
+
l
2
(
q
-
2
*
-
1
)
×
(
W
11
(
2
)
(
-
1
)
e
-
i
τ
0
ω
0
+
1
2
W
20
(
2
)
(
-
1
)
e
i
τ
0
ω
0
+
1
2
q
-
2
W
20
(
1
)
(
-
1
)
e
i
τ
0
ω
0
+
q
2
W
11
(
1
)
(
-
1
)
e
-
i
τ
0
ω
0
1
2
)
+
a
1
(
q
-
4
*
-
1
)
×
(
W
11
(
4
)
(
0
)
+
1
2
W
20
(
4
)
(
0
)
+
1
2
q
-
4
W
20
(
1
)
(
0
)
+
q
4
W
11
(
1
)
(
0
)
)
+
a
2
(
q
-
4
*
-
q
-
2
*
)
×
(
q
2
W
11
(
4
)
(
0
)
+
1
2
q
-
2
W
20
(
4
)
(
0
)
+
1
2
q
-
4
W
20
(
2
)
(
0
)
+
q
4
W
11
(
2
)
(
0
)
)
]
×
z
2
z
-
+
⋯
1
2
}
.
Comparing the coefficients in (57) with those in (60), it follows that
(61)
g
20
=
2
τ
0
ρ
-
[
l
1
(
1
-
q
-
2
*
)
q
2
2
e
-
2
i
τ
0
ω
0
+
l
2
(
q
-
2
*
-
1
)
q
2
e
-
2
i
τ
0
ω
0
+
a
1
(
q
-
4
*
-
1
)
q
4
+
a
2
(
q
-
4
*
-
q
-
2
*
)
q
2
q
4
q
2
2
e
-
2
i
τ
0
ω
0
]
,
g
02
=
2
τ
0
ρ
-
[
l
1
(
1
-
q
-
2
*
)
q
-
2
2
e
2
i
τ
0
ω
0
+
l
2
(
q
-
2
*
-
1
)
q
-
2
e
2
i
τ
0
ω
0
+
a
1
(
q
-
4
*
-
1
)
q
-
4
+
a
2
(
q
-
4
*
-
q
-
2
*
)
q
-
2
q
-
4
q
2
2
e
-
2
i
τ
0
ω
0
]
,
g
11
=
τ
0
ρ
-
[
l
1
(
1
-
q
-
2
*
)
2
q
2
q
-
2
+
l
2
(
q
-
2
*
-
1
)
(
q
2
+
q
-
2
)
+
a
1
(
q
-
4
*
-
1
)
(
q
4
+
q
-
4
)
+
a
2
(
q
-
4
*
-
q
-
2
*
)
(
q
2
q
-
4
+
q
-
2
q
4
)
]
,
g
21
=
2
τ
0
ρ
-
[
1
2
q
-
4
W
20
(
2
)
l
1
(
1
-
q
-
2
*
)
×
(
2
q
2
W
11
(
2
)
(
-
1
)
e
-
i
τ
0
ω
0
+
q
-
2
W
20
(
2
)
(
-
1
)
e
i
τ
0
ω
0
)
+
l
2
(
q
-
2
*
-
1
)
×
(
W
11
(
2
)
(
-
1
)
e
-
i
τ
0
ω
0
+
1
2
W
20
(
2
)
(
-
1
)
e
i
τ
0
ω
0
+
1
2
q
-
2
W
20
(
1
)
(
-
1
)
e
i
τ
0
ω
0
+
q
2
W
11
(
1
)
(
-
1
)
e
-
i
τ
0
ω
0
)
+
a
1
(
q
-
4
*
-
1
)
×
(
W
11
(
4
)
(
0
)
+
1
2
W
20
(
4
)
(
0
)
+
1
2
q
-
4
W
20
(
1
)
(
0
)
+
q
4
W
11
(
1
)
(
0
)
)
+
a
2
(
q
-
4
*
-
q
-
2
*
)
×
(
q
2
W
11
(
4
)
(
0
)
+
1
2
q
-
2
W
20
(
4
)
(
0
)
+
1
2
q
-
4
W
20
(
2
)
(
0
)
+
q
4
W
11
(
2
)
(
0
)
)
]
.Since theW
20
(
θ
) and W
11
(
θ
) exist in (61), we need to compute them. From (44) and (54), we have
(62)
W
˙
=
x
˙
t
-
z
˙
q
-
z
-
˙
·
q
-
=
A
x
t
+
R
x
t
-
[
i
τ
0
ω
0
z
+
q
-
*
(
0
)
f
0
(
z
,
z
-
)
]
q
-
[
-
i
τ
0
ω
0
z
-
+
q
*
(
0
)
f
-
0
(
z
,
z
-
)
]
q
-
=
A
(
W
+
2
Re
(
z
q
)
)
+
R
x
t
-
2
Re
[
q
-
*
(
0
)
f
0
(
z
,
z
-
)
q
]
-
2
Re
[
i
τ
0
ω
0
z
q
]
=
{
A
W
-
2
Re
[
q
-
*
(
0
)
f
0
(
z
,
z
-
)
q
(
θ
)
]
,
-
1
≤
θ
<
0
,
A
W
-
2
Re
[
q
-
*
(
0
)
f
0
(
z
,
z
-
)
q
(
θ
)
]
+
f
,
θ
=
0
,
≜
A
W
+
H
(
z
,
z
-
,
θ
)
,
where
(63)
H
(
z
,
z
-
,
θ
)
=
H
20
(
θ
)
z
2
2
+
H
11
(
θ
)
z
z
-
+
H
02
(
θ
)
z
-
2
2
+
⋯
.
Substituting (55) and (63) into (62), we get
(64)
W
˙
=
(
A
W
20
+
H
20
)
z
2
2
+
(
A
W
11
+
H
11
)
z
z
-
+
(
A
W
02
+
H
02
)
z
-
2
2
+
⋯
.
Taking the derivate of W with respect to t in (55), we have
(65)
W
˙
=
W
z
z
˙
+
W
z
-
z
-
˙
=
(
W
20
z
+
W
11
z
-
+
W
30
z
2
2
+
⋯
)
(
i
τ
0
ω
0
z
+
g
)
+
(
W
11
z
+
W
02
z
-
+
⋯
)
(
-
i
τ
0
ω
0
z
-
+
g
-
)
=
i
τ
0
ω
0
W
20
z
2
+
z
z
-
(
i
τ
0
ω
0
W
11
-
i
τ
0
ω
0
W
11
)
-
i
τ
0
ω
0
W
02
z
-
2
+
⋯
.
Then, together with the two above equations, we obtain
(66)
A
W
20
+
H
20
=
2
i
τ
0
ω
0
W
20
,
A
W
11
+
H
11
=
0
.
By (62), we know that
(67)
H
(
z
,
z
-
,
θ
)
=
-
g
q
(
θ
)
-
g
-
·
q
-
(
θ
)
+
R
x
t
=
-
(
g
20
z
2
2
+
g
11
z
z
-
+
g
02
z
-
2
2
+
⋯
)
q
(
θ
)
-
(
g
-
20
z
-
2
2
+
g
-
11
z
z
-
+
g
-
02
z
2
2
+
⋯
)
q
-
(
θ
)
+
R
x
t
.
We know that when -
1
≤
θ
<
0  R
x
t
=
0, comparing the coefficients of the z
2, z
z
- in (63) and (67), respectively, we can find
(68)
H
20
(
θ
)
=
-
g
20
q
(
θ
)
-
g
-
02
q
-
(
θ
)
,
H
11
(
θ
)
=
-
g
11
q
(
θ
)
-
g
-
11
q
-
(
θ
)
.
Combining (66) with (68), it follows that
(69)
W
˙
20
(
θ
)
=
2
i
τ
0
ω
0
W
20
(
θ
)
+
g
20
q
(
θ
)
+
g
-
02
q
-
(
θ
)
,
W
˙
11
(
θ
)
=
g
11
q
(
θ
)
+
g
-
11
q
-
(
θ
)
,
iiiiiiiiiiiiiiiiiiiiiiiiii
-
1
≤
θ
<
0
.
We solve (69) and get the solutions in the following form:
(70)
W
20
(
θ
)
=
i
g
20
τ
0
ω
0
q
(
0
)
e
i
τ
0
ω
0
θ
-
g
-
02
3
i
τ
0
ω
0
q
-
(
0
)
e
-
i
τ
0
ω
0
θ
+
E
1
e
2
i
τ
0
ω
0
θ
,
W
11
(
θ
)
=
g
11
i
τ
0
ω
0
q
(
0
)
e
i
τ
0
ω
0
θ
-
g
-
11
i
τ
0
ω
0
q
-
(
0
)
e
-
i
τ
0
ω
0
θ
+
E
2
.
Next, we consider (66) again when θ
=
0; we can see that
(71)
A
W
20
(
0
)
=
∫
-
1
0
d
η
(
θ
)
W
20
(
θ
)
=
2
i
τ
0
ω
0
W
20
(
0
)
-
H
20
(
0
)
,
A
W
11
(
0
)
=
∫
-
1
0
d
η
(
θ
)
W
11
(
θ
)
=
-
H
11
(
0
)
.
Furthermore, from (62) and (63), we have
(72)
H
20
(
0
)
=
-
g
20
q
(
0
)
-
g
-
02
q
-
(
0
)
+
2
τ
0
(
l
1
q
2
2
e
-
2
i
τ
0
ω
0
-
l
2
q
2
e
-
2
i
τ
0
ω
0
-
a
1
q
4
-
l
1
q
2
2
e
-
2
i
τ
0
ω
0
+
l
2
q
2
e
-
2
i
τ
0
ω
0
-
a
2
q
2
q
4
0
a
1
q
4
+
a
2
q
2
q
4
)
,
(73)
H
11
(
0
)
=
-
g
11
q
(
0
)
-
g
-
11
q
-
(
0
)
+
τ
0
(
2
l
1
|
q
2
|
2
-
l
2
(
q
2
+
q
-
2
)
-
a
1
(
q
4
+
q
-
4
)
-
2
l
1
|
q
2
|
2
+
l
2
(
q
2
+
q
-
2
)
-
a
2
(
q
2
q
-
4
+
q
-
2
q
4
)
0
a
1
(
q
4
+
q
-
4
)
+
a
2
(
q
2
q
-
4
+
q
-
2
q
4
)
)
.Substitute (70) and (72) into (71) and notice that
(74)
(
i
τ
0
ω
0
I
-
∫
-
1
0
e
i
τ
0
ω
0
θ
d
η
(
θ
)
)
q
(
0
)
=
0
,
(
-
i
τ
0
ω
0
I
-
∫
-
1
0
e
-
i
τ
0
ω
0
θ
d
η
(
θ
)
)
q
-
(
0
)
=
0
.
It is easy to yield
(75)
∫
-
1
0
d
η
(
θ
)
E
1
e
2
i
τ
0
ω
0
θ
=
2
i
τ
0
ω
0
E
1
-
2
τ
0
(
l
1
q
2
2
e
-
2
i
τ
0
ω
0
-
l
2
q
2
e
-
2
i
τ
0
ω
0
-
a
1
q
4
-
l
1
q
2
2
e
-
2
i
τ
0
ω
0
+
l
2
q
2
e
-
2
i
τ
0
ω
0
-
a
2
q
2
q
4
0
a
1
q
4
+
a
2
q
2
q
4
)
,
which leads to(76)
E
1
=
2
(
2
i
ω
0
+
c
3
+
c
1
e
-
2
i
τ
0
ω
0
c
2
e
-
2
i
τ
0
ω
0
-
σ
a
1
S
-
c
1
e
-
2
i
τ
0
ω
0
2
i
ω
0
+
c
4
-
c
2
e
-
2
i
τ
0
ω
0
0
a
2
I
0
-
δ
2
i
ω
0
+
σ
+
μ
0
-
a
1
A
-
a
2
A
0
2
i
ω
0
)
-
1
·
(
l
1
q
2
2
e
-
2
i
τ
0
ω
0
-
l
2
q
2
e
-
2
i
τ
0
ω
0
-
a
1
q
4
-
l
1
q
2
2
e
-
2
i
τ
0
ω
0
+
l
2
q
2
e
-
2
i
τ
0
ω
0
-
a
2
q
2
q
4
0
a
1
q
4
+
a
2
q
2
q
4
)
≜
2
M
1
-
1
M
2
.By the similar way we can get theE
2 as in the following:
(77)
E
2
=
(
c
3
+
c
1
c
2
-
σ
a
1
S
-
c
1
c
4
-
c
2
0
a
2
I
0
-
δ
σ
+
μ
0
-
a
1
A
-
a
2
A
0
0
)
-
1
×
(
2
l
1
|
q
2
|
2
-
l
2
(
q
2
+
q
-
2
)
-
a
1
(
q
4
+
q
-
4
)
-
2
l
1
|
q
2
|
2
+
l
2
(
q
2
+
q
-
2
)
-
a
2
(
q
2
q
-
4
+
q
-
2
q
4
)
0
a
1
(
q
4
+
q
-
4
)
+
a
2
(
q
2
q
-
4
+
q
-
2
q
4
)
)
.Thus, we can calculate theW
20
(
θ
) and W
11
(
θ
) from (70). Consequently, g
21 in (61) can be received. Then, we need to compute the following parameters [19]:
(78)
C
1
(
0
)
=
i
2
τ
0
ω
0
(
g
20
g
11
-
2
|
g
11
|
2
-
|
g
02
|
2
3
)
+
g
21
2
,
(79)
μ
2
=
-
Re
{
C
1
(
0
)
}
Re
{
λ
′
(
τ
0
)
}
,
β
2
=
2
Re
{
C
1
(
0
)
}
,
T
2
=
-
Im
⁡
{
C
1
(
0
)
}
+
μ
2
Im
⁡
{
λ
′
(
τ
0
)
}
τ
0
ω
0
,
which determine the characteristic of bifurcating periodic solutions in the center manifold at the critical value τ
0. More details are given in the following theorem.Theorem 10.
Under the conditions of Theorem6 one has the following.(1)
u
=
0 is Hopf bifurcation value of system (39).
(2)
The direction of Hopf bifurcation is determined by the sign ofμ
2: if μ
2
>
0, the Hopf bifurcation is supercritical; if μ
2
<
0, the Hopf bifurcation is subcritical.
(3)
The stability of bifurcating periodic solutions is determined byβ
2: if β
2
<
0, the periodic solutions are stable; if β
2
>
0, they are unstable.
(4)
The sign ofT
2 determines the period of the bifurcating periodic solutions: if T
2
>
0, the period increases; if T
2
<
0, the period decreases.


## 3.1. Virus-Free Equilibrium Point

Under the condition of virus-free, namely, we assumeI
=
0, there is no need to equip the computers with the antivirus programs, which means A
=
0. Then, bringing the equilibrium point E
*
=
[
S
*
,
I
*
,
R
*
,
A
*
] into (2), we get
(3)
N
-
a
1
S
*
A
*
-
β
S
*
I
*
1
+
α
I
*
-
μ
S
*
+
σ
R
*
=
0
,
β
S
*
I
*
1
+
α
I
*
-
a
2
A
*
I
*
-
δ
I
*
-
μ
I
*
=
0
,
δ
I
*
-
σ
R
*
-
μ
R
*
=
0
,
a
1
S
*
A
*
+
a
2
A
*
I
*
-
μ
A
*
=
0
.Based on (3), the virus-free equilibrium point can be calculated:
(4)
E
1
=
(
S
1
,
I
1
,
R
1
,
A
1
)
=
(
N
μ
,
0,0
,
0
)
.The characteristic equation of (2) at E
1 is given by the following:(5)
det
⁡
(
λ
+
μ
N
μ
β
e
-
λ
τ
-
σ
a
1
N
μ
0
λ
-
N
μ
β
e
-
λ
τ
+
(
δ
+
μ
)
0
0
0
-
δ
λ
+
σ
+
μ
0
0
0
0
λ
-
N
μ
a
1
+
μ
)
=
0
,which equals
(6)
(
λ
+
μ
)
(
λ
-
N
μ
β
e
-
λ
τ
+
(
δ
+
μ
)
)
×
(
λ
+
σ
+
μ
)
(
λ
-
N
μ
a
1
+
μ
)
=
(
λ
+
μ
)
(
λ
+
σ
+
μ
)
(
λ
-
N
μ
a
1
+
μ
)
×
[
λ
+
(
δ
+
μ
)
(
1
-
K
0
e
-
λ
τ
)
]
=
0
,
where K
0
=
β
N
/
μ
(
δ
+
μ
).Clearly, (6) always has two negative eigenvalues λ
1
=
-
μ, λ
2
=
-
σ
-
μ and two indefinite eigenvalues λ
3
=
(
N
/
μ
)
a
1
-
μ and λ
4
=
-
(
δ
+
μ
)
(
1
-
K
0
e
-
λ
τ
). We let K
1
=
a
1
N
/
μ
2.It is clear to see that ifK
1
<
1 and K
0
<
1, the other two eigenvalues must be negative too. So, the following theorem can be acquired.Theorem 1.
For the system (2), if K
0
=
β
N
/
μ
(
δ
+
μ
)
<
1 and K
1
=
a
1
N
/
μ
2
<
1 are satisfied, then the virus-free equilibrium point E
1
=
(
N
/
μ
,
0,0
,
0
) is locally asymptotically stable. Besides, if K
1
>
1, the equilibrium point E
1 is unstable.Remark 2.
Theorem1 investigates the local stability of the virus-free equilibrium point by analyzing the eigenvalues of the corresponding characteristic equation. We can see that when viruses do not appear, all computers in the network are subjected to possible infection.


## 3.2. Endemic Equilibrium Points

Endemic equilibrium points are characterized by the existence of infected ones in the network; that is,I
≠
0. First, we consider the case when the network has no antidotal node; namely, A
=
0. Then, it is not difficult to solve (3) when A
=
0, I
≠
0 and the solution is
(7)
E
2
=
(
S
2
,
I
2
,
R
2
,
A
2
)
=
(
S
2
,
I
2
,
R
2
,
0
)
,
where
(8)
I
2
=
(
σ
+
μ
)
[
β
N
-
μ
(
δ
+
μ
)
]
(
σ
+
μ
)
(
δ
+
μ
)
(
β
+
μ
α
)
-
β
σ
δ
,
S
2
=
(
δ
+
μ
)
(
1
+
α
I
2
)
β
,
R
2
=
δ
I
2
σ
+
μ
.It indicates that whenK
0
<
1, I
2
<
0. Consequently, the condition K
0
<
1 can avoid the existence of the equilibrium point E
2.Another more important case isA
≠
0 when I
≠
0. Again, calculating (3) and the endemic point is given by
(9)
E
3
=
(
S
3
,
I
3
,
R
3
,
A
3
)
=
(
μ
-
a
2
I
3
a
1
,
I
3
,
δ
I
3
σ
+
μ
,
β
(
μ
-
a
2
I
3
)
a
1
a
2
(
1
+
α
I
3
)
-
δ
+
μ
a
2
)
.Bringing theE
3 point into the first equation of (3), this leads to
(10)
N
-
a
1
S
3
A
3
-
β
S
3
I
3
1
+
α
I
3
-
μ
S
3
+
σ
R
3
=
N
-
a
1
μ
-
a
2
I
3
a
1
(
β
(
μ
-
a
2
I
3
)
a
1
a
2
(
1
+
α
I
3
)
-
δ
+
μ
a
2
)
-
β
I
3
1
+
α
I
3
·
μ
-
a
2
I
3
a
1
-
μ
μ
-
a
2
I
3
a
1
+
σ
δ
I
3
σ
+
μ
=
0
.Then, we expand (10) and merge the similar terms; we let I
=
I
3 without confusion:
(11)
p
1
(
1
+
α
I
)
+
p
2
I
(
1
+
α
I
)
+
p
3
(
μ
-
a
2
I
)
+
p
4
(
μ
-
a
2
I
)
(
1
+
α
I
)
=
0
,
where
(12)
p
1
=
a
1
a
2
N
(
σ
+
μ
)
,
p
2
=
a
1
a
2
σ
δ
,
p
3
=
-
μ
β
(
σ
+
μ
)
,
p
4
=
(
σ
+
μ
)
(
μ
a
1
+
δ
a
1
-
μ
a
2
)
.Furthermore, a quadratic equation of the variableI is obtained:
(13)
b
1
I
2
+
b
2
I
+
b
3
=
0
,
where
(14)
b
1
=
p
2
α
-
p
4
α
a
2
,
b
2
=
p
1
α
+
p
2
-
p
3
a
2
+
p
4
α
μ
-
p
4
a
2
,
b
3
=
p
1
+
p
3
μ
+
p
4
μ
.From (13), we have the following result.Theorem 3.
Suppose thatS, I, R, A are all positive, and K
2
>
1, K
3
>
1. Then (1)
ifb
1
=
0 and b
3
/
b
2
<
0 hold, there exists unique positive solution (9) for (3), and I
3
=
-
b
3
/
b
2;
(2)
whenΔ
>
0 always holds, then if b
2
/
b
1
<
0 and b
3
/
b
1
>
0, there exist two positive solutions, where I
3
(
1
)
=
I
+ and I
3
(
2
)
=
I
-; if b
3
/
b
1
<
0, there is only one positive solution I
3
=
I
+ with b
1
>
0 or I
3
=
I
- with b
1
<
0; if b
3
=
0 and b
2
/
b
1
<
0, there is only one positive solution I
3
=
-
b
2
/
b
1;
(3)
ifΔ
=
0 and b
2
/
b
1
<
0, there is only one positive solution I
3
=
-
b
2
/
2
b
1,where K
2
=
(
β
(
μ
-
a
2
I
)
)
/
(
a
1
(
1
+
α
I
)
(
δ
+
μ
)
), K
3
=
μ
/
a
2
I, Δ
=
b
2
2
-
4
b
1
b
3, I
+
=
(
-
b
2
+
b
2
2
-
4
b
1
b
3
)
/
2
b
1, and I
-
=
(
-
b
2
-
b
2
2
-
4
b
1
b
3
)
/
2
b
1.Remark 4.
Theorem3 mainly focuses on the existence of the positive equilibrium point of the system. We can see that all conditions in the theorem are easily verified.
When the equilibrium pointE
3 exists, the characteristic equation of (2) at the point is(15)
det
⁡
(
λ
+
a
1
A
3
+
β
I
3
e
-
λ
τ
1
+
α
I
3
+
μ
β
S
3
e
-
λ
τ
(
1
+
α
I
3
)
2
-
σ
a
1
S
3
-
β
I
3
e
-
λ
τ
1
+
α
I
3
λ
-
β
S
3
e
-
λ
τ
(
1
+
α
I
3
)
2
+
a
2
A
3
+
δ
+
μ
0
a
2
I
3
0
-
δ
λ
+
σ
+
μ
0
-
a
1
A
3
-
a
2
A
3
0
λ
)
=
0
.We setc
1
=
β
I
3
/
(
1
+
α
I
3
), c
2
=
β
S
3
/
(
1
+
α
I
3
)
2, c
3
=
a
1
A
3
+
μ, and c
4
=
a
2
A
3
+
δ
+
μ. For notational simplicity, we use the (
S
,
I
,
R
,
A
) in place of E
3. Then, the following four-degree exponential polynomial equation is obtained:
(16)
λ
4
+
d
1
λ
3
+
d
2
λ
2
+
d
3
λ
+
d
4
+
e
-
λ
τ
(
d
5
λ
3
+
d
6
λ
2
+
d
7
λ
+
d
8
)
=
0
,
where
(17)
d
1
=
c
3
+
c
4
+
σ
+
μ
,
d
2
=
a
1
2
A
S
+
a
2
2
A
I
+
c
3
c
4
+
(
σ
+
μ
)
(
c
3
+
c
4
)
,
d
3
=
a
1
2
A
S
(
c
4
+
σ
+
μ
)
+
a
2
2
A
I
(
c
3
+
σ
+
μ
)
+
c
3
c
4
(
σ
+
μ
)
,
d
4
=
(
σ
+
μ
)
(
a
1
2
c
4
A
S
+
a
2
2
c
3
A
I
)
+
a
1
a
2
σ
δ
A
I
,
d
5
=
c
1
-
c
2
,
d
6
=
c
1
(
c
4
+
σ
+
μ
)
-
c
2
(
c
3
+
σ
+
μ
)
,
d
7
=
A
(
a
2
I
+
a
1
S
)
(
c
1
a
2
-
c
2
a
1
)
-
c
1
σ
δ
+
(
c
1
c
4
-
c
2
c
3
)
(
σ
+
μ
)
,
d
8
=
A
(
a
2
I
+
a
1
S
)
(
c
1
a
2
-
c
2
a
1
)
(
σ
+
μ
)
.Multiplyinge
λ
τ on both sides of (16), it is obvious to get
(18)
J
=
(
λ
4
+
d
1
λ
3
+
d
2
λ
2
+
d
3
λ
+
d
4
)
e
λ
τ
+
(
d
5
λ
3
+
d
6
λ
2
+
d
7
λ
+
d
8
)
=
0
.Letλ
=
i
ω
0, τ
=
τ
0, and substituting this into (18), for the sake of simplicity, denote ω
0 and τ
0 by ω, τ, respectively; then (18) becomes
(19)
(
ω
4
-
i
d
1
ω
3
-
d
2
ω
2
+
i
d
3
ω
+
d
4
)
(
cos
⁡
ω
τ
+
i
sin
ω
τ
)
+
(
-
i
d
5
ω
3
-
d
6
ω
2
+
i
d
7
ω
+
d
8
)
=
0
.Separating the real and imaginary parts, we have(20)
(
ω
4
-
d
2
ω
2
+
d
4
)
cos
⁡
ω
τ
+
(
d
1
ω
3
-
d
3
ω
)
sin
ω
τ
=
d
6
ω
2
-
d
8
,
(
-
d
1
ω
3
+
d
3
ω
)
cos
⁡
ω
τ
+
(
ω
4
-
d
2
ω
2
+
d
4
)
sin
ω
τ
=
d
5
ω
3
-
d
7
ω
.By simple calculation, the following equations are obtained:(21)
cos
⁡
ω
τ
=
e
5
ω
6
+
e
6
ω
4
+
e
7
ω
2
+
e
8
ω
8
+
e
1
ω
6
+
e
2
ω
4
+
e
3
ω
2
+
e
4
,
(22)
sin
ω
τ
=
e
9
ω
7
+
e
10
ω
5
+
e
11
ω
3
+
e
12
ω
ω
8
+
e
1
ω
6
+
e
2
ω
4
+
e
3
ω
2
+
e
4
,
where
(23)
e
1
=
d
1
2
-
2
d
2
,
e
2
=
-
2
d
1
d
3
+
d
2
2
+
2
d
4
,
e
3
=
-
2
d
2
d
4
+
d
3
2
,
e
4
=
d
4
2
,
e
5
=
d
6
-
d
1
d
5
,
e
6
=
d
1
d
7
-
d
2
d
6
+
d
3
d
5
-
d
8
,
e
7
=
d
2
d
8
-
d
3
d
7
+
d
4
d
6
,
e
8
=
-
d
4
d
8
,
e
9
=
d
5
,
e
10
=
d
1
d
6
-
d
2
d
5
-
d
7
,
e
11
=
-
d
1
d
8
+
d
2
d
7
-
d
3
d
6
+
d
4
d
5
,
e
12
=
d
3
d
8
-
d
4
d
7
.As is known to all thatsi
n
2
ω
τ
+
cos
⁡
2
ω
τ
=
1, we get
(24)
ω
16
+
f
7
ω
14
+
f
6
ω
12
+
f
5
ω
10
+
f
4
ω
8
+
f
3
ω
6
+
f
2
ω
4
+
f
1
ω
2
+
f
0
=
0
,
where
(25)
f
7
=
2
e
1
-
e
9
2
,
f
6
=
e
1
2
+
2
e
2
-
e
5
2
-
2
e
9
e
10
,
f
5
=
2
e
1
e
2
+
2
e
3
-
2
e
5
e
6
-
2
e
9
e
11
-
e
10
2
,
f
4
=
2
e
1
e
3
+
e
2
2
+
2
e
4
-
2
e
5
e
7
-
e
6
2
-
2
e
9
e
12
-
2
e
10
e
11
,
f
3
=
2
e
1
e
4
+
2
e
2
e
3
-
2
e
5
e
8
-
2
e
6
e
7
-
2
e
10
e
12
-
e
11
2
,
f
2
=
2
e
2
e
4
+
e
3
2
-
2
e
6
e
8
-
e
7
2
-
2
e
11
e
12
,
f
1
=
2
e
3
e
4
-
2
e
7
e
8
-
e
12
2
,
f
0
=
e
4
2
-
e
8
2
.
Denote z
=
ω
2; (24) can be rewritten as
(26)
z
8
+
f
7
z
7
+
f
6
z
6
+
f
5
z
5
+
f
4
z
4
+
f
3
z
3
+
f
2
z
2
+
f
1
z
+
f
0
=
0
.
We suppose that(H1)
(26) has at least one positive real root.Without loss of generality, we can assume the equation hasl (1
≤
l
≤
8) positive real roots, which are represented as z
i (1
≤
i
≤
l); then ω
i
=
z
i.By (21), we get
(27)
τ
k
j
=
1
ω
k
{
arccos
(
e
5
ω
k
6
+
e
6
ω
k
4
+
e
7
ω
k
2
+
e
8
ω
k
8
+
e
1
ω
k
6
+
e
2
ω
k
4
+
e
3
ω
k
2
+
e
4
)
+
2
j
π
}
,
i
i
i
i
i
i
i
i
i
i
i
i
i
i
i
i
i
i
i
i
i
i
i
i
i
i
i
i
i
i
i
i
i
i
i
i
i
i
i
i
i
i
i
i
i
i
i
i
i
i
i
i
i
i
i
i
i
i
i
i
i
i
i
i
i
i
i
j
=
0,1
,
…
.From the early discussions, we know that the±
i
ω
k are a pair of purely imaginary roots of (16) with τ
k
j. Define
(28)
τ
0
=
τ
k
0
0
=
min
⁡
k
∈
{
1
,
…
,
l
}
{
τ
k
0
}
,
ω
0
=
ω
k
0
.It is noted that whenτ
=
0, (16) becomes
(29)
λ
4
+
d
1
λ
3
+
d
2
λ
2
+
d
3
λ
+
d
4
+
(
d
5
λ
3
+
d
6
λ
2
+
d
7
λ
+
d
8
)
=
λ
4
+
(
d
1
+
d
5
)
λ
3
+
(
d
2
+
d
6
)
λ
2
+
(
d
3
+
d
7
)
λ
+
(
d
4
+
d
8
)
=
0
.By virtue of the well-known Routh-Hurwitz criteria, a set of necessary and sufficient conditions for all roots of (29) to have the negative real part is given in the following form:
(30)
D
1
=
d
1
+
d
5
>
0
,
(31)
D
2
=
|
d
1
+
d
5
d
3
+
d
7
1
d
2
+
d
6
|
=
(
d
1
+
d
5
)
(
d
2
+
d
6
)
-
(
d
3
+
d
7
)
>
0
,
(32)
D
3
=
|
d
1
+
d
5
d
3
+
d
7
0
1
d
2
+
d
6
d
4
+
d
8
0
d
1
+
d
5
d
3
+
d
7
|
=
(
d
1
+
d
5
)
[
(
d
2
+
d
6
)
(
d
3
+
d
7
)
-
(
d
1
+
d
5
)
(
d
4
+
d
8
)
]
-
(
d
3
+
d
7
)
2
>
0
,
(33)
D
4
=
|
d
1
+
d
5
d
3
+
d
7
0
0
1
d
2
+
d
6
d
4
+
d
8
0
0
d
1
+
d
5
d
3
+
d
7
0
0
1
d
2
+
d
6
d
4
+
d
8
|
=
(
d
4
+
d
8
)
D
3
>
0
.If (30)–(33) hold, (29) has four roots with negative real parts, and therefore when τ
=
0, system (2) is stable near the equilibrium point E
3.In order to give the main results, it is necessary to make the following assumption:(H2)
Re
(
d
λ
d
τ
)
|
τ
=
τ
0
≠
0.In order to calculate the derivative ofλ with respect to τ in (18), it is followed by
(34)
d
λ
d
τ
=
-
∂
J
/
∂
τ
∂
J
/
∂
λ
=
-
(
λ
(
λ
4
+
d
1
λ
3
+
d
2
λ
2
+
d
3
λ
+
d
4
)
e
λ
τ
)
×
(
(
4
λ
3
+
3
d
1
λ
2
+
2
d
2
λ
+
d
3
)
e
λ
τ
+
τ
e
λ
τ
×
(
λ
4
+
d
1
λ
3
+
d
2
λ
2
+
d
3
λ
+
d
4
)
+
(
3
d
5
λ
2
+
2
d
6
λ
+
d
7
)
(
4
λ
3
+
3
d
1
λ
2
+
2
d
2
λ
+
d
3
)
e
λ
τ
+
τ
e
λ
τ
)
-
1
;
thus
(35)
(
d
λ
d
τ
)
-
1
=
(
(
4
λ
3
+
3
d
1
λ
2
+
2
d
2
λ
+
d
3
)
e
λ
τ
+
(
3
d
5
λ
2
+
2
d
6
λ
+
d
7
)
)
×
(
-
λ
e
λ
τ
(
λ
4
+
d
1
λ
3
+
d
2
λ
2
+
d
3
λ
+
d
4
)
)
-
1
-
τ
λ
=
(
(
4
λ
3
+
3
d
1
λ
2
+
2
d
2
λ
+
d
3
)
e
λ
τ
+
(
3
d
5
λ
2
+
2
d
6
λ
+
d
7
)
)
×
(
d
5
λ
4
+
d
6
λ
3
+
d
7
λ
2
+
d
8
λ
)
-
1
-
τ
λ
.Whenτ
=
τ
0, ±
i
ω
0 are a pair of purely imaginary roots of (16). Substituting the τ
0, i
ω
0 into (35) and denoting τ
0 and ω
0 by τ, ω for simplicity, then
(36)
(
d
λ
d
τ
)
-
1
|
λ
=
i
ω
0
,
τ
=
τ
0
=
(
(
-
3
d
1
ω
2
+
d
3
)
cos
⁡
ω
τ
+
(
4
ω
3
-
2
d
2
ω
)
sin
ω
τ
-
3
d
5
ω
2
+
d
7
)
×
(
(
d
5
ω
4
-
d
7
ω
2
)
+
i
(
-
d
6
ω
3
+
d
8
ω
)
)
-
1
+
i
(
(
-
4
ω
3
+
2
d
2
ω
)
cos
⁡
ω
τ
+
(
-
3
d
1
ω
2
+
d
3
)
sin
ω
τ
+
2
d
6
ω
)
×
(
(
d
5
ω
4
-
d
7
ω
2
)
+
i
(
-
d
6
ω
3
+
d
8
ω
)
)
-
1
-
τ
i
ω
.We setQ
=
(
d
5
ω
4
-
d
7
ω
2
)
2
+
(
-
d
6
ω
3
+
d
8
ω
)
2,
(37)
Q
·
Re
(
d
λ
d
τ
)
-
1
|
τ
=
τ
0
=
[
(
-
3
d
1
ω
2
+
d
3
)
cos
⁡
ω
τ
+
(
4
ω
3
-
2
d
2
ω
)
sin
ω
τ
-
3
d
5
ω
2
+
d
7
]
×
(
d
5
ω
4
-
d
7
ω
2
)
+
[
(
-
4
ω
3
+
2
d
2
ω
)
cos
⁡
ω
τ
+
(
-
3
d
1
ω
2
+
d
3
)
sin
ω
τ
+
2
d
6
ω
]
×
(
-
d
6
ω
3
+
d
8
ω
)
,
sgn
⁡
{
Re
(
d
λ
d
τ
)
|
τ
=
τ
0
}
=
sgn
⁡
{
Q
·
Re
(
d
λ
d
τ
)
-
1
|
τ
=
τ
0
}
.On the basis of (H2), we can know that the roots of characteristic equation (16) cross the imaginary axis as τ continuously varies from a number less than τ
0 to one greater than τ
0 by Rouche’s theorem [17]. Therefore, the transversality condition holds and the conditions for Hopf bifurcation are then satisfied at τ
=
τ
0.Lemma 5 (see [18]).
Consider the exponential polynomial(38)
P
(
λ
,
e
-
λ
τ
1
,
…
,
e
-
λ
τ
m
)
=
λ
n
+
p
1
(
0
)
λ
n
-
1
+
⋯
+
p
n
-
1
(
0
)
λ
+
p
n
(
0
)
+
[
p
1
(
1
)
λ
n
-
1
+
⋯
+
p
n
-
1
(
1
)
λ
+
p
n
(
1
)
]
e
-
λ
τ
1
+
⋯
+
[
p
1
(
m
)
λ
n
-
1
+
⋯
+
p
n
-
1
(
m
)
λ
+
p
n
(
m
)
]
e
-
λ
τ
m
,
where τ
i
≥
0
(
i
=
1,2
,
…
,
m
) and p
j
(
i
)
(
i
=
0,1
,
…
,
m
;
j
=
1,2
,
…
,
n
) are constants. As (
τ
1
,
τ
2
,
…
,
τ
m
) vary, the sum of the order of the zeros of P
(
λ
,
e
-
λ
τ
1
,
…
,
e
-
λ
τ
m
) on the open right plane can change only if a zero appears on or crosses the imaginary axis.
From Lemma5, it is easy to get the theorem.Theorem 6.
Suppose that (H1) and (H2) hold; then (1)
for system (2), the equilibrium point E
3 is asymptotically stable for τ
∈
[
0
,
τ
0
);
(2)
system (2) undergoes a Hopf bifurcation at the E
3 when τ
=
τ
0. That is to say, it has a branch of periodic solutions bifurcating from the E
3 near τ
=
τ
0, where τ
0 can be calculated by (28).
Furthermore, we can investigate the stability and direction of bifurcating periodic solutions by analyzing higher order terms according to Hassard et al. [19].Remark 7.
The analyses above study the existence of the Hopf bifurcation and obtain the critical value of the parameterτ. We can apply the theorem into the system; thus, through changing some parameters, some bad performances of the model can be avoided.Remark 8.
When we setN
=
0, μ
=
0, α
=
0, and τ
=
0, the present model can be transformed to the model in [14]. However, [14] only calculated the disease-free/endemic equilibrium points and analyzed the stability. This paper expands the model and makes a detailed analysis about the bifurcations.Remark 9.
In this paper, we propose a new group called antidotal population, that is,A
(
t
), which is more reasonable when we study the computer viruses. It is well known that when dealing with Hopf bifurcation, the complexity of computation increases significantly as the dimension of system increases. Sometimes, it even cannot be calculated, especially when nonlinear terms β
S
(
t
-
τ
)
I
(
t
-
τ
)
/
(
1
+
α
I
(
t
-
τ
)
) exist in our system. Partly because of this, the majority of the literatures published use the traditional SIR three-dimensional model to study the epidemic model or virus [7, 10, 20], which has some limitations.


## 3.3. Stability and Direction of the Hopf Bifurcation

We have obtained the conditions under which a family of periodic solutions bifurcate from the positive equilibriumE
3 at the critical value of τ
0. In this subsection, the formulae for determining the direction of Hopf bifurcation and stability of bifurcating periodic solutions of system at τ
0 will be presented by employing the normal form theory and the center manifold reduction [19, 21–24].For convenience, letx
1
=
S
-
S
3, x
2
=
I
-
I
3, x
3
=
R
-
R
3, x
4
=
A
-
A
3, x
-
i
(
t
)
=
x
i
(
τ
t
), and τ
=
τ
0
+
u; we drop the bars for simplification of notations. In the light of multivariate Taylor expansion, system (2) can be transformed into an FDE in C
=
C
(
[
-
1,0
]
,
R
3
) as
(39)
x
˙
(
t
)
=
L
u
(
x
t
)
+
f
(
u
,
x
t
)
,
where x
(
t
)
=
(
x
1
(
t
)
,
x
2
(
t
)
,
x
3
(
t
)
,
x
4
(
t
)
)
T
∈
R
4, and L
u
:
C
→
R, f
:
R
×
C
→
R are given, respectively, by
(40)
L
u
(
ϕ
)
=
(
τ
0
+
u
)
(
-
c
3
0
σ
-
a
1
S
0
-
c
4
0
-
a
2
I
0
δ
-
σ
-
μ
0
a
1
A
a
2
A
0
0
)
(
ϕ
1
(
0
)
ϕ
2
(
0
)
ϕ
3
(
0
)
ϕ
4
(
0
)
)
+
(
τ
0
+
u
)
(
-
c
1
-
c
2
0
0
c
1
c
2
0
0
0
0
0
0
0
0
0
0
)
(
ϕ
1
(
-
1
)
ϕ
2
(
-
1
)
ϕ
3
(
-
1
)
ϕ
4
(
-
1
)
)
=
(
τ
0
+
u
)
B
1
ϕ
(
0
)
+
(
τ
0
+
u
)
B
2
ϕ
(
-
1
)
,
f
(
u
,
ϕ
)
=
(
τ
0
+
u
)
×
(
l
1
ϕ
2
2
(
-
1
)
-
l
2
ϕ
1
(
-
1
)
ϕ
2
(
-
1
)
-
a
1
ϕ
1
(
0
)
ϕ
4
(
0
)
+
h
.
o
.
t
-
l
1
ϕ
2
2
(
-
1
)
+
l
2
ϕ
1
(
-
1
)
ϕ
2
(
-
1
)
-
a
2
ϕ
2
(
0
)
ϕ
4
(
0
)
+
h
.
o
.
t
0
a
1
ϕ
1
(
0
)
ϕ
4
(
0
)
+
a
2
ϕ
2
(
0
)
ϕ
4
(
0
)
)
,
where l
1
=
α
c
2
/
(
1
+
α
I
), l
2
=
c
2
/
S, and L
u is a one-parameter family of bounded linear operators in C
[
-
1,0
]
→
R
4.By the Riesz representation theorem, there exists a functionη
(
θ
,
u
) of bounded variation for θ
∈
[
-
1,0
], such that
(41)
L
u
ϕ
=
∫
-
1
0
d
η
(
θ
,
u
)
ϕ
(
θ
)
for
ϕ
∈
C
1
(
[
-
1,0
]
,
R
4
)
.
In fact, we can choose
(42)
η
(
θ
,
u
)
=
(
τ
0
+
u
)
B
1
δ
(
θ
)
-
(
τ
0
+
u
)
B
2
δ
(
θ
+
1
)
,
where δ
(
θ
) is the Dirac delta function.Forϕ
∈
C
1
(
[
-
1,0
]
,
R
4
), define
(43)
A
(
u
)
ϕ
=
{
d
ϕ
(
θ
)
d
θ
,
-
1
≤
θ
<
0
,
∫
-
1
0
d
η
(
θ
,
u
)
ϕ
(
θ
)
,
θ
=
0
,
R
(
u
)
ϕ
=
{
0
,
θ
∈
[
-
1,0
)
,
f
(
u
,
ϕ
)
,
θ
=
0
.Then, system (39) is equivalent to
(44)
x
˙
t
=
A
(
u
)
x
t
+
R
(
u
)
x
t
,
where x
t
(
θ
)
=
x
(
t
+
θ
), for θ
∈
[
-
1,0
), and x
=
(
x
1
,
x
2
,
x
3
,
x
4
)
T.Forψ
∈
C
1
(
[
0,1
]
,
(
R
4
)
*
), define
(45)
A
*
(
u
)
ψ
=
{
-
d
ψ
(
s
)
d
s
,
0
<
s
≤
1
,
∫
-
1
0
d
η
T
(
t
,
u
)
ϕ
(
-
t
)
,
s
=
0
.In order to normalize the eigenvectors of operatorA and adjoint operator A
*, the following bilinear inner product is needed to introduce
(46)
〈
ψ
(
s
)
,
ϕ
(
θ
)
〉
=
ψ
-
(
0
)
ϕ
(
0
)
-
∫
-
1
0
∫
ξ
=
0
θ
ψ
-
(
ξ
-
θ
)
d
η
(
θ
)
ϕ
(
ξ
)
d
ξ
,
where η
(
θ
)
=
η
(
θ
,
0
).From the early discussions, we know that±
i
τ
0
ω
0 are eigenvalues of A
(
0
) and other eigenvalues have strictly negative real parts; thus, they are also eigenvalues of A
*. Define that
(47)
q
(
θ
)
=
(
1
,
q
2
,
q
3
,
q
4
)
T
e
i
τ
0
ω
0
θ
,
-
1
<
θ
≤
0
,
which is the eigenvector of A
(
0
) belonging to the eigenvalue i
τ
0
ω
0; namely, A
(
0
)
q
(
θ
)
=
i
τ
0
ω
0
q
(
θ
). Then, we can easily obtain
(48)
τ
0
(
i
ω
0
+
c
3
+
c
1
e
-
i
τ
0
ω
0
c
2
e
-
i
τ
0
ω
0
-
σ
a
1
S
-
c
1
e
-
i
τ
0
ω
0
i
ω
0
+
c
4
-
c
2
e
-
i
τ
0
ω
0
0
a
2
I
0
-
δ
i
ω
0
+
σ
+
μ
0
-
a
1
A
-
a
2
A
0
i
ω
0
)
×
q
(
0
)
=
(
0
0
0
0
)
.
Hence, we obtain
(49)
q
2
=
i
ω
0
c
1
e
-
i
τ
0
ω
0
-
a
1
a
2
I
A
i
ω
0
(
i
ω
0
+
c
4
-
c
2
e
-
i
τ
0
ω
0
)
+
a
2
2
I
A
,
q
3
=
δ
(
i
ω
0
c
1
e
-
i
τ
0
ω
0
-
a
1
a
2
I
A
)
(
i
ω
0
+
σ
+
μ
)
[
i
ω
0
(
i
ω
0
+
c
4
-
c
2
e
-
i
τ
0
ω
0
)
+
a
2
2
I
A
]
,
q
4
=
a
1
A
i
ω
0
+
a
2
A
(
i
ω
0
c
1
e
-
i
τ
0
ω
0
-
a
1
a
2
I
A
)
-
ω
0
2
(
i
ω
0
+
c
4
-
c
2
e
-
i
τ
0
ω
0
)
+
i
ω
0
a
2
2
I
A
.
Suppose that the eigenvector q
* of A
* belonging to the eigenvalue -
i
τ
0
ω
0 is
(50)
q
*
(
s
)
=
1
ρ
(
1
,
q
2
*
,
q
3
*
,
q
4
*
)
e
i
τ
0
ω
0
s
,
0
≤
s
<
1
.
Similar to the calculation of (49), we can get
(51)
q
2
*
=
i
ω
0
(
i
ω
0
-
c
3
-
c
1
e
i
τ
0
ω
0
)
+
a
1
2
A
S
-
i
ω
0
c
1
e
i
τ
0
ω
0
-
a
1
a
2
I
A
,
q
3
*
=
σ
-
i
ω
0
+
σ
+
μ
,
q
4
*
=
a
1
S
i
ω
0
+
a
2
I
i
ω
0
(
i
ω
0
-
c
3
-
c
1
e
i
τ
0
ω
0
)
+
a
1
2
a
2
I
A
S
ω
0
2
c
1
e
i
τ
0
ω
0
-
i
ω
0
a
1
a
2
I
A
.
Let 〈
q
*
,
q
〉
=
1; then
(52)
〈
q
*
(
s
)
,
q
(
θ
)
〉
=
1
ρ
-
(
1
+
q
2
q
2
*
+
q
3
q
3
*
+
q
4
q
4
*
)
-
∫
θ
=
-
1
0
∫
ξ
=
0
θ
1
ρ
-
(
1
,
q
-
2
*
,
q
-
3
*
,
q
-
4
*
)
e
-
i
τ
0
ω
0
(
ξ
-
θ
)
d
η
(
θ
)
×
(
1
q
2
q
3
q
4
)
e
i
τ
0
ω
0
ξ
d
ξ
=
1
ρ
-
(
1
+
q
2
q
2
*
+
q
3
q
3
*
+
q
4
q
4
*
)
-
∫
-
1
0
1
ρ
-
(
1
,
q
-
2
*
,
q
-
3
*
,
q
-
4
*
)
θ
e
i
τ
0
ω
0
θ
d
η
(
θ
)
(
1
q
2
q
3
q
4
)
=
1
ρ
-
{
τ
0
e
-
i
τ
0
ω
0
(
1
+
q
2
q
2
*
+
q
3
q
3
*
+
q
4
q
4
*
)
+
τ
0
e
-
i
τ
0
ω
0
[
(
q
-
2
*
-
1
)
(
c
1
+
α
c
2
)
]
}
,
which leads to
(53)
ρ
-
=
(
1
+
q
2
q
2
*
+
q
3
q
3
*
+
q
4
q
4
*
)
+
τ
0
e
-
i
τ
0
ω
0
[
(
q
-
2
*
-
1
)
(
c
1
+
α
c
2
)
]
.In the following, we apply the method in [19] to compute the coordinates describing the center manifold Ω
0 near u
=
0. Let x
t be the solution of (39) when u
=
0. We define
(54)
z
(
t
)
=
〈
q
*
,
x
t
〉
,
W
(
t
,
θ
)
=
x
t
(
θ
)
-
2
Re
[
z
(
t
)
q
(
θ
)
]
.On the center manifoldΩ
0, we have W
(
t
,
θ
)
=
W
(
z
(
t
)
,
z
-
(
t
)
,
θ
), where
(55)
W
(
z
(
t
)
,
z
-
(
t
)
,
θ
)
=
W
20
(
θ
)
z
2
2
+
W
11
(
θ
)
z
z
-
+
W
02
(
θ
)
z
-
2
2
+
W
30
(
θ
)
z
3
6
+
⋯
.In fact,z and z
- are local coordinates of center manifold Ω
0 in the direction of q and q
*, respectively. For solution x
t
∈
Ω
0 of (44), since u
=
0, we get
(56)
z
˙
(
t
)
=
〈
q
*
,
x
˙
t
〉
=
i
τ
0
ω
0
z
(
t
)
+
q
-
*
(
0
)
×
f
(
0
,
W
(
t
,
0
)
+
2
Re
[
z
(
t
)
q
(
0
)
]
)
=
i
τ
0
ω
0
z
(
t
)
+
q
-
*
(
0
)
f
0
(
z
,
z
-
)
≜
i
τ
0
ω
0
z
(
t
)
+
g
(
z
,
z
-
)
.
It is noted that
(57)
g
(
z
,
z
-
)
=
q
-
*
(
0
)
f
0
(
z
,
z
-
)
=
g
20
z
2
2
+
g
11
z
z
-
+
g
02
z
-
2
2
+
g
21
z
2
z
-
2
+
⋯
.By (54), we know that
(58)
x
t
(
θ
)
=
(
x
1
t
(
θ
)
,
x
2
t
(
θ
)
,
x
3
t
(
θ
)
,
x
4
t
(
θ
)
)
=
W
(
t
,
θ
)
+
z
q
(
θ
)
+
z
-
·
q
-
(
θ
)
.Considering (47) and (55), we have
(59)
x
1
t
(
0
)
=
z
+
z
-
+
W
20
(
1
)
(
0
)
z
2
2
+
W
11
(
1
)
(
0
)
z
z
-
+
W
02
(
1
)
(
0
)
z
-
2
2
+
o
(
|
(
z
,
z
-
)
|
3
)
,
x
2
t
(
0
)
=
q
2
z
+
q
-
2
z
-
+
W
20
(
2
)
(
0
)
z
2
2
+
W
11
(
2
)
(
0
)
z
z
-
+
W
02
(
2
)
(
0
)
z
-
2
2
+
o
(
|
(
z
,
z
-
)
|
3
)
,
x
4
t
(
0
)
=
q
4
z
+
q
-
4
z
-
+
W
20
(
4
)
(
0
)
z
2
2
+
W
11
(
4
)
(
0
)
z
z
-
+
W
02
(
4
)
(
0
)
z
-
2
2
+
o
(
|
(
z
,
z
-
)
|
3
)
,
x
1
t
(
-
1
)
=
z
e
-
i
τ
0
ω
0
+
z
-
e
i
τ
0
ω
0
+
W
20
(
1
)
(
-
1
)
z
2
2
+
W
11
(
1
)
(
-
1
)
z
z
-
+
W
02
(
1
)
(
-
1
)
z
-
2
2
+
o
(
|
(
z
,
z
-
)
|
3
)
,
x
2
t
(
-
1
)
=
z
q
2
e
-
i
τ
0
ω
0
+
z
-
·
q
-
2
e
i
τ
0
ω
0
+
W
20
(
2
)
(
-
1
)
z
2
2
+
W
11
(
2
)
(
-
1
)
z
z
-
+
W
02
(
2
)
(
-
1
)
z
-
2
2
+
o
(
|
(
z
,
z
-
)
|
3
)
.It follows that(60)
g
(
z
,
z
-
)
=
q
-
*
(
0
)
f
0
(
z
,
z
-
)
=
τ
0
ρ
-
(
1
,
q
-
2
*
,
q
-
3
*
,
q
-
4
*
)
×
(
l
1
x
2
t
2
(
-
1
)
-
l
2
x
1
t
(
-
1
)
x
2
t
(
-
1
)
-
a
1
x
1
t
(
0
)
x
4
t
(
0
)
+
h
.
o
.
t
-
l
1
x
2
t
2
(
-
1
)
+
l
2
x
1
t
(
-
1
)
x
2
t
(
-
1
)
-
a
2
x
2
t
(
0
)
x
4
t
(
0
)
+
h
.
o
.
t
0
a
1
x
1
t
(
0
)
x
4
t
(
0
)
+
a
2
x
2
t
(
0
)
x
4
t
(
0
)
)
=
τ
0
ρ
-
{
1
2
[
l
1
(
1
-
q
-
2
*
)
q
2
2
e
-
2
i
τ
0
ω
0
+
l
2
(
q
-
2
*
-
1
)
q
2
e
-
2
i
τ
0
ω
0
+
a
1
(
q
-
4
*
-
1
)
q
4
+
a
2
(
q
-
4
*
-
q
-
2
*
)
q
2
q
4
l
2
(
q
-
2
*
-
1
)
q
2
e
-
2
i
τ
0
ω
0
]
z
2
+
[
l
1
(
1
-
q
-
2
*
)
q
-
2
2
e
2
i
τ
0
ω
0
+
l
2
(
q
-
2
*
-
1
)
q
-
2
e
2
i
τ
0
ω
0
+
a
1
(
q
-
4
*
-
1
)
q
-
4
+
a
2
(
q
-
4
*
-
q
-
2
*
)
q
-
2
q
-
4
l
2
(
q
-
2
*
-
1
)
q
2
e
-
2
i
τ
0
ω
0
]
z
-
2
+
[
l
1
(
1
-
q
-
2
*
)
2
q
2
q
-
2
+
l
2
(
q
-
2
*
-
1
)
(
q
2
+
q
-
2
)
+
a
1
(
q
-
4
*
-
1
)
(
q
4
+
q
-
4
)
+
a
2
(
q
-
4
*
-
q
-
2
*
)
(
q
2
q
-
4
+
q
-
2
q
4
)
]
z
z
-
+
[
1
2
q
-
2
W
20
(
4
)
(
0
)
l
1
(
1
-
q
-
2
*
)
×
(
2
q
2
W
11
(
2
)
(
-
1
)
e
-
i
τ
0
ω
0
+
q
-
2
W
20
(
2
)
(
-
1
)
e
i
τ
0
ω
0
)
+
l
2
(
q
-
2
*
-
1
)
×
(
W
11
(
2
)
(
-
1
)
e
-
i
τ
0
ω
0
+
1
2
W
20
(
2
)
(
-
1
)
e
i
τ
0
ω
0
+
1
2
q
-
2
W
20
(
1
)
(
-
1
)
e
i
τ
0
ω
0
+
q
2
W
11
(
1
)
(
-
1
)
e
-
i
τ
0
ω
0
1
2
)
+
a
1
(
q
-
4
*
-
1
)
×
(
W
11
(
4
)
(
0
)
+
1
2
W
20
(
4
)
(
0
)
+
1
2
q
-
4
W
20
(
1
)
(
0
)
+
q
4
W
11
(
1
)
(
0
)
)
+
a
2
(
q
-
4
*
-
q
-
2
*
)
×
(
q
2
W
11
(
4
)
(
0
)
+
1
2
q
-
2
W
20
(
4
)
(
0
)
+
1
2
q
-
4
W
20
(
2
)
(
0
)
+
q
4
W
11
(
2
)
(
0
)
)
]
×
z
2
z
-
+
⋯
1
2
}
.
Comparing the coefficients in (57) with those in (60), it follows that
(61)
g
20
=
2
τ
0
ρ
-
[
l
1
(
1
-
q
-
2
*
)
q
2
2
e
-
2
i
τ
0
ω
0
+
l
2
(
q
-
2
*
-
1
)
q
2
e
-
2
i
τ
0
ω
0
+
a
1
(
q
-
4
*
-
1
)
q
4
+
a
2
(
q
-
4
*
-
q
-
2
*
)
q
2
q
4
q
2
2
e
-
2
i
τ
0
ω
0
]
,
g
02
=
2
τ
0
ρ
-
[
l
1
(
1
-
q
-
2
*
)
q
-
2
2
e
2
i
τ
0
ω
0
+
l
2
(
q
-
2
*
-
1
)
q
-
2
e
2
i
τ
0
ω
0
+
a
1
(
q
-
4
*
-
1
)
q
-
4
+
a
2
(
q
-
4
*
-
q
-
2
*
)
q
-
2
q
-
4
q
2
2
e
-
2
i
τ
0
ω
0
]
,
g
11
=
τ
0
ρ
-
[
l
1
(
1
-
q
-
2
*
)
2
q
2
q
-
2
+
l
2
(
q
-
2
*
-
1
)
(
q
2
+
q
-
2
)
+
a
1
(
q
-
4
*
-
1
)
(
q
4
+
q
-
4
)
+
a
2
(
q
-
4
*
-
q
-
2
*
)
(
q
2
q
-
4
+
q
-
2
q
4
)
]
,
g
21
=
2
τ
0
ρ
-
[
1
2
q
-
4
W
20
(
2
)
l
1
(
1
-
q
-
2
*
)
×
(
2
q
2
W
11
(
2
)
(
-
1
)
e
-
i
τ
0
ω
0
+
q
-
2
W
20
(
2
)
(
-
1
)
e
i
τ
0
ω
0
)
+
l
2
(
q
-
2
*
-
1
)
×
(
W
11
(
2
)
(
-
1
)
e
-
i
τ
0
ω
0
+
1
2
W
20
(
2
)
(
-
1
)
e
i
τ
0
ω
0
+
1
2
q
-
2
W
20
(
1
)
(
-
1
)
e
i
τ
0
ω
0
+
q
2
W
11
(
1
)
(
-
1
)
e
-
i
τ
0
ω
0
)
+
a
1
(
q
-
4
*
-
1
)
×
(
W
11
(
4
)
(
0
)
+
1
2
W
20
(
4
)
(
0
)
+
1
2
q
-
4
W
20
(
1
)
(
0
)
+
q
4
W
11
(
1
)
(
0
)
)
+
a
2
(
q
-
4
*
-
q
-
2
*
)
×
(
q
2
W
11
(
4
)
(
0
)
+
1
2
q
-
2
W
20
(
4
)
(
0
)
+
1
2
q
-
4
W
20
(
2
)
(
0
)
+
q
4
W
11
(
2
)
(
0
)
)
]
.Since theW
20
(
θ
) and W
11
(
θ
) exist in (61), we need to compute them. From (44) and (54), we have
(62)
W
˙
=
x
˙
t
-
z
˙
q
-
z
-
˙
·
q
-
=
A
x
t
+
R
x
t
-
[
i
τ
0
ω
0
z
+
q
-
*
(
0
)
f
0
(
z
,
z
-
)
]
q
-
[
-
i
τ
0
ω
0
z
-
+
q
*
(
0
)
f
-
0
(
z
,
z
-
)
]
q
-
=
A
(
W
+
2
Re
(
z
q
)
)
+
R
x
t
-
2
Re
[
q
-
*
(
0
)
f
0
(
z
,
z
-
)
q
]
-
2
Re
[
i
τ
0
ω
0
z
q
]
=
{
A
W
-
2
Re
[
q
-
*
(
0
)
f
0
(
z
,
z
-
)
q
(
θ
)
]
,
-
1
≤
θ
<
0
,
A
W
-
2
Re
[
q
-
*
(
0
)
f
0
(
z
,
z
-
)
q
(
θ
)
]
+
f
,
θ
=
0
,
≜
A
W
+
H
(
z
,
z
-
,
θ
)
,
where
(63)
H
(
z
,
z
-
,
θ
)
=
H
20
(
θ
)
z
2
2
+
H
11
(
θ
)
z
z
-
+
H
02
(
θ
)
z
-
2
2
+
⋯
.
Substituting (55) and (63) into (62), we get
(64)
W
˙
=
(
A
W
20
+
H
20
)
z
2
2
+
(
A
W
11
+
H
11
)
z
z
-
+
(
A
W
02
+
H
02
)
z
-
2
2
+
⋯
.
Taking the derivate of W with respect to t in (55), we have
(65)
W
˙
=
W
z
z
˙
+
W
z
-
z
-
˙
=
(
W
20
z
+
W
11
z
-
+
W
30
z
2
2
+
⋯
)
(
i
τ
0
ω
0
z
+
g
)
+
(
W
11
z
+
W
02
z
-
+
⋯
)
(
-
i
τ
0
ω
0
z
-
+
g
-
)
=
i
τ
0
ω
0
W
20
z
2
+
z
z
-
(
i
τ
0
ω
0
W
11
-
i
τ
0
ω
0
W
11
)
-
i
τ
0
ω
0
W
02
z
-
2
+
⋯
.
Then, together with the two above equations, we obtain
(66)
A
W
20
+
H
20
=
2
i
τ
0
ω
0
W
20
,
A
W
11
+
H
11
=
0
.
By (62), we know that
(67)
H
(
z
,
z
-
,
θ
)
=
-
g
q
(
θ
)
-
g
-
·
q
-
(
θ
)
+
R
x
t
=
-
(
g
20
z
2
2
+
g
11
z
z
-
+
g
02
z
-
2
2
+
⋯
)
q
(
θ
)
-
(
g
-
20
z
-
2
2
+
g
-
11
z
z
-
+
g
-
02
z
2
2
+
⋯
)
q
-
(
θ
)
+
R
x
t
.
We know that when -
1
≤
θ
<
0  R
x
t
=
0, comparing the coefficients of the z
2, z
z
- in (63) and (67), respectively, we can find
(68)
H
20
(
θ
)
=
-
g
20
q
(
θ
)
-
g
-
02
q
-
(
θ
)
,
H
11
(
θ
)
=
-
g
11
q
(
θ
)
-
g
-
11
q
-
(
θ
)
.
Combining (66) with (68), it follows that
(69)
W
˙
20
(
θ
)
=
2
i
τ
0
ω
0
W
20
(
θ
)
+
g
20
q
(
θ
)
+
g
-
02
q
-
(
θ
)
,
W
˙
11
(
θ
)
=
g
11
q
(
θ
)
+
g
-
11
q
-
(
θ
)
,
iiiiiiiiiiiiiiiiiiiiiiiiii
-
1
≤
θ
<
0
.
We solve (69) and get the solutions in the following form:
(70)
W
20
(
θ
)
=
i
g
20
τ
0
ω
0
q
(
0
)
e
i
τ
0
ω
0
θ
-
g
-
02
3
i
τ
0
ω
0
q
-
(
0
)
e
-
i
τ
0
ω
0
θ
+
E
1
e
2
i
τ
0
ω
0
θ
,
W
11
(
θ
)
=
g
11
i
τ
0
ω
0
q
(
0
)
e
i
τ
0
ω
0
θ
-
g
-
11
i
τ
0
ω
0
q
-
(
0
)
e
-
i
τ
0
ω
0
θ
+
E
2
.
Next, we consider (66) again when θ
=
0; we can see that
(71)
A
W
20
(
0
)
=
∫
-
1
0
d
η
(
θ
)
W
20
(
θ
)
=
2
i
τ
0
ω
0
W
20
(
0
)
-
H
20
(
0
)
,
A
W
11
(
0
)
=
∫
-
1
0
d
η
(
θ
)
W
11
(
θ
)
=
-
H
11
(
0
)
.
Furthermore, from (62) and (63), we have
(72)
H
20
(
0
)
=
-
g
20
q
(
0
)
-
g
-
02
q
-
(
0
)
+
2
τ
0
(
l
1
q
2
2
e
-
2
i
τ
0
ω
0
-
l
2
q
2
e
-
2
i
τ
0
ω
0
-
a
1
q
4
-
l
1
q
2
2
e
-
2
i
τ
0
ω
0
+
l
2
q
2
e
-
2
i
τ
0
ω
0
-
a
2
q
2
q
4
0
a
1
q
4
+
a
2
q
2
q
4
)
,
(73)
H
11
(
0
)
=
-
g
11
q
(
0
)
-
g
-
11
q
-
(
0
)
+
τ
0
(
2
l
1
|
q
2
|
2
-
l
2
(
q
2
+
q
-
2
)
-
a
1
(
q
4
+
q
-
4
)
-
2
l
1
|
q
2
|
2
+
l
2
(
q
2
+
q
-
2
)
-
a
2
(
q
2
q
-
4
+
q
-
2
q
4
)
0
a
1
(
q
4
+
q
-
4
)
+
a
2
(
q
2
q
-
4
+
q
-
2
q
4
)
)
.Substitute (70) and (72) into (71) and notice that
(74)
(
i
τ
0
ω
0
I
-
∫
-
1
0
e
i
τ
0
ω
0
θ
d
η
(
θ
)
)
q
(
0
)
=
0
,
(
-
i
τ
0
ω
0
I
-
∫
-
1
0
e
-
i
τ
0
ω
0
θ
d
η
(
θ
)
)
q
-
(
0
)
=
0
.
It is easy to yield
(75)
∫
-
1
0
d
η
(
θ
)
E
1
e
2
i
τ
0
ω
0
θ
=
2
i
τ
0
ω
0
E
1
-
2
τ
0
(
l
1
q
2
2
e
-
2
i
τ
0
ω
0
-
l
2
q
2
e
-
2
i
τ
0
ω
0
-
a
1
q
4
-
l
1
q
2
2
e
-
2
i
τ
0
ω
0
+
l
2
q
2
e
-
2
i
τ
0
ω
0
-
a
2
q
2
q
4
0
a
1
q
4
+
a
2
q
2
q
4
)
,
which leads to(76)
E
1
=
2
(
2
i
ω
0
+
c
3
+
c
1
e
-
2
i
τ
0
ω
0
c
2
e
-
2
i
τ
0
ω
0
-
σ
a
1
S
-
c
1
e
-
2
i
τ
0
ω
0
2
i
ω
0
+
c
4
-
c
2
e
-
2
i
τ
0
ω
0
0
a
2
I
0
-
δ
2
i
ω
0
+
σ
+
μ
0
-
a
1
A
-
a
2
A
0
2
i
ω
0
)
-
1
·
(
l
1
q
2
2
e
-
2
i
τ
0
ω
0
-
l
2
q
2
e
-
2
i
τ
0
ω
0
-
a
1
q
4
-
l
1
q
2
2
e
-
2
i
τ
0
ω
0
+
l
2
q
2
e
-
2
i
τ
0
ω
0
-
a
2
q
2
q
4
0
a
1
q
4
+
a
2
q
2
q
4
)
≜
2
M
1
-
1
M
2
.By the similar way we can get theE
2 as in the following:
(77)
E
2
=
(
c
3
+
c
1
c
2
-
σ
a
1
S
-
c
1
c
4
-
c
2
0
a
2
I
0
-
δ
σ
+
μ
0
-
a
1
A
-
a
2
A
0
0
)
-
1
×
(
2
l
1
|
q
2
|
2
-
l
2
(
q
2
+
q
-
2
)
-
a
1
(
q
4
+
q
-
4
)
-
2
l
1
|
q
2
|
2
+
l
2
(
q
2
+
q
-
2
)
-
a
2
(
q
2
q
-
4
+
q
-
2
q
4
)
0
a
1
(
q
4
+
q
-
4
)
+
a
2
(
q
2
q
-
4
+
q
-
2
q
4
)
)
.Thus, we can calculate theW
20
(
θ
) and W
11
(
θ
) from (70). Consequently, g
21 in (61) can be received. Then, we need to compute the following parameters [19]:
(78)
C
1
(
0
)
=
i
2
τ
0
ω
0
(
g
20
g
11
-
2
|
g
11
|
2
-
|
g
02
|
2
3
)
+
g
21
2
,
(79)
μ
2
=
-
Re
{
C
1
(
0
)
}
Re
{
λ
′
(
τ
0
)
}
,
β
2
=
2
Re
{
C
1
(
0
)
}
,
T
2
=
-
Im
⁡
{
C
1
(
0
)
}
+
μ
2
Im
⁡
{
λ
′
(
τ
0
)
}
τ
0
ω
0
,
which determine the characteristic of bifurcating periodic solutions in the center manifold at the critical value τ
0. More details are given in the following theorem.Theorem 10.
Under the conditions of Theorem6 one has the following.(1)
u
=
0 is Hopf bifurcation value of system (39).
(2)
The direction of Hopf bifurcation is determined by the sign ofμ
2: if μ
2
>
0, the Hopf bifurcation is supercritical; if μ
2
<
0, the Hopf bifurcation is subcritical.
(3)
The stability of bifurcating periodic solutions is determined byβ
2: if β
2
<
0, the periodic solutions are stable; if β
2
>
0, they are unstable.
(4)
The sign ofT
2 determines the period of the bifurcating periodic solutions: if T
2
>
0, the period increases; if T
2
<
0, the period decreases.


## 4.  Numerical Example

In this section, we will present some numerical simulations for verifying our theoretical analysis. Our example involves 9 parameters, including the delayτ.Case 1.
Considera
1
=
0.025, a
2
=
0.25, σ
=
0.8, α
=
0.5, β
=
0.5, δ
=
0.4, μ
=
0.2, and N
=
0.2. In this example, according to Theorems 1 and 3, the system in Case 1 only has one reasonable equilibrium point [
1,0
,
0,0
] and this equilibrium is locally asymptotically stable with arbitrarily τ, because K
0
=
0.8333
<
1, K
1
=
0.1250
<
1, K
2
=
0.6009
<
1, and K
3
=
0.5621
<
1. When we let τ
=
0.5, with the randomly chosen initial value [
0.6145
0.5077
1.6924
0.5913
], the time response curves are shown in Figures 1 and 2.Figure 1
All components of the system converge to the equilibrium.Figure 2
The first component of the system converges to 1.Case 2.
Considera
1
=
0.025, a
2
=
0.25, σ
=
0.8, α
=
1, β
=
4, δ
=
0.3, μ
=
0.1, and N
=
0.2. We can calculate that K
2
=
1.7838
>
1, K
3
=
1.0653
>
1, D
1
=
2.3180
>
0, D
2
=
4.3667
>
0, D
3
=
2.1135
>
0, D
4
=
0.0711
>
0, and τ
0
=
2.3581; we can see that all conditions of Theorem 6 are satisfied; thus, it is easy to obtain the following results.
Whenτ
∈
[
0
,
τ
0
), there is only one positive equilibrium point [
0.2454
0.3755
0.1252
1.2540
], which is asymptotically stable. In simulation, we choose [
0.255
0.38
0.11
1.24
] as initial values; when τ
=
2
<
τ
0, the simulation results are Figures 3 and 4; the system in Case 2 undergoes a Hopf bifurcation at the equilibrium point when τ
=
2.3581
=
τ
0, as shown in Figures 5 and 6.
Furthermore, when adopting Theorem10, it is easy to acquire more details about the bifurcating periodic solutions. By (79), we can compute that C
1
(
0
)
=
-
4.1125
-
4.2861
i, μ
2
=
43.8461
>
0, β
2
=
-
8.2251
<
0, and T
2
=
7.8990
>
0. Hence, by Theorem 10, we know that the bifurcating point is supercritical, the periodic solutions are stable, and the period increases.Figure 3
The first component of the system converges to 0.2454.Figure 4
Other three components of the system converge to[
0.3755
0.1252
1.2540
].Figure 5
The first component oscillates as a periodic solution.Figure 6
The system undergoes a Hopf bifurcation.


## 5. Conclusions

In this paper, a modified SIRA model with time delay and nonlinear terms has been proposed, and sufficient conditions on the existence of the virus-free and endemic equilibrium points have been derived. We also obtained several results guaranteeing the stability of the equilibrium and the occurrence of the Hopf bifurcation at the critical value. By using normal form and center manifold theory, the explicit formulae which determine the stability, direction, and other properties of bifurcating periodic solutions have been established. Numerical simulations have been presented to verify the accuracy of our results. The present model can be extended to formulate the more general network for computer virus spread.



---

*Source: 784684-2014-06-05.xml*