**1. Introduction**

A small piece of software code or program in a computer system that works on a system without the consent of the user may cause damage or steal information for the exploitation of the desired targets. In strategic conflicting environments, as well as in the financial market, computer viruses can be used in a network operation as a digital weapon against the desired targets, e.g., a computer spyware program used as an information collection platform in the Syrian war [1], or Shamoon and Stuxnet viruses for cyber incidents [2]. The tools used for cyberwar vary from a tiny code that exhibits annoying messages on the console to a complicated routine that physically damages the system, such as Stuxnet [3]. Stuxnet was discovered at Natanz, Iran, a nuclear enrichment facility, in June 2010 [4]. The name of the Stuxnet virus was derived from two keywords in its source code, .stub and mrxnet.sys. The Stuxnet virus is a sophisticated piece of code that mainly targets the supervisory control and data acquisition systems (SCADA), exploits zero-day vulnerabilities/bugs to attack the targeted hosts, and uses advanced technology to hide from guard programs. The Stuxnet virus exploits different services, such as a print spooler (MS 10-061), the zero-day vulnerability of the windows system, network shares, file-sharing and server message block (SMB), etc. Stuxnet virus monitors the frequency of motors operating centrifuge machines before modification, which must be in the range of from 807 Hertz to 1210 Hertz. Stuxnet virus controls the running frequency

**Citation:** Masood, Z.; Raja, M.A.Z.; Chaudhary, N.I.; Cheema, K.M.; Milyani, A.H. Fractional Dynamics of Stuxnet Virus Propagation in Industrial Control Systems. *Mathematics* **2021**, *9*, 2160. https:// doi.org/10.3390/math9172160

Academic Editors: Mikhail Posypkin, Andrey Gorshenin and Vladimir Titarev

Received: 20 June 2021 Accepted: 31 August 2021 Published: 4 September 2021

**Publisher's Note:** MDPI stays neutral with regard to jurisdictional claims in published maps and institutional affiliations.

