Dynamical Behavior Analysis of a Time-Delay SIRS-L Model in Rechargeable Wireless Sensor Networks

: The traditional SIRS virus propagation model is used to analyze the malware propagation behavior of wireless rechargeable sensor networks (WRSNs) by adding a new concept: the low-energy status nodes. The SIRS-L model has been developed in this article. Furthermore, the inﬂuence of time delay during the charging behavior of the low-energy status nodes needs to be considered. Hopf bifurcation is studied by discussing the time delay that is chosen as the bifurcation parameter. Finally, the properties of the Hopf bifurcation are explored by applying the normal form theory and the center manifold theorem.


Introduction
With the development and application of communication technology and sensor technology, wireless sensor networks (WSNs) have come into being. WSNs are widely used globally, such as in the military industry, agricultural production, ecological monitoring, disaster warning, and infrastructure status monitoring [1]. Due to the self-organized network of WSNs, sensor nodes are vulnerable to various attacks in the process of information transmission [2][3][4][5][6][7][8][9][10]. Thus, it is of practical significance to prevent and deal with malicious attacks in WSNs.
The characteristics of WSNs can be listed as follows: adapting wireless communication, the topology being able to be changed dynamically, the security mechanism being imperfect, and the energy of nodes being limited, etc. Thus, the security of wireless sensor networks is facing great challenges. The attacks against sensor network nodes can be roughly divided into the following types: network worms, bots, and Trojans, among others [11][12][13][14]. Most malware attacks in WSNs mainly aim at paralyzing them or losing their power. In order to overcome the defect of short life cycle of WSNs, wireless rechargeable sensor networks (WRSNs) have been developed and have caught more attention. Consequently, a relatively special attack has appeared to against the nodes' power, namely, the denial of charge (DOC) [15]. The prevention and control of malware spread in WRSNs have become a hot issue.
The propagation mechanism of malware is similar to that of viruses in biological groups. Many scholars have made great contributions to malware propagation dynamics. In the beginning, the spread of malware mainly focuses on the traditional Internet. For example, Zou et al. put forward the worm propagation models represented by SI (susceptible-infected), SIS (susceptible-infected-susceptible), and SIR (susceptibleinfected-recovered) [16][17][18][19]. These studies lay a foundation for the later studies on the propagation behavior of malware in WSNs. However, it is worth noting that there are some differences between WSNs and the traditional Internet: (1) Because of their topological structure, the communication range of WSNs has physical limitations. (2) The limited energy characteristics lead to uneven communication capabilities. (3) The link-layer access conflict and avoidance mechanism of WSNs sharing wireless channel also introduces spacetime correlation for the propagation dynamics of WSNs. (4) WSNs have a high degree of self-organization. Therefore, it is of great significance to study the malware propagation model in WSNs, which has been proposed in recent years. For example, Tanachaiwiwat began to study the propagation behavior of network worms in 2009 [20]. Khayam et al. used signal processing technology to study the propagation model of worms in 2006 [21]. Furthermore, the SEIR (susceptible-exposed-infected-recovered) model is adapted to study the spread of malware [22], while the SEIRS-V (susceptible-exposed-infected-recoveredsusceptible and vaccinated) model explores the application of Hopf bifurcation during the malware spreading [23,24]. With the development, the research on WRSNs was pushed forward recently [25][26][27][28][29][30][31], considering the new concept of L (low-energy status) nodes. It is worth noting that the transition from the low-energy status nodes to the normal status nodes is called charging. At present, there are many charging methods for WRSNs, such as unmanned aerial vehicle (UAV) charging [32].
The related techniques and modeling basis of this paper mainly refer to the papers in Table 1: Table 1. Related research on the modeling and analysis of this paper.

Model Characteristics Reference Content
Zhu et al. [33] SIRS (susceptible-infectedrecovered-susceptible) The authors put forward the time delay of the immune validity; the SIRS model is applied to WSNs analysis.
SIRS model is taken as the premise of modeling; the SIRS-L model considering the low-energy state nodes is proposed in this paper.
It provides the analysis reference of Hopf bifurcation and the corresponding mathematical processing method.
Liu et al. [25] SIS-L (susceptible-infected-susceptiblelow-energy status) It first proposes the low-energy state nodes and combines them into the research of WSRNs.
It provides a theoretical basis for the low-energy status modeling.
Liu et al. [26] SIAS-L (susceptible-infected-antimalware-susceptible-lowenergy status) The status of anti-malware is proposed and the optimal strategy is considered.
It provides a theoretical basis for the low-energy status modeling.
Liu et al. [29] SIS-L (susceptible-infected-susceptiblelow-energy status) Time delay is considered for the first time in the model with the low-energy state nodes. However, the bifurcation is not discussed.
It provides a theoretical reference for time-delay analysis and the feasibility in modeling.
In order to promote the study of virus dynamics in WRSNs, the modeling in this paper is mainly based on the SIRS (susceptible-infected-recovered-susceptible) infectious disease model. Combined with the charging delay in WRSNs, the bifurcation in the process of charging needs to be considered. Therefore, the general contributions of this paper are as follows: 1.

