Malware attacks pose a serious security risk, which threatens our ever-expanding wireless sensor networks (WSNs). The characteristic of malware propagation in WSNs is similar to how an epidemic spreads among humans. The malware at each infective node may seek to contract more susceptible nodes by amplifying the transmission range and the media scanning rate, thereby accelerating its spread. Kephart et al [1] originally applied the mean-field theory to study the modeling of malware propagation.

Since then, a number of mathematical epidemic models, ones have looked into dynamic behaviors of malware propagation. In the past few decades, many researchers focused on the Susceptible-Infectious-Recovered (SIR) epidemic model to describe behaviors of malware propagation, such as Youssef et al. [2], Feng et al. [3], and Rey et al [4]. In their works, all nodes are divided into the susceptible node, infectious node, and recovered node. The interaction of these three classes is governed by the following model, consisting of three non-linear differential equations: d S t d t = − γ S t I t , d I t d t = γ S t I t + μ I t , d R t d t = − μ I t , (1) where S(t), I(t), and R(t) are the densities of the susceptible node, infectious node, and recovered node at time t, respectively, and γ and μ are the multiplication rate and the recovered rate, respectively.

In recent years, many malware propagation models based on the classical SIR model have been developed. For example, Xiao et al. [5] established a SIRS model of malware propagation with diffusion and time delays. They forecasted the occurrence of tipping points and studied the tipping dynamics due to the Turing instability and Hopf bifurcation. Dong et al. [6] proposed a fractional network SIRS epidemic model with fuzzy transmission and saturated treatment function. Carnier et al. [7] generalized the methodology of derivation of the exact Markov chain for any malware model based on the simplest compartmental model SIR model. In those proposed models, the susceptible, exposed, infected, and recovered (SEIR) model is the most adopted to characterize the spreading of malware. In addition, calculation of the basic reproduction numbers and local stability of a non-trivial equilibrium and an endemic equilibrium are widely concerned in mathematical epidemiology. For example, Prajapati et al. [8] proposed an epidemic model to describe the spread of malicious objects in the network due to removable devices. Equilibrium points, both endemic and malware free, were obtained, and they also formulated a reproduction number. Liu et al. [9] proposed a delayed e-epidemic SEIRS malware propagation model with a generalized non-monotone incidence rate. Shakya et al. [10] proposed the correlation-based SIR model, which takes into account the spatial correlation characteristics of WSNs. The basic reproduction number was derived, and the local stability and Hopf bifurcation analysis were performed in their works. Moreover, Liu et al. [11] proposed a distributed continuous-time model in which two competing viruses spread over a network. The unique equilibrium and global stability were performed in their work. Zhang et al. [12] proposed an e-epidemic time-delay epidemic model to study the appearance of delay dynamics and performed non-linear stability analysis, Hopf bifurcation analysis, and an analysis of its stability. Due to the exposed nodes in WSNs, they may have different infection rates during malware propagation, Yu et al. [13] proposed the improved SEIR model with two infectious rates in the cyber-physical systems to explore the transition mechanism. Dmitriy et al. [14] divided infectious nodes into susceptible-infected by strain 1 or by strain 2-susceptible with duty cycles. The basic reproduction number and the stability were analyzed in their work. Dong et al. [15, 16] proposed the improved SEIR model with two distinct compartment exposed nodes. The dynamics of the network-based fractional order epidemic model were studied in their work. Nwokoye et al. [17] combined the strength of the above models and proposed a multi-group model, which represents multiple exposed infections due to worms and viruses. In addition to the SEIR model, some researchers have tried to divide nodes into more classes. For example, Ojha et al. [18] proposed an improved epidemic model that aggregates quarantine and vaccination techniques. The basic reproduction number and equilibrium points were analyzed in their work. Hosseini et al. [19–21] proposed a new dynamic model of malware propagation in heterogeneous networks based on the rumor diffusion model. In their works, all the nodes were divided into six classes: susceptible nodes, exposed nodes, infectious nodes, recovered nodes, vaccinated nodes, and quarantine nodes. In addition to their characteristic of the spreading over the network, the optimal strategies are also a concern for many researchers, such as Muthukrishnan et al. [22], who investigated an optimal control strategy to reduce malware propagation in WSNs. Liu et al. [23] extended the traditional SIR model by adding another delitescent compartment to address the behavior of malware. The optimal control theory was employed to study malware immunization strategies. Moreover, Nwokoye et al. [24] reviewed the epidemic models of malware propagation and control in WSNs. Jain et al. [25] introduced an optimal control of rumor spreading in a homogeneously mixed population to minimize the density of rumor adopters and control cost.

