Dynamics Analysis of a Wireless Rechargeable Sensor Network for Virus Mutation Spreading

Virus spreading problems in wireless rechargeable sensor networks (WSNs) are becoming a hot topic, and the problem has been studied and discussed in recent years. Many epidemic spreading models have been introduced for revealing how a virus spreads and how a virus is suppressed. However, most of them assumed the sensors are not rechargeable sensors. In addition, most of existing works do not consider virus mutation problems. This paper proposes a novel epidemic model, including susceptible, infected, variant, low-energy and dead states, which considers the rechargeable sensors and the virus mutation factor. The stability of the proposed model is first analyzed by adopting the characteristic equation and constructing Lyapunov functions methods. Then, an optimal control problem is formulated to control the virus spread and decrease the cost of the networks by applying Pontryagin’s maximum principle. Finally, all of the theoretical results are confirmed by numerical simulation.


Research Background
In recent years, wireless sensor networks (WSNs) have played a crucial role in the internet of things (IoT), and WSNs has been widely developed and applied in many fields, for example: industrial, military and healthcare applications. Benefiting from the recent breakthrough in WSNs, they have attracted increasingly attention [1,2]. WSNs encompass numerous sensor nodes and nodes connect with each other with the radio signal, but it is difficult to create a complex security protective structure. Due to the vulnerability, the network is always destroyed by malware and this leads to information leakage and even network paralysis. Thus, security is an essential problem in WSNs to ensure accuracy and efficiency. Scholars have done a lot of research on the security of WSNs [3][4][5][6][7].
WSNs suffer from the issue of the vulnerability of the network and energy limitation. Wireless rechargeable sensor networks (WRSNs) [8] are considered wireless power transfer (WPT) technology which greatly improve WSNs. The security of WSNs has been pushed to a new level by optimizing the charging scheduling and analyzing the results of denial of charge attacks [9,10]. It is popular to apply epidemic dynamics analysis in WRSNs attacked by malwares. The classical mathematical models have been researched by scholars, including susceptible, infected (SI), susceptible, infected, susceptible (SIS), susceptible, infected, removed (SIR), susceptible, exposed, infected, removed (SEIR), etc. However, the malware can be changed when the malware attacks the network [11], and the same is true for a virus mutation. Once the virus mutation happens in WRSNs, it is devastating, and the system of the security strategy is destroyed. However, it has seldom been researched by scholars before.
Thus, this paper includes virus mutation to establish a new model, which is very effective to reduce the spreading of the virus, reduce the harm of the virus mutation to WRSNs and greatly improve the security of WRSNs.

Related Work
In order to improve the security in WRSNs, the efficient scheme and energy efficient secured ring routing (E2SR2) protocol was first proposed by Shafie et al. and Bhushan et al. [12,13]. For the problem of physical node capture attacks, a response strategy was proposed by Bonaci et al. [6] to ensure security and stability of the network connectivity when the network is attacked. Considering detection and correction, Sing et al. [3] also proposed a selective forwarding attacks technology which increases the QOS and provides better data transmission.
Applying epidemic dynamics analysis to the study of WSNs, Kephart et al. [14,15] firstly proposed a model to study and predict the spread of the virus. Since then, many new models have been proposed to study the problem. Cao et al. [16] considered the recovered factor to construct the model of WSNs, and the authors patched the dissemination of security and immunized or healed the nodes to use the security strategy; Zheng et al. [17] considered the vaccination strategies with temporary immunity and a quarantined strategy to construct an SEIQR model; Han et al. [18] considered the system of the nonlinear stochastic system to construct an SEIR model. Liu et al. [19] considered the low-energy factor to construct a model of WRSNs, and the stability of the model was proved. According to Liu et al. [19][20][21][22], the system is a good combination of an epidemic dynamic model and wireless sensor network. Thus, the system skillfully blends the information of the WRSNs with biological characteristics. However, the all of above models do not consider the problem of virus mutation. Above, Table 1 lists the recent relevant studies.