2.
The equilibrium solutions of the SIRS-L model are obtained, and the basic reproductive number R0 is defined through the regeneration matrix [34].

3.
Revealing of the stability of the SIRS-L model when the charging delay is ignored.

4.
The variation of the solutions of the characteristic equation are discussed if the charging delay is considered through the theory in [35], and the occurrence conditions of Hopf bifurcation are figured out.
The properties of the Hopf bifurcation are explored by applying the normal form theory and the center manifold theorem [36].
The content of this paper is organized as follows: In Section 2, we establish the SIRS-L (susceptible-infected-recovered-susceptible-low-energy) model combined with the concept of low-energy nodes and the charging delay; in Section 3, the equilibrium points analysis and the corresponding Hopf bifurcation conditions are presented; the properties of the Hopf bifurcation are given in Section 4; the simulation analysis are revealed in Section 5; the conclusions are given in Section 6.

Modeling
On the basis of the classical SIRS model, we obtained the SIRS-L model as follows (1), with the parameters being explained in Table 2: (1) Table 2. Interpretation of the model parameters.

S(t)
The Susceptible Nodes The infected nodes The recovered nodes

LS(t)
The low-energy status susceptible nodes

LI(t)
The low-energy status infected nodes

LR(t)
The low-energy status recovered nodes Immune failure rate of the recovered nodes; the recovered nodes will be re-exposed to the malware and may be infected again β 1 Low-energy node conversion rate, which is to describe the process of the general nodes dropping to the low-energy nodes Node deactivation rate The following formula is obtained through the injection rate and the deactivation rate: . Therefore, we can obtain N(∞) = ∧ b if the total number of nodes is stable. In order to make the analysis more clearer, the system (1) is simplified as follows (2) with the meaning x(t) = X(t) N(∞) , where x(t) = (s(t), i(t), r(t), ls(t), li(t), lr(t)) T and X(t) = (S(t), I(t), R(t), LS(t), LI(t), LR(t)) T . Thus, x(t) represents the proportion of nodes in the system (1).

Local Stability and Analysis of Hopf Bifurcation
The disease-free solution e 0 and endemic solution e + of the system (2) can be obtained as follows: where ε = β 1 b+β 2 . The basic reproductive number R0 is obtained through the next generation matrix method [34]. Set: and Thus, It can be found that if R0 > 1, the model (2) has the unique endemic solution e + . By linearization technique, the characteristic matrix of system (2) about e + can be obtained as follows. where The characteristic equation of the model (2) at e + is shown as: The detailed representations of X n (X = A, B, C, D; n = 1, 2, 3 . . .) are given in Appendix A. Firstly, we study the local stability of system (2) at e + when τ = 0. If τ = 0, Equation (8) can be rebuilt as: where The following formulas can be obtained: If hypothesis (H1) E 5 , E 4 , E 3 , E 2 , E 1 , E 0 > 0, 1 , 2 , 3 , 4 , 5 , 6 > 0 and τ = 0, in Equation (10) holds, and all the characteristic solutions of Equation (9) are negative. Thus, e + is locally asymptotically stable under the theory of the Hurwitz criterion.
Then, we begin to consider the change of the characteristic solutions of Equation (9) when the charging time delay τ is introduced.
If τ > 0, multiplying e λτ on both sides of Equation (8), we can obtain the following equation: Mathematics 2021, 9, 2007 6 of 19 Referring to the processing method of the paper [35], let λ = iω (ω > 0) be the root of Equation (11). After separating the imaginary and real parts, we obtain where Squaring and adding the two Equations in Equation (12), we obtain Because of sin τω = ± √ 1 − cos 2 τω, the following two cases are discussed: (13) could be represented as follows: Equation (14) can be calculated as follows: where Let cos τω = y, and the following equation can be obtained: Set where The following solutions can be calculated: where Then, the expression of cos τω can be obtained as f 1 (ω) = cos τω; combined with Equation (13), we can obtain f 2 (ω) = sin τω. A function concerning ω can be obtained by Thus, (H2) Equation (20) has finite positive roots ω 1i (i = 1, 2, . . . , n). The corresponding time delay can be obtained as follows: (13) could be represented as follows: Similar to Case 1, a function concerning ω can be obtained by Thus, (H3) Equation (23) has finite positive roots ω 2i (i = 1, 2, . . . , n). The corresponding time delay can be obtained as follows: Let Then, if τ = τ 0 , Equation (13) has a pair of purely imaginary roots ±iω 0 . Next, the conditions for bifurcation are obtained in the following analysis. From Equation (11), we obtain dλ dτ where following Theorem 1 can be obtained by the Hopf bifurcation theorem [36].

