Abstract
The malware attacks targeting the industrial control network are gradually increasing, and the nonlinear phenomenon makes it difficult to predict the propagation behavior of malware. Once the dynamic system becomes unstable, the propagation of malware will be out of control, which will seriously threaten the security of the industrial control network. So, it is necessary to model and study the propagation of malware in the industrial control network. In this paper, a SIDQR model with dual delay is proposed by fully considering the characteristics of the industrial control network. By analyzing the nonlinear dynamics of the model, the Hopf bifurcation is discussed in detail when the value of dual delay is greater than zero, and the expression for the threshold is also provided. The results of the experiments indicate that the system may have multiple bifurcation points. By comparing different immune and quarantine rates, it is found that the immune rate can be appropriately increased and the isolation rate can be appropriately reduced in the industrial control network, which can suppress the spread of malware without interrupting the industrial production.
1. Introduction
The industrial control network is the key foundation for realizing digital transformation. It is an emerging business form and application model formed by the deep integration of new information technology and the industrial economy. ICS (industrial control system) has become an important component of many national infrastructures. With the continuous upgrading of the ICS, the connection between the industrial control network and the Internet is increasingly close, which leads to the further increase of security risks.
The number of the industrial control network security accidents is increasing year by year. In 2010, the Bushehr nuclear power plant in Iran was attacked by the Stuxnet worm [1]. Since then, many malwares targeting the industrial control network have been discovered, such as Night Dragon [2], Flame [3], Duqu/Duqu2.0 [4], Blaster [5], and BlackEnergy [6]. The typical PLC (programmable logic controller) worm Blaster spread through the control system of Siemens SIMATIC S7-1200 without any PCs. PLC Blaster can scan the ICS networks to find new targets and then attack the PLC and complete self-replication in the infected PLC. Some new types of malware can spread not only in PC networks but also in PLCs in the industrial control network. In 2012, Russian security experts discovered the flame virus spreading widely in the energy industry in the Middle East [7]. In December 2015, Ukraine was attacked by hackers which led to a large-scale power outage [8]. The malware Industroyer [9], which was found in 2017, is aimed at key ICS and can lead to power outages. In December 2017, due to the zero-day vulnerability of Schneider’s Triconex SIS, a power plant in the Middle East was attacked and ultimately had to shut down [10]. The ransomware WannaCry can spread crazily globally in the same year and attack critical infrastructure [11]. In 2018, a chlorine gas station in Ukraine was attacked by VPNFilter virus [12]. In 2019, the power grid of Venezuela was attacked, and it led to large-scale power outages across the country [13]. The industrial control network has its own characteristics, for example, industrial control protocols are lack of built-in security mechanisms. Also, the processing capacity of ICS is weak, and the system update is lagging behind. Based on the above situation, in order to address the security issues faced by the industrial control network, it is necessary to understand the propagation patterns of malware in the industrial control network and propose appropriate containment strategies. Therefore, it is particularly necessary to model and analyze the propagation behavior of malicious software in industrial control networks.
Researchers have proposed some important epidemic models to explore the dynamic behavior of malware propagation. For example, in the SIS model [14–17] and the SIR model [18–22], the early studies were focused on ordinary networks. On the basis of these traditional models, researchers have proposed new models to study the propagation of different types of malware. Chen et al. analyzed the propagation behavior of malware in Bluetooth and mobile applications and proposed a malicious software propagation model in mobile networks [23]. Inspired by the SEIR model [24], Xiao et al. introduced a new state (i.e., quarantined state) in the epidemic model, which is a malware propagation model in WiFi environments [25]. Wang et al. proposed a microscopic mathematical model to describe the propagation behavior of malware in a sensor network and designed a LDS (local defense strategy) using mobile “patches” (mobile components that can distribute patches) [26]. In the malware propagation model, if time delay (such as patch release and quarantine) is not considered, the model is an ordinary differential dynamical system, and the above models all belong to this category. If such delay factors are considered, the resulting model is a delay differential dynamic system. Yao et al. considered the delay caused by IDS (intrusion detection systems), analyzed nonlinear phenomena, and proposed a threshold for bifurcation [27]. Ren et al. considered the delay factor on the basis of the SIR model and analyzed the stability conditions of the dynamic system [28]. Subsequently, they proposed an epidemic model with time-varying latency and analyzed the bifurcation phenomenon [29]. On this basis, Wang et al. analyzed and discussed the chaos phenomenon of time-delay models [30]. Feng et al. considered the situation of dual delay and antimalware measures and studied the Hopf bifurcation phenomenon in malware propagation [31]. Wang et al. studied the threshold problem of stability in dynamic systems when the propagation rate varies linearly [32]. Khan et al. investigated a discretized two-dimensional model, and the results for the existence and uniqueness, and conditions for local stability with topological classifications of the equilibrium solutions are determined [33]. Wang et al. investigated the selection mechanism of the minimal wave speed for traveling waves to an epidemic model, and a threshold is defined by the eigenvalue problem of the linearized system [34].
Currently, there is limited research on the spread of malware in the industrial control network. The network of real industrial control systems is relatively complex, and the nonlinear phenomena (such as bifurcation and chaos) in the propagation process of malware also make it difficult to predict its propagation behavior, which can also lead to the failure of containment strategies, and then, it will lead to system instability and hinder the normal operation of industrial production. So we will propose a malware propagation model for the industrial control network. The industrial control network is different from the Internet, the industrial control equipment generally do not have powerful processors like computers, and the bandwidth requirements of the industrial control networks are much lower than those of the Internet. At the same time, the quarantine of industrial control equipment also needs to take into account production activities. Therefore, when dealing with infected industrial control equipment, it is difficult to repair them as quickly as computers. So we introduce immune delay and quarantine delay into the model, and it can describe the influence of malware containment strategies on the propagation in the industrial control network more accurately. Currently, some researchers have conducted some research on the problem of dual delay. Zhang et al. proposed the conditions for the asymptotic stability of Hopf bifurcation with dual delay [35]. Fan et al. introduced the stability and bifurcation of a coupled HR model with dual delay [36]. He et al. proposed a neural network model with unidirectional coupling delay and discovered double Hopf bifurcations in this model [37].
Based on the above works, we consider the dual delay in the industrial control network and propose a new malware propagation model and the propagation behavior of malware, and the bifurcation phenomena is analyzed under different cases. The innovation of this model lies in the inclusion of two different delays, which makes it more suitable for the actual situation of the industrial control network. This model can provide a security defense strategy against the spread of malware for the industrial control networks without affecting industrial production as much as possible. In addition, in the industrial control networks, our research results demonstrate how to suppress the spread of malware while maintaining the stability of industrial control systems. The organization of the paper is as follows: Section 2 explains how the model is established in the industrial control network, Section 3 analyzes the stability of the equilibrium of the dynamic system, Section 4 presents the experimental results, and Section 5 is a conclusion.
2. Model Formulation
In actual industrial control networks, the propagation delay of malware objectively exists and can take various forms. For example, the delay caused by malware latency, and the delay caused by upgrading and patching the software and hardware of susceptible equipment. During the detection process, the time window mechanism is used for quarantine, and it will cause quarantine delay. Based on these characteristics, we propose a malware propagation model with dual delay. Namely, the delay is caused by upgrading and patching the software and hardware of susceptible equipment, and it is called immune delay. Another delay is caused by using the time window mechanism to quarantine infected equipment, which is referred to as quarantine delay. The time window mechanism can improve the accuracy of detection, so as not to affect the normal production of the factory due to false alarms. That is to say, when abnormal behavior is detected, an alarm will not be triggered immediately. Only when this abnormal behavior reaches a preset threshold, it will be considered an intrusion behavior, and an alarm will be issued. Time window mechanism will cause a delay before quarantine, and it will bring complex dynamic changes to the spread of malware. We will use the stability switching principle [38] to study the stability of the dynamic system with dual delay in the next section, and the assumptions for model formulation are listed as follows:(a)In our model, the industrial control network is assumed as a homogeneous network(b)It is assumed that the total number of all equipment remains unchanged, and the number is N(c)The equipment in the industrial control network has functions such as software upgrade and patching, and a time window mechanism is adopted
In our proposed SIDQR dual-delay model, each industrial control equipment may have six states: susceptible (S) state, infected (I) state, immune delay (D1) state, quarantine delay (D2) state, quarantine (Q) state, and recovery (R) state. The recovery rate of the susceptible equipment is , the quarantine rate of the infected equipment is , the infection rate of the susceptible equipment is , and the recovered rate of the quarantined equipment is . When facing new malware, there is a probability that recovered equipment will return to susceptible equipment. The transition diagram among the different states is shown in Figure 1. In summary, assuming the total number of all equipment is N, the differential equation system (1) of the SIDQR model with dual delay can be obtained, and the immune delay is , and the quarantine delay is .

