# Fractional Dynamics of Computer Virus Propagation

Authors: Carla M. A. Pinto; J. A. Tenreiro Machado
Journal: Mathematical Problems in Engineering
 (2014)

---

## Abstract

We propose a fractional model for computer virus propagation. The model includes the interaction between computers and removable devices. We simulate numerically the model for distinct values of the order of the fractional derivative and for two sets of initial conditions adopted in the literature. We conclude that fractional order systems reveal richer dynamics than the classical integer order counterpart. Therefore, fractional dynamics leads to time responses with super-fast transients and super-slow evolutions towards the steady-state, effects not easily captured by the integer order models.


---

## Body

## 1. Introduction

The global information world in what we live has brought numerous advantages into our lives. Daily commodities like making traveling arrangements, checking bank accounts, buying online books, groceries, and souvenirs, and other important tasks in our lives, such as working, are readily available at the click of a button. Nevertheless, this easy way into the internet has also some disadvantages, namely, in what concerns computer virus propagation. Computer viruses spread and cause huge money losses to companies and customers. Annually, millions of dollars are lost due to virus infections [1].During the last few years, authors have applied mathematical models for epidemics to computer virus propagation. Computer virus and biological virus have very similar behavioral patterns [2–4]. Based on this observation, some biological epidemic models, such as the susceptible-infected-susceptible (SI/SIS) or the susceptible-infected-recovered (SIR) models, were applied in modeling the epidemic behavior of computer virus patterns [5–12].In [5, 8], authors study the propagation of viruses and worms in distinct network topologies. In 2005, Zou et al. [6] propose an Internet worm monitoring system that detects a worm in its early propagation stage, using Kalman filter estimation. Moreover, the model can also predict the vulnerable susceptible computer population and estimate the infected computers for uniform-scan worms, such as Code Red. Feng et al. [9] propose a modified SIRS model for computer virus propagation, with dual delays and multistate antivirus measures. Some of the measures may be cleaning, patching, and filtering. Authors use the central manifold theorem and normal form theory to establish explicit formulas for determining the stability and direction of periodic solutions produced by the model. The appearance of periodic solutions is an aggressive state in a computer network, since it means the virus prevalence will not be constant (i.e., will not be easier to combat). A way to control this phenomenon is to combine bifurcation control and reduce virus prevalence strategies, such as using distinct topological network structures. Mishra and Jha [7] develop an SIRS model that includes temporary immunity. The latter is observed when an antivirus software is run in a computer network, after a node gets affected by a malicious object. The authors observe that the endemic equilibrium of the model might be unstable, leading to new outcomes of an SIRS model. Zhu et al. [10] use optimal control methods to fight computer virus propagation. They consider a controlled delayed model and then apply an optimal control strategy, assuming a tradeoff between the control costs and the effects. Ren et al. [11] study a computer virus propagation model, where the effect of antivirus ability is considered. The ability of an antivirus software can be measured by the number of computers recovered, from infected computers, per unit time, due to running the antivirus software. The ability of an antivirus software is usually proportional to its cost. Backward and Hopf bifurcations can be found in the model, for variation of the antivirus ability. Authors also show the existence of bistable states. Zhu et al. [12] propose a model for computer virus propagation, where the interaction between computers and external removable devices is taken into account. Authors assume reasonable assumptions for this interaction and prove that the state without virus is globally asymptotically stable for a value of the reproduction number, R
0
<
1. They also prove that the epidemic state is globally asymptotically stable for values of R
0
>
1. Thus, a good strategy to control computer virus transmission is to develop measures such that R
0
<
1.In this paper, we analyze the fractional order version of the integer order model proposed by Zhu et al. [12], for computer and removable devices virus propagation. We simulate numerically the model for different values of the order of the fractional derivative α. In this line of thought, the paper is organized as follows. In Section 3, we describe the model proposed for computer viruses propagation. In Section 4, we analyze several simulations of the model, for distinct values of the fractional derivative, and discuss implications of the results. In the last section, we present the main conclusions and outline some future research topics.


## 2. Fractional Calculus: A Review

