Analysis of Time-Delay Epidemic Model in Rechargeable Wireless Sensor Networks

With the development of wireless rechargeable sensor networks (WRSNs), many scholars began to attach attention to network security under the spread of viruses. This paper mainly studies a novel low-energy-status-based model SISL (Susceptible, Infected, Susceptible, Low-Energy). The conversion process from low-energy nodes to susceptible nodes is called charging. It is noted that the time delay of the charging process in WRSNs should be considered. However, the charging process and its time delay have not been investigated in traditional epidemic models in WRSNs. Thus, the model SISL is proposed. The basic reproduction number, the disease-free equilibrium point, and the endemic equilibrium point are discussed here. Meanwhile, local stability and global stability of the disease-free equilibrium point and the endemic equilibrium point are analyzed. The addition of the time-delay term needs to be analyzed to determine whether it affects the stability. The intervention treatment strategy under the optimal control is obtained through the establishment of the Hamiltonian function and the application of the Pontryagin principle. Finally, the theoretical results are verified by simulations.


Introduction
With the development of science and technology, the research of wireless sensor networks (WSNs) has attracted scholars' attention in recent years. By organizing and combining sensor nodes freely through wireless communication technology, WSNs are established [1,2]. A WSN is a distributed sensor network, and its end is a sensor that can sense and check the outside environment. Sensors in WSNs communicate wirelessly, so the network setting is flexible, and the device position can be changed at any time. WSNs are widely applied in various fields such as military facilities, home automation, and the interconnection of transportation systems [3]. Wireless rechargeable sensor network (WRSN) technology came into being under the wireless power transfer (WPT) technology. Thus, the research of wireless sensors has been pushed to a new level.
However, WRSNs are usually vulnerable to malware attacks [4][5][6], such as denial of service, malicious code distribution, and information exfiltration [7]. WRSNs especially suffer from the Denial of Charge (DOC) [8]. It is of great significance to study the network security of WRSNs.
The security of sensor networks has become a hot research field [9,10]. The research methods here are as follows: 1.
Establish a mathematical epidemics model.

2.
Analyze the equilibrium point in the system, and study its stability. 3.
Achieve the optimal control to make the control of malware propagation more effective. According to previous studies, the classical mathematical models established by predecessors mainly include the following: SIS (Susceptible, Infected, Susceptible), SIR (Susceptible, Infected, Removed), and SIRS (Susceptible, Infected, Removed, Susceptible). Consideration of immune validity period and the analysis of the optimal control.
Zhu, L. et al. [20] SEIR (Susceptible-Exposed-Infected-Recovered) Time delay of the incubation 1 Time-delay analysis of infection process based on SEIR system.
Al-Darabsah, I. et al. [21] SEIRS-V (Susceptible-Exposed-Infected-Recovered-Susceptible and Vaccinated) Time delay of the exposure 1 Time-delay analysis considering specific effective contact rate. However, we know that the time delay of the charging process in WRSNs has not been investigated previously. This paper mainly provides a kind of influence on time delay based on the SIS model combined with low-energy nodes (L) [23][24][25]. Thus, a novel epidemic model (1) of virus spreading in WRSNs is developed.
The virus-propagation rules and the optimal control strategy of the model (1) is proposed as follows:

1.
A novel low-energy-status-based model is introduced to describe the propagation process of malicious software (virus) in WRSNs.

2.
The basic reproduction number is defined by the disease-free equilibrium solution and the epidemic equilibrium solution. The Routh criterion is applied to prove the local stability, and the Lyapunov universal function is constructed to prove the global attraction.

3.
Based on Pontryagin's minimum principle, the optimal control variables satisfying the minimization of the objective function are obtained. 4.
The stability problem under time delay is specially revealed.
The rest of the paper is as follows: the introduction and analysis of the model are presented in Section 2; the dynamic stability analysis is given in Section 3; the optimal control theory is presented in Section 4; the analysis of the influence of time delay is given in Section 5;. the simulation results are given in Section 6; and the conclusions are summarized in Section 7.