The majority of previous studies have focused on modeling malware propagation but ignored the characteristics of WSNs. How to design efficient control strategies to reduce the spread of malware is little discussed. For example, the malware transmission rate is a parameter related to the communication radius of the node communication volume; the node recovery rate can be realized by a node system update. Therefore, combining node characteristics can optimize the wireless sensor network by reducing the propagation radius and increasing the update frequency. To overcome such weaknesses, an improved malware propagation model is introduced in this article. The basic reproduction number is derived in detail, which is the main contribution we make. The novelty of this work is to explore malware propagation mechanisms in networks through the proposed SUIQR model and finding the optimized strategies for controlling the spread of malware based on the reproduction number. The rest of the article is organized as follows. The main contribution of this article is briefly summarized in Section 2. A new difference model is proposed, called a SUIQR model. Then, the basic reproduction number of the proposed model is derived, and local stability is analyzed in Section 3. In Section 4, optimal strategies are discussed to control the spread of malware propagation. Mathematical results are illustrated by numerical simulations, and some control strategies are given in Section 5. Finally, the conclusions and further questions are presented in Section 6.

In this section, a new differential equation epidemic model is built to describe the dynamic behaviors of malware propagation. The total nodes in WSNs are divided into five categories: susceptible, unexposed, infectious, isolation, and removed, called SUIQR. The relationship of these categories is shown in Figure 1, and the meaning of all the symbols in Figure 1 is shown in Table 1.

The susceptible nodes will be infected by the unexposed nodes and infectious nodes at the rate of γ. There are γS(U + I) nodes, which will remove from class S into class U. Due to the system being updated at a certain period, the nodes in WSNs will enter the class R at the rate of μ. The recovered nodes will enter class S at the rate of τ due to protection failure [26, 27]. Therefore, the transfer relationships between S and other classes can be expressed as follows: d S d t = − γ S U + I + τ R − μ S . (2)

The density of the unexposed nodes will increase when the susceptible nodes are infected. Some unexposed nodes will be found as the isolated class Q at the rate of η through regular checking. Some unexposed nodes will move into class I before being recovered at the rate of β. Therefore, the transfer relationships between U and other classes can be expressed as follows: d U d t = γ S U + I − μ U − β U − η U . (3)

The density of the infectious nodes will increase because some unexposed nodes will move into class I. Some nodes in class I will remove into class R due to the periodic system updates at the rate of μ. Therefore, the transfer relationships between I and other classes can be expressed as the following equations: d I d t = β U − μ I , (4)

The density of the isolation nodes will increase when the unexposed nodes are founded. Some nodes in class Q will remove into class R due to the periodic system updates at the rate of μ. Therefore, the transfer relationships between Q and other classes can be expressed as the following equations: d Q d t = η U − μ Q . (5)

Assuming the total number of the nodes in WSNs is a constant. Therefore, the transfer relationships between R and other classes can be obtained as follows: d R d t = μ S + μ U + μ I + μ Q − τ R . (6)

Combining the above ideas, the SUIQR model governing the transmission of malware is described by the following system of non-linear differential equations: ∂ S ∂ t = − γ S U + I + τ R − μ S , ∂ U ∂ t = γ S U + I − μ U − β U − η U , ∂ I ∂ t = β U − μ I , ∂ Q ∂ t = η U − μ Q , ∂ R ∂ t = μ S + μ U + μ I + μ Q − τ R . (7)