The generalization of the derivative operatorD
α
f
(
x
) to fractional values of α, the order of the derivative, started with the theory of differential calculus, namely, when Leibniz wrote about D
1
/
2
f
(
x
). The development of the fractional calculus (FC) is due to many contributions of mathematicians such as Euler, Liouville, Riemann, and Letnikov [13–16]. In the fields of physics and engineering, FC is presently associated with long term memory effects [17–21]. This research is still giving its first steps and new areas of application of FC, such as the modelling of dynamical systems, are emerging [22].Several definitions of fractional derivatives were proposed. The most used definitions of a fractional derivative of orderα are the Riemann-Liouville (RL), Grünwald-Letnikov (GL), and Caputo (C) formulations. GL is defined as
(1)
a
GL
D
t
α
f
(
t
)
=
lim
⁡
h
→
0
1
h
α
∑
k
=
0
[
(
t
-
a
)
/
h
]
‍
(
-
1
)
k
(
α
k
)
f
(
t
-
k
h
)
,
t
>
a
,
α
>
0
,
where Γ
(
·
) is Euler’s gamma function, [
x
] means the integer part of x, and h is the step time increment.These expressions capture the history of the past dynamics, contrary to the integer counterpart that is a “local” operator. This property was recognized in several phenomena and their modelling becomes easier using the FC formalism, while integer order models are often much more complicated.The GL definition inspired a discrete-time calculation algorithm, based on the approximation of the time incrementh by means of the sampling period T, yielding the equation in the z domain:
(2)
Z
{
D
α
f
(
t
)
}
Z
{
f
(
t
)
}
=
1
T
α
∑
k
=
0
∞
‍
(
-
1
)
k
Γ
(
α
+
1
)
k
!
Γ
(
α
-
k
+
1
)
z
-
k
=
(
1
-
z
-
1
T
)
α
,
where Z denotes the Z-transform operator.One implementation of (2) is accomplished by means of a r-term truncated series:
(3)
Z
{
D
α
f
(
t
)
}
Z
{
f
(
t
)
}
=
1
T
α
∑
k
=
0
r
‍
(
-
1
)
k
Γ
(
α
+
1
)
k
!
Γ
(
α
-
k
+
1
)
z
-
k
,
where, in order to have good approximations, a large r and a small value of T are required.Expression (3) represents the Euler, or first backward difference, approximation in the so-called s
→
z conversion scheme. Another possibility consists in the Tustin conversion rule. The Euler and Tustin rational expressions, ψ
0
(
z
-
1
)
=
(
(
1
-
z
-
1
)
/
T
) and ψ
1
(
z
-
1
)
=
(
2
/
T
)
(
(
1
-
z
-
1
)
/
(
1
+
z
-
1
)
), are called generating approximants of zero and first order, respectively. Therefore, the generalization of these conversion methods leads to the noninteger order α results [23]:
(4)
s
α
≈
(
1
-
z
-
1
T
)
α
,
s
α
≈
(
2
T
1
-
z
-
1
1
+
z
-
1
)
α
.The expressionψ
0
α
(
z
-
1
)
=
[
ψ
0
(
z
-
1
)
]
α and ψ
1
α
(
z
-
1
)
=
[
ψ
1
(
z
-
1
)
]
α weighted by the factors p and 1
-
p generate a family of fractional differentiator:
(5)
ψ
a
v
α
(
z
-
1
)
=
p
ψ
0
α
(
z
-
1
)
+
(
1
-
p
)
ψ
1
α
(
z
-
1
)
.In order to get a rational expression, the final approximation corresponds to a truncated Taylor series or a rational fraction expansion. The arithmetic mean (5) motivates the study of an averaging method [22] based on the generalized formula of averages (often called average of order q
∈
R):
(6)
ψ
a
v
α
(
z
-
1
)
=
{
p
[
ψ
0
α
(
z
-
1
)
]
q
+
(
1
-
p
)
[
ψ
1
α
(
z
-
1
)
]
q
}
1
/
q
,
where (
p
,
q
) are two tuning degrees of freedom, corresponding q to the order of the averaging expression and p to the weighting factor. For example, when q
=
{
-
1,0
,
1
}, in expression (6), we get the well-known expressions for the {harmonic, geometric, arithmetic} averages.


## 3. The Model

