Dynamical Analysis and Optimal Control for a SEIR Model Based on Virus Mutation in WSNs

: With the rapid development of science and technology, the application of wireless sensor networks (WSNs) is more and more widely. It has been widely concerned by scholars. Viruses are one of the main threats to WSNs. In this paper, based on the principle of epidemic dynamics, we build a SEIR propagation model with the mutated virus in WSNs, where E nodes are infectious and cannot be repaired to S nodes or R nodes. Subsequently, the basic reproduction number R 0 , the local stability and global stability of the system are analyzed. The cost function and Hamiltonian function are constructed by taking the repair ratio of infected nodes and the repair ratio of mutated infected nodes as optimization control variables. Based on the Pontryagin maximum principle, an optimal control strategy is designed to effectively control the spread of the virus and minimize the total cost. The simulation results show that the model has a guiding signiﬁcance to curb the spread of mutated virus in WSNs.


Research Background
WSN is a self-organizing network system composed of a large number of sensor nodes deployed in a specific area through wireless communication [1,2]. It can achieve dynamic and real-time information monitoring, sensing and collection of monitoring objects in the network coverage area through the cooperation between nodes. WSNs are widely used in various fields, such as manufacturing, security monitoring and even military fields [3,4].
However, due to the strong openness of nodes, it is easy to be attacked by various types of viruses [5,6]. The virus against wireless devices can be transmitted directly from the device to the device using wireless communication technologies, such as Bluetooth [7][8][9]. Once the virus infects too many nodes in the network, it will lead to network interruption and paralysis. The virus may mutate in the computer network [10,11], its complexity and unpredictability are unmatched by the non-mutated virus, such as CIH virus. In addition, viruses may mutate in wireless sensor networks because of the similarities between computers network and WSNs. However, most of the research on the mutated virus focuses on infectious diseases [12,13].
Therefore, it is urgent and important to study the transmission of mutated virus in WSNs based on the principle of epidemic dynamics [14,15]. Through the establishment of mutated virus propagation model, we can have a deeper understanding of the mutated virus propagation behavior in the network. Then, we can take effective control and response strategies before the large-scale spread of the mutated virus, greatly reduce the harm of the mutated virus to the network.

Related Work
The establishment and dynamic analysis of virus propagation model based on WSNs has been concerned by many scholars. Kephart, J.O. and White, S.R. [16,17] first used the epidemiological model to study and predict the spread of virus in the network. Subsequently, many epidemiological models [18][19][20] are proposed and applied to WSNs. Due to the different characteristics and situations of infection, the models are constantly updated, and targeted epidemic models are proposed. Different authors proposed different epidemic models to study the dynamics of malicious objects propagation and control WSNs.
In the past, Mishra B.K. et al. [21] considered the effects of exposure status, vaccination and reinfection process and proposed the SEIRS-V model. However, this paper only proves the equilibrium of disease-free equilibrium. In Reference [22], Mishra B.K. and Srivastava S.K. proposed a quarantine model of worm propagation behavior and studied the stability of the model based on the basic reproduction number. The quarantine is a way to isolate infected nodes from the network so as not to infect susceptible nodes. In Reference [23], the authors studied an SLBRS computer virus model with two delays and obtain sufficient conditions for the existence of local Hopf bifurcation and the local stability. Subsequently, in Reference [24], Srivastava P.K. et al. proposed a worm propagation model (SEIQR). The model described the influence of quarantined state and recovery state on worm propagation. Moreover, they found the equilibrium point in disease-free and endemic cases, and proved the local stability and global stability, respectively. Besides, on the basis of the SEIQR model, Mishra B.K. and Tyagi I [25] added the vaccination and proposed SEIQRS-V model. Based on the SIR model, Zhu L. et al. [26] studied the nonlinear propagation model of malware and obtained the sufficient conditions for the existence of Hopf bifurcation and the local stability.
In recent years, the authors of Reference [27] proposed the SEIRS-D agent-based model and the characteristics of the model can be adjusted to more realistic values based on the environment by integrating new elements. Huang D.W. et al. [28] proposed a model with patch injection mechanism (SIPS) to evaluate the performance of patch injection rate in inhibiting computer infection. In Reference [29], the authors considered the communication radius and distribution density of nodes and proposed a delayed SEIQRS-V model. In addition, they obtained the sufficient conditions for the local stability and the existence of Hopf bifurcation. In addition, a variety of models have been proposed, such as SILRD [30], I2S2R [31], SEIQRV [32], VCQPS [33], etc.
In reality, after a period of time, the virus will be included in the database of the anti-virus program. Then, it will be intercepted by the anti-virus program and lose its infectivity. However, in order to continue spreading the virus, the virus maker will also update the virus. Therefore, it is necessary to study and control the mutated virus in the network. However, most of the researches on mutated viruses focus on how to detect them. As early as in [34], the authors proposed a computer virus detection method based on neural network, which can detect and destroy the mutated virus. In Reference [11], EVs and ECs are proposed to detect and analyze computer viruses for the first time. This method can effectively detect all mutations of corresponding viruses with one signature. In Reference [35], Rad B.B. developed a Hidden Markov Model that can recognize and detect other mutations in the same metamorphic virus family. It is rare to study the propagation behavior of mutated virus in the network. This paper enriches the content of this aspect.
In this paper, we will apply the theories of epidemic dynamics and optimal control to study the security of WSNs. Optimal control [36][37][38], as a method of studying the optimal dynamic strategy, has also been widely applied in wireless sensor networks. In Reference [39], the attack behavior of malicious programs is studied by combining the epidemic model and the loss equation. The maximum attack behavior of virus is discussed by taking the node communication radius and the media scanning rate as the optimal control variables. In Reference [40], the best defense strategy model between malware and WSNs is constructed and it can achieve the best defense effect with less resource consumption. In Reference [41], Huang Y.H. et al. proposed an anti-virus weight adaptive strategy, which effectively alleviated the spread of the virus in the network. In Reference [26], the optimal dynamic strategy for the system and malware based on Pontryagin Maximum Principle is given and it can reduce the total cost and inhibit the spread of malware to a certain extent.