The basic reproductive number is a key parameter to represent the infected numbers in an average infection period. Firstly, the global stability of the malware-free equilibrium (MFE) is introduced [28]. We can get the SUIQR model 7) to have a malware-free equilibrium point P f * and the malware-endemic equilibrium point P e * by solving the following Eq. 8: − γ S × U + I + τ R − μ S = 0 , γ S × U + I − μ U − β U − η U = 0 , β U − μ I = 0 , η U − μ Q = 0 , S + U + R + I + Q = 1 . (8)

The malware-free equilibrium point P f * and the malware-endemic equilibrium point P e * are shown in Eq.9 and Eq.(10), respectively: P f * S f * , U f * , I f * , Q f * , R f * = τ μ + τ , 0,0,0 , μ μ + τ . (9) P e ∗ S e ∗ , U e ∗ , I e ∗ , Q e ∗ , R e ∗ = β μ + η μ + μ 2 β γ + γ μ , − μ β μ 2 + η μ 2 + μ 2 τ + μ 3 − β γ τ + β μ τ + η μ τ − γ μ τ μ + τ γ β 2 + 2 γ β μ + η γ β + γ μ 2 + η γ μ , − β β μ 2 + η μ 2 + μ 2 τ + μ 3 − β γ τ + β μ τ + η μ τ − γ μ τ μ + τ γ β 2 + 2 γ β μ + η γ β + γ μ 2 + η γ μ , − η μ 3 + η 2 μ 2 + η μ 2 τ + η 2 μ τ + β η μ 2 − β η γ τ + β η μ τ − η γ μ τ μ + 1 γ β 2 + 2 γ β μ + η γ β + γ μ 2 + η γ μ , μ μ + τ (10)

Then, the next-generation method (NGM) is applied to calculate the basic reproductive number R ₀ [29]. The main advantage of the NGM is that it allows the research to ignore any uninfected classes and focus only on the infected classes. There are three infected classes in the proposed SUIQR model. Let X = (U,I,Q) ^T , the model 7) equals to the following form: d X d t = F − V . (11)

Where F = γ S ( U + I ) 0 0 and V = μ U + β U + η U − β U + μ I − η U + μ Q .

We define f and v as the Jacobian matrices of F and V evaluated at the malware-free equilibrium point P f * : f = ∂ F ∂ X | P f * = γ S f * γ S f * 0 0 0 0 0 0 0 = τ γ μ + τ τ γ μ + τ 0 0 0 0 0 0 0 . (12) v = ∂ V ∂ X | P f * = μ + β + η 0 0 − β μ 0 − η 0 μ . (13)

The basic reproductive number R ₀ is the largest eigenvalue of the matrix fv ⁻¹ given by Eq. 14: f v − 1 = β + μ τ γ μ μ + τ β + η + μ γ τ μ μ + τ 0 0 0 0 0 0 0 . (14)

Then, we can obtain the basic reproductive number R ₀. R 0 = ρ f v − 1 = β + μ τ γ μ μ + τ β + η + μ . (15)

When R ₀ < 1, it means that the SUIQR model will reach a malware-free situation eventually; when R ₀ > 1, the SUIQR model always has malware and has an endemic equilibrium.

Finally, we analyze the sensitivity of each parameter (μ, η, β, λ) about the basic reproductive number and obtain normed forward sensitivity indexes as follows: ∂ R 0 R 0 ∂ λ λ = 1 > 0 , ∂ R 0 R 0 ∂ β β = η β β + μ β + μ + η > 0 , ∂ R 0 R 0 ∂ μ μ = τ γ μ 2 + μ τ β + η + μ − τ γ β + μ 3 μ 2 + 2 μ β + η + τ + τ β + τ η μ + τ β + η + μ β + μ τ γ < 0 , ∂ R 0 R 0 ∂ η η = − η β + μ + η < 0 . (16)