The fractional model considered here for computer viruses propagation consists of one SIR (susceptible-infected-recovered) and one SI (susceptible-infected) coupled models. It is derived as follows. We consider that the entire computer and removable devices population is divided into five classes, the susceptible computers,S, the infected computers, I, the immune (protected with effective antivirus) computers, R, the susceptible removable devices, R
S, and the infected removable devices, R
I. The first system SIR is for computer virus modeling and the second system SI is for removable devices viruses propagation.Computers are connected to the network at a rateλ
1. Susceptible computers can be infected by other computers, at a rate β
1
S
I, or by a removable device, at a rate β
2
(
R
I
/
R
N
)
S, and, after that, move to the infected class, I. Parameters β
1 and β
2 are the contact infective forces between susceptible and infective computers and between computers and removable devices, respectively. Infected computers can be repaired at a rate σ
1
I and move to the immune class, R, where σ
1 is the recovery rate of infected computers, due to antivirus software. Every computer can be disconnected from the network, at a rate μ
1, in all classes.The removable devices are recruited, at a rateλ
2, to the susceptible class, R
S. They are infected by the computers at a rate β
2
(
I
/
N
)
R
S and then move to the infected removable class, R
I. These infected devices may recover at a rate σ
2
(
R
/
N
)
R
I and after they move back to the susceptible class, R
S. Parameter σ
2 is the recovery rate of removable devices, due to antivirus software. All removable devices may break down at a rate μ
2 in all classes.The model assumes that the antivirus software provides protection of computers, for all times. This is not a realistic assumption, but it simplifies the model at this stage. Future work will consider adding temporary immunity to computers; that is, the model will assume, more realistically, that antivirus software must be updated (i.e., immunity is temporary).In Figure1, we depict the schematic diagram of the model.Figure 1
Schematic diagram of the fractional order model (7) for computer viruses propagation [12].The system of fractional nonlinear ordinary differential equations, for the proposed model, is given by(7)
d
S
α
d
t
=
λ
1
-
β
1
S
I
-
β
2
S
R
I
R
N
-
μ
1
S
,
d
I
α
d
t
=
β
1
S
I
+
β
2
S
R
I
R
N
-
(
μ
1
+
σ
1
)
I
,
d
R
α
d
t
=
σ
1
I
-
μ
1
R
,
d
R
I
α
d
t
=
λ
2
-
β
2
R
S
I
N
+
σ
2
R
I
R
N
-
μ
2
R
S
,
d
R
S
α
d
t
=
β
2
R
S
I
N
-
σ
2
R
I
R
N
-
μ
2
R
I
,
where α
∈
]
0,1
] is the order of the fractional derivative.The derivation of model (7) including fractional orders, in a pure mathematical way, embodies some criticism. In fact, it derives from a classical integer order expression and assumes, by including a new degree of freedom, that we can have a better fit between real-world data and theoretical formulation. Furthermore, possibly, the definition of different fractional orders can lead to slight better models but at cost of an increased complexity. Therefore, the proposed methodology follows recent studies [24–26] where fractional models constitute merely an abstract conjecture that needs to be validated in the future, in the presence of real data. In what concerns the meaning of α, while the abstract interpretation can be related solely to “anomalous” or “super-slow” dynamics, its clinical significance remains to be investigated. Possible research directions point to DNA relationship [27, 28].The basic reproduction number,R
0, of the integer order model (α
=
1) is computed in [12] to be
(8)
R
0
=
β
2
2
+
μ
2
β
1
(
λ
1
/
μ
1
)
μ
2
(
μ
1
+
σ
1
)
.
The basic reproduction number, R
0, is defined as the number of secondary infections due to a single infection in a completely susceptible population. For R
0
<
1 the disease-free equilibrium is globally asymptotically stable. If R
0
>
1, the endemic equilibrium is globally asymptotically stable [12].


## 4. Numerical Results