Research Field Model Content
Yang et al. [23] Dynamic analysis of the virus mutation model SIS Proof of the local and global stability of the system Dong-Mei et al. [24] Model analysis of disease viruses mutated in the process of transmission SEIR Proof of the local and global stability of the system that considers the exposed Gao et al. [25] An SEIR epidemic model analysis with logistic death rate of virus mutation SEIR Proof of the local and global stability of the system that considers the logistic death rate of virus mutation Tong et al. [26] Dynamic model analysis with delay of the virus mutation SIS Proof of the local and global stability of the system that considers the time delay Dong-Mei et al. [27] SIR model analysis with delay of the virus mutation SIR Proof of the local and global stability of the system that considers recovered factor and the time delay De-gang et al. [28] A variation epidemic model's propagation and analysis in complex networks SIVR Proof of the local and global stability of the system Cai et al. [29] Model analysis of spread of the pathogen with mutant strain and vaccination SIVR Proof of the local and global stability and analysis of the Hopf bifurcation of the system Xu et al. [30] Optimal control of the SIVRS epidemic spreading model with virus mutation in complex networks SIVRS Considers the optimal strategy and calculates the optimal control results of the system As shown in Table 1, based on the virus mutation, scholars have done some studies. However, as we all know, the theories of the epidemic dynamics with virus mutation applied in WRSNs are rarely studied. The infected nodes infect other nodes to cause mutations because of coding errors, decryption problems, etc. [11], hence, the idea of virus mutation of artificial life was introduced, and a new way needs to be proposed to analyze and solve the security of the mutation virus model. Thus, in this paper, a new scheme is put forward of a model combined with virus mutation and is analyzed to solve the WRSN spreading problem of a mutated virus.

Contributions
In the previous studies, the charging behavior and mutation virus activity of wireless sensor networks has rarely been investigated. In this paper, the main goal is to introduce a model that considers the low-energy states and virus mutation. Combined with a practical application, the system can be applied when the stability of the system is proven and the other conditions are not taken into account, and the cost can be controlled by an optimal strategy. It is divided into five variable states. Our contributions are as follows: 1.
An epidemic model suitable for WRSNs is established to describe the propagation process of malwares (the virus and mutated virus).

2.
To analyze and calculate the basic reproductive numbers R 1 and R 2 . Then, considering the existence of equilibrium of the system, the local and global asymptotic stability is proved by adopting the characteristic equation and the Lyapunov principle.
Numerical simulation is carried out to confirm the results.

3.
By constructing the objective function and applying Pontryagin's maximum principle, we can obtain the optimal control variable which satisfies the optimal control objective of the security problem.
The rest of the paper is as follows: In Section 2, the introduction and analysis of the model is presented; in Section 3, the local and global asymptotic stability is proved; in Section 4, the optimal strategy is presented; in Section 5, the numerical simulation verifies the proposed theoretical results; in Section 6, the relevant conclusion of the model is presented.

Model Analysis
In this section, the epidemic model of virus mutation in WRSNs is introduced. Sensors are divided into five sensor node states. The five sensor node states include: S, I 1 , I 2 , L and D. The S state is susceptible to virus attacks. I 1 is the infected node state, I 2 is the mutant virus node state, L is the low-energy state of all nodes, D represents death nodes (the sensor nodes cannot work). N stands for the sum of all nodes, and it is a certain number. Thus, N(t) = S(t) + I 1 (t) + I 2 (t) + L(t). The state transition diagram is shown in Figure 1.
Obviously, adding Equation (1), we obtain: As shown in Figure 1, A is the number of injected nodes; β 1 is the conversion rate of susceptible to infected; β 2 is the conversion rate of the sensor nodes becoming infected with the virus when the virus mutates; γ 1 is the rate of cleaning I 1 virus; γ 2 is the rate of cleaning I 2 virus; d is the conversion rate when sensors are unable to work; c is the rate of nodes going into a low-energy state; u is the rate of recovery from L energy, and during the low-energy state, I 1 and I 2 have removed the virus; ε is the rate of virus mutation. In this paper, when infected nodes and mutant virus nodes are converted to nodes of the low-energy state, a virus remover is created to eliminate the virus and mutant virus. According to the parameters and nodes of states, the dynamics equation of the model (1) is given as the following system (1): Obviously, adding Equation (1), we obtain: As Then, considering the limit of (3), we can obtain N(t) ≤ A d , so the feasible region for (1) is given by: Above, Ψ is a positively invariant.