Contributions
This paper is the first time to apply mutation model to WSNs to study the spread and control of virus, and proposes an SEIR epidemic model based on virus mutation in WSNs. Our contributions are as follows: 1.
An SEIR model based on virus mutation is established to describe the propagation process of mutated virus in WSNs.

2.
Calculating the basic reproduction number R 0 of the improved model by the next generation matrix method. Besides, the local and global stability of the two equilibria are proved and simulated by the Routh criterion and the Lyapunov stability method. Moreover, the influence of the repair rate γ 1 and γ 2 on the basic reproduction number is also revealed in the simulation part.

3.
Based on the Pontryagin maximum principle, the optimal control variable pairs of the repair ratio of infected nodes and the repair ratio of mutated infected nodes are obtained. The simulation results show that the optimal control strategy ensures the security of wireless sensor networks and minimizes the maintenance cost.
The rest of the paper is as follows. The compositions of the model are introduced in Section 2. The local and global stability is proved and the optimal strategy is designed in Section 3. The simulation results, which support our theoretical predictions, are shown in Section 4. The relevant conclusions are given in Section 5.

Dynamic Equation
The proposed model is global and deterministic, and the sensor nodes are divided into six parts: susceptible (S), exposed (E), infected (I 1 ), mutated infected (I 2 ), recovered (R) and death (D). The relationship between the six compartments is depicted in Figure 1. S nodes are vulnerable to viruses; E nodes are compromised with the attacker, but they work normally; I 1 nodes transform from E nodes, but they are unable to work normally; I 2 nodes are obtained by virus mutation from I 1 nodes, they are also unable to work normally; R nodes are immune to both pre mutated and post mutated viruses; D nodes are completely damaged. Specifically, E, I 1 and I 2 nodes all can infect other nodes. Now, let us make the following assumptions. The number of nodes increases at rate b, they all belong to S nodes. Meanwhile, new E nodes are generated with λ1I1(t)S(t), λ2I2(t)S(t) and λ3E(t)S(t). The parameters λ1, λ2 and λ3 represent the transmission rate of I1, I2 and E nodes, respectively. Because E nodes are mainly latent and the mutated virus is more infectious, λ2 is greater than λ1 and λ1 is greater than λ3. E nodes are converted to I1 nodes at rate ε, and I1 nodes mutate to I2 nodes The number of nodes increases at rate b, they all belong to S nodes. Meanwhile, new E nodes are generated with λ 1 I 1 (t)S(t), λ 2 I 2 (t)S(t) and λ 3 E(t)S(t). The parameters λ 1 , λ 2 and λ 3 represent the transmission rate of I 1 , I 2 and E nodes, respectively. Because E nodes are mainly latent and the mutated virus is more infectious, λ 2 is greater than λ 1 and λ 1 is greater than λ 3 . E nodes are converted to I 1 nodes at rate ε, and I 1 nodes mutate to I 2 nodes at rate µ. Because I 1 and I 2 nodes are unable to work normally, they can be detected by the anti-virus program, and repaired at rate γ 1 and γ 2 , respectively. In addition, the five compartments S, E, I 1 , I 2 and R have the same mortality b. All parameters are greater than 0. The parameters are summarized in Table 1.