**Copyright:** © 2021 by the authors. Licensee MDPI, Basel, Switzerland. This article is an open access article distributed under the terms and conditions of the Creative Commons Attribution (CC BY) license (https:// creativecommons.org/licenses/by/ 4.0/).

of centrifuge machines for a short interval of time to 1410 Hertz and then decreases to 2 Hz and increases to 1064 Hertz. A change in the output frequency of the motors essentially sabotages the working of machines [5]. Due to the attack of the Stuxnet virus, approximately 1000 centrifuge machines were out of order, of a total of 5000 machines operating in the Iran nuclear facility at Natanz [6]. The purpose of the virus was not just to infect the computers, but to cause real-world physical damage.

A theoretical study of the Stuxnet's malicious code behavior was conducted through the strength of epidemic modeling of virus spread [7–9]. The control scheme of these malicious codes is very challenging because they often hide, and may exploit zero-day vulnerabilities, gain administrative rights and execute code as an authenticated program. The development in technologies creates new issues regarding the safety and security of the critical infrastructure of the countries in the presence of these vulnerabilities and smart viruses. The desire to manufacture an automated process immensely increases software dependencies, which ultimately require lengthy and complex routines.

These complex codes are challenging to screen out completely using software testing mechanisms, and leftover vulnerabilities in these codes can compromise the whole system [10]. Therefore, the comprehensive and dynamic study of these codes is a promising domain for research communities to investigate.

The spread of the virus in a computer network is closely related to the spread of biological viruses in the population. Mathematical and statistical models are often based on concepts and methods borrowed from physics. Models play an important role in infection control by quickly predicting and understanding disease outbreaks. In recent decades, new infectious diseases have been observed, together with the development of eliminated technologies.

The ability to quickly measure the unfolding of outbreaks, communications, and movements is key to capturing the spread of a virus. The inherent complexity of such methods limits the study of these processes. However, developments in technology are helping to lift these limitations [11]. Classical approaches and linear thinking are unable to effectively mitigate the problem due to the lack of equilibrium and non-linear nature of the problems. A complex system, its counter-intuitive behavior, and other macrolevel changes can be addressed by applying complex sciences. The usual models did not provide an in-depth picture of real system dynamics because these systems neglect feedback scenarios, cascade effects, and instabilities. To predict the global-scale spread of disease dynamics, several factors, such as demographic disparity, mobility scenarios which include air-flow system, commuter movement in the area, disease-specific information, and control mechanisms, should be acccounted for. There has long been work on the development of mathematical models for use in the analysis of infectious disease behavior. The mathematical model of Daniel Bernoulli against smallpox disease was published in 1766. Mathematical models of these types were designed to elaborate the behavior of an epidemic over the course of time, in which every single population of the virus is assumed to interact with the individual of other populations. The ability to monitor hidden outbreaks, as well as contact and communication, are key to the portrayal of disease-spreading [12]. It is known that immunizing a large fraction of the population or a computer network, the epidemic that spreads upon contact between infected nodes or individuals can be stopped.

Some diseases require 80–90% immunization (measles requires 95%), and the same is true for the computer, where 100% immunization from the Internet may stop viruses in connected networks [13]. Mathematical modeling of infectious disease or viruses in biology or in computer systems gives us a thorough understanding of the problem and helps us to devise a reliable, viable, and robust control strategy [14]. It was observed that the state of the various biological organisms at a certain time depends on its past states and fractional derivatives that also contains those characteristics. Thus, a fractional derivative is a natural approach to the solution of these biological systems. Mathematical modeling is used in numerous disciplines of science and engineering problems [15,16]. Kermack and McKendrick founded mathematical modeling at the beginning of twentieth century with a series of publications, and introduced a susceptible, infected and recovered epidemic model [17]. In this field, several other scientists, biologists, computer engineers and mathematicians have worked on epidemic modeling and published work in this area, such as time delay virus models [18], a fractional epidemiological model [19], antivirus strategy for computer virus model [20], modified susceptible, infected and susceptible models [21] and epidemic models with two control mechanisms, quarantine and immunity [22], and models that highlight the topological facets of the network [23]. Besides these, the role of fundamental concepts and underlying theories of fractional calculus was effectively applied in modeling complex systems in diversified fields with rich dynamics compared to its integer counterparts [24–27]. Considering these facts, the current study aims to exploit the rich heritage of fractional dynamics for the development of the fractional Stuxnet virus model by using the features of the Stuxnet model to illustrate the virus spread in SCADA systems [28]. In this study, a fractional-order mathematical model of the Stuxnet virus is presented to study the ultra-fast transient and slow evolutions of the virus spread dynamic and attack pattern on isolated critical infrastructures, managed by industrial control computers. The contribution of the proposed fractional Stuxnet virus model is briefly described as:


#### **2. Fractional Calculus Fundamentals**

## *2.1. Preliminaries*

Fractional calculus is a branch of mathematics and a generalization of the calculus theory of integrals and derivatives of a real number or complex number power. The discussion of fractional calculus was started 300 years ago, and the idea of fractional calculus was first listed in the literature with a letter from Leibniz to L'Hostital in 1696. In this letter, a half-derivative term was introduced, i.e., the generalization of the derivative operator *D<sup>α</sup> f*(), where *α*, representing the order of a fractional derivative. The history of the fractional derivative is as long as the classical differential operators in calculus, but the inherent strength of the fractional operator is relatively less exploited in engineering domains until the early 1980s. The physical interpretation of the fractional derivative outcomes is still ambiguous, and remained an open debate for clarity in the research community. However, the fractional derivative is an inspiring operator to describe the physics of many modeling phenomena, which are difficult to realize through integer-order derivatives. Recently, the kernel function of a fractional derivative is referred to as a memory function, and fractional-order derivative is proposed as a memory index [29,30] with different types of kernel [31–36]. The theory development of fractional calculus belonged to the efforts of several scientists, such as Letnikov, Liouville, Euler, and Riemann [37,38]. Different definitions of fractional order derivatives have existed; the most-used definitions are those of Riemann–Liouville (RL), Caputo (CP), and Grunwald–Letnikov (GL) [39]. The GL definition of fractional derivative is as follows:

$$\sideset{\_{a}}{}{\mathop{L}}{}{\mathop{L}}^{a}\_{t}f}(t) = \lim\_{h \to 0} \frac{1}{h^{\mathfrak{a}}} \sum\_{m=0}^{(t-a)/h} (-1)^{m} \binom{a}{m} f(t-mh), \ t > a, a > 0. \tag{1}$$