Computing the Equilibrium Points and Basic Reproductive Number
In order to discuss the existence of equilibria, let the left-hand side of Equation (1) be zero. Obviously, there are three solutions that satisfy this situation. The disease-free equilibrium is (S 0 , 0, 0, L 0 ), the individual plant virus equilibrium is (S, 0, I 2 , L) and the endemic equilibrium is (S * , I * 1 , I * 2 , L * ). It is noted that the basic reproductive number determines the existence of different equilibria. In this paper, there are two basic reproductive numbers (R 1 and R 2 ), where R 1 = S 0 S * , R 2 = S 0 S . S 0 , S and S * are given as: Thus, the basics reproductive numbers are given as: According to the existence of the equilibrium point, if R 1 ≤ 1, R 2 ≤ 1, the system (1) has only a disease-free equilibrium point E 0 = (S 0 , 0, 0, L 0 ) given as: Entropy 2021, 23, 572

of 16
If R 2 > 1, the system has a disease-free point and an individual plant virus equilibrium point, and the individual plant virus equilibrium E = (S, 0, I 2 , L) is given as: If R 1 > 1, R 1 > R 2 , the system has a disease-free equilibrium point and an endemic equilibrium point. The endemic equilibrium E * = S * , I * 1 , I * 2 , L * is given as: Thus, we have the following Theorem 1: there are two equilibriums including E 0 and E * .

Dynamic Stability Analysis
In this section, the epidemic model of virus mutation in WRSNs is divided into local and global stability analysis.
System (1) is a nonlinear system, and we must turn it into a linear system so as to prove the local stability of this system. In order to reduce the calculated amount of the system, L in the first equation of the system (1) can be replaced with N(t) − (S(t) + I 1 (t) + I 2 (t)). Thus, in system (1), the last equation is independent of the other 3 equations, and the linear equations can be given as the system (15): Then, the Jacobian matrix of (15) is given by:

Local Stability
Theorem 2. If R 1 ≤ 1 and R 2 ≤ 1, the disease-free equilibrium is locally asymptotically stable.
Proof. According to the Lyapunov criterion [31], if all of the characteristic values are negative, the equilibrium is locally asymptotically stable. Thus, the Jacobian matrix is given by: Hence, we have the following characteristic Equation (18): All of the characteristic values are negative, so the disease-free equilibrium is locally asymptotically stable. Theorem 3. If R 2 > 1 and R 2 > R 1 , the individual plant virus equilibrium is locally asymptotically stable.
Proof. The Jacobian matrix is given as: Hence, the characteristic equation is given by: Then, the characteristic value is Obviously, if R 2 > R 1 , λ 2 is always a negative number. Hence, the individual plant virus equilibrium is locally asymptotically stable. Theorem 4. If R 1 > 1 and R 1 > R 2 , the endemic equilibrium is locally asymptotically stable.
Proof. The Jacobian matrix is given by: By the transformation of (21), we obtain We have the following characteristic equation: Obviously, λ 1 = −(d + c + u), and the other characteristic values can be calculated by the following Equation (24): According to Vieta's formulas [32], the zero solution of (24) is given by , thus, if R 1 > 1 and R 1 > R 2 , the value of λ 2 and λ 3 can be proved to be λ 2 + λ 3 < 0 and λ 2 * λ 3 > 0. Hence, the endemic equilibrium is locally asymptotically stable.