Symbol Description b
Birth rate or Death rate λ 1 Transmission rate of I 1 nodes λ 2 Transmission rate of I 2 nodes λ 3 Transmission rate of E nodes ε Probability at which E nodes are converted to I 1 nodes µ Probability of virus mutation γ 1 Repair rate of I 1 nodes γ 2 Repair rate of I 2 nodes S(t), E(t), I 1 (t), I 2 (t), and R(t) are the ratio of susceptible, exposed, infected, mutated infected and recovered nodes at t. The novel dynamical system is given as follows: (1) Because R(t) and D(t) are independent of the other four equations, the following system is considered: (2) For the system (2), we suppose that N(t) = S(t) + E(t) + I 1 (t) + I 2 (t) and it satisfies which implies N(t) < 0. Therefore, we get 0 < N ≤ 1. we get a feasible region as follow: Therefore, the solution of system (2) is bounded and independent of the initial condition. Therefore, D is an invariant set. In addition, 0 < N ≤ 1, D also is a positive set. Therefore, we will consider the stability of system (2) on the set D

Calculation of the Equilibrium Point and the Basic Reproduction Number
In this subsection, two equilibrium points of the system (2) on the set D are calculated. The first point is the disease-free equilibrium point N(t) = S(t) + E(t) + I 1 (t) + I 2 (t) P 0 = (S 0 , E 0 , I 1 0 , I 2 0 ), and The second point is the endemic equilibrium point P * = (S * , E * , I 1 * , I 2 * ), and Furthermore, by the next generation matrix method, the basic reproduction number R 0 is obtained. Set and The next generation matrix is FV −1 , the basic reproduction number R 0 is its spectral radius.

Dynamic Analysis and Optimal Strategy
In this section, the local and global stability of two equilibrium points are proved by using the Routh criterion [42] and the Lyapunov stability method [43]. The existence of local stability and global stability indicates that the development of this infectious disease will not appear large-scale repeated infection, and will eventually maintain a static equilibrium with the passage of time. Moreover, based on Pontryagin maximum principle [44], an optimal strategy is proposed to control the spread of virus with minimum cost.

Subs Stability Analysis of P 0
Theorem 1. When R 0 < 1, the disease-free equilibrium point P 0 is locally asymptotically stable.
Proof. First, we can get the disease-free equilibrium point P 0 = (1, 0, 0, 0) according to the system (2). Thus, the Jacobian matrix of the disease-free equilibrium point P 0 is: Thus, −b is one of the eigenvalues of J(P 0 ), so just consider the following Jacobian matrix: The characteristic polynomial of (15) is where and Moreover, a simple calculation shows a 1 a 2 − a 0 > 0. Thus, according to the Routh criterion, P 0 is locally asymptotically stable if R 0 < 1. Theorem 2. When R 0 < 1, the disease-free equilibrium point P 0 is globally asymptotically stable.
Proof. We can prove it by the Lyapunov stability method. Consider the following Lyapunov function: We have: According to mean inequality n √ a 1 a 2 · · · a n ≤ a 1 +a 2 +···+a n n , S(t) + 1 S(t) ≥ 2. Then, According to the Lyapunov stability theorem, P 0 is globally asymptotically stable.

Stability Analysis of P*
Theorem 3. When R 0 > 1, the epidemic equilibrium point P * is locally asymptotically stable.
Theorem 4. When R 0 > 1, the epidemic equilibrium point P* is globally asymptotically stable.