The definition of Caputos fractional derivatives can be written as:

$$\, \_a^{CP}D\_t^a f(t) = \frac{1}{\Gamma(n-a)} \int\_a^t \frac{f^n(x)}{(t-x)^{n-n+1}} dx,\tag{2}$$

for (*n* − 1 < *α* < *n*) and where Γ(·) is a gamma function.

The RL definition is given as:

$$\, \_a^{RL}D\_t^a f(t) = \frac{1}{\Gamma(n-a)} \frac{d^n}{dt} \int\_a^t \frac{f(\mathbf{x})}{(t-\mathbf{x})^{n-n+1}} d\mathbf{x}.\tag{3}$$

For (*<sup>n</sup>* <sup>−</sup> <sup>1</sup> <sup>&</sup>lt; *<sup>α</sup>* <sup>&</sup>lt; *<sup>n</sup>*), while *<sup>a</sup>* and *<sup>t</sup>* are the bounds of the operation for *aD<sup>α</sup> <sup>t</sup>* , the Laplace transform method is normally used with CP, GL and-RL fractional derivatives under zero initial conditions, as: [40]

$$\mathcal{E}\{\_{a}D\_{t}^{\pm \mathfrak{a}}f(t);s\} = s^{\pm \mathfrak{a}}F(s),\tag{4}$$

while the analytical expressions are represented by Mittag–Leffler (ML)-type functions [41] introduced by Agarwal and Humbert [42] and are given mathematically as:

$$E\_{a,\emptyset}(z) = \sum\_{k=0}^{\infty} \frac{z^k}{\Gamma(\beta + ak)},\tag{5}$$

$$
\alpha, \beta, z \in \mathbb{C}, \text{ } \Re(\alpha) > 0, \Re(\beta) > 0, \text{ } \cdot
$$

where *C* represents the set of complex numbers and *Eα*,*<sup>β</sup>* is a two-parameter-based ML function.

#### *2.2. Grunwald–Letnikov-Based Numerical Solver for FDEs*

Analytical solution to the fractional differential equations (FDEs) generally determined through the Laplace transform method (4), and these expressions are commonly represented by the ML function (5), while, for the numerical solutions, the most commonly used method is based on GL definition.

To introduce the numerical solver based on GL [43] for FDEs, let a general from of an FDE, along with its initial conditions, is given as follows:

$$\begin{aligned} \, \_aD\_t^a f(t) &= f(y(t), t), \\ \, \_y^{(k)}(0) &= y\_0^{(k)}, k = 0, 1, 2, \dots, n-1, \end{aligned} \tag{6}$$

where (*n* − 1 < *α* < *n*) , using Equation (1), Ivo Petras [44] provided the final recursive expression of a GL-based solver is as follows:

$$\frac{1}{h^a} \sum\_{j=0}^{\left[\left(t-a\right)/h\right]} (-1)^j \binom{a}{k} y(t-jh) \approx f(y(t), t),$$

simplifying above relation, we have

$$y(t) + \sum\_{j=1}^{\left[\left(t-a\right)/h\right]} (-1)^j \binom{a}{k} y(t-jh) \approx h^{-a} f(y(t), t).$$

In case of discrete input grids between interval *t* ∈ [0, *T*]=[0, *h*, 2*h*, ... , *Mh* = *T*], where *h* represents the step size parameter, so [0, *T*]=[*t*<sup>0</sup> = 0, *t*1, ... , *tM* = *T*] and any grid

points in the interval are represented as *tm* = *mh* for *m* = 0, 1, 2, ... , *M*. Thus, in a discrete form, the above equation is written as:

$$y(t\_m) + \sum\_{j=1}^{m} (-1)^j \binom{a}{j} y(t\_m - jh) = h^{-a} f(y(t\_m), t\_m), m = 0, 1, 2, \dots, M.$$

In simple usage, the above term is written as:

$$y(t\_m) + \sum\_{j=1}^{m} c\_j^a y(t\_m - jh) = h^{-a} f(y(t\_m), t\_m), \\ m = 0, 1, 2, \dots, M,$$

where *c* (*α*) *<sup>j</sup>* is defined as:

$$c^a\_j = (-1)^j \binom{a}{j} \, , \tag{7}$$

or equivalently with *c<sup>α</sup>* <sup>0</sup> = 1,

$$c\_j^{\alpha} = \left(1 - \frac{1+\alpha}{j}\right) c\_{j-1}^{\alpha}, j = 0, 1, \dots$$

GL numerical solver in the recursive form is written as:

$$y(t\_m) = f(y(t\_m), t\_m)h^{-a} - \sum\_{j=1}^{k} c\_j^a y(t\_{m-j}), m = 0, 1, 2, \dots, M. \tag{8}$$

A further elaboration of the Grunwald–Letnikov (GL)-based numerical solver can be seen in [45].

#### **3. Model Formulation of Fractional Order Stuxnet Virus**

The formulation of a fractional-order Stuxnet virus model (FO-SVM) is presented here. A detailed workflow of the proposed FO-SVM is shown in Figure 1. The entire FO-SVM is segmented into five classes: three for computer population, i.e., susceptible *S*(*t*), infected *I*(*t*), and damaged *M*(*t*), and two for removable storage media, i.e., susceptible storage media *Us*(*t*) and infected storage media *Us*(*t*). However, *N*(*t*) represents the total population, i.e., *N*(*t*) = *S*(*t*) + *I*(*t*) + *M*(*t*), and total removable devices *U*(*t*), i.e., *U*(*t*) = *Us*(*t*) + *UI*(*t*). In the rest of the article, the variables with respect to time *t*, *S*(*t*), *I*(*t*), *M*(*t*), *Us*(*t*), *Us*(*t*), *N*(*t*), and *U*(*t*) are denoted by *S*, *I*, *M*, *Us*, *UI*, *N*, and *U*, respectively. Let *A*<sup>1</sup> and *A*<sup>2</sup> represent the arrival of new computing nodes and removable storage media, respectively, damage rate caused to PLC's due to virus infection is represented by *ρ* , *β*<sup>1</sup> is the infectious contact rate of susceptible nodes with infected nodes during the network scan, and *β*<sup>2</sup> denotes the contact rate of infectious-removable storage media with susceptible computer nodes, *r*<sup>1</sup> and *r*<sup>2</sup> represent the natural removal (death) of computer nodes and removable devices from the network, respectively. The number of nodes in Internet protocol version 4 (IPv4) is 232, and the probability of finding susceptible nodes in IPv4 scheme is *S*/232. Susceptible nodes can be infected at the rate *β*1*SI* or at *β*2*SUI* " *N*, while the removable storage media could be infected at a rate of *β*2*UsI* " *N*. Removable storage media is a common source of virus spread in critical industrial air-gapped networks, which are isolated from normal networks. The removable storage devices facilitate the flow of information to and from the networks that make them as an easy prey for intruders [46]. In this study, fractional-order virus model is used to explain the spread of the virus, especially Stuxnet [47,48] in industrial networks through removable storage media. A proposed flow chart diagram of the Stuxnet virus model is shown in Figure 2, and the fundamental mathematical equations of the model are given as:

$$\begin{aligned} D^a S &= A\_1 - \frac{\beta\_1 SI}{2^{32}} - \frac{\beta\_2 SI\_I}{N} - r\_1 S\_\prime \\ D^a I &= \frac{\beta\_1 SI}{2^{32}} + \frac{\beta\_2 SI\_I}{N} - \rho I - r\_1 I\_\prime \\ D^a M &= \rho I - r\_1 M\_\prime \\ D^a U\_s &= A\_2 - \frac{\beta\_2 Ul\_s I}{N} - r\_2 U\_{s\prime} \\ D^a U\_I &= \frac{\beta\_2 Ul\_s I}{N} - r\_2 U\_I \end{aligned} \tag{9}$$

where *<sup>α</sup>* <sup>∈</sup> [0, 1] represents the order of the fractional derivative term *<sup>D</sup><sup>α</sup>* <sup>=</sup> *<sup>d</sup>α*" *dtα*. For the value of *α* = 1, the above-mentioned FO-SVM system provided in a set of Equation (9) will be converted into a first-order system. From the differential equations mentioned in (9), solving the equations by taking the value of *α* = 1, we get

$$\begin{aligned} \frac{dN}{dt} &= A\_1 - r\_1 N\_\prime\\ \frac{dII}{dt} &= A\_2 - r\_2 II. \end{aligned} \tag{10}$$

The change in population is given by *c*<sup>1</sup> = *A*<sup>1</sup> − *r*<sup>1</sup> and *c*<sup>2</sup> = *A*<sup>2</sup> − *r*2, and the values of these constants may be negative, positive or zero.

**Figure 1.** FO-SVM model proposed graphical overview.

**Figure 2.** FO-SVM model schematic flow diagram.

Solving the set of Equation (10), we get

$$\begin{aligned} N(t) &\to \frac{A\_1}{r\_1} \stackrel{\Delta}{=} N^\*, \; t \to \infty, \\\downarrow \mathcal{U}(t) &\to \frac{A\_2}{r\_2} \stackrel{\Delta}{=} \mathcal{U}^\*, \; t \to \infty. \end{aligned} \tag{11}$$

The system given in Equation (9) can be simplified by incorporating *N* and *U* variables, as in:

$$\begin{aligned} D^a I &= \frac{\beta\_1 (N - I - M) I}{2^{32}} + \frac{\beta\_2 (N - I - M) I I\_I}{N} - \rho I - r\_1 I\_\prime \\ D^a M &= \rho I - r\_1 M\_\prime \\ D^a U\_I &= \frac{\beta\_2 (U - U\_I) I}{N} - r\_2 U\_{I\prime} \end{aligned} \tag{12}$$

where

$$\begin{aligned} N(t) &= N^\* + (N(0) - N^\*)e^{-r\_1 t}, \\ \mathcal{U}(t) &= \mathcal{U}^\* + (\mathcal{U}(0) - \mathcal{U}^\*)e^{-r\_2 t}. \end{aligned} \tag{13}$$

Using Equation (11) in system (12), one may obtain a limit system (*IMUI*), as in [49,50]:

$$\begin{split} D^a I &= \frac{\beta\_1 (N^\*-I-M)I}{2^{32}} + \frac{\beta\_2 (N^\*-I-M) \mathcal{U}\_I}{N^\*} - \rho I - r\_1 I, \\ D^a M &= \rho I - r\_1 \mathcal{M}, \\ D^a \mathcal{U}\_I &= \frac{\beta\_2 (\mathcal{U}^\* - \mathcal{U}\_I) I}{N^\*} - r\_2 \mathcal{U}\_I. \end{split} \tag{14}$$

The equations in system (14), are the reduced version of (9), and will be used in further investigations.
