**5. 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].

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



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.

**Figure 3.** Simulation of Stuxnet virus spread with available data of parameters *A*<sup>1</sup> = 0.042, *A*<sup>2</sup> = 0.042, *<sup>β</sup>*<sup>1</sup> = 0.366, *<sup>β</sup>*<sup>2</sup> = 0.6, *<sup>ρ</sup>* = 0.00265, *<sup>r</sup>*<sup>1</sup> = 0.1126, *<sup>r</sup>*<sup>2</sup> = 0.0088, *<sup>S</sup>* = 2.3 <sup>×</sup> <sup>10</sup>6, *<sup>I</sup>* = 10,000, *<sup>M</sup>* = 10, *Us* = 50,000, *UI* = 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 *Us* and infected removable-storage media *UI*, respectively, for nine model cases.

**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.

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 *β*<sup>1</sup> from 36.6% to 60% is highlighted in case 1 of Equation (9) (value of *β*<sup>1</sup> 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*<sup>1</sup> and *r*<sup>2</sup> (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.

**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.

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 *ρ*.

**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.

In case 7, the values of *β*<sup>1</sup> and *β*<sup>2</sup> of case 6 are swapped to observe the behavior of the model. In case 7, the value of *β*<sup>1</sup> = 0.745, as compared to 0.4 in case 6, and the value of *β*<sup>2</sup> = 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 *β*<sup>2</sup> infectious contact rate of removable storage media (Figure 5e).

However, by increasing *β*<sup>2</sup> value (removable-storage media's infectious contact rate with susceptible computers) and *A*<sup>2</sup> (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 *β*<sup>2</sup> 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 *Us* and infected-removable-storage media *UI* 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 *Us* and infected-removablestorage media *UI* 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 *Us*, and infected-removable-storage media *UI*. 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.

**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.** 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.** 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.** 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 *Us* and infected-removable-storage media *UI*.

**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 *Us* and infected-removable-storage media *UI*.

## **6. 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 *β*<sup>2</sup> = 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*<sup>0</sup> for the value of *α* = 1. Equilibrium points of the model are globally and asymptomatically stable for *R*<sup>0</sup> < 1 and *R*<sup>0</sup> > 1, respectively.

In the future, one may exploit the strength of stochastic numerical solvers [56–61] based on fractional evolutionary and swarming techniques [62–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.

**Author Contributions:** Conceptualization, Z.M. and M.A.Z.R.; methodology, Z.M. and M.A.Z.R.; software, Z.M.; validation, M.A.Z.R.; formal analysis, N.I.C. and M.A.Z.R.; investigation, Z.M.; writing—original draft preparation, Z.M.; writing—review and editing, N.I.C. and M.A.Z.R.; visualization, N.I.C. and M.A.Z.R.; project administration, K.M.C. and A.H.M.; funding acquisition, K.M.C. and A.H.M. All authors have read and agreed to the published version of the manuscript.

**Funding:** This research received no external funding.

**Institutional Review Board Statement:** Not applicable.

**Informed Consent Statement:** Not applicable.

**Data Availability Statement:** Not applicable.

**Conflicts of Interest:** The authors declare no conflict of interest.

#### **References**