Optimal Strategy
In this paper, our goal is not only to effectively control the spread of the virus, but also to make the cost of sensor node repair as low as possible. Selecting the repair rate of I 1 nodes γ 1 and the repair rate of I 2 nodes γ 2 as the control variable, the cost function is as follows: In this formula, c 1 and c 2 indicate the repair cost parameter of I 1 nodes and I 2 nodes, respectively. Here, t f represents the terminal moment. The terms c 1 γ 1 2 (t)I 1 2 (t) and c 2 γ 2 2 (t)I 2 2 (t) describe the repair cost of I 1 nodes and I 2 nodes at time t. The Equation (34) shows that the sum of the number of E nodes, I 1 nodes and I 2 nodes at the terminal moment and the cost from the start time to the terminal moment is the least by controlling γ 1 (t) and γ 2 (t). We apply the Pontryagin maximum principle to solve the optimization problem and consider the following Hamiltonian: .

Simulation
This section will be divided into three parts, the first two parts are to further verify the stability of the system, and the latter part is to show the effectiveness of the optimal control. Simulation parameters are obtained by simulating different scenes.

Stable Analysis of Disease-Free Equilibrium
In this subsection, we suppose λ 1 = 0.1, λ 2 = 0.2, λ 3 = 0.05, γ 1 = 0.1, γ 2 = 0.05, µ = 0.05, ε = 0.05, and b = 0.1. According to the Formula (11), we can get R 0 = 0.56 < 1. Theoretically, S(∞) = 1, E(∞) = 0, I 1 (∞) = 0, and I 2 (∞) = 0 according to the Formula (3). The simulation results are illustrated in Figure 2, and it is consistent with Theorems 1 and 2. In the three-dimensional graph with I2 nodes scale as X-axis, E nodes scale as Y-axis and S nodes scale as Z-axis, the curves finally converge to (0, 0, 1). In Figure 2c, in the beginning, there contain only I2 and I1 nodes when the curves start from the X-axis; there contain only E and I1 nodes when the curves start from the Y-axis; and there contain only S and I1 nodes when the curves start from the Zaxis. In Figure 2d, the curves contain only I2 and E nodes in the X-Y plane; the curves contain only I2 and S nodes in the X-Z plane; and the curves contain only E and S nodes in the Y-Z plane. In summary, those curves will gather together, and then finally converge to (0, 0, 1) in a variety of cases of R0 < 1, as shown in Figure 2. These results are according to Theorems 1 and 2.

Stable Analysis of Epidemic Equilibrium
In this subsection,  In Figure 2a,c, the curves all start from the X, Y and Z axes, and the curves in Figure 2b,d all start from the hypotenuse of the X-Y, X-Z and Y-Z planes.
As shown in Figure 2a,b in the three-dimensional graph with I 1 nodes scale as X axis, E nodes scale as Y axis and S nodes scale as Z axis, the curves finally converge to (0, 0, 1). In Figure 2a, there are only I 1 and I 2 nodes in the beginning if the curves start from the X-axis; there are only E and I 2 nodes in the beginning if the curves start from the Y-axis; and there are only S and I 2 nodes in the beginning if the curves start from the Z-axis. These assumptions ensure that the virus exists at the beginning. In Figure 2b, the ratio sum of I 1 nodes and E nodes is 1 in the X-Y plane at the initial moment; the ratio sum of I 1 nodes and S nodes is 1 in the X-Z plane at the initial moment; and the ratio sum of E nodes and S nodes is 1 in the Y-Z plane at the initial moment. Figure 2c,d are similar to Figure 2a,b. In the three-dimensional graph with I 2 nodes scale as X-axis, E nodes scale as Y-axis and S nodes scale as Z-axis, the curves finally converge to (0, 0, 1). In Figure 2c, in the beginning, there contain only I 2 and I 1 nodes when the curves start from the X-axis; there contain only E and I 1 nodes when the curves start from the Y-axis; and there contain only S and I 1 nodes when the curves start from the Z-axis. In Figure 2d, the curves contain only I 2 and E nodes in the X-Y plane; the curves contain only I 2 and S nodes in the X-Z plane; and the curves contain only E and S nodes in the Y-Z plane. In summary, those curves will gather together, and then finally converge to (0, 0, 1) in a variety of cases of R 0 < 1, as shown in Figure 2. These results are according to Theorems 1 and 2.