Properties of the Hopf bifurcation
The properties of the Hopf bifurcation are explored in this section by applying the normal form theory and the center manifold theorem [36]. Let τ = τ 0 + µ, µ R. The transformations are given as u τ . Thus, the model (2) can be written as: . where According to the application of Riesz representation theorem, there exists a 6 × 6 matrix function η(θ, µ) : [−1, 0] → R 6×6 such as: The following equation is chosen: with δ is the Dirac delta function.
Following the computation introduced in [36], we obtain with Thus, the following values can be obtained: Theorem 2 is given as: , the bifurcating periodic solutions are stable (unstable). If T 2 > 0 (T 2 < 0), the bifurcating periodic solutions increase (decrease).

Parameter Dependence of R 0
In this section, prevention methods are analyzed in system (2). The basic reproduction number R 0 is generally used to describe the infectivity of virus. Therefore, it is also applicable to the spread of malware in WRSNs. On the assumption that the system is stable, if R 0 ≤ 1, malware will be cleared eventually. If R 0 > 1, malware will always exist.
The size of R 0 is affected by different parameters. What is worth concerning about is the influence of the parameter (α 2 , α 3 , α 4 or α 5 ) on the malware spread under different diffusion rates α 1 . Here, we first assume that only one parameter (α 2 , α 3 , α 4 or α 5 ) and α 1 are changed while the other parameters are unchanged. We observe the spread of malware in the system.
(1) The parameters are as follows in Figure 1a: 1], and α 2 ∈ [0, 1]. Figure 1a shows that with the increase of the diffusion rate α 1 , the malware is more easier prevailing. However, the self-healing rate of the infected nodes α 2 has a inhibitory effect on the spread of malware. Figure 1a also shows that the charging behavior (β 2 = 0.2) makes the spread of malware easier. (2) The parameters are as follows in Figure 1b: 1], and α 3 ∈ [0, 1]. Figure 1b shows that with the increase of the recovery rate of the infected nodes α 3 , the spread of malware can be effectively suppressed, which will provide us with the reference value in the control of malware. Similarly, Figure  1b also shows that without the charging behavior (β 2 = 0), the control of malware will become easier.  1], and α 5 ∈ [0, 1]. It is shown that the increase of the immune prevention rate α 4 can suppress the spread of malware but the increase of the immune failure rate α 5 helps the malware to spread. (4) Figure 1a-d together reflect that the charging behavior (β 2 = 0.2) will encourage the spread of malware, which will provide us with the data reference for the control of the malware spread in the SIRS-L model.
Next, the influence on the malware spread under different low-energy-node conversion rate β 1 and charging success rate β 2 are also discussed.
(5) The parameter settings are as follows in Figure 1e: 1], and β 2 ∈ [0, 1]. Figure 1e reflects the influence of the low-energy node conversion rate β 1 and the charging success rate β 2 on malware propagation. It is shown that if the low-energy node conversion rate β 1 is less than 0.2, the increase of the charging success rate β 2 can easily lead to the prevalence of malware. On the other side, if the low-energy node conversion rate β 1 is larger than 0.4, we can appropriately reduce the charging success rate β 2 to suppress the spread of malware. 0.7, ∈ [0, 1], and ∈ [0, 1]. Figure 1b shows that with the increase of the recovery rate of the infected nodes , the spread of malware can be effectively suppressed, which will provide us with the reference value in the control of malware. Similarly, Figure 1b also shows that without the charging behavior ( = 0), the control of malware will become easier. . It is shown that the increase of the immune prevention rate can suppress the spread of malware but the increase of the immune failure rate helps the malware to spread. (4) Figure 1a-d together reflect that the charging behavior ( = 0.2) will encourage the spread of malware, which will provide us with the data reference for the control of the malware spread in the SIRS-L model.
Next, the influence on the malware spread under different low-energy-node conversion rate and charging success rate are also discussed.  Figure 1e reflects the influence of the low-energy node conversion rate and the charging success rate on malware propagation. It is shown that if the low-energy node conversion rate is less than 0.2, the increase of the charging success rate can easily lead to the prevalence of malware. On the other side, if the low-energy node conversion rate is larger than 0.4, we can appropriately reduce the charging success rate to suppress the spread of malware.

