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

Abstract

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₀, 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 significance to curb the spread of mutated virus in WSNs.

1. Introduction

1.1. 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.

1.2. 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.

1.3. 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:

• An SEIR model based on virus mutation is established to describe the propagation process of mutated virus in WSNs.
• Calculating the basic reproduction number R₀ 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 γ₁ and γ₂ on the basic reproduction number is also revealed in the simulation part.
• 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.

2. Modeling

2.1. Dynamic Equation

The proposed model is global and deterministic, and the sensor nodes are divided into six parts: susceptible (S), exposed (E), infected (I₁), mutated infected (I₂), 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₁ nodes transform from E nodes, but they are unable to work normally; I₂ nodes are obtained by virus mutation from I₁ 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₁ and I₂ 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 λ₁I₁(t)S(t), λ₂I₂(t)S(t) and λ₃E(t)S(t). The parameters λ₁, λ₂ and λ₃ represent the transmission rate of I₁, I₂ and E nodes, respectively. Because E nodes are mainly latent and the mutated virus is more infectious, λ₂ is greater than λ₁ and λ₁ is greater than λ₃. E nodes are converted to I₁ nodes at rate ε, and I₁ nodes mutate to I₂ nodes at rate μ. Because I₁ and I₂ nodes are unable to work normally, they can be detected by the anti-virus program, and repaired at rate γ₁ and γ₂, respectively. In addition, the five compartments S, E, I₁, I₂ and R have the same mortality b. All parameters are greater than 0. The parameters are summarized in

Table 1. The parameters of the model.

Symbol	Description
b	Birth rate or Death rate
λ1	Transmission rate of I1 nodes
λ2	Transmission rate of I2 nodes
λ3	Transmission rate of E nodes
ε	Probability at which E nodes are converted to I1 nodes
μ	Probability of virus mutation
γ1	Repair rate of I1 nodes
γ2	Repair rate of I2 nodes

.