Stable Analysis of Epidemic Equilibrium
In this subsection, λ 1 = 0.2, λ 2 = 0.3, λ 3 = 0.05, γ 1 = 0.1, γ 2 = 0.05, µ = 0.05, ε = 0.1, and b = 0.05. It is easily obtained that R 0 = 1.5 > 1. According to the Formulas (4) to (6), we know S(∞) = 0.667, E(∞) = 0.111, I 1 (∞) = 0.056, and I 2 (∞) = 0.028. The simulation results conform to Theorems 3 and 4, as shown in Figure 3. The assumptions are the same as in Figure 2. In Figure 3a,b, I1 nodes ratio is the Xaxis; the E nodes ratio is the Y-axis; and the S node ratio is the Z-axis. The simulation results show that the curve converges to (0.056, 0.111, 0.667) from the axis and the hypotenuse, respectively. In Figure 3c,d, I2 nodes ratio is the X-axis; E nodes ratio is the Y-axis; and S node ratio is the Z-axis. The simulation results show that the curve converges to (0.028, 0.111, 0.667) from the axis and the hypotenuse, respectively.

Optimal Control
In this subsection, four parts will be considered: the evolution of sensor nodes, the total cost, the evolution of control variables, and the influence of control variables on R0. The assumptions are the same as in Figure 2. In Figure 3a,b, I 1 nodes ratio is the X-axis; the E nodes ratio is the Y-axis; and the S node ratio is the Z-axis. The simulation results show that the curve converges to (0.056, 0.111, 0.667) from the axis and the hypotenuse, respectively. In Figure 3c,d, I 2 nodes ratio is the X-axis; E nodes ratio is the Y-axis; and S node ratio is the Z-axis. The simulation results show that the curve converges to (0.028, 0.111, 0.667) from the axis and the hypotenuse, respectively.

Optimal Control
In this subsection, four parts will be considered: the evolution of sensor nodes, the total cost, the evolution of control variables, and the influence of control variables on R 0 . Suppose that except for b = 0.1, the other parameters remain unchanged, as described in Section 4.2. Besides, set The optimal control calculation process is shown in Figure 4.  Firstly, the control variables γ1, γ2 are be given an initial guess, and then Runge-Kutta method is used to solve the numerical solution of the state variables of the system (2).
With the values of the obtained state variables, the initial values of the optimization variables and transversal conditions, the numerical solutions of the costate differential equations (36) are obtained by inverse integral. In order to better curb the spread of the virus, the value of control is always one when calculating the covariate. Finally, the values of the optimization control variables are obtained according to the Equation (39). Repeat the above process until the curves of the control variable are no longer changed. Firstly, the control variables γ 1 , γ 2 are be given an initial guess, and then Runge-Kutta method is used to solve the numerical solution of the state variables of the system (2). With the values of the obtained state variables, the initial values of the optimization variables and transversal conditions, the numerical solutions of the costate differential Equation (36) are obtained by inverse integral. In order to better curb the spread of the virus, the value of control is always one when calculating the covariate. Finally, the values of the optimization control variables are obtained according to the Equation (39). Repeat the above process until the curves of the control variable are no longer changed.

Evolution of Sensor Nodes
In this part, we will discuss the evolution of sensor nodes under non-optimal control and optimal control, as shown in Figure 5. One of our optimization objectives is to effectively control the ratio of nodes with virus at the terminal moment. By comparing Figure 5b and Figure 5a, it can be seen that the nodes with virus under the optimal control are eliminated faster than those with fixed repair rate.

Total Cost and Control Variables
Another control goal is to reduce the repair cost on the basis of curbing the spread of the virus. The actual total cost of repair equation is In this subsection, the fixed repair rates of 1 was selected as the non-optimal control group. By applying for approximate method, the total cost as shown in Figure 6.
When γ1 and γ2 are always 1, the ratio of nodes with virus decreases fastest, the evolution of sensor nodes as shown in Figure 7 is similar to the evolution of sensor nodes under the optimal control as shown in Figure 4a. It further shows that the goal of effectively controlling the ratio of nodes with virus is achieved. In addition, as shown in Figure  6, it can be seen that the total cost under the optimal control is lower than the total cost under non-optimal control. Specifically, the cost finally reaches 0.801 under the non-optimal control and the cost finally reaches 0.674 under optimal control. Therefore, the optimal strategy not only effectively controls the spread of the virus, but also minimizes the cost of the security strategy. One of our optimization objectives is to effectively control the ratio of nodes with virus at the terminal moment. By comparing Figures 5b and 5a, it can be seen that the nodes with virus under the optimal control are eliminated faster than those with fixed repair rate.