3. Stability of Equilibrium
The stability of the equilibrium of system (1) is studied in this section, we focus on discussing the situation of , and we provide an expression for the threshold. For system (1), the following theorem can be obtained.
Theorem 1. System (1) has an unique positive equilibrium point when , where R0 is the basic reproduction number, and it means that the basic reproduction number is positive as an initial condition.
Proof. When system (1) is stable, that is, the left side of the differential equation system is equal to zero, and thus, the equilibrium point can be obtained:Since the total number of all equipment is N, the equation with as the root can be obtained asObviously, equation (3) has a unique positive real root and a unique positive equilibrium point . Then, , and system (1) can be simplified into the following form:The Jacobi matrix of the equilibrium point isThen, the following equation can be obtained:where .
Due to the existence of two delays in this model, in the process of solving characteristic equations, it is necessary to reduce the order and obtain the algebraic cofactors. The first algebraic cofactor isthe second algebraic cofactor isand the third algebraic cofactor isThen, we calculate the three cofactors separately and add them together to obtain the characteristic equation of the system:whereFor system (1), the stability of the equilibrium point needs to be discussed in different cases. The cases , , and are essentially the same as the single delay model discussed in paper [27–29, 32], and we will not repeat the proof here. Our main analysis here is the stability of the equilibrium in the case of . In this case, the stability analysis of the system requires fixing the value range of one delay within the threshold, that is, or ; at this point, the variable can be considered as one of the delays ( or ). Let us take the case of as an example to discuss here. The root of the system characteristic equation is . By substituting it into equation (10) and separating the real and imaginary parts, we can obtainBy combining equation (12) and (13), it can be obtained thatwhereAssuming that equation (15) has finite positive roots , using the Routh–Hurwitz criterion, for each value of k, the corresponding threshold of delay iswhere k and j are both positive real numbers, let and the transversality condition holds, and the following conclusion can be obtained based on Rouche’s theorem:
Theorem 2. Assuming that , then the following can be obtained, and it is quoted in reference [27, 39].(1)When , for system (1), the positive equilibrium point is locally asymptotically stable. When , the positive equilibrium point is unstable.(2)When system (1) satisfies , the positive equilibrium point of the system will undergo a Hopf bifurcation at , and the system will lose stability.(3)In equation (16), k and j are both positive real numbers, and the value of is also affected by . Therefore, when , the system may have multiple bifurcation points.
4. Experiments
In order to demonstrate the impact of dual delay on the propagation of malware, numerical experiments analysis is conducted in this section. The infection rate is assumed that , the recovered rate of the susceptible equipment is assumed that , the quarantine rate of the infected equipment is assumed that , the recovered rate of the quarantined equipment is assumed that , and the probability of the recovered equipment becoming susceptible to infection is . At the initial stage, the total number of equipment is 10000, assuming that the number of infected equipment (I) is 50 and the other equipment were susceptible (S). Due to the existence of a double delay, this section first presents the overall nonlinear phenomenon of the model through the bifurcation diagram.
4.1. Case 1 with
The first to discuss is the bifurcation phenomenon of the system when , as shown in Figure 2. The dynamic system in this case is equivalent to a single delay situation, where Hopf bifurcation occurs when the delay value exceeds the threshold. That is to say, in this situation, it is necessary to control the quarantine delay in the ICS to ensure that the malware will not get out of control. The quarantine delay needs to be less than the threshold, so that the dynamic system will eventually reach equilibrium after oscillation. Also, the curves of the equipment in different states are shown in Figure 3 () and Figure 4 (), which also validate Theorem 5 of the single delay model in paper [27].