S(t), E(t), I₁(t), I₂(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:

$$\{\begin{array}{l}\frac{dS(t)}{dt}=b-bS(t)-{\lambda }_{1}S(t)I_1(t)-{\lambda }_{2}S(t)I_2(t)-{\lambda }_{3}S(t)E(t),\\ \frac{dE(t)}{dt}={\lambda }_{1}S(t)I_1(t)+{\lambda }_{2}S(t)I_2(t)+{\lambda }_{3}S(t)E(t)-(\epsilon +b)E(t),\\ \frac{d{I}_1(t)}{dt}=\epsilon E(t)-({\gamma }_1+b+\mu )I_1\left(t\right),\\ \frac{d{I}_2(t)}{dt}=\mu I_1(t)-({\gamma }_2+b)I_2(t),\\ \frac{dR(t)}{dt}={\gamma }_1{I}_1(t)+{\gamma }_2{I}_2(t)-bR(t),\\ \frac{dD(t)}{dt}=bS(t)+bE(t)+b{I}_1(t)+b{I}_2(t)+bR(t).\end{array}$$

(1)

Because R(t) and D(t) are independent of the other four equations, the following system is considered:

$$\{\begin{array}{l}\frac{dS(t)}{dt}=b-bS(t)-{\lambda }_{1}S(t)I_1(t)-{\lambda }_{2}S(t)I_2(t)-{\lambda }_{3}S(t)E(t),\\ \frac{dE(t)}{dt}={\lambda }_{1}S(t)I_1(t)+{\lambda }_{2}S(t)I_2(t)+{\lambda }_{3}S(t)E(t)-(\epsilon +b)E(t),\\ \frac{d{I}_1(t)}{dt}=\epsilon E(t)-({\gamma }_1+b+\mu )I_1\left(t\right),\\ \frac{d{I}_2(t)}{dt}=\mu I_1(t)-({\gamma }_2+b)I_2(t).\end{array}$$

(2)

For the system (2), we suppose that N(t) = S(t) + E(t) + I₁(t) + I₂(t) and it satisfies

$$\frac{dN(t)}{dt}=b-bN(t)-{\gamma }_1{I}_1(t)-{\gamma }_2{I}_2(t),$$

(3)

which implies $\frac{dN(\mathrm{t})}{dt}\le b-bN(t)$, thus $\lim {\sup }_{t\to \infty }N(t)\le 1$. Also, if $N(t)>1$ then $\stackrel{.}{N}(t)<0$. Therefore, we get $0<N\le 1$. we get a feasible region as follow:

$$D=\left\{\left(S(t),E(t),I_1(t),I_2(t)\right)\in {\mathrm{R}}_{+}^4:0<S(t)+E(t)+I_1(t)+I_2(t)\le 1\right\}.$$

(4)

Therefore, the solution of system (2) is bounded and independent of the initial condition. Therefore, D is an invariant set. In addition, $0<N\le 1$, D also is a positive set. Therefore, we will consider the stability of system (2) on the set D

2.2. 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₁(t) + I₂(t) $P^0=(S^0,E^0,I_1{}^0,I_2{}^0)$, and

$$S^0=1,E^0=0,I_1{}^0=0,I_2{}^0=0.$$

(5)

The second point is the endemic equilibrium point $P^{*}=(S^{*},E^{*},I_1{}^{*},I_2{}^{*})$, and

$$S^{*}=\frac{\left({\gamma }_1+b+\mu \right)\left({\gamma }_2+b\right)\left(\epsilon +b\right)}{{\lambda }_1\epsilon \left({\gamma }_2+b\right)+{\lambda }_2\epsilon \mu +{\lambda }_3\left({\gamma }_1+b+\mu \right)\left({\gamma }_2+b\right)},$$

(6)

$$E^{*}=\frac{{\lambda }_{1}b\epsilon \left({\gamma }_2+b\right)+{\lambda }_{2}b\epsilon \mu +{\lambda }_{3}b\left({\gamma }_1+b+\mu \right)\left({\gamma }_2+b\right)-b\left({\gamma }_1+b+\mu \right)\left({\gamma }_2+b\right)\left(\epsilon +b\right)}{\epsilon \left({\gamma }_2+b\right)\left(\epsilon +b\right)\left[{\lambda }_1+\frac{{\lambda }_2\mu }{{\gamma }_2+b}+\frac{{\lambda }_3\left({\gamma }_1+b+\mu \right)}{\epsilon }\right]},$$

(7)

$$I_1{}^{*}=\frac{{\lambda }_{1}b\epsilon \left({\gamma }_2+b\right)+{\lambda }_{2}b\epsilon \mu +{\lambda }_{3}b\left({\gamma }_1+b+\mu \right)\left({\gamma }_2+b\right)-b\left({\gamma }_1+b+\mu \right)\left({\gamma }_2+b\right)\left(\epsilon +b\right)}{\left({\gamma }_1+b+\mu \right)\left({\gamma }_2+b\right)\left(\epsilon +b\right)\left[{\lambda }_1+\frac{{\lambda }_2\mu }{{\gamma }_2+b}+\frac{{\lambda }_3\left({\gamma }_1+b+\mu \right)}{\epsilon }\right]},$$

(8)

$$I_2{}^{*}=\frac{\mu \left[{\lambda }_{1}b\epsilon \left({\gamma }_2+b\right)+{\lambda }_{2}b\epsilon \mu +{\lambda }_{3}b\left({\gamma }_1+b+\mu \right)\left({\gamma }_2+b\right)-b\left({\gamma }_1+b+\mu \right)\left({\gamma }_2+b\right)\left(\epsilon +b\right)\right]}{\left({\gamma }_1+b+\mu \right){\left({\gamma }_2+b\right)}^2\left(\epsilon +b\right)\left[{\lambda }_1+\frac{{\lambda }_2\mu }{{\gamma }_2+b}+\frac{{\lambda }_3\left({\gamma }_1+b+\mu \right)}{\epsilon }\right]}.$$

(9)

Furthermore, by the next generation matrix method, the basic reproduction number R₀ is obtained.

Set

$$F=\left[\begin{array}{ccc}{\lambda }_3& {\lambda }_1& {\lambda }_2\\ 0& 0& 0\\ 0& 0& 0\end{array}\right]$$

(10)

and

$$V=\left[\begin{array}{ccc}\epsilon +b& 0& 0\\ -\epsilon & {\gamma }_1+b+\mu & 0\\ 0& -\mu & {\gamma }_2+b\end{array}\right].$$

(11)

The next generation matrix is FV⁻¹, the basic reproduction number R₀ is its spectral radius.

$$F{V}^{-1}=\left[\begin{array}{ccc}\frac{{\lambda }_1\epsilon \left({\gamma }_2+b\right)+{\lambda }_2\epsilon \mu +{\lambda }_3\left({\gamma }_1+b+\mu \right)\left({\gamma }_2+b\right)}{\left({\gamma }_1+b+\mu \right)\left({\gamma }_2+b\right)\left(\epsilon +b\right)}& \frac{{\lambda }_1\left({\gamma }_2+b\right)+{\lambda }_2\mu }{\left({\gamma }_1+b+\mu \right)\left({\gamma }_2+b\right)}& \frac{{\lambda }_2}{\left({\gamma }_2+b\right)}\\ 0& 0& 0\\ 0& 0& 0\end{array}\right],$$

(12)

$$R_0=\frac{{\lambda }_1\epsilon \left({\gamma }_2+b\right)+{\lambda }_2\epsilon \mu +{\lambda }_3\left({\gamma }_1+b+\mu \right)\left({\gamma }_2+b\right)}{\left({\gamma }_1+b+\mu \right)\left({\gamma }_2+b\right)\left(\epsilon +b\right)}.$$

(13)

3. 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.

3.1. Subs Stability Analysis of P⁰

Theorem 1.

When R₀ < 1, the disease-free equilibrium point P⁰ is locally asymptotically stable.

Proof. 

First, we can get the disease-free equilibrium point P⁰ = (1, 0, 0, 0) according to the system (2). Thus, the Jacobian matrix of the disease-free equilibrium point P⁰ is:

$$J\left(P^0\right)=\left[\begin{array}{cccc}-b& -{\lambda }_3& -{\lambda }_1& -{\lambda }_2\\ 0& -\epsilon -b& {\lambda }_1& {\lambda }_2\\ 0& \epsilon & -{\gamma }_1-b-\mu & 0\\ 0& 0& \mu & -{\gamma }_2-b\end{array}\right]$$

(14)

Thus, −b is one of the eigenvalues of J(P⁰), so just consider the following Jacobian matrix:

$$J^{\prime }\left(P^0\right)=\left[\begin{array}{ccc}-\epsilon -b& {\lambda }_1& {\lambda }_2\\ \epsilon & -{\gamma }_1-b-\mu & 0\\ 0& \mu & -{\gamma }_2-b\end{array}\right]$$

(15)

The characteristic polynomial of (15) is

$$P\left(\lambda \right)=a_3{\lambda }^3+a_2{\lambda }^2+a_1\lambda +a_0$$

(16)

where

$$a_3=1>0,$$

(17)

$$a_2=3b+{\gamma }_1+{\gamma }_2+\mu +\epsilon >0,$$

(18)

$$a_1=\frac{{\lambda }_1\epsilon \left({\gamma }_2+b\right)\left(1-R_0\right)+{\lambda }_2\epsilon \mu +{\lambda }_3\left({\gamma }_1+b+\mu \right)\left({\gamma }_2+b\right)}{R_0\left({\gamma }_2+b\right)}+\left({\gamma }_1+b+\mu \right)\left({\gamma }_2+b\right)+\left(\epsilon +b\right)\left({\gamma }_2+b\right)>0,$$

(19)

and

$$a_0=\frac{\left(1-R_0\right)\left[{\lambda }_1\epsilon \left({\gamma }_2+b\right)+{\lambda }_2\epsilon \mu \right]+{\lambda }_3\left({\gamma }_1+b+\mu \right)\left({\gamma }_2+b\right)}{R_0}>0.$$

(20)

Moreover, a simple calculation shows $a_1{a}_2-a_0>0$. Thus, according to the Routh criterion, P⁰ is locally asymptotically stable if R₀ < 1. □

Theorem 2.

When R₀ < 1, the disease-free equilibrium point P⁰ is globally asymptotically stable.

Proof. 

We can prove it by the Lyapunov stability method. Consider the following Lyapunov function:

$$V\left(t\right)=S\left(t\right)-1-\ln S\left(t\right)+E\left(t\right)+\frac{{\lambda }_1\left({\gamma }_2+b\right)+{\lambda }_2\mu }{\left({\gamma }_1+b+\mu \right)\left({\gamma }_2+b\right)}{I}_1(t)+\frac{{\lambda }_2}{{\gamma }_2+b}{I}_2\left(t\right)>0.$$

(21)

We have:

$$\begin{array}{l}\frac{dV(t)}{dt}=-b\left[S(t)+\frac{1}{S(t)}-2\right]+{\lambda }_1{I}_1\left(t\right)+{\lambda }_2{I}_2\left(t\right)+{\lambda }_{3}E\left(t\right)-\left(\epsilon +b\right)E\left(t\right)+\\ \begin{array}{cc}& \end{array}\begin{array}{c}\end{array}\frac{{\lambda }_1\left({\gamma }_2+b\right)+{\lambda }_2\mu +\frac{{\lambda }_3\left({\gamma }_1+b+\mu \right)\left({\gamma }_2+b\right)}{\epsilon }-\frac{{\lambda }_3\left({\gamma }_1+b+\mu \right)\left({\gamma }_2+b\right)}{\epsilon }}{\left({\gamma }_1+b+\mu \right)\left({\gamma }_2+b\right)}\left[\epsilon E\left(t\right)-\left({\gamma }_1+b+\mu \right)I_1\left(t\right)\right]+\\ \begin{array}{cc}& \end{array}\begin{array}{c}\end{array}\frac{{\lambda }_2}{\left({\gamma }_2+b\right)}\left[\mu I_1\left(t\right)-\left({\gamma }_2+b\right)I_2\left(t\right)\right]\\ \begin{array}{cc}& \end{array}=-b\left[S(t)+\frac{1}{S(t)}-2\right]-\left(\epsilon +b\right)E\left(t\right)+\frac{{\lambda }_1\left({\gamma }_2+b\right)+{\lambda }_2\mu +\frac{{\lambda }_3\left({\gamma }_1+b+\mu \right)\left({\gamma }_2+b\right)}{\epsilon }}{\left({\gamma }_1+b+\mu \right)\left({\gamma }_2+b\right)}\epsilon E\left(t\right)\\ \begin{array}{cc}\begin{array}{c}\end{array}& \begin{array}{c}\end{array}\end{array}\begin{array}{c}=\end{array}-b\left[S(t)+\frac{1}{S(t)}-2\right]-\left(\epsilon +b\right)\left[1-\frac{{\lambda }_1\left({\gamma }_2+b\right)+{\lambda }_2\mu +\frac{{\lambda }_3\left({\gamma }_1+b+\mu \right)\left({\gamma }_2+b\right)}{\epsilon }}{\left({\gamma }_1+b+\mu \right)\left({\gamma }_2+b\right)\left(\epsilon +b\right)}\epsilon \right]E\left(t\right)\\ \begin{array}{cc}& \end{array}=-b\left[S(t)+\frac{1}{S(t)}-2\right]-\left(\epsilon +b\right)\left(1-R_0\right)E\left(t\right)\end{array}$$

(22)

According to mean inequality $\sqrt[n]{a_1{a}_2\cdots a_n}\le \frac{a_1+a_2+\cdots +a_n}{n}$, $S\left(t\right)+\frac{1}{S\left(t\right)}\ge 2$. Then, $\frac{dV(t)}{dt}\le 0$ if $R_0\le 1$. It is noted that $\frac{dV\left(t\right)}{dt}|{}_{P^0}=0$. According to the Lyapunov stability theorem, P⁰ is globally asymptotically stable. □

3.2. Stability Analysis of P*

Theorem 3.

When R₀ > 1, the epidemic equilibrium point P * is locally asymptotically stable.

Proof. 

First of all, S*, E*, I₁* and I₂* are given in Equations (4)–(7). Thus, the Jacobian matrix at the epidemic equilibrium P* is:

$$J\left(P^{*}\right)=\left[\begin{array}{cccc}-{\lambda }_1{I}_1{}^{*}-{\lambda }_2{I}_2{}^{*}-{\lambda }_3{E}^{*}-b& -{\lambda }_3{S}^{*}& -{\lambda }_1{S}^{*}& -{\lambda }_2{S}^{*}\\ {\lambda }_1{I}_1{}^{*}+{\lambda }_2{I}_2{}^{*}+{\lambda }_3{E}^{*}& -\epsilon -b& {\lambda }_1{S}^{*}& {\lambda }_2{S}^{*}\\ 0& \epsilon & -\left({\gamma }_1+b+\mu \right)& 0\\ 0& 0& \mu & -\left({\gamma }_2+b\right)\end{array}\right].$$

(23)

The characteristic polynomial of Equation (13) is

$$P^{*}\left(\lambda \right)=a_4{}^{*}{\lambda }^4+a_3{}^{*}{\lambda }^3+a_2{}^{*}{\lambda }^2+a_1{}^{*}\lambda +a_0{}^{*}.$$

(24)

where

$$a_4{}^{*}=1>0,$$

(25)

$$a_3{}^{*}=b{R}_0+3b+{\gamma }_1+{\gamma }_2+\epsilon +\mu >0,$$

(26)

$$a_1{}^{*}=b{R}_0\left({\delta }_1+{\delta }_2\right)+\frac{{\lambda }_3\left({\gamma }_1+b+\mu \right)\left({\gamma }_2+b\right)}{R_0}+\frac{b\left(R_0-1\right)\left[{\lambda }_3\left({\gamma }_1+{\gamma }_2+2b+\mu \right)+{\lambda }_1\epsilon \right]}{R_0}>0,$$

(27)

and

$$a_0{}^{*}=b\left(R_0-1\right)\left({\gamma }_1+b+\mu \right)\left({\gamma }_2+b\right)\left(\epsilon +b\right)+{\lambda }_{3}b\left({\gamma }_1+b+\mu \right)\left({\gamma }_2+b\right)>0.$$

(28)

where

$${\delta }_1=\left({\gamma }_2+b\right)\left(\epsilon +b\right)+\left({\gamma }_1+b+\mu \right)\left({\gamma }_2+b\right),$$

(29)

$${\delta }_2=\frac{{\lambda }_2\epsilon \mu +{\lambda }_3\left({\gamma }_1+b+\mu \right)\left({\gamma }_2+b\right)}{R_0\left({\gamma }_2+b\right)}.$$

(30)

In addition, a simple calculation shows $a_1{}^{*}{a}_2{}^{*}{a}_3{}^{*}-a_1{{}^{*}}^2-a_0{}^{*}{a}_3{{}^{*}}^2>0$ and $a_2{}^{*}{a}_3{}^{*}-a_1{}^{*}>0$. Thus, if R₀ > 1, according to the Routh criterion, P* is locally asymptotically stable.

Theorem 4.

When R₀ > 1, the epidemic equilibrium point P* is globally asymptotically stable.

Proof. 

In this subsection, S, E, I₁, I₂, R and V denote S(t), E(t), I₁(t), I₂(t), R(t) and V(t), respectively. Consider the following Lyapunov function:

$$\begin{array}{l}V=S-S^{*}-S^{*}\ln \frac{S}{S^{*}}+E-E^{*}-E^{*}\ln \frac{E}{E^{*}}+\frac{\epsilon +b-{\lambda }_3{S}^{*}}{\epsilon }\left(I_1-I_1{}^{*}-I_1{}^{*}\ln \frac{I_1}{I_1{}^{*}}\right)\\ \begin{array}{c}\end{array}\begin{array}{c}\end{array}\begin{array}{c}\end{array}+\frac{{\lambda }_2{S}^{*}}{{\gamma }_2+b}\left(I_2-I_2{}^{*}-I_2{}^{*}\ln \frac{I_2}{I_2{}^{*}}\right)\end{array}$$

(31)

For simplicity, let $t_0=\frac{S}{S^{*}}$, $t_1=\frac{E}{E^{*}}$, $t_2=\frac{I_1}{I_1^{*}}$ and $t_3=\frac{I_2}{I_2^{*}}$, we obtain:

$$\begin{array}{l}\frac{dV}{dt}=b+b{S}^{*}+(\epsilon +b)E^{*}+\frac{\epsilon +b-{\lambda }_3{S}^{*}}{\epsilon }\left({\gamma }_1+b+\mu \right)I_1^{*}+{\lambda }_2{S}^{*}{I}_2^{*}\\ \begin{array}{c}\end{array}\begin{array}{c}\end{array}\begin{array}{c}\end{array}\begin{array}{c}\end{array}-(b{S}^{*}+{\lambda }_3{S}^{*}{E}^{\mathrm{*}})t_0-b\frac{1}{t_0}-(\epsilon +b-{\lambda }_3{S}^{*})E^{*}\frac{t_1}{t_2}-\frac{{\lambda }_2\mu S^{*}{I}_1^{*}}{{\gamma }_2+b}\frac{t_2}{t_3}\\ \begin{array}{c}\end{array}\begin{array}{c}\end{array}\begin{array}{c}\end{array}\begin{array}{c}\end{array}-{\lambda }_1{S}^{*}{I}_1{}^{*}\frac{t_0{t}_2}{t_1}-{\lambda }_2{S}^{*}{I}_2{}^{*}\frac{t_0{t}_3}{t_1}\end{array}$$

(32)

By substituting Equations (4)–(7) in (32), we can get

$$\begin{array}{l}\frac{dV}{dt}=2\left(b{S}^{*}+{\lambda }_3{S}^{*}{E}^{*}\right)+3{\lambda }_1{S}^{*}{I}_1{}^{*}+4{\lambda }_2{S}^{*}{I}_2{}^{*}-(b{S}^{*}+{\lambda }_3{S}^{*}{E}^{*})t_0\\ \begin{array}{c}\end{array}\begin{array}{c}\end{array}\begin{array}{c}\end{array}\begin{array}{c}\end{array}-\left(b{S}^{*}+{\lambda }_1{S}^{*}{I}_1{}^{*}+{\lambda }_2{S}^{*}{I}_2{}^{*}+{\lambda }_3{S}^{*}{E}^{*}\right)\frac{1}{t_0}-\left({\lambda }_1{S}^{*}{I}_1{}^{*}+{\lambda }_2{S}^{*}{I}_2{}^{*}\right)\frac{t_1}{t_2}\\ \begin{array}{c}\end{array}\begin{array}{c}\end{array}\begin{array}{c}\end{array}\begin{array}{c}\end{array}-{\lambda }_2{S}^{*}{I}_2{}^{*}\frac{t_2}{t_3}-{\lambda }_1{S}^{*}{I}_1{}^{*}\frac{t_0{t}_2}{t_1}-{\lambda }_2{S}^{*}{I}_2{}^{*}\frac{t_0{t}_3}{t_1}\\ \begin{array}{c}\end{array}\begin{array}{c}\end{array}\begin{array}{c}\end{array}=-\left(b{S}^{*}+{\lambda }_3{S}^{*}{E}^{*}\right)\left(t_0+\frac{1}{t_0}-2\right)-{\lambda }_1{S}^{*}{I}_1{}^{*}\left(\frac{1}{t_0}+\frac{t_0{t}_2}{t_1}+\frac{t_1}{t_2}-3\right)\\ \begin{array}{c}\end{array}\begin{array}{c}\end{array}\begin{array}{c}\end{array}\begin{array}{c}\end{array}-{\lambda }_2{S}^{*}{I}_2{}^{*}\left(\frac{1}{t_0}+\frac{t_1}{t_2}+\frac{t_2}{t_3}+\frac{t_0{t}_3}{t_1}-4\right)\end{array}$$

(33)

According to mean inequality $\sqrt[n]{a_1{a}_2\cdots a_n}\le \frac{a_1+a_2+\cdots +a_n}{n}$, $t_0+\frac{1}{t_0}\ge 2$, $\frac{1}{t_0}+\frac{t_0{t}_2}{t_1}+\frac{t_1}{t_2}\ge 3$ and $\frac{1}{t_0}+\frac{t_1}{t_2}+\frac{t_2}{t_3}+\frac{t_0{t}_3}{t_1}\ge 4$ are obtained. Then $\frac{dV}{dt}\le 0$ if $R_0\ge 1$. According to the Lyapunov stability theorem, the epidemic equilibrium is globally asymptotically stable. □

3.3. 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₁ nodes γ₁ and the repair rate of I₂ nodes γ₂ as the control variable, the cost function is as follows:

$$J\left({\gamma }_1,{\gamma }_2\right)=\underset{{\gamma }_1,{\gamma }_2}{\min }\left\{E\left(t_f\right)+I_1\left(t_f\right)+I_2\left(t_f\right)+{\displaystyle \int {}_0^{t_f}\left(c_1{\gamma }_1{}^2\left(t\right)I_1{}^2\left(t\right)+c_2{\gamma }_2{}^2\left(t\right)I_2{}^2\left(t\right)\right)}\begin{array}{c}\end{array}dt\right\},$$

(34)

In this formula, c₁ and c₂ indicate the repair cost parameter of I₁ nodes and I₂ nodes, respectively. Here, t_f represents the terminal moment. The terms $c_1{\gamma }_1{}^2\left(t\right)I_1{}^2\left(t\right)$ and $c_2{\gamma }_2{}^2\left(t\right)I_2{}^2\left(t\right)$ describe the repair cost of I₁ nodes and I₂ nodes at time t. The Equation (34) shows that the sum of the number of E nodes, I₁ nodes and I₂ nodes at the terminal moment and the cost from the start time to the terminal moment is the least by controlling γ₁(t) and γ₂(t). We apply the Pontryagin maximum principle to solve the optimization problem and consider the following Hamiltonian:

$$H=c_1{\gamma }_1{}^2\left(t\right)I_1{}^2\left(t\right)+c_2{\gamma }_2{}^2\left(t\right)I_2{}^2\left(t\right)+{\beta }_1\left(t\right)\stackrel{.}{S}\left(t\right)+{\beta }_2\left(t\right)\stackrel{.}{E}\left(t\right)+{\beta }_3\left(t\right){\stackrel{.}{I}}_1\left(t\right)+{\beta }_4\left(t\right){\stackrel{.}{I}}_2\left(t\right),$$

(35)

where $\beta \left(t\right)=\left\{{\beta }_1\left(t\right),{\beta }_2\left(t\right),{\beta }_3\left(t\right),{\beta }_4\left(t\right)\right\}$ is a set of covariates. According to the Pontryagin Maximum Principle, the differential equations of covariates can be obtained as follow:

$$\{\begin{array}{l}\stackrel{.}{{\beta }_1}\left(t\right)=\left[{\beta }_1\left(t\right)-{\beta }_2\left(t\right)\right]\left[{\lambda }_1{I}_1\left(t\right)+{\lambda }_2{I}_2\left(t\right)+{\lambda }_{3}E\left(t\right)\right]+b{\beta }_1\left(t\right)\\ \stackrel{.}{{\beta }_2}\left(t\right)=\left[{\beta }_1\left(t\right)-{\beta }_2\left(t\right)\right]{\lambda }_{3}S\left(t\right)+b{\beta }_2\left(t\right)\\ \stackrel{.}{{\beta }_3}\left(t\right)=-2c_1{I}_1\left(t\right){\gamma }_1{}^2\left(t\right)+\left[{\beta }_1\left(t\right)-{\beta }_2\left(t\right)\right]{\lambda }_{1}S\left(t\right)+\left({\gamma }_1+b+\mu \right){\beta }_3\left(t\right)-\mu {\beta }_4\left(t\right)\\ \stackrel{.}{{\beta }_4}\left(t\right)=-2c_2{I}_2\left(t\right){\gamma }_2{}^2\left(t\right)+\left[{\beta }_1\left(t\right)-{\beta }_2\left(t\right)\right]{\lambda }_{2}S\left(t\right)+\left({\gamma }_2+b\right){\beta }_4\left(t\right)\end{array},$$

(36)

where ${\beta }_1\left(t_f\right)=0$ and ${\beta }_2\left(t_f\right)={\beta }_3\left(t_f\right)={\beta }_4\left(t_f\right)=1$.

According to the Pontryagin Maximum Principle, we get

$$\{\begin{array}{l}\frac{\partial H}{\partial {\gamma }_1\left(t\right)}=2c_1{\gamma }_1\left(t\right)I_1{}^2\left(t\right)-{\beta }_3\left(t\right)I_1\left(t\right)=0,\\ \frac{\partial H}{\partial {\gamma }_2\left(t\right)}=2c_2{\gamma }_2\left(t\right)I_2{}^2\left(t\right)-{\beta }_4\left(t\right)I_2\left(t\right)=0.\end{array}$$

(37)

From the Equation (38), it can be concluded that:

$${\gamma }_1\left(t\right)=\frac{{\beta }_3\left(t\right)}{2c_1{I}_1\left(t\right)},{\gamma }_2\left(t\right)=\frac{{\beta }_4\left(t\right)}{2c_2{I}_2\left(t\right)}.$$

(38)

And ${\gamma }_1\left(t\right),{\gamma }_2\left(t\right)\in \left[0,1\right]$, so the optimal control pairs are as follows:

$${\gamma }_1*\left(t\right)=\min \left\{\max \left(0,\frac{{\beta }_3\left(t\right)}{2c_1{I}_1\left(t\right)}\right),\begin{array}{c}\end{array}1\right\},{\gamma }_2*\left(t\right)=\min \left\{\max \left(0,\frac{{\beta }_4\left(t\right)}{2c_2{I}_2\left(t\right)}\right),\begin{array}{c}\end{array}1\right\}.$$

(39)

4. 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.

4.1. Stable Analysis of Disease-Free Equilibrium

In this subsection, we suppose ${\lambda }_1=0.1$, ${\lambda }_2=0.2$, ${\lambda }_3=0.05$, ${\gamma }_1=0.1$, ${\gamma }_2=0.05$, $\mu =0.05$, $\epsilon =0.05$, and $b=0.1$. According to the Formula (11), we can get $R_0=0.56<1$. Theoretically, $S(\infty )=1$, $E(\infty )=0$, $I_1(\infty )=0$, and $I_2(\infty )=0$ according to the Formula (3). The simulation results are illustrated in Figure 2, and it is consistent with Theorems 1 and 2.

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₁ 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₁ and I₂ nodes in the beginning if the curves start from the X-axis; there are only E and I₂ nodes in the beginning if the curves start from the Y-axis; and there are only S and I₂ 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₁ nodes and E nodes is 1 in the X-Y plane at the initial moment; the ratio sum of I₁ 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₂ 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₂ and I₁ nodes when the curves start from the X-axis; there contain only E and I₁ nodes when the curves start from the Y-axis; and there contain only S and I₁ nodes when the curves start from the Z-axis. In Figure 2d, the curves contain only I₂ and E nodes in the X-Y plane; the curves contain only I₂ 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₀ < 1, as shown in Figure 2. These results are according to Theorems 1 and 2.

4.2. Stable Analysis of Epidemic Equilibrium

In this subsection, ${\lambda }_1=0.2$, ${\lambda }_2=0.3$, ${\lambda }_3=0.05$, ${\gamma }_1=0.1$, ${\gamma }_2=0.05$, $\mu =0.05$, $\epsilon =0.1$, and $b=0.05$. It is easily obtained that R₀ = 1.5 > 1. According to the Formulas (4) to (6), we know $S(\infty )=0.667$, $E(\infty )=0.111$, $I_1(\infty )=0.056$, and $I_2(\infty )=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, I₁ 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₂ 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.

4.3. 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₀. Suppose that except for b = 0.1, the other parameters remain unchanged, as described in Section 4.2. Besides, set $\gamma \left(t\right)=\left\{{\gamma }_1\left(t\right),{\gamma }_2\left(t\right)\right\}$,$X\left(t\right)=\left\{S\left(t\right),E\left(t\right),I_1\left(t\right),I_2\left(t\right),R\left(t\right)\right\}$, $\beta \left(t\right)=\left\{{\beta }_1\left(t\right),{\beta }_2\left(t\right),{\beta }_3\left(t\right),{\beta }_4\left(t\right)\right\}$, $t_f=200$, $c_1=1$, $c_2=2$, $S\left(0\right)=0.4$, $E\left(0\right)=0.3$, $I_1\left(0\right)=0.2$, $I_2\left(0\right)=0.1$, and $R\left(0\right)=0$. The optimal control calculation process is shown in Figure 4.

Firstly, the control variables γ₁, γ₂ 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.

4.3.1. 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.

4.3.2. 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={\displaystyle {\int }_0^{{\mathrm{t}}_{\mathrm{f}}}\left[c_1{\gamma }_1\left(t\right)+c_2{\gamma }_2\left(t\right)\right]}\begin{array}{c}\end{array}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.

When γ₁ and γ₂ 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 I₁ nodes is higher and the repair cost c₂ of I₂ nodes is larger than the repair cost c₁ of I₁ nodes, so γ₂ < γ₁. Then, due to the decreasing ratio of I₁ nodes and the higher transmission rate λ₂ of I₂ nodes, I₂ nodes become the biggest threat, so γ₂ > γ₁. 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. Figure 8 provides the optimal combination of I₁ nodes repair rate γ₁ and I₂ nodes repair rate γ₂ at all times.

4.3.3. Influence of Control Variables on the Basic Reproduction Number R₀

The influence of control variables on R₀ is shown in Figure 9. The simulation results show that R₀ decreases with the increase of the repair rate γ₁ and γ₂. The lower the R₀, the faster the virus is eliminated. When the control variables converge to 1, then R₀ < 1, as shown in Figure 8. Thus, the system will eventually achieve the disease-free equilibrium under optimal control. In addition, compared to γ₂, γ₁ has a greater impact on R₀ under these parameters as described in Section 4.3.

5. 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 γ₁ and γ₂. 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 γ₁ and γ₂, 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.