In this section we develop several numerical simulations of the new fractional order model (7). The dynamical behavior of the model is studied numerically for variation of the noninteger order derivative α. Initial conditions used in the simulations are found in Table 1 and parameter values in Table 2. For integrating the differential equations of the type D
α
x
i
=
φ
i
(
x
1
,
…
,
x
n
) a standard integration of order 1 followed by a differentiation of order 1
-
α, approximated by means of (3), was implemented. Furthermore, the numerical calculations adopt the values r
=
1
0
3 and T
=
2.5
·
1
0
-
4 year.Table 1
Initial conditions used in the numerical simulations of model (7).
Variable
S
I
R
R
S
R
I
Case  1
0.0
4.0
3.0
0.1
0.5
Case  2
5.0
1.0
0.0
0.5
0.1Table 2
Parameters used in the numerical simulations of model (7).
Parameter
Case  1
Case  2
λ
1
1
1
λ
2
0.1
0.1
β
1
0.01
0.035
β
2
0.01
0.035
σ
1
0.02
0.02
σ
2
0.005
0.005
μ
1
0.1
0.1
μ
2
0.1
0.1In Figure2, we plot the dynamics of the number of susceptible computers, S, versus time, t, of model (7), for different values of order of the fractional derivative; namely, α
∈
{
0.1
,
…
,
1
}. Initial conditions and parameter values are those of Case 1. For these parameter values, the integer order system has reproduction number R
0
≃
0.84
<
1 [12].Figure 2
Dynamics of the susceptible computersS versus time t, of system (7), for α
∈
{
0
.
1
,
…
,
1
}. Parameter values and initial conditions are those of Case 1.In Figures3 and 4 we plot, respectively, the dynamics of the number of infected computers, I, and immune computers, R, versus time t, of model (7), for α
∈
{
0.1
,
…
,
1
}. Initial conditions and parameters are those of Case 1.Figure 3
Dynamics of the infected computersI versus time t, of system (7), for α
∈
{
0.1
,
…
,
1
}. Parameter values and initial conditions are those of Case 1.Figure 4
Dynamics of the immune computersR versus time t, of system (7), for α
∈
{
0.1
,
…
,
1
}. Parameter values and initial conditions are those of Case 1.In Figures5 and 6, we plot, respectively, the values of the susceptible removable devices, R
S, and of the infected removable devices, R
I, versus time, t, respectively, for distinct values of the fractional derivative α. Initial conditions are those of Case 1.Figure 5
Dynamics of the susceptible removable devicesR
S versus time t, of system (7), for α
∈
{
0.1
,
…
,
1
}. Parameter values and initial conditions are those of Case 1.Figure 6
Dynamics of the infected removable devicesR
I versus time t, of system (7), for α
∈
{
0.1
,
…
,
1
}. Parameter values and initial conditions are those of Case 1.From Figures2–6, we observe that by varying α, we get a smooth variation of the dynamics. It is important to note that we can vary not only the “velocity” of the initial transient but also the evolution towards the steady-state. This effect is typical in fractional dynamical systems, where we can have a super-fast initial transient, followed by a super-slow convergence to the final value. Moreover, we verify also that for α
=
0.1 we have the slowest evolution, while α
=
0.7 yields the fastest transient.In Figures7–11, we consider initial conditions and parameter values as in Case 2 and vary the fractional derivative α. We increase the values of the contact infective forces to be β
1
=
β
2
=
0.035. For these parameter values, the integer order system has reproduction number R
0
≃
3.02
>
1 [12].Figure 7
Dynamics of the susceptible computersS versus time t, of system (7), α
∈
{
0.1
,
…
,
1
}. Parameter values and initial conditions are as in Case 2, except for β
1
=
β
2
=
0.035.Figure 8
Dynamics of the infected computersI versus time t, of system (7), for α
∈
{
0.1
,
…
,
1
}. Parameter values and initial conditions are as in Case 2, except for β
1
=
β
2
=
0.035.Figure 9
Dynamics of the immune computersR versus time t, of system (7), for α
∈
{
0.1
,
…
,
1
}. Parameter values and initial conditions are as in Case 2, except for β
1
=
β
2
=
0.035.Figure 10
Dynamics of the susceptible removable devicesR
S versus time t, of system (7), α
∈
{
0.1
,
…
,
1
}. Parameter values and initial conditions are as in Case 2, except for β
1
=
β
2
=
0.035.Figure 11
Dynamics of the infected removable devicesR
I versus time t, of system (7), for α
∈
{
0.1
,
…
,
1
}. Parameter values and initial conditions are as in Case 2, except for β
1
=
β
2
=
0.035.In Figure7, we plot the dynamics of the number of susceptible computers, S, versus time, t, of model (7), for α
∈
{
0.1
,
…
,
1
}.In Figures8 and 9, we plot, respectively, the dynamics of the number of infected computers, I, and of immune computers, R, versus time t, of model (7), for α
∈
{
0.1
,
…
,
1
}.In Figures10 and 11, we plot, respectively, the values of the susceptible removable devices, R
S, and of the infected removable devices, R
I, versus time t, for α
∈
{
0.1
,
…
,
1
}.These charts reveal that Case 2 leads to much richer dynamics than Case 1. From Figures7–11, we observe that by varying α, we can tune the dynamics of the transient. Furthermore, we verify again that for α
=
0.1, we have the slowest evolution, whereas for α
=
0.7, the transient is the fastest. Moreover, from Figure 7, we conclude that the overshoots for different values of α remain globally almost invariant. Even so, for α
=
0.3, we get the smallest amplitude of the overshoots. As for the transients of the other charts, the fastest dynamics is at α
=
0.7; therefore, the smaller peak time occurs for α
=
0.7 and the larger one for α
=
0.1. It is interesting to note a slight overshoot for the number of susceptible removable devices, R
S, for α
=
{
0.1,0.2,0.3
}.


## 5. Conclusions

We proposed a fractional order model for computer virus propagation, where the interaction of computers and removable devices is included. We simulated the system for distinct values of the orderα of the fractional derivative and two sets of initial conditions. We conclude that we can tackle a family of distinct dynamical responses, including very fast transients and super-slow responses, usual in systems with long range memory. The set of initial conditions denoted as Case 1 produces a set of monotonous family of responses, while Case 2 presents more complex patterns. In both cases, the fractional orders α
=
0.1 and α
=
0.7 are those that lead to the slowest final evolution and the fastest transients, respectively. Globally, we note that fractional calculus is a mathematical tool well suited to model dynamical systems including important memory phenomena.



---

*Source: 476502-2014-06-03.xml*