Total Cost and Control Variables
Another control goal is to reduce the repair cost on the basis of curbing the spread of the virus. The actual total cost of repair equation is C = t f 0 [c 1 γ 1 (t) + c 2 γ 2 (t)] dt. In this subsection, the fixed repair rates of 1 was selected as the non-optimal control group. By applying for approximate method, the total cost as shown in Figure 6. One of our optimization objectives is to effectively control the ratio of nodes with virus at the terminal moment. By comparing Figure 5b and Figure 5a, it can be seen that the nodes with virus under the optimal control are eliminated faster than those with fixed repair rate.

Total Cost and Control Variables
Another control goal is to reduce the repair cost on the basis of curbing the spread of the virus. The actual total cost of repair equation is In this subsection, the fixed repair rates of 1 was selected as the non-optimal control group. By applying for approximate method, the total cost as shown in Figure 6.
When γ1 and γ2 are always 1, the ratio of nodes with virus decreases fastest, the evolution of sensor nodes as shown in Figure 7 is similar to the evolution of sensor nodes under the optimal control as shown in Figure 4a. It further shows that the goal of effectively controlling the ratio of nodes with virus is achieved. In addition, as shown in Figure  6, it can be seen that the total cost under the optimal control is lower than the total cost under non-optimal control. Specifically, the cost finally reaches 0.801 under the non-optimal control and the cost finally reaches 0.674 under optimal control. Therefore, the optimal strategy not only effectively controls the spread of the virus, but also minimizes the cost of the security strategy.  When γ 1 and γ 2 are always 1, the ratio of nodes with virus decreases fastest, the evolution of sensor nodes as shown in Figure 7 is similar to the evolution of sensor nodes under the optimal control as shown in Figure 4a. It further shows that the goal of effectively controlling the ratio of nodes with virus is achieved. In addition, as shown in Figure 6, it can be seen that the total cost under the optimal control is lower than the total cost under non-optimal control. Specifically, the cost finally reaches 0.801 under the non-optimal control and the cost finally reaches 0.674 under optimal control. Therefore, the optimal strategy not only effectively controls the spread of the virus, but also minimizes the cost of the security strategy. The results of optimizing variables are shown in Figure 8. In the early stage of security measures, because the ratio of I1 nodes is higher and the repair cost c2 of I2 nodes is larger than the repair cost c1 of I1 nodes, so γ2 < γ1. Then, due to the decreasing ratio of I1 nodes and the higher transmission rate λ2 of I2 nodes, I2 nodes become the biggest threat, so γ2 > γ1. In the later stage, our goal is to diminish the number of nodes with virus, so the removal intensity of nodes with virus should be improved to make γ1 = γ2 = 1. Figure 8 provides the optimal combination of I1 nodes repair rate γ1 and I2 nodes repair rate γ2 at all times.

Influence of Control Variables on the Basic Reproduction Number R0
The influence of control variables on R0 is shown in Figure 9. The simulation results show that R0 decreases with the increase of the repair rate γ1 and γ2. The lower the R0, the faster the virus is eliminated. When the control variables converge to 1, then R0 < 1, as shown in Figure 8. Thus, the system will eventually achieve the disease-free equilibrium under optimal control. In addition, compared to γ2, γ1 has a greater impact on R0 under The results of optimizing variables are shown in Figure 8. In the early stage of security measures, because the ratio of I 1 nodes is higher and the repair cost c 2 of I 2 nodes is larger than the repair cost c 1 of I 1 nodes, so γ 2 < γ 1 . Then, due to the decreasing ratio of I 1 nodes and the higher transmission rate λ 2 of I 2 nodes, I 2 nodes become the biggest threat, so γ 2 > γ 1 . In the later stage, our goal is to diminish the number of nodes with virus, so the removal intensity of nodes with virus should be improved to make γ 1 = γ 2 = 1. Figure 8 provides the optimal combination of I 1 nodes repair rate γ 1 and I 2 nodes repair rate γ 2 at all times. The results of optimizing variables are shown in Figure 8. In the early stage of secu rity measures, because the ratio of I1 nodes is higher and the repair cost c2 of I2 nodes i larger than the repair cost c1 of I1 nodes, so γ2 < γ1. Then, due to the decreasing ratio of I nodes and the higher transmission rate λ2 of I2 nodes, I2 nodes become the biggest threat so γ2 > γ1. In the later stage, our goal is to diminish the number of nodes with virus, so th removal intensity of nodes with virus should be improved to make γ1 = γ2 = 1. Figure 8 provides the optimal combination of I1 nodes repair rate γ1 and I2 nodes repair rate γ2 a all times.