Analysis and Display of Equilibrium Solutions
In this section, the distribution of nodes' states under different parameters of R 0 and τ are discussed to verify the stability of the system and to discuss the bifurcation.
(1) The parameter settings are as follows in Figure 2 (2), it will gradually disappear if R 0 < 1. It also can be shown that if β 2 = 0, the total proportion of the low-energy nodes (ls, li, lr) is larger than that of the general nodes (s, i, r). (2) The parameter settings are as follows in Figure 3 It shows that if the malware appears in the model (2), it will gradually disappear if R 0 < 1. It also can be shown that if β 2 = 0.2, the total proportion of the low-energy nodes (ls, li, lr) is smaller than the proportion of the low-energy nodes in Figure 2 if the charging operation is not performed β 2 = 0. (4) The parameter settings are as follows in Figure 5: It shows that the malware will prevail if R 0 > 1. Similarly, it also can be shown that if β 2 = 0.2, the total proportion of the low-energy nodes (ls, li, lr) is smaller than the proportion of the low-energy nodes in Figure 2 if the charging operation is not performed β 2 = 0. The total proportion of the low-energy nodes is similar to the corresponding proportion in Figure 4, which represents the proportion of low-energy nodes is only related to β 1 and β 2 .
Mathematics 2021, 9, 2007 22 of 30 In this section, the distribution of nodes' states under different parameters of and are discussed to verify the stability of the system and to discuss the bifurcation.  (2), it will gradually disappear if < 1. It also can be shown that if = 0, the total proportion of the low-energy nodes (ls, li, lr) is larger than that of the general nodes (s, i, r).     can be shown that if = 0.2, the total proportion of the low-energy nodes (ls, li, lr) is smaller than the proportion of the low-energy nodes in Figure 2 if the charging operation is not performed ( = 0).     can be shown that if = 0.2, the total proportion of the low-energy nodes (ls, li, lr) is smaller than the proportion of the low-energy nodes in Figure 2 if the charging operation is not performed ( = 0).   The above results (1)-(4) can be summarized as follows: the system is locally stable on the premise that the parameters meet the hypothesis (H1). In addition, if R 0 ≥ 1, the malicious software will persist in the system, while if R 0 < 1, the malicious software will be finally cleared. Next, we verify the bifurcation of the system. . It can be calculated that τ 0 = 11.5370 and ω 0 = 0.2254. In Figure 6, it can be seen that the model (2) is asymptotically stable if τ = 10.25 < τ 0 = 11.5370. What's more, in Figure 7, the model (2) undergoes a Hopf bifurcation if τ = 11.74 > τ 0 = 11.5370. According to Equation (41) and Theorem 2, the following parameters can be obtained: λ (τ 0 ) = 0.0292 + 0.0328i, C 1 (0) = −1.4520 + 0.7622i, µ 2 = 49.7260 > 0, β * 2 = −2.9040 < 0, and T 2 = −0.9203 < 0 . The result can be concluded that the Hopf bifurcation is supercritical, the bifurcating periodic solutions are stable and the period of the periodic solutions decreases.
It shows that the malware will prevail if > 1. Similarly, it also can be shown that if = 0.2, the total proportion of the low-energy nodes (ls, li, lr) is smaller than the proportion of the low-energy nodes in Figure 2 if the charging operation is not performed ( = 0). The total proportion of the low-energy nodes is similar to the corresponding proportion in Figure 4, which represents the proportion of low-energy nodes is only related to and . The above results (1)-(4) can be summarized as follows: the system is locally stable on the premise that the parameters meet the hypothesis (H1). In addition, if ≥ 1, the malicious software will persist in the system, while if < 1, the malicious software will be finally cleared. Next, we verify the bifurcation of the system.

Conclusions
This article puts forward a novel SISR-L model that considers the low-energy status nodes. Combined with the reality, the charging delay is introduced in the model to study the bifurcation. We have proved that the system is locally stable under certain conditions. In addition, the influence of different parameters on the spread of malware is displayed in the simulation. It is worth noting that in the SISR-L system, the charging behavior will make it more difficult to control the spread of malware. On the analysis of bifurcation, we reveal that the SISR-L model is locally asymptotically stable if the charging delay is less than , and the model undergoes a Hopf bifurcation at the endemic solution if the charging delay is larger than . Finally, the properties of Hopf bifurcation were explored through applying the normal form method and the center mani-

Conclusions
This article puts forward a novel SISR-L model that considers the low-energy status nodes. Combined with the reality, the charging delay τ is introduced in the model to study the bifurcation. We have proved that the system is locally stable under certain conditions. In addition, the influence of different parameters on the spread of malware is displayed in the simulation. It is worth noting that in the SISR-L system, the charging behavior will make it more difficult to control the spread of malware. On the analysis of bifurcation, we reveal that the SISR-L model is locally asymptotically stable if the charging delay is less than τ 0 , and the model undergoes a Hopf bifurcation at the endemic solution e + if the charging delay is larger than τ 0 . Finally, the properties of Hopf bifurcation were explored through applying the normal form method and the center manifold theorem. All the analyses can provide theoretical reference for the control of malware propagation in the SISR-L model.