Global Stability
In this part, three Lyapunov functions are constructed to verify the global stability of the disease-free equilibrium E 0 , individual plant virus equilibrium E and endemic equilibrium E * . Theorem 5. If R 1 ≤ 1,R 2 ≤ 1 and R 1 + ε γ 1 +c+ε+d < 1, the disease-free equilibrium is globally stable.
Proof. According to the second Equation (1), the number of infected nodes is I 0 at t 0 (t > t 0 ). We can obtain.
When t → ∞ , I 1 (t) will be stable at zero. Thus, the limit Equation (27) is given as: As I 1 (t) converges to a limit value zero, the Lyapunov function can be given as: The derivative of V(S, I 2 ) of (28) is given by: Obviously, if R 2 > 1, R 1 < 1, the individual plant virus equilibrium is globally stable.
Proof. It is noted that the rate of infection is higher than the rate of virus mutation. If the rate of virus mutation is enough small and is equal to be zero, it is assumed that ε = 0. The Equation (1) is given as: The Lyapunov function V(S, I 1 , I 2 ) is given as follows: V(S, I 1 , Obviously, V(S, I 1 , I 2 ) ≥ 0. Then, the derivative of V(S, I 1 , I 2 ) of the solution of (31) is given by: If R 1 > R 2 , the endemic equilibrium is globally stable if ε = 0. Theorem 6 is proved.

Optimal Strategy
In this subsection, the cost of cleaning the virus and the cost of charging low-energy nodes are considered. For this purpose, the control U(·) is introduced as the objective function, γ 1 (0 < γ 1 < 1) is the rate of cleaning I 1 virus; γ 2 (0 < γ 2 < 1) is the rate of cleaning I 2 virus and u(0 < u < 1) is the rate of recovery energy. Q I 1 is the treatment cost coefficient of I 1 , Q I 2 is the treatment cost coefficient of I 2 and Q L is the charging cost coefficient of low-energy nodes. Thus, the optimal control problem is to minimize the objective function as follows: It is clear that the feasible region of the control variable set U of the system is [0, 1]. Hence, the optimization goal is to diminish I 1 and I 2 during the time interval [0, t f ]. According to the Pontryagin maximum principle, the Hamiltonian function is constructed as follows: X is the state variable, U is the control variable set, α i (i = 1, , 2, 3, 4) is the adjoint variables. The adjoint variables are determined by solving the following equations.
At the moment, the transversality condition of optimal control satisfies the following equation: According to differential equations of the covariant variable and transversal condition of optimal control, the optimization condition is given by: Finally, Pontryagin's maximum principle is applied to obtain the optimal control variable set. The result is given as:

Numerical Simulation
In this section, all of the theoretical analyses are proved and the numerical results of the system (1) are presented to support the analytic results. The result of the numerical simulation is given as follows.

Stability Simulation
The parameters are given so that A = 10, d = 0.1, c = 0.5, u = 0.2, ε = 0.2. It can be seen that rates of β 1 and γ 1 have significant impacts on the basic reproduction number R 1 , shown in Figure 2a. With the increasing of β 1 , R 1 increases rapidly, and with the increasing of γ 1 , R 1 increases slowly. In the same way, the rates of β 2 and γ 2 also have significant impacts on the basic reproduction number R 2 , shown in Figure 2b. With the increasing of β 2 , R 2 increases rapidly, and with the increasing of γ 2 , R 2 increases slowly. The relationship between 1 , 2 and equilibrium is shown in Figure 3. The parameters of the system (1) are set as = 10, = 0.5, = 0.2. It is noted that the individual plant virus equilibrium ̅ 2 is determined by and the basic reproduction number 2 . The rate of death has significant impacts on the ̅ 2 if is small enough, as shown in Figure 3a. It can be seen that the relationship between the rate of virus mutation and 2 * is a linear, as shown in Figure 3b. The relationships among 1 , 2 and 1 , 2 are shown in Figure 3c,d when the system (1) satisfies the conditions of the existence of endemic equilibrium ( 1 > 1, 1 > 2 ). The relationship between R 1 , R 2 and equilibrium is shown in Figure 3. The parameters of the system (1) are set as A = 10, c = 0.5, u = 0.2. It is noted that the individual plant virus equilibrium I 2 is determined by d and the basic reproduction number R 2 . The rate of death d has significant impacts on the I 2 if d is small enough, as shown in Figure 3a. It can be seen that the relationship between the rate of virus mutation ε and I * 2 is a linear, as shown in Figure 3b. The relationships among R 1 , R 2 and I 1 , I 2 are shown in Figure 3c,d when the system (1) satisfies the conditions of the existence of endemic equilibrium (R 1 > 1, R 1 > R 2 ). plant virus equilibrium ̅ 2 is determined by and the basic reproduction number 2 . The rate of death has significant impacts on the ̅ 2 if is small enough, as shown in Figure 3a. It can be seen that the relationship between the rate of virus mutation and 2 * is a linear, as shown in Figure 3b. The relationships among 1 , 2 and 1 , 2 are shown in Figure 3c,d when the system (1) satisfies the conditions of the existence of endemic equilibrium ( 1 > 1, 1 > 2 ). Figure 3. Relationship between I 1 , I 2 , the basic reproductive number: (a) Relationship between I 2 and d, R 2 ; (b) relationship between I * 2 and ε; (c) relationship between I 1 and R 1 , R 2 ; (d) relationship between I 2 and R 1 , R 2 If R 1 < 1 and R 2 < 1, the parameters are set as β 1 = 0.003, β 2 = 0.002, γ 1 = 0.2, γ 2 = 0.4, c = 0.5, u = 0.2, d = 0.1, ε = 0.2 and A = 10. Hence, the basic reproductive numbers are calculated as R 1 = 0.1125 < 1 and R 2 = 0.075 < 1. N = S + I 1 + I 2 + L = 100. The disease-free equilibrium E 0 (S 0 , 0, 0, L 0 ) can been calculated as S 0 = 37.5, L 0 = 62.5, I 1 = 0, I 2 = 0, as shown in Figure 4.
As shown in Figure 4a-c, according to Theorem 2, the initial values S(0), I 1 (0), I 2 (0) and L(0) do not have any influence on system stability and the system (1) will be stable at the disease-free equilibrium E 0 eventually.
According to Theorem 3, it is noted that the system (1) will be stable at the individual plant virus equilibrium eventually. From Figure 5a-c, the initial values do not have any influence on system stability and the system (1) will be stable at the individual plant virus equilibrium E ̅ eventually. According to Theorem 3, it is noted that the system (1) will be stable at the individual plant virus equilibrium eventually. From Figure 5a-c, the initial values do not have any influence on system stability and the system (1) will be stable at the individual plant virus equilibrium E eventually.
Obviously, as shown in Figure 6a-c, it is can be seen that R 1 = 1.125 > 1 and R 1 > R 2 = 0.75, and the initial values do not have any influence on system stability and the system (1) will be stable at the endemic equilibrium eventually.
Briefly, the relationship between the system and parameters is shown as a simulation diagram. The rate of ε has a great influence on the stable solution of the system, and other parameters have a great influence on the basic reproduction number to affect the equilibrium solution. Thus, a good system can be constructed by setting reasonable parameters. With the endemic equilibrium * ( * , 1 * , 2 * , * ), we can obtain * =26.7, 1 * = 13.3, 2 * , = 10, * = 50. The simulation diagram is shown in Figure 6.
Obviously, as shown in Figure 6a-c, it is can be seen that 1 = 1.125 > 1 and 1 > 2 = 0.75, and the initial values do not have any influence on system stability and the system (1) will be stable at the endemic equilibrium eventually.
Briefly, the relationship between the system and parameters is shown as a simulation diagram. The rate of has a great influence on the stable solution of the system, and other parameters have a great influence on the basic reproduction number to affect the equilibrium solution. Thus, a good system can be constructed by setting reasonable parameters.

Optimal Strategy Simulation
As for the optimal control strategy, the control variable set U(γ 1 , γ 2 , u) is calculated by Pontryagin's maximum principle, and the numerical result is calculated by implementing a fourth order Runge-Kutta method [33]. Firstly, for the system (1), the nodes' number of each state can be obtained by initializing control variable values, the time interval is [0, t f ] and the transversality condition Equation (35) is satisfied. It is noted that the numerical result of adjoint variables can be calculated by the system (34). Then, by updating U(γ 1 , γ 2 , u) iteratively, the control cost converges to a limited value.

Optimal Strategy Simulation
As for the optimal control strategy, the control variable set ( 1 , 2 , ) is calculated by Pontryagin's maximum principle, and the numerical result is calculated by implementing a fourth order Runge-Kutta method [33]. Firstly, for the system (1), the nodes' number of each state can be obtained by initializing control variable values, the time interval is  As shown in Figure 8 considering the optimal control, the rate of cleaning I 1 virus (i.e., γ 1 ) increases to 1 gradually; the rate of cleaning I 2 virus (i.e., γ 2 ) is increased to 0.0226 and the rate of recovery energy (i.e., u) increases to 0.0067.
It is noted that the cost under the optimal strategy is greatly reduced, as shown in Figure 9, which is equal to 86.86 at the terminal time t f , shown in Figure 9a. As shown in Figure 9b, the values of I 1 and I 2 are smaller than that without optimization control. It is obvious that the performance under the optimal strategy is superior to that without optimization control. As shown in Figure 8 considering the optimal control, the rate of cleaning 1 virus (i.e., 1 ) increases to 1 gradually; the rate of cleaning 2 virus (i.e., 2 ) is increased to 0.0226 and the rate of recovery energy (i.e., ) increases to 0.0067. It is noted that the cost under the optimal strategy is greatly reduced, as shown in Figure 9, which is equal to 86.86 at the terminal time , shown in Figure 9a. As shown in Figure 9b, the values of 1 and 2 are smaller than that without optimization control. It is obvious that the performance under the optimal strategy is superior to that without optimization control. As shown in Figure 8 considering the optimal control, the rate of cleaning 1 virus (i.e., 1 ) increases to 1 gradually; the rate of cleaning 2 virus (i.e., 2 ) is increased to 0.0226 and the rate of recovery energy (i.e., ) increases to 0.0067. It is noted that the cost under the optimal strategy is greatly reduced, as shown in Figure 9, which is equal to 86.86 at the terminal time , shown in Figure 9a. As shown in Figure 9b, the values of 1 and 2 are smaller than that without optimization control. It is obvious that the performance under the optimal strategy is superior to that without optimization control.

Conclusions and Future Work
In this paper, a novel epidemic model in WRSNs is proposed, including susceptible, infected, variant and low-energy states. The model considers rechargeable sensors and the

Conclusions and Future Work
In this paper, a novel epidemic model in WRSNs is proposed, including susceptible, infected, variant and low-energy states. The model considers rechargeable sensors and the virus mutation factor. By analyzing the basic reproductive number, the existence of equilibriums is first proved, and the local stability and global asymptotic stability are proved by the Lyapunov stability criterion. Meanwhile, the influence of the rate of virus mutation and the number of mutated nodes on the endemic equilibrium is revealed. In addition, an optimal strategy is proposed to minimize the numbers of the infected nodes and the virus mutation nodes, the cost of cleaning the virus and the cost of charging low-energy nodes. Finally, the numerical simulation validates the theoretical results.
Future research will attempt to consider the time delay to move closer to practical application. For WRSNs, it is a key topic of future work to consider the time delay of the virus mutation and the time delay of charging at the same time. Additionally, the model can become closer to reality when the time delay is proposed.