Modeling Analysis
As mentioned above, a classical model called SIS can be applied in wireless sensor networks. However, the SIS model can be modified by increasing the low-energy status L(t) in WRSNs [11]. The modified model here is shown in Figure 1. Different from the model in [11], in order to control malware more efficiently, charging low-energy-infected nodes is neglected. The process of the nodes L(t) converting into S(t) is called charging. The feature of the charging process is its latency. First, without considering the time delay, the system (1) is discussed here (with the initial values (0) = Φ S , . Parameters are shown in Table 2. However, we know that the time delay of the charging process in WRSNs has not been investigated previously. This paper mainly provides a kind of influence on time delay based on the SIS model combined with low-energy nodes (L) [23][24][25]. Thus, a novel epidemic model (1) of virus spreading in WRSNs is developed.
The virus-propagation rules and the optimal control strategy of the model (1) is proposed as follows: 1. A novel low-energy-status-based model is introduced to describe the propagation process of malicious software (virus) in WRSNs. 2. The basic reproduction number is defined by the disease-free equilibrium solution and the epidemic equilibrium solution. The Routh criterion is applied to prove the local stability, and the Lyapunov universal function is constructed to prove the global attraction. 3. Based on Pontryagin's minimum principle, the optimal control variables satisfying the minimization of the objective function are obtained. 4. The stability problem under time delay is specially revealed.
The rest of the paper is as follows: the introduction and analysis of the model are presented in Section 2; the dynamic stability analysis is given in Section 3; the optimal control theory is presented in Section 4; the analysis of the influence of time delay is given in Section 5;. the simulation results are given in Section 6; and the conclusions are summarized in Section 7.

Modeling Analysis
As mentioned above, a classical model called SIS can be applied in wireless sensor networks. However, the SIS model can be modified by increasing the low-energy status L(t) in WRSNs [11]. The modified model here is shown in Figure 1. Different from the model in [11], in order to control malware more efficiently, charging low-energy-infected nodes is neglected. The process of the nodes L(t) converting into S(t) is called charging. The feature of the charging process is its latency. First, without considering the time delay, the system (1) is discussed here (with the initial values (0) = ， (0) = ， (0) = 0， (0) = 0, 0 ≤ , ≤ ∧ ). Parameters are shown in Table 2.

S(t)
The quantity of the susceptible node I(t) The quantity of the infected node L(t) The quantity of the low-energy node ∧ Injection rate of system α 1 Transmission rate of malware α 2 Self-disinfection rate of the infected nodes c Charging success rate µ Low-energy-node conversion rate d Node deactivation rate In order to further study the dynamic balance of the system, and the propagation law of malicious programs, the equilibrium solution of the system needs to be found first. By studying the stability and transformation behavior of each equilibrium solution, the general law of malware control could be drawn, and intervention for controlling the spreading of malicious programs could be easily put forward.
From the numerical relationship, the system (1) has two equilibrium solutions (the disease-free equilibrium solution E 0 and the endemic equilibrium solution E + ). The solutions E 0 and E + are depicted as follows: The basic reproduction number R 0 is given as If (d + µ)(d + c) − cµ > 0, we can put forth the following theorem: Theorem 1. If R 0 < 1, the system only has a unique disease-free equilibrium solution E 0 . While R 0 ≥ 1, the system (1) has an endemic equilibrium solution E + and a disease-free equilibrium solution E 0 .

Dynamic-Stability Analysis
In this section, we focus on analyzing the local stability of the system (1) by using the linearization technique and the Routh-Herwitz criterion. However, the global stability is analyzed by constructing a Lyapunov universal function.

Local Stability Analysis
Theorem 2. If R 0 < 1, the disease-free equilibrium point E 0 is locally stable.
Proof. First of all, the system is linearized as follows: J(E * ) is the corresponding Jacobian matrix, in which "*" represents "0" or "+". The characteristic equation is constructed using Formula (6): which gives the following equation: Theorem 2 is valid under the theory of the Hurwitz criterion.
Theorem 3. If R 0 ≥ 1, the endemic equilibrium point E + is locally stable.
Proof. The analysis method is the same as above. The characteristic equation of J(E + ) is found as follows: which donates the following equation: Set: Mathematics 2021, 9, 978 6 of 19 The following results could be obtained: Theorem 3 is proved under the theory of the Hurwitz criterion.

Global-Stability Analysis
In this section, we discuss the global stability by constructing a suitable Lyapunov universal function, which is positive definite except for the stable point.
Proof. The Lyapunov function is defined as follows: If "t" tends to infinity, the disease-free equilibrium solution satisfies the following condition: S ≤ S 0 . Thus: Equation (15) shows Theorem 5. If R 0 ≥ 1, the endemic equilibrium point E + is globally asymptotically stable.
Proof. Considering the conditions: the Lyapunov function is defined as follows: Thus: Equation (18) shows It is noted that V E + (t) = 0 and V E + (t) = 0 if and only if S = S + , I = I + and L = L + . Thus, Theorem 5 is confirmed. This means that the spread of malware will continue to exist in the system (1).

Optimal-Control Analysis
In the practical application of malware control, intervention treatment is usually added to ensure that the spread of malware will not be epidemic. However, there is a certain cost due to intervention treatment. The Pontryagin minimum principle is adopted to solve the cost-minimization problem. Due to the characteristics of rechargeable systems, our optimization objectives are as follows: First, we must effectively control the spread of malware. The second objective is to ensure that the number of susceptible nodes in the system is not too small, and the low-energy nodes should be controlled within a certain number to ensure the normal operation of the system. Lastly, the cost of control needs to be feasible economically. Therefore, the system (1) is modified as follows (19): We establish an objective function: where c 1 is the cost coefficient of implementing intervention control, which mainly comes from the identification and removal of malicious software.; and u(t) is the intervention ratio, which satisfies 0 ≤ u(t) ≤ 1; tf is the end time of intervention. It is obvious that the number of low-energy nodes and infected nodes at tf should be as few as possible, and the economic cost should be as low as possible. The Pontryagin minimum principle is applied to find the optimal proportion of intervention. The corresponding Hamiltonian function is as follows: in which β i (t) (i = 1, 2, 3) are covariates. The differential equation of covariant is as follows: The terminal constraints of the co-state variables are as follows: The optimal strategies are established by: Then, the optimization variable can be obtained:

Time-Delay System Analysis
If the number of system nodes is relatively large, the phenomenon of charging delay cannot be ignored. The time-delay item (t − τ) is considered in this section. In order to better study the influence of time delay, we need to explore whether the delay term will affect the stability of the system. The system is established as follows:

Local-Stability Analysis
In this section, the characteristic solutions of the transcendental equations are explored by the theories in [26]. Theorem 6. If R 0 < 1, the local stability at E 0 of the system (26) is consistent with the stability of the system (1).
Suppose that for the characteristic equation of the system (26) at E 0 , it has a pure imaginary characteristic solution iω 0 .
The characteristic equation at E 0 is: The following formula can be obtained: The Square and add the two Equations (30) and (31): Set: Then, We can conclude that: It can be obtained that Hypothesis 1 is not valid. Thus, Theorem 6 is proved.
Theorem 7. If R 0 ≥ 1, the local stability at E + of the system (26) is consistent with the stability of the system (1).
Suppose that for the characteristic equation of the system (26) at E + , it has a pure imaginary characteristic solution iω 0 .
For E + , the solution is the same. Its characteristic equation at E + is: We get: Similarly, the pure imaginary is substituted into the equation: Separate the real and imaginary parts of Equation (38): Let ω 2 0 = v 0 > 0. Square and add the two Equations (39) and (40): These terms satisfy the following equations: and Thus, the solutions of v 0 are negative. This contradicts the condition v 0 > 0. Therefore, Hypothesis 2 does not hold, and Theorem 7 can be proved.

Global-Stability Analysis
In this section, Lyapunov functions are constructed to verify the stability of the equilibriums of the system (26). Theorem 8. If R 0 < 1, the disease-free equilibrium point E 0 is globally asymptotically stable.
Proof. The Lyapunov function is defined as follows: What is more, V E 0 (t) = 0 and V E 0 (t) = 0 if and only if S = S 0 and I = 0. Thus, Theorem 8 is valid. The spread of malware will not prevail.
However, for the epidemic-spreading equilibrium point, the Lyapunov function is constructed as follows: Equation (47) shows V E + (t) ≤ 0 if |µ − c| ≤ 2d. Therefore, the following Theorem 9 is proved. Theorem 9. If R 0 ≥ 1, the endemic equilibrium point E + is globally asymptotically stable.

Parameter Dependence of R 0 and I(∞)
In this section, some possible prevention schemes are provided to analyze the influence of different parameters in the system (1) on R 0 and I(∞). R 0 determines whether the spread of malware is prevalent or not. Meanwhile, I(∞) reflects the seriousness of malware spread.
The parameter settings are given as follows in Figure 2a: ∧ = 10, d = 0.05, µ = c = 0.02 and α 1 ∈ [0, 0.1], α 2 ∈ [0, 1]. The phenomenon could be found that if the infection rate α 1 is at a very low level, the spread of malware will not prevail. If α 2 increases, the virus will break out. The self-recovery rate can hardly inhibit the spread of the virus if α 1 is larger than 0.006.
The parameter settings are given as follows in Figure Figure 2c shows that if d is small enough, ∧ is more likely to cause a virus outbreak. As shown in Figure 2d, to keep the virus from spreading, the larger the node mortality rate d is, the wider the allowable range of injection rate ∧ is.
The parameter settings are given as follows in Figure  The parameter settings are given as follows in Figure 2e: ∧= 2, = 0.05, = 0.002, = 0.03, ∈ 0, 1 , ∈ 0, 1 . The charging rate and the conversion rate of low-energy nodes have no decisive effect on whether the virus breaks out or not.

Stability of Equilibrium Point
The parameter settings are given as follows in Figure

Stability of Equilibrium Point
The parameter settings are given as follows in Figure 3: ∧ = 10, d = 0.12, α 1 = 0.002, α 2 = 0.03, µ = c = 0.02 and S(0) ∈ [30, 84], I(0) ∈ [30, 39], L(0) = 0. It can be found that the initial values of states will not change the convergence of the system. It is noted that S(∞), I(∞), L(∞) = 72.9167, 0, 10.4167. Even if the number of initial infected nodes is relatively large, the system virus will eventually disappear as long as R 0 ≤ 1. The results shown in Figure 3 confirm Theorem 4. found that the initial values of states will not change the convergence of the system. It is noted that S(∞), I(∞), L(∞)= 72.9167, 0, 10.4167. Even if the number of initial infected nodes is relatively large, the system virus will eventually disappear as long as ≤ 1.
The results shown in Figure 3 confirm Theorem 4.

Influence of Time Delay on Malware Propagation
In this section, the influence of time delay on the convergence is discussed. found that the initial values of states will not change the convergence of the system. It is noted that S(∞), I(∞), L(∞)= 72.9167, 0, 10.4167. Even if the number of initial infected nodes is relatively large, the system virus will eventually disappear as long as ≤ 1.
The results shown in Figure 3 confirm Theorem 4.

Influence of Time Delay on Malware Propagation
In this section, the influence of time delay on the convergence is discussed.

Influence of Time Delay on Malware Propagation
In this section, the influence of time delay on the convergence is discussed.
The parameter settings are given as follows in Figure 5: ∧ = 10, d = 0.05, α 1 = 0.002, α 2 = 0.03, µ = c = 0.02 and S(0) ∈ [30, 84], I(0) ∈ [1, 10], L(0) = 0, τ = 10. The delay term will not affect the convergence of the system in Figure 5b. Because of the introduction of τ, the convergence values will have a very small wave, as shown in Figures 4a and 5a. Two groups of convergence values from the system (1) and the system (26) will get closer and closer as time passes. As shown in Figure 5b, different delay values only affect the convergence speed.
The parameter settings are given as follows in Figure 5: ∧ = 10, = 0.05, = 0.002, = 0.03, = = 0.02 and (0) ∈ 30, 84 , (0) ∈ 1, 10 , (0) = 0, = 10. The delay term will not affect the convergence of the system in Figure 5b. Because of the introduction of τ, the convergence values will have a very small wave, as shown in Figures  4a and 5a. Two groups of convergence values from the system (1) and the system (26) will get closer and closer as time passes. As shown in Figure 5b, different delay values only affect the convergence speed. As shown in Figure 6, with the small delay item (0. 1-0.8), the delay of the S nodes' convergence is not obvious. However, it is concerning that if < 1 in Figure 6a, the S nodes will be in the zero-quantity state for a period under the effect of large time delay item (100-800), which should be avoided by increasing charging power to reduce the time delay. If > 1 in Figure 6b, the delay term will only introduce the speed influence of convergence. As shown in Figure 7a, the different time-delay terms have little effect on the convergence of the I node if < 1. As shown in Figure 7b, I nodes will go through a sud- As shown in Figure 6, with the small delay item τ(0. 1-0.8), the delay of the S nodes' convergence is not obvious. However, it is concerning that if R 0 < 1 in Figure 6a, the S nodes will be in the zero-quantity state for a period under the effect of large time delay item τ(100-800), which should be avoided by increasing charging power to reduce the time delay. If R 0 > 1 in Figure 6b, the delay term will only introduce the speed influence of convergence.
The parameter settings are given as follows in Figure 5: ∧ = 10, = 0.05, = 0.002, = 0.03, = = 0.02 and (0) ∈ 30, 84 , (0) ∈ 1, 10 , (0) = 0, = 10. The delay term will not affect the convergence of the system in Figure 5b. Because of the introduction of τ, the convergence values will have a very small wave, as shown in Figures  4a and 5a. Two groups of convergence values from the system (1) and the system (26) will get closer and closer as time passes. As shown in Figure 5b, different delay values only affect the convergence speed. As shown in Figure 6, with the small delay item (0. 1-0.8), the delay of the S nodes' convergence is not obvious. However, it is concerning that if < 1 in Figure 6a, the S nodes will be in the zero-quantity state for a period under the effect of large time delay item (100-800), which should be avoided by increasing charging power to reduce the time delay. If > 1 in Figure 6b, the delay term will only introduce the speed influence of convergence. As shown in Figure 8, similarly, the delay term does not affect the convergence value but delays the convergence process obviously. It is worth noticing that the time-delay term will bring a fluctuation of the convergence process in Figure 8b if > 1.

Realization of Optimal Control
The flowchart of simulation implementation is shown in Figure 9. The intervention treatment ratio (t) is obtained. As shown in Figure 8, similarly, the delay term does not affect the convergence value but delays the convergence process obviously. It is worth noticing that the time-delay term will bring a fluctuation of the convergence process in Figure 8b if R 0 > 1. As shown in Figure 8, similarly, the delay term does not affect the convergence value but delays the convergence process obviously. It is worth noticing that the time-delay term will bring a fluctuation of the convergence process in Figure 8b if > 1.

Realization of Optimal Control
The flowchart of simulation implementation is shown in Figure 9. The intervention treatment ratio (t) is obtained.

Realization of Optimal Control
The flowchart of simulation implementation is shown in Figure 9. It shows that the objective function converges through iterations: (G(u) = J u k , The intervention treatment ratio u k (t) is obtained. Mathematics 2021, 9, x FOR PEER REVIEW  Figure 10 shows the convergence of the system without intervention treatm can be seen that systemic viruses are always prevalent. Figure 11 presents the o intervention-treatment proportion. The optimal intervention-treatment curve is z the beginning, and gradually increases to a value of 1.   Figure 10 shows the convergence of the system without intervention treatment. It can be seen that systemic viruses are always prevalent. Figure 11 presents the optimal intervention-treatment proportion. The optimal intervention-treatment curve is zero at the beginning, and gradually increases to a value of 1.
To better compare the economic utility brought by the optimal control, state variable values at the terminal time (tf = 200) in Figure 12a are adjusted to be close to the state variable values in Figure 12b. The fixed treatment coefficient in Figure 12a is u(t) = 0.3. Comparing the size of J calculated by the objective function (20), the cost J in Figure 13b is 32.0104. However, the cost J in Figure 13a is 1901.5. Combined with the state variable values at the terminal time of Figure 12b, the strategy under optimal control not only stops the spread of the virus, but also greatly reduces the cost required by the security policy.  Figure 10 shows the convergence of the system without intervention treatment. It can be seen that systemic viruses are always prevalent. Figure 11 presents the optimal intervention-treatment proportion. The optimal intervention-treatment curve is zero at the beginning, and gradually increases to a value of 1.  To better compare the economic utility brought by the optimal control, state variable values at the terminal time (tf = 200) in Figure 12a are adjusted to be close to the state variable values in Figure 12b. The fixed treatment coefficient in Figure 12a is u(t) = 0.3. Comparing the size of J calculated by the objective function (20), the cost J in Figure 13b is 32.0104. However, the cost J in Figure 13a is 1901.5. Combined with the state variable values at the terminal time of Figure 12b, the strategy under optimal control not only stops the spread of the virus, but also greatly reduces the cost required by the security policy.  To better compare the economic utility brought by the optimal control, state variable values at the terminal time (tf = 200) in Figure 12a are adjusted to be close to the state variable values in Figure 12b. The fixed treatment coefficient in Figure 12a is u(t) = 0.3. Comparing the size of J calculated by the objective function (20), the cost J in Figure 13b is 32.0104. However, the cost J in Figure 13a is 1901.5. Combined with the state variable values at the terminal time of Figure 12b, the strategy under optimal control not only stops the spread of the virus, but also greatly reduces the cost required by the security policy.

Conclusions
To better inhibit the spread of viruses, a novel epidemic model of virus spreading in WRSNs was developed after taking the time delay of the charging process into consideration. The local stability and global stability of the disease-free equilibrium point and the epidemic equilibrium point were proved by defining the basic regeneration number and the application of the Routh criterion. In addition, we concluded that the time-delay term would not change the convergence of the system, but would only affect the convergence speed. It is worth noting that Figure 6a shows that if R < 1, the time delays

Conclusions
To better inhibit the spread of viruses, a novel epidemic model of virus spreading in WRSNs was developed after taking the time delay of the charging process into considera-tion. The local stability and global stability of the disease-free equilibrium point and the epidemic equilibrium point were proved by defining the basic regeneration number and the application of the Routh criterion. In addition, we concluded that the time-delay term would not change the convergence of the system, but would only affect the convergence speed. It is worth noting that Figure 6a shows that if R 0 < 1, the time delays will make the S nodes zero for a while, which will greatly affect the normal operation of the network. In terms of parameter relationships, they can reflect the benefits of the modified SILS model in Figure 2e. It shows that the increase of charging power will not make the virus spread out of control, which is different from the SILS model. The asymptotic stability of the SISL model was strictly verified in simulation. At the same time, we verified the effectiveness of the proposed optimal control strategy, which not only effectively controlled the spread of the virus, but also greatly reduced the maintenance cost. The SISL model considers the influence of time delay and provides a new research perspective. We hope the results of this paper can provide some reference for related researchers.