4.2. Case 2 with
What needs to be discussed next is the bifurcation phenomenon of the system when . As shown in Figure 5, the Hopf bifurcation diagram in this case is different from the ordinary single delay bifurcation. When the value of delay exceeds the threshold, a curve is in a fluctuating state, but it has no effect on the system bifurcation, and the entire dynamic system is still in a bifurcation state. This can indicate that in this dual-delay model, has a greater impact on the dynamic system, which means that immune delay will make the nonlinear phenomena of the system more complex. Figures 6 and 7 show the curves of the number of equipment at and , respectively. The propagation process shown in Figures 6 and 7 can verify the results of Figure 5.



When , that is, the quarantine delay is zero, the dynamic system may be in a bifurcation state. When the immune delay is less than the threshold value ( < 1100), the system will be stable, and the curves finally reach equilibrium and no longer fluctuate, as it is shown in Figure 6. In Figure 7, the number of equipments in different states fluctuates continuously over time, indicating that the system cannot be controlled. In this situation, the number of infected equipment and the spread trend of malware become difficult to predict. This will pose a serious threat to the equipment in the ICS.
4.3. Case 3 with
When , it can be seen from the two situations mentioned above that immune delay has a significant impact on the stability of the system. Therefore, let us take to further discuss the impact on the system. Figure 8 shows the Hopf bifurcation diagram of the system in this case. It can be seen that when neither delay is zero and the value of is fixed, then the system will have multiple bifurcation points. This indicates that the system will experience bifurcation and then return to a stable state, and then, as the value of increases, the system will experience bifurcation again, which is consistent with the content of Theorem 2 in Section 3. To verify the result, the value of will be taken within the red circle range in Figure 8, and the curves of different state equipment are also shown in Figures 9–11.