Influence of Control Variables on the Basic Reproduction Number R0
The influence of control variables on R0 is shown in Figure 9. The simulation result show that R0 decreases with the increase of the repair rate γ1 and γ2. The lower the R0, th faster the virus is eliminated. When the control variables converge to 1, then R0 < 1, a shown in Figure 8. Thus, the system will eventually achieve the disease-free equilibrium under optimal control. In addition, compared to γ2, γ1 has a greater impact on R0 unde these parameters as described in Section 4.3. The influence of control variables on R 0 is shown in Figure 9. The simulation results show that R 0 decreases with the increase of the repair rate γ 1 and γ 2 . The lower the R 0 , the faster the virus is eliminated. When the control variables converge to 1, then R 0 < 1, as shown in Figure 8. Thus, the system will eventually achieve the disease-free equilibrium under optimal control. In addition, compared to γ 2 , γ 1 has a greater impact on R 0 under these parameters as described in Section 4.3.
Mathematics 2020, 17, x FOR PEER REVIEW 17 of 19 Figure 9. Influence of control variables on the basic reproduction number R0.

Conclusions
In this paper, we use epidemiology to propose an improved SEIR model based on virus mutation, where E nodes are infectious and cannot be repaired to S nodes or R nodes, describing the propagation of mutated virus in WRSNs. Meanwhile, the basic reproduction number of the improved model is calculated by the next generation matrix method and the local and global stability of the two equilibrium points are proved by analyzing the stability of the model. Besides, in order to minimize the ratio of nodes with virus and the total cost associated with the measure, an optimal control strategy is proposed by optimizing the combination of the repair rates γ1 and γ2. In addition, the influence of control variables on the basic reproduction number is analyzed.
The simulation results verify the correctness of the stability theory. Further, the optimal control strategy is proved to be effective in controlling the ratio of nodes with virus and minimizing the maintenance cost. Figure 9 shows that the virus can be eliminated by adjusting the repair rate γ1 and γ2, but it needs to rise to a certain threshold. The SEIR propagation model with the mutated virus provides new insights into the virus propagation in WSNs. However, our model only considers one case of mutated virus, which is not comprehensive enough. In the future work, we will consider the case of time delay, and apply stochastic modeling and advanced mathematical theory when the ability permits. We hope the work of this paper can give some enlightenment to the relevant researchers.

Conclusions
In this paper, we use epidemiology to propose an improved SEIR model based on virus mutation, where E nodes are infectious and cannot be repaired to S nodes or R nodes, describing the propagation of mutated virus in WRSNs. Meanwhile, the basic reproduction number of the improved model is calculated by the next generation matrix method and the local and global stability of the two equilibrium points are proved by analyzing the stability of the model. Besides, in order to minimize the ratio of nodes with virus and the total cost associated with the measure, an optimal control strategy is proposed by optimizing the combination of the repair rates γ 1 and γ 2 . In addition, the influence of control variables on the basic reproduction number is analyzed.
The simulation results verify the correctness of the stability theory. Further, the optimal control strategy is proved to be effective in controlling the ratio of nodes with virus and minimizing the maintenance cost. Figure 9 shows that the virus can be eliminated by adjusting the repair rate γ 1 and γ 2 , but it needs to rise to a certain threshold. The SEIR propagation model with the mutated virus provides new insights into the virus propagation in WSNs. However, our model only considers one case of mutated virus, which is not comprehensive enough. In the future work, we will consider the case of time delay, and apply stochastic modeling and advanced mathematical theory when the ability permits. We hope the work of this paper can give some enlightenment to the relevant researchers.