It is obvious that λ, β and R ₀ are proportional, η, μ, and R ₀ are inverse proportional. The increase of η or μ may result in the decrease of R ₀. The decrease of λ and β may result in the decrease of R ₀. Above all, we can take some measures to control the spread of malware.

• Reduce the period of system updates. This can improve the defense capabilities of nodes, that is, to increase μ.

The Jacobian matrix of the SUIQR model at P e * is as follows: J = − γ U e ∗ + I e ∗ − μ − γ S e ∗ − γ S e ∗ 0 τ γ U e ∗ + I e ∗ γ S e − μ − β − η γ S e 0 0 0 β − μ 0 0 0 η 0 − μ 0 μ μ μ μ − τ . (17)

The corresponding characteristic equation is shown in Eq. 18: λ 5 + A 4 λ 4 + A 3 λ 3 + A 2 λ 2 + A 1 λ + A 0 = 0 . (18) A 4 = β + η + 4 μ − τ + γ I e ∗ − γ S e ∗ + γ U e ∗ , A 3 = 3 β μ + 3 η μ − β τ − η τ − 5 μ τ + 6 μ 2 + β γ I e ∗ + η γ I e ∗ − β γ S e ∗ + 3 γ μ I e ∗ + β γ U e ∗ + η γ U e ∗ − γ τ I e ∗ − 3 γ μ S e ∗ + 3 γ μ U e ∗ + γ τ S e ∗ − γ τ U e ∗ , A 2 = 3 β μ 2 + 3 η μ 2 − 9 μ 2 τ + 4 μ 3 − 4 β μ τ − 4 η μ τ + 3 γ μ 2 I e ∗ − 3 γ μ 2 S e ∗ + 3 γ μ 2 U e ∗ + 2 β γ μ I e ∗ + 2 η γ μ I e ∗ − β γ τ I e ∗ − 2 β γ μ S e ∗ − η γ τ I e ∗ + 2 β γ μ U e ∗ + 2 η γ μ U e ∗ + β γ τ S e ∗ − 4 γ μ τ I e ∗ − β γ τ U e ∗ − η γ τ U e ∗ + 4 γ μ τ S e ∗ − 4 γ μ τ U e ∗ , A 1 = β μ 3 + η μ 3 − 7 μ 3 τ + μ 4 − 5 β μ 2 τ − 5 η μ 2 τ + γ μ 3 I e ∗ − γ μ 3 S e ∗ + γ μ 3 U e ∗ + β γ μ 2 I e ∗ + η γ μ 2 I e ∗ − β γ μ 2 S e ∗ + β γ μ 2 U e ∗ + η γ μ 2 U e ∗ − 5 γ μ 2 τ I e ∗ + 5 γ μ 2 τ S e ∗ − 5 γ μ 2 τ U e ∗ − 3 β γ μ τ I e ∗ − 3 η γ μ τ I e ∗ + 3 β γ μ τ S e ∗ − 3 β γ μ τ U e ∗ − 3 η γ μ τ U e ∗ , A 0 = 2 γ μ 3 τ S e ∗ − 2 β μ 3 τ − 2 η μ 3 τ − 2 γ μ 3 τ I e ∗ − 2 μ 4 τ − 2 γ μ 3 τ U e ∗ − 2 β γ μ 2 τ I e ∗ − 2 η γ μ 2 τ I e ∗ + 2 β γ μ 2 τ S e ∗ − 2 β γ μ 2 τ U e ∗ − 2 η γ μ 2 τ U e ∗ . (19)

According to the Routh-Hurwitz criterion [27], when R ₀ > 1, the endemic equilibrium P e * S e * , U e * , I e * , Q e * , R e * is asymptotically stable.