In Figures 9–11, the quarantine delay of the system is fixed at , and the value of is 255, 300, and 400, respectively. It can be seen that as the value of increases, the system undergoes a process of stability, bifurcation, and then stability. This is completely consistent with the results in the red circle in Figure 8. It indicates that in a dual-delay system, it is necessary to clarify the impact of different delays on system stability. For example, in Case 3, we cannot simply demand that the immune delay be as small as possible, because it may also lead to the bifurcation. It requires us to specifically analyze the impact of immune delay on the dynamic system, accurately control it to ensure system stability, and suppress the spread of malware.
The experiments above demonstrate how to maintain the stability of the dynamic system by controlling the value of immune delay and quarantine delay. In addition, from equation (16), it can be seen that the stability of the dynamic system can also be maintained by adjusting the immune rate and quarantine rate , that is, by changing the parameters so that the delay is less than the threshold. Figure 12 shows the curve of infected equipment under different immune rates with . It can be seen that as the immune rate increases, the system changes from a bifurcation state to a stable state, and the number and peak value of infected equipment also decrease. It indicates that increasing the immune rate can maintain system stability and control the spread of malware.

Figure 13 shows the curve of infected equipment under different quarantine rates with . From Figure 13, it can be observed that when the quarantine rate , the curve of infected equipment continues to fluctuate, and the system is in a bifurcation state. As the value of the quarantine rate decreases, the number of infected equipment will gradually reach a stable state. From the results, it can be seen that the spread of malware cannot be controlled solely through a large number of quarantine equipment, which may lead to instability and also have a significant impact on industrial production. In summary, in the protection of the ICS, the value of the immune rate can be appropriately increased and the value of the isolation rate can be appropriately reduced, which can ensure the stability of the system and suppress the spread of malware.

5. Conclusion
We fully consider the characteristics of the ICS and equipment in the industrial control network and propose a malware propagation model with dual delay. On this basis, we study the propagation of malware in the industrial control network and the stability and Hopf bifurcation of the dynamic system. Moreover, the containment strategy for malware in the industrial control network is proposed. In particular, the following conclusions can be obtained:(1)In the industrial control network, the characteristics of immunity and quarantine in actual industrial production are considered, and the SIDQR model with dual delay is established. The model includes six states of industrial equipment: susceptible (S) state, infected (I) state, immune delay (D1) state, quarantine delay (D2) state, quarantine (Q) state, and recovery (R) state.(2)The positive equilibrium point of the system is proven with Jacobian matrix and reduced order, the stability of the dynamic system at is discussed in detail, and an expression for the threshold is provided. When the delay exceeds the threshold, the system becomes unstable and Hopf bifurcation occurs. Also, the system may have multiple bifurcation points at .(3)Under different cases, the experiments demonstrate the propagation of malware in the industrial control networks, and the possible bifurcation points of the system in different cases are shown. In addition, after comparing different immune and quarantine rates, the experimental results show that the immune rate can be appropriately increased and the quarantine rate can be appropriately reduced, which can ensure the stability of the ICS and suppress the spread of malware.
However, when the industrial control network encounters cross-network attacks, the model is not applicable. In addition, the importance of different equipment in the industrial control system also varies, and the key equipment such as SCADA servers is likely to play the role of key nodes. If malware prioritizes attacking these critical nodes, the propagation trend will undergo significant changes. How will the effectiveness of the containment strategy change when prioritizing the control of these key nodes during the defense process? These works need to be improved in subsequent research.
Data Availability
No underlying data were collected or produced in this study.
Conflicts of Interest
The authors declare that there are no conflicts of interest regarding the publication of this paper.
Acknowledgments
This paper was supported by National Key R&D Program of China under Grant no. 2021YFB3101700, Applied Basic Research Program of Liaoning Province under Grant no. 2022JH2/101300240, Fundamental Research Funds for the Central Universities under Grant no. N2324004-12, and Scientific Research Project of Liaoning Province Education Department under Grant no. LJKQZ20222457.