Fractional Dynamics of Stuxnet Virus Propagation in Industrial Control Systems

: The designed fractional order Stuxnet, the virus model, is analyzed to investigate the spread of the virus in the regime of isolated industrial networks environment by bridging the air-gap between the traditional and the critical control network infrastructures. Removable storage devices are commonly used to exploit the vulnerability of individual nodes, as well as the associated networks, by transferring data and viruses in the isolated industrial control system. A mathematical model of an arbitrary order system is constructed and analyzed numerically to depict the control mechanism. A local and global stability analysis of the system is performed on the equilibrium points derived for the value of α = 1. To understand the depth of fractional model behavior, numerical simulations are carried out for the distinct order of the fractional derivative system, and the results show that fractional order models provide rich dynamics by means of fast transient and super-slow evolution of the model’s steady-state behavior, which are seldom perceived in integer-order counterparts.


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, .stuband 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 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][8][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][25][26][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:

•
A novel fractional-order Stuxnet virus model is proposed by exploiting the rich heritage of fractional calculus in an isolated and air-gapped network environment.

•
Stability analysis of Stuxnet virus model for both local and global equilibrium points when disease-free, and endemic spread is performed.

•
Correctness of the proposed Grunwald-Letnikov-based fractional numerical solver is ascertained, with close results to the state-of-the-art Runge-Kutta numerical solver for integer-order variants of the model.

•
Numerical simulation with Grunwald-Letnikov-based fractional numerical solver for a distinct order of the fractional derivative terms in the system shows that fractionalorder models offer rich characteristics by way of ultrafast transience and ultra-slow advancements of steady-state.

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 α 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][32][33][34][35][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: The definition of Caputos fractional derivatives can be written as: for (n − 1 < α < n) and where Γ(•) is a gamma function.
The RL definition is given as: (3) For (n − 1 < α < n), while a and t are the bounds of the operation for a D α t , the Laplace transform method is normally used with CP, GL and-RL fractional derivatives under zero initial conditions, as: [40] £{ a D ±α t f (t); s} = s ±α F(s), while the analytical expressions are represented by Mittag-Leffler (ML)-type functions [41] introduced by Agarwal and Humbert [42] and are given mathematically as: where C represents the set of complex numbers and E α,β is a two-parameter-based ML function.

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: where (n − 1 < α < n) , using Equation (1), Ivo Petras [44] provided the final recursive expression of a GL-based solver is as follows: simplifying above relation, we have In case of discrete input grids between interval t ∈ [0, T] = [0, h, 2h, . . ., Mh = T], where h represents the step size parameter, so [0, T] = [t 0 = 0, t 1 , . . ., t M = T] and any grid points in the interval are represented as t m = mh for m = 0, 1, 2, . . ., M. Thus, in a discrete form, the above equation is written as: In simple usage, the above term is written as: where c (α) j is defined as: or equivalently with GL numerical solver in the recursive form is written as: A further elaboration of the Grunwald-Letnikov (GL)-based numerical solver can be seen in [45].

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 U s (t) and infected storage media U s (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) = U s (t) + U I (t).In the rest of the article, the variables with respect to time t, S(t), I(t), M(t), U s (t), U s (t), N(t), and U(t) are denoted by S, I, M, U s , U I , N, and U, respectively.Let A 1 and A 2 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 ρ , β 1 is the infectious contact rate of susceptible nodes with infected nodes during the network scan, and β 2 denotes the contact rate of infectious-removable storage media with susceptible computer nodes, r 1 and r 2 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 2 32 , and the probability of finding susceptible nodes in IPv4 scheme is S/2 32 .Susceptible nodes can be infected at the rate β 1 SI or at β 2 SU I N, while the removable storage media could be infected at a rate of β 2 U s I 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: where α ∈ [0, 1] represents the order of the fractional derivative term D α = d α 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 The change in population is given by c and the values of these constants may be negative, positive or zero.Solving the set of Equation ( 10), we get The system given in Equation ( 9) can be simplified by incorporating N and U variables, as in: where Using Equation (11) in system ( 12), one may obtain a limit system (I MU I ), as in [49,50]: The equations in system (14), are the reduced version of ( 9), and will be used in further investigations.

Model Analysis
In this unit, stability analysis of the model is performed through the derivation of basic reproduction number, R 0 .The endemic and disease-free equilibrium points of the system are investigated for a local as well as global stability analysis.

Basic Reproduction Number (R 0 )
In epidemiology modeling, a basic reproduction number is defined as the advent of a new infection in an entirely susceptible population due to an infected individual, and is usually represented by R 0 .The value of R 0 determines the spread of infection; for R 0 > 1 infection will spread in the population, and for R 0 < 1 infection will soon end [51].
To simplify the derivation process, a reduced model ( 14) has been utilized for further investigation of R 0 .The calculation of R 0 is based on the value of α = 1.The necessary condition of a disease epidemic is based on the increase in the infected individuals, with the supposition that, initially, the entire population is susceptible.
For the case of D α I > 0, we have and, accordingly, in case of D α U I > 0, we have With the assumption that all the population is susceptible at the start, the above expressions may be written as: Simplifying the above relations, we have Accordingly, Equation ( 20) represents the basic reproduction number derived for the model.

Equilibria Studies
In this subsection, we study the equilibrium points of FO-SVM model Equation ( 14).The FO-SVM model has virus-free equilibrium and endemic equilibrium points.In the endemic equilibrium point, the spread of infection is observed.
For equilibria studies, we have equilibrium points of system (14) for virus-free and endemic are as: K 0 = (I, M, U I ) = (0, 0, 0) and K * = (I * , M * , U * I ) for R 0 > 1.The analysis for the endemic equilibria of model ( 14) is written as: Solving the equations in set (21), we obtain the expressions for the endemic equilibrium point (I * , M * , U * I ) as: where It is evident from Equation ( 22) that the possibility of infection spread, i.e., I * > 0, is only verified for the value of R 0 > 1.

Disease Free Equilibrium
Theorem 1. Disease-free equilibrium (DFE) point of a system is locally and asymptotically stable at K 0 , if R 0 < 1.
Proof.The DFE point of a system is locally asymptotically stable at K 0 = (I, M, U I ) = (0, 0, 0).The Jacobian matrix of function f : R 3 → R 3 with components: Thus, the Jacobian matrix at K 0 , DFE point of integer-order model ( 14) is given as: and simplify as: The corresponding Eigen values of the above relation are Simplifying the above expression to find the remaining Eigenvalues and rearranging the above expression and, for R 0 < 1, Equation ( 9) can be written as: Using the expression (31) in Section 4.3, make the coefficient positive for R 0 < 1, which shows that system Section 4.3 eigenvalues are in a stable region; this confirms that the system is asymptotically stable for point K 0 when R 0 < 1.If system is stable for the value of α = 1, it will be stable for the value of α < 1, as reported in [52].This completes the proof.Theorem 2. If R 0 < 1, then point K 0 is globally asymptotically stable, and otherwise unstable.

Proof. Considering the Lyapunov function mentioned below,
The function in R 3 is positive, for R 3 = (I, M, U I) and (I > 0, M > 0, U I > 0).For α = 1, the derivative of Lyapunov function (32) is For R 0 < 1, this implies that D α L ≤ 0 and K 0 is the only invariant set of system (21).According to the LaSalle Invariance Principle, K 0 is proven to be globally asymptotically stable.Hence, equilibrium point K 0 is globally asymptotically stable for R 0 < 1.Additionally, if the system is stable for the value of α = 1, then the system will be stable for α < 1, as described in [52].

Endemic Stability
The endemic stability of equilibrium point K * = (I * , M * , U * I ) is investigated in this section for the values of R 0 > 1 and I * ≥ 0.
Proof.Consider the function f : R 3 → R 3 with components and the Jacobian matrix of the system (14) as: The endemic equilibrium of system ( 14) is K * = (I * , M * , U * I ), for the value of α = 1, the Jacobian matrix at endemic point is mentioned below. where For stability analysis, Hurwitz criteria may be used, as reported in [53,54] for system (36).Equating the Equation ( 36 Determinants (D 1 , D 2 and D 3 ) of the Equation ( 36) are stated in Hurwitz as: using the value of Equation ( 20) for R 0 > 1 as: Here, D 2−1 represent the positive terms in D 2 , while, for the remaining terms, represented with D 2−2 , we have using the value of R 0 > 1, and after simplification, the above expression becomes as a result Thus, all the values of D 1 , D 2 and D 3 are positive, so all the eigenvalues of the Equation ( 36) are negative, for R 0 > 1.This proves that the endemic equilibrium point K * is locally asymptotically stable.The proof of theorem is completed.

Simulation and Results
In this section, the results of numerical simulations for FO-SVM are presented to understand the dynamics of virus spread in a critical network infrastructure in the presence of removable storage connectivity, which may compromise the air-gap between the networks.Numerical experimentation is conducted for the designed FO-SVM as given in Equation ( 9) for a different variation in parameters and initial start-up scenarios, as given in Tables 1 and 2, respectively.The dynamic behavior of the fractional order (FO) model is studied by varying the non-integer order derivative α.Most FO differential systems lack exact analytical solutions, so the numerical solver based on Grunwald-Letnikov (GL) procedure, as described in Section 2 is exploited for an approximate solution to the model.The security firms, including Symantec, tracked 100,000 infected computers as of 29 September 2010, in the world.Additionally, available real data are used to validate the accuracy and convergence of the model for the Stuxnet virus spread.The virus infects approximately 100,000 users from 155 different countries, and 63% were only in Iran.Due to this attack, the number of hosts that lost functionality (hardware connected to these hosts was damaged due to sudden increase in frequency of up to 1410 Hz, which then decreased to 2 Hz and increased to 1064 Hz in spite of the normal working range of from 807 Hz to 1210 Hz) due to virus attack.A virus operates the machines connected with the hosts at an extreme range of frequencies dictated by Stuxnet and caused physical damage to 1500 centrifuge machines (approximately 1200 in Iran only).Approximately 3280 unique samples and variants of the Stuxnet virus were recorded by Symantec and other security corporations [3,6,55].In order to establish the working accuracy of GL-based numerical solvers, the results of the scheme are compared with state-of-the-art numerical solvers based on the Runge-Kutta (RK) method for the value of α = 1.The results are determined for nine cases of integer order models (9) by a GL-based computing technique for inputs t ∈ [0, 60] with step size h = 0.001 (time t represents months).Numerical solutions to the model for the same inputs are also calculated by the RK method for each variation.Figure 3 highlights the comparison of model behavior with Stuxnet virus real-world data.FO-SVM model results shown in Figure 3 are calculated using the RK method to assume the value of fractional order α = 1.β 1 = 0.366, β 2 = 0.6, ρ = 0.00265, r 1 = 0.1126, r 2 = 0.0088, S = 2.3 × 10 6 , I = 10,000, M = 10, U s = 50,000, U I = 10,000.
In Figure 3, the number of hosts versus time in months is plotted, which shows the effect of the Stuxnet attack on the number of hosts as time passes.The number of infected hosts is 96,760 (real infected host number was 100,000), and the number of damaged hosts is 1500 (real damaged host number was 1500) in 23 months time, which shows the model accuracy for real-world virus data, as shown in Figure 3, with red and blue dots, respectively.In this case, removable media are considered to be 60,000, and, after increasing the number of removable-storage media, infection in the host nodes also increases (96,760 after 23 months).
The number of infected removable-storage devices is 43,740 in 23 months, and in 24 months, the time number of infected devices increases to 44,920.An increase in the number of damaged hosts is observed after the increase in infected hosts in 24 months' time.This highlights the role of removable-storage media in spreading the infection in isolated critical networks in the absence of any remedial strategy in the model.Stuxnet is an advanced, persistent threat (APT) type of malicious code that penetrates in the remote system in a quasi-autonomous fashion.Then, a 23-month decline in the number of infected hosts is observed due to the availability of remedial technique and other controlling mechanisms.However, the Stuxnet virus was carried by removable-storage media spreads in other•networks.
In Figure 4, the solutions to the RK method with GL solver is compared with an error analysis of susceptible hosts S: a and b for cases 2 to 4, c and d for cases 5 to 7, and e and f for cases 8 to 10. Comparisons of results from both the RK numerical solver and GL-based method (for fractional-order α = 1) are presented for susceptible hosts S in nine cases.The error analysis, based on the absolute difference between the two approaches, is also plotted in Figure 4 to assess closeness.The results show a matching of both solutions of up to three decimal places of accuracy.The small errors in these plots show that the results of the GL method are in good agreement with the standard RK numerical technique, which establishes the working accuracy of the GL-based solver.In Figure 5, the solution of the RK method with the GL solver is compared in the case of infected hosts I and damaged hosts M: a and b for cases 1 to 3, c and d for cases 4 to 6, and e and f for cases 7 to 9. Figure 4 compares solutions for the RK method with GL solver in case of susceptible and infected removable-storage media: a and b for cases 1 to 3, c and d for cases 4 to 6, and e and f for cases 7 to 9. In Figures 5 and 6, the solution of the RK method with a GL solver are compared and presented for infected nodes I, damaged node M, susceptible removable-storage media U s and infected removable-storage media U I , respectively, for nine model cases.These nine cases also explain the virus spread behavior in different scenarios.Considering Figures 4-6, and the different cases simulated, we have the following comments.
The effect of changing the infectious contact rate β 1 from 36.6% to 60% is highlighted in case 1 of Equation ( 9) (value of β 1 in Figure 3 is 36.6%).It is observed that the number of infected hosts in 24 months is 96,760, as shown in Figure 5a (in Figure 3, the number of infected hosts in 24 months is 96,270), which shows a slight increase in infected hosts due to β 1 .In case 2, the number of initially infected hosts is assumed to be 30,000.Increasing the contact rate of infectious removable media (in case 2) reduces the number of susceptible hosts rapidly as compared to case 1 (Figure 4a).However, the number of infected hosts is reduced (Figure 5a) due to an increase in the natural removal rate of hosts and removable storage r 1 and r 2 (hosts are removed to save them from the Stuxnet attack).In case 3, we reduce the damage rate and the quantity of initial susceptible removable-storage media, which reduces infected removable-storage media number (Figure 6b) and increases the infected hosts, as in Figure 5a).A decrease in damaged hosts is observed in case 3, despite the increase in the number of infected hosts.In case 4, FO-SVM model dynamics are observed by increasing the arrival rate of new nodes and the arrival rate of new removable-storage devices, as mentioned in Tables 1  and 2. The results show that increasing the arrival rate of new hosts and arrival rate of new removable-storage media will not spread the infection faster without the presence of a sufficient number of infected removable-storage devices, as shown in Figure 5c.In cases 5 and 6, we further increase the values of the arrival rate of new nodes as well as removable-storage devices for an in-depth behavior analysis of the model.Both cases have similar parameters, except case 6, which represents a higher damage rate (especially for zero-day vulnerability or for a new virus attack) that increases the number of damaged computers and reduces the number of infected computers (removed due to high damage rate) in the networks as compared to case 5. Case 5 shows the high number of infected nodes (Figure 5c) because the Stuxnet virus only destroys the machines with specific hardware (Siemens specific PLCs) and remains dormant till it finds the target.In case 6 (Figure 5c,d), the number of infected hosts decreases; however an increase in the number of damaged hosts is observed due to an increase in damage rate ρ.(a) Comparison of RK method with GL solver for susceptible removable storage media in cases 1 to 3, (b) Comparison of RK method with GL solver for infected removable storage media in cases 1 to 3, (c) Comparison of RK method with GL solver for susceptible removable storage media in cases 4 to 6, (d) Comparison of RK method with GL solver for infected removable storage media in cases 4 to 6, (e) Comparison of RK method with GL solver for susceptible removable storage media in cases 7 to 9, (f) Comparison of RK method with GL solver for infected removable storage media in cases 7 to 9.
In case 7, the values of β 1 and β 2 of case 6 are swapped to observe the behavior of the model.In case 7, the value of β 1 = 0.745, as compared to 0.4 in case 6, and the value of β 2 = 0.4, as compared to 0.745 in case 6.These swaps are carried out to observe the devastation effect of infected removable storage media as compared to the effect of infected nodes in the model, because infected removable media have a greater devastation effect.Simulation results show that the number of damaged nodes in case 6 is 35,000, whereas, in case 7, it is 5000, due to a decrease in the value of β 2 infectious contact rate of removable storage media (Figure 5e).
However, by increasing β 2 value (removable-storage media's infectious contact rate with susceptible computers) and A 2 (the arrival of removable-storage media) for case 8 will also increase the infection in the network.This outlines the importance of removablestorage media in spreading the virus in air-gapped networks (Figure 5e).In case 9, the contact rate of susceptible computer nodes with infectious removable-storage media β 2 is reduced, which results in a reduction in damaged nodes (Figure 5f) and infected nodes (Figure 5e), and an increase in the number of susceptible storage devices (Figure 6e).Case 9 further elaborates the scenario presented in case 8.
The derivative order α = 1 is presented in Figures 4-6.The effect of change in fractional order α is presented in Figures 7-11.A detailed analysis of the FO-SVM model is conducted by changing the fractional order α in the system ( 9), such that one may observe fastchanging as well as super-slow growth in the model dynamics.The fractional order solution of the FO-SVM model for different values of the fractional order α is solved using a GL-based solver.The solutions are determined for nine cases of FO-SVM by a GL-based numerical procedure for different fractional orders, i.e., α = [0.5, 0.75, 0.8, 0.85, 0.9, 0.95, 1], for the inputs t ∈ [0, 60] with step size h = 0.001.Results for the dynamics of the FO-SVM model in terms of susceptible S, infected I, and damaged M computers are plotted in Figures 7-9 for cases 1-3, 4-6, and 7-9, respectively.Susceptible removable-storage media U s and infected-removable-storage media U I are plotted in Figures 10 and 11 for cases 1-4 and 5-9, respectively, for different values of the fractional order α.
Figure 7 shows a simulation of fractional order, i.e., α = [0.5, 0.75, 0.8, 0.85, 0.9, 0.95, 1] for the FO-SVM model for different values of fractional order α for case 1-3 of susceptible S, infected I and damaged hosts M. In Figure 7, the number of susceptible, infected and damaged hosts is plotted versus time for cases 1-3 for different values of α = [0.5, 0.75, 0.8, 0.85, 0.9, 0.95, 1].A consistent pattern is observed in the evolution of curves with the value of α.The value of infected hosts in case 1 with α = 1 is 96,760, and for α = 0.95, the value of infected hosts is approximately 56,000 for t = 24 months, as shown in Figure 7b.In Figure 7c, the number of damaged hosts (hosts that were connected with specific models of Siemens PLCs) for the value of α = 0.95 are 1000 for t = 30.Adjusting the value of α to 0.98 may adjust the number of damaged hosts to 1500, which matches the real published data of the Stuxnet virus.This illustrates the controllability feature of α for tuning the model.Despite the rapid spreadability of the Stuxnet virus, it causes little or no harm to the systems that do not have specific hardware.
Figure 8 shows the simulation of fractional order dynamics of the FO-SVM model for different values of fractional order α for cases 4-6, and Figure 9 depicts the simulation of fractional-order dynamics of the FO-SVM model for case 7-9.Figures 8 and 9 highlight the results for variation in fractional order α, which shows that variation in α gives smooth variations in the dynamics of the model.For α = 0.1, we have the slowest evolution.Simulation of fractional order dynamics of FO-SVM model for different values of fractional order α for case 1-5 of susceptible removable-storage media U s and infected-removablestorage media U I are illustrated in Figure 10. Figure 11 shows the simulation of fractional order dynamics of FO-SVM model for different values of fractional order α, for cases 5-9 of susceptible removable-storage media U s , and infected-removable-storage media U I .In Figures 10 and 11, the number of susceptible storage media and infected storage media is plotted for case 1-9 against the time variation for different values of fractional order α = [0.5, 0.75, 0.8, 0.85, 0.9, 0.95, 1].It is observed that tuning the values of α tunes the dynamics of transients, as shown in Figure 10a.The value of susceptible storage media for t = 1 and α = 0.95 is 35,000, which is effectively reduced to 10,000 by a slight change in fractional order α from 0.95 to 0.8.In contrast, a slow change is observed in the dynamics of the FO-SVM model for α = 0.1.Increasing the fractional order α increases the rate of change of the variables.Fractional-order virus models provide extra tunable parameters in the form of α, which highlight more minute changes in the model dynamics.

Conclusions
A detailed analysis of the novel design of the fractional order Stuxnet virus model is presented, with richer dynamics for the transmission of virus spread in an isolated critical network through removable-storage media.The fractional-order Stuxnet-virusbased mathematical models are found to be at least as stable as integer-order models.The fractional order value α of the proposed fractional Stuxnet virus model more effectively controls the solution reachability towards a steady state point.Additionally, the fractional order system of the Stuxnet virus model can tackle the different responses, including super-slow evolutions and very fast transients; these responses are found to have long memory characteristics in the system.Taking the value of α = 0.98, one may adjust the number of damaged hosts to 1500 in case 1, which matches the damage caused by the Stuxnet virus.The transformation process of the classical model into a fractional model is very sensitive to the value of the order of differentiation α, and can be converted to a simple SIR model if we choose the values of the infectious contact rate β 2 = 0.A theoretical analysis of the model capturing the Stuxnet virus-spreading characteristics is determined by a mathematical derivation of the basic reproduction number R 0 for the value of α = 1.Equilibrium points of the model are globally and asymptomatically stable for R 0 < 1 and R 0 > 1, respectively.
In the future, one may exploit the strength of stochastic numerical solvers [56][57][58][59][60][61] based on fractional evolutionary and swarming techniques [62][63][64][65][66][67] for a detailed analysis of the designed fractional-order Stuxnet virus model.Additionally, new definitions of the fractional operator, such as Yang-Machado [35] and Yang-Abdel-Aty-Cattani [36] fractional derivatives looks promising for the development of new computing solvers for the numerical solution of the fractional-order Stuxnet virus model and other fractionalorder systems with better theoretical justifications, a better applicability domain, proof of the accuracy, convergence, stability, and robustness.

Figure 4 .
Figure 4. Solution comparison of the RK method with GL solver and error analysis with susceptible S hosts: a and b for cases 2 to 4, c and d for cases 5 to 7, and e and f for cases 8 to 10.(a) Solution comparison of the RK method with GL solver for cases 2 to 4, (b) error analysis for cases 2 to 4, (c) Solution comparison of the RK method with GL solver for cases 5 to 7, (d) error analysis for cases 5 to 7, (e) Solution comparison of the RK method with GL solver for cases 8 to 10, (f) error analysis for cases 8 to 10.

Figure 5 .
Figure 5. Solution comparison of RK method with GL solver for infected hosts I and damaged hosts M; a and b for cases 1 to 3, c and d for cases 4 to 6 while e and f for cases 7 to 9. (a) Comparison of RK method with GL solver for infected hosts in cases 1 to 3, (b) Comparison of RK method with GL solver for damaged hosts in cases 1 to 3, (c) Comparison of RK method with GL solver for infected hosts in cases 4 to 6, (d) Comparison of RK method with GL solver for damaged hosts in cases 4 to 6, (e) Comparison of RK method with GL solver for infected hosts in cases 7 to 9, (f) Comparison of RK method with GL solver for damaged hosts in cases 7 to 9.

Figure 6 .
Figure 6.Solution comparison RK method with GL solver for susceptible and infected-removablestorage media: a and b for cases 1 to 3, c and d for cases 4 to 6, and e and f for cases 7 to 9.(a) Comparison of RK method with GL solver for susceptible removable storage media in cases 1 to 3, (b) Comparison of RK method with GL solver for infected removable storage media in cases 1 to 3, (c) Comparison of RK method with GL solver for susceptible removable storage media in cases 4 to 6, (d) Comparison of RK method with GL solver for infected removable storage media in cases 4 to 6, (e) Comparison of RK method with GL solver for susceptible removable storage media in cases 7 to 9, (f) Comparison of RK method with GL solver for infected removable storage media in cases 7 to 9.

Figure 7 .
Figure 7. Simulation of fractional order dynamics of FO-SVM model for different values of fractional order α for cases 1 (a-c), 2 (b-f) and 3 (g-i) of susceptible S, infected I and damaged hosts M.

Figure 8 .
Figure 8. Simulation of fractional order dynamics of FO-SVM model for different values of fractional order α for cases 4 (a-c), 5 (b-f) and b (g-i) of susceptible S, infected I and damaged hosts M.

Figure 9 .
Figure 9. Simulation of fractional order dynamics of FO-SVM model for different values of fractional order α for cases 7 (a-c), 8 (b-f) and 9 (g-i) of susceptible S, infected I and damaged hosts M.

Figure 10 .
Figure 10.Simulation of fractional order dynamics of FO-SVM model for different values of fractional order α for cases 1 (a,b), 2 (c,d) and 3 (e,f), 4 (g,h) and 5 (i) of susceptible removable-storage media U s and infected-removable-storage media U I .

Figure 11 .
Figure 11.Simulation of fractional order dynamics of FO-SVM model for different values of fractional order α for cases 5 (a), 6 (b,c) and 7 (d,e), 8 (f,g) and 9 (h,i) of susceptible removable-storage media U s and infected-removable-storage media U I .

Table 1 .
Values of parameters used in model simulation for different scenarios.

Table 2 .
Starting values of variables used in the simulation of different scenarios.