We assume that the distribution of the nodes follows the Poisson point process in WSN. The probability of there being k nodes in the communication range can be calculated by the following equation (30): P X = k = λ π R 1 2 k e − λ π R 1 2 k ! . (20) Where λ is the density of the nodes in WSNs, and R ₁ is the communication range of each node, which is always an adjustable parameter.

Therefore, the conversion rate that susceptible nodes contract with the other nodes can be calculated by the following equation: γ = E X = λ π R 1 2 . (21)

The sensitivity of the communication range about the basic reproductive number and obtaining a normed forward sensitivity index is shown as Eq. 22: ∂ R 0 R 0 ∂ R 1 R 1 = 2 R 1 > 0 . (22)

According to the analysis above, we need to adjust the parameters (μ or R ₁) to keep the basic reproductive number at less than 1. By substituting Eq. 21 into Eq. 15, the relationship between μ and R ₁ is given by the following inequality: β + μ τ λ π R 1 2 μ μ + τ β + η + μ < 1 . (23)

•Case 1. R ₁ is the adjustable parameter and μ is the fixed constant:

Based on the above inequality, we set the threshold of the communication range. When the range is less than the threshold, the SUIQR model reaches a malware-free situation eventually.

•Case 2. μ is the adjustable parameter and R ₁ is the fixed constant:

In this section, some numerical simulations are performed to illustrate and complement our analytical results. There are 100 nodes in the networks with the area of 1000 × 1000. We fix the parameters β = 0.6, λ = 0.0001, η = 0.3, μ = 0.1, R ₁ = 50, τ = 0.3. The initial density of each class is S(0) = 0.6, U(0) = 0.1, I(0) = 0.1, Q(0) = 0.1, R(0) = 0.1. The evolutions of each class along with time are shown in Figure 2.

In this example, we can calculate the basic reproductive number R ₀ = 4.1233 based on Eq. 15 As can be seen from Figure 2, the reproduction number is determined by R ₀ > 1. Here, we can see the densities of class S and class U decrease rapidly over time and eventually stabilize the equilibrium point P e * , with a high rate of conversion and recovery. The densities of class I and class R increase rapidly over time and eventually stabilize the equilibrium point. The malware-endemic equilibrium point P e * = [ 0.1819 , 0.0568 , 0.3409 , 0.2500 ] . The malware will exist in this WSN. The corresponding characteristic equation is shown as follows: λ 5 + 2.5406 λ 4 + 2.0997 λ 3 + 0.6419 λ 2 + 0.0692 λ + 0.0024 = 0 (26)

The roots of the characteristic equation are [-1.1874, −0.8467, - 0.3370, −0.1000, −0.0695]. The Routh-Hurwitz table is given by: H = 1.0000 2.0997 0.0692 2.5406 0.6419 0.0024 1.8470 0.0683 0 0.5480 0.0024 0 0.0603 0 0 0.0024 0 0 .

According to the Routh-Hurwitz criterion, this system is stable based on the given parameters.

Then, we analyze the influence of the optimal control strategy on the evolutions of each class. In the above example, the malware will not be eliminated. When only one parameter is adjustable, R ₁ can be calculated by Eq.24, R ₁ < 24.6233 if R ₀ < 1. We fixed R ₁ = 24, and the evolutions of each class are shown in Figure 3.

In this example, we can calculate the basic reproductive number R ₀ = 0.9500 based on Eq. 15 As can be seen from Figure 3, the reproduction number is determined by R ₀ < 1. Here, the densities of class S, class U, and class Q can be seen to increase rapidly and decay slowly to 0, with a high recovery rate. The density of class S decreases rapidly and then slowly increases to eventually stabilize the equilibrium point. The malware-free equilibrium point P f * = [ 0.75 , 0,0,0,0.25 ] . The corresponding characteristic equation is shown as follows: λ 5 + 1.6188 λ 4 + 0.5099 λ 3 + 0.2729 λ 2 − 0.0012 λ − 0.00003 = 0 (27)

The roots of the characteristic equation are [-1.2197, −0.3179, −0.1000, 0.0391, −0.0202]. The Routh-Hurwitz table is given by: H = 1.0000 0.5099 − 0.0012 1.6188 0.0273 − 0.0000 0.4930 − 0.0011 0 0.0310 − 0.0000 0 − 0.0007 0 0 − 0.0000 0 0 .

According to the Routh-Hurwitz criterion, this system is stable based on the given parameters.

When only one parameter is adjustable, μ can be calculated by Eq. 24, μ > 0.2961. We fixed the communication range of each R ₁ = 0.3, and the evolutions of each class are shown in Figure 4.

In this example, we can calculate the basic reproductive number R ₀ = 0.9817 based on Eq. 15 As can be seen from Figure 4, it is observed that the reproduction number is determined by R ₀ < 1. Here, the densities of class S, class U, and class Q can be seen to increase rapidly and decay slowly to 0, with a high recovery rate. The density of class S decreases rapidly and then slowly increases to eventually stabilize the equilibrium point. The SUIQR model will reach a malware-free situation eventually, and the malware-free equilibrium point P f * = [ 0.5 , 0,0,0,0.5 ] . The corresponding characteristic equation is shown as follows: λ 5 + 3.2799 λ 4 + 3.3250 λ 3 + 1.2538 λ 2 − 0.1903 λ − 0.0099 = 0 (28)

The roots of the characteristic equation are [-1.7649, −0.9278, −0.3000, −0.1658, −0.1214]. The Routh-Hurwitz table is given by: H = 1.0000 3.3250 0.1903 3.2799 1.2538 0.0099 2.9427 0.1873 0 1.0450 0.0099 0 0.1594 0 0 0.0099 0 0 . According to the Routh-Hurwitz criterion, this system is stable based on the given parameters.

We set the update period of the system μ, which satisfies the above inequality to ensure the SUIQR model will reach a malware-free situation eventually. As shown in Figures 3, 4, the system is optimal and stable when the R ₁ or μ is set. We can obtain the control strategy to reduce the spread of malware by setting reasonable parameters for nodes in WSNs.

Sometimes, the results of adjusting the single parameter μ or R ₁ are not satisfactory. In these scenes, we need to adjust μ and R ₁ simultaneously; the relationship between μ and R ₁ is given in Eq. 25 Based on the given parameters, the relationship between the two parameters is shown in Figure 5.

As shown in Figure 5, the region below is feasible, while the region above is not. We can adjust these two parameters in the curve to ensure the SUIQR model reaches a malware-free situation eventually. A schematic of the results of controlling both the communication radius and the frequency of system updates to the system is shown in Figure 5. The previous control method highlights that malware propagation in networks can be simulated regardless of whether the radius of communication is decreased or the frequency of updates to the system is increased. In the curve shown in Figure 5, the entire space is divided into two parts. The upper part of the space is not feasible, that is, the control strategy in the upper part of the space cannot achieve the purpose of suppressing the spread of malware. The second half of the space is the feasible part, that is, the control strategy in the second half of the space can effectively suppress the malware. The points on the curve represent the critical values and are the most up-to-date regulatory solution. If the current control strategy can fail to achieve the goal of suppressing malware propagation, control can be achieved by adjusting the radius of communication at a triangle checkpoint or by increasing the frequency of system updates at the circle checkpoint.

This article proposed a SUIQR epidemic model to describe the spreading of malware. The basic reproductive number was derived by using the next-generation method. Finally, we carried out numerical simulations to observe the malware spreading in WSNs to verify the obtained theoretical results. Furthermore, we also investigated the communication range of the nodes to make the results more complete. As a result, we can adjust the communication range and the period of system updating to reduce the spreading of malware in WSNs. Due to the lack of experimental data, only numerical simulations were performed in this article, which is also a limitation of the work. In our further work, we will continue to investigate and perform relevant experiments and compare them with other infectious disease models to verify that the method is feasible.