Analysis and Control of Malware Mutation Model in Wireless Rechargeable Sensor Network with Charging Delay

: In wireless rechargeable sensors (WRSNs), the system is vulnerable to be attacked by malware. Because of the distributed network structure of WRSNs, the malware attack has great inﬂuence on the security system of WRSNs. With the variability in malware, the problem of decryption and coding errors will lead to the malware mutating. In this paper, there are two problems to be solved, including the malware mutation and the charging delay in WRSNs. The malware mutation state and the low-energy state are introduced. Firstly, three different equilibrium solutions of the mutation model are given. Then, the local stability is proven by the characteristic equation, and the system will be stabilized at different equilibrium solutions when the base reproductive number is different. With the condition of charging delay, the bifurcation phenomenon is investigated by using the Hopf bifurcation theory. Furthermore, to improve the security of WRSNs and decrease the control cost, the Pontryagin’s Maximum principle is applied to obtain an optimal control scheme under mutation and charging delay. Finally, the numerical simulation is applied by Matlab to conﬁrm this model. The simulation results show that the mutation malware can be controlled when the delay is less than the maximum threshold.


Introduction
With the development of various network applications, there are many network security issues for networks [1][2][3][4]. Thus, the security of wireless sensor network (WSN) is a hot topic [5][6][7]. Many authors have studied WSNs from different aspects. A WSN is a distributed network structure. It has a great number of sensor nodes connected by a common wireless channel. Thus, it is applied in many fields, such as the military, environmental monitoring, medical treatment, disaster management, automotive applications, and so on [8,9]. However, WSN is constrained by limited resource, processing speed, and battery life. Because of the deployment of the nodes, it is hard to create a complicated safeguard ability for a WSN. It is also easy to lead to information transmission failure or even network paralysis when the WSN is attacked by malware. Therefore, it is urgent and significant to take measures to ensure accuracy and efficiency in WSN.
More and more scholars have studied the propagation of malware. At first, Kephart et al. [10,11] studied the propagation of the virus by incorporating concepts of epidemic dynamics. Since then, this method has opened a new way to solve the security problems of networks. Liu et al. [12] studied the malware propagation on complex networks which follow the power-law degree distribution. A SEIR model is proposed by Zhu et al. [13], and the authors analyzed the threshold, time delay, and global stability of a rumor spreading on social networks. Then, authors [14,15] also studied the spread of rumors in complex networks. Huang et al. [16] proposed spreading dynamics of social behaviors in multiplex networks. With the mathematical model, a number of researchers have made attempts using different models for investigating the malware disseminating dynamics in WSNs. When considering the spatial and temproal parameters, Khayam et al. [17] proposed a worm propagation model by using signal processing techniques in WSN. A worm propagation model in a WSN is proposed in Refs. [18,19], and the authors took into account the communication radius and the distributed density of the nodes. In [20], Ojha et al. also proposed an improved model which aggregates quarantine and vaccination techniques by Feng et al. Muthukrishnan et al. [21] established the node-based malware trace and patching model. They studied the topological structure and proposed an optimal control strategy. In Ref. [22], the authors proposed the SIR1R2 with secondary immunity to analyze the feature of malware propagation in WSNs. In Ref. [18], a mathematical equation is presented as follows: where, S, I, and R are the state nodes. They mean the numbers of the susceptible, the infected, and the recovered nodes, respectively. As shown in Table 1, the parameters are explained. Table 1. Notation.

N
The number of all nodes in the WSN µ The proportion of birth and death sensor nodes ξ The proportion of susceptible nodes becoming infected nodes ω The proportion of susceptible nodes becoming recovered nodes γ The proportion of infected nodes becoming recovered nodes The proportion of recovered nodes becoming susceptible nodes To study the malware spread specifically, other factors are taken into account in the mawlare spreading model, including time delay, virus mutation, the battery life, etc. Because of the battery capacity limitations, many scholars also pay attention to rechargeable WSN schemes [23][24][25] to solve the battery life problem. Li et al. [26] proposed an energyefficient patching strategy to solve the problems of virus attacks and energy constraints. Liu et al. [27,28] considered the impact of pulse charging and charging delay on the system. Keshri et al. [29] proposed an SEIR model with two time delays (including the period exposed and the period of temporary immunity). They studied the effect of delay in the WSN during multiple attacks. In Ref. [30], the authors proposed the propagation of a worm model with three delays and analyzed the influence of multiple delays on a system. For these time delay models [31][32][33], the authors investigated the condition of Hopf bifurcation. However, the energy constraint is the most important problem to be solved for all delay problems. It is necessary to study the low-energy state and the influence of charging delay on a WSN's efficiency. Liu et al. [34] proposed the malware propagation model when considering the low-energy state. They revealed the relationship between the charging rate and the basic reproductive number. However, they did not analyze the effect of a charging delay in WRSNs.
At the moment, it is worth noting that virus mutation [35] means that a virus is likely to convert to another virus when transmitting the virus on networks. In most networks, the malware mutation is real [36]. The condition of mutation needs to be considered, but few people have studied dynamic models of malware mutation. Li et al. [37,38] proposed an epidemic model with virus mutation and considered the virus mutation during transmission. Afterwards, they considered the virus's spontaneous variation, and they only established a delay model with virus mutation; they did not analyze the effect of delay on the mutation model. Hao [39] proposed the virus variation model on complex networks. They mainly analyzed the stability. In Ref. [40], the authors proposed a SIVR model to investigate the dynamics of epidemic propogation with virus mutation on complex networks, and they only discussed different epidemic immunization strategies to control the propagation of viruses and their mutations. Regardless of mutation, there are some control methods, including dissemination of security patches or recovery in Ref. [41]. Zhang et al. [42] considered the control strategy of the SIQRS model and provided three measures of vaccination, quarantine, and treatment. Zhu et al. [43] established an optimal control model with delay, and four different control strategies were investigated. Many excellent works have also been described in Refs. [44][45][46]. Chen et al. [47] considered the virus mutation with drug resistance and proposed an SIVRS epidemic spreading model on scale-free networks. They analyzed the mutation under heterogeneous networks. The model has only one reproductive number, and the model could not reflect the influence of the variation and infection node well. Liu et al. [48,49] analyzed the propagation of malware with virus mutation in wireless sensor networks. They merely studied the stability and control. Based on the above survey, there are many problems with mutation models. However, few people have studied charging delays and mutation models simultaneously. Moreover, the problem of optimal control needs to be solved.
In this paper, a new model is proposed to consider charging delay and mutation. A complex SIIRL model is proposed in WSN. S, I 1 , I 2 , R, and L mean the numbers of the susceptible, the infected, the mutation malware, the recovered, and the low-energy nodes, respectively. Our contributions are as follows: 1. For existing mutation models, the WSN has an energy constraint. As a solution, charging needs to be completed within a certain period of time. To address this problem, the mutation model with charging delay is proposed to analyse the process of spreading malware. The effect of charging delay on mutation model is also analyzed. 2. Through the analysis of the mutation model with charging delay, two different basic reproductive numbers, R 1 and R 2 , are obtained. Then, the locally asymptotically stability is proven. Moreover, the bifurcation phenomenon and periodic oscillation of the mutation model are analyzed. 3. Based on the mutation model with charging delay, the optimal charging delay control strategy is discussed. According to the Pontryagin's Maximun Principle, the optimal variable u * (t) can be derived. It is found that the charging delay has great influence on optimal control. Meanwhile, no treatment control, average control, and full treatment are also analyzed. By comparison with the three control methods, the optimal control is the best control method of them.
The follow of this paper is organized as follows. In Section 2, the SIIRL model is proposed. The mean-field equations are analyzed. In Section 3, the equilibrium points are calculated, and the local stability and the bifurcation phenomenon are investigated. In Section 4, the direction of Hopf bifurcation is analyzed. In Section 5, the charging delay control strategy is proposed for mutation malware spreading in WRSNs. In Section 6, some simulations are described to support the theoretical analysis. Finally, in Section 7, the conclusion of the model is described.

Modeling
In this section, the model of malware spreading is introduced. In general, the condition of normal and abnormal depends on if the sensors are infected or not. The normal working sensors are not infected by malware such as S and R. For example, the energy consumption of S is normal. R will be not be attacked. Thus, the energy consumption of R is normal.
The abnormal working sensors are infected by malware such as I 1 and I 2 . The nodes will consume a lot of energy when the sensors run the malware. The assumption is proposed that the conversion ratio of S or R to L is the same c 1 and the conversion ratio of I 1 or I 2 to L is the same c 2 . As depicted above, the model of equation is given as follows: According to Equation (2), A is the quantity of the injected nodes; β 1 is the conversion ratio of susceptible to infected; β 2 is the conversion ratio of the susceptible to virus mutation; ε is the conversion ratio of the infected to virus mutation; c 1 is the conversion ratio of the susceptible or recovered to low energy; c 2 is the conversion ratio of the infected or virus mutation to low energy; d is the conversion ratio of the sensor death; γ 1 is the conversion ratio of the infected to the recovered; γ 2 is the conversion ratio of the virus mutation to the recovered; µ is the conversion ratio of recovery from low energy; α is the failure ratio of recovery; and τ is the time delay.

Local Stability and Existence of Hopf Bifurcation
According to the system (2), the solutions of the system, including the disease-free equilibrium, the boundary equilibrium, and the endemic equilibrium, are given by: The disease-free equilibrium is E 0 (S 0 , 0, 0, R 0 , L 0 ). According to the next generation, the basic reproductive numbers can be calculated. The matrix is given by: Thus, the basic reproductive numbers are given by: If R 1 < 1 and R 2 < 1, the system has a disease-free equilibrium solution; if R 2 > 1 and R 2 > R 1 , the system has a disease-free equilibrium solution and a boundary equilibrium solution; if R 1 > 1 and R 1 > R 2 , the system has a disease-free equilibrium solution and an endemic equilibrium solution. The equilibrium is given as: where, The characteristic matrix is given as: where The characteristic equation is given by: where A 5 = 1, and the others are given by: Proof. If τ = 0, Equation (24) can be depicted: Thus, through putting E 0 into Equation (22), we can obtain b 1 = b 2 = b 3 = 0. Equation (31) can be given as follows: where Assumption 1. For the system (2), when τ = 0, if the conditions (R 1 < 1, is locally asymptotically stable in accordance with the Routh-Hurwitz criterion. Theorem 2. If R 1 < 1 and R 2 > 1,Ē(S, 0,Ī,R,L) is locally asymptotically stable when τ = 0 .
The characteristic equation can be obtained: If the root of Equation (37) is λ = iω, Equation (37) can be depicted as follows: Hence, the imaginary component and real component can be given as: Thus, ( 1 2 + 2 2 ) can be obtained as follows: (40) is given by: Obviously, if H 01 < 0, Equation (42) has one positive root at least. Thus, the assumption can be obtained as follows.

Assumption 2.
We assume Equation (42) has five solutions and each solution has a different τ.
The results can be given by the following theorem.

Direction and Stability of the Hopf Bifurcation
In this section, the Hopf bifurcation's properties of the system (2) will be analyzed. Let τ = τ 0 + µ, µ ∈ R (In this part, µ is different from above). µ = 0 is value of the bifurcation of the system (2). It can be obtained that µ 1 (t) = S(t) − S * , µ 2 (t) = I 1 (t) − I * 1 , µ 3 (t) = I 2 (t) − I * 2 , µ 4 (t) = R(t) − R * , µ 5 (t) = L(t) − L * . The system is written as: where With The adjoint operator A * (ϕ) of A(0) is defined as: The formula can be defined as: ) T e iω 0 τs is the eigenvector of A * (0). The eigenvector of A * (0) and A(0) are −iω 0 τ 0 and iω 0 τ 0 , and the matrix operations are given by: The results can be obtained as follows: and where h * (s), h(θ) = 1. From Equation (68), the equation can be given by Then, we can choosē According to the center manifold theory, the coordinate is calculated to express the manifold C 0 at µ = 0. When µ = 0, define: With C 0 , it is given by: Then, z andz are local coordinates for C 0 in the direction of h * andh * , and the real solution is considered. Thus: that is: With g(z,z) =h * (0) f 0 (z,z) = g 20 z 2 2 + g 11 zz + g 02z According to Equation (84), Equation (86), and f 0 (z,z) = F(0, µ(t)), it is given by The coefficients of the Hopf bifurcation's properties are given by: However, there are W 11 (θ) and W 20 (θ) in g 21 . The formula can be obtained as follows: and From Equation (89), it is given by: Thus, substituting (90), H 20 (0) and H 11 (0), we can obtain the matrix: Consequently, the values can be calculated: Theorem 6. For the system (2), µ 2 , ψ 2 , and T 2 explain the bifurcation properties. µ 2 determines the direction of the Hopf bifurcation (subcritical if µ 2 < 0, supercritical if µ 2 > 0); ψ 2 determines the stability of the periodic solutions of bifurcation (stable if ψ 2 < 0, unstable if ψ 2 > 0); and T 2 determines the period of the periodic solutions of bifurcation (decrease if T 2 < 0, increase if T 2 > 0).

Optimal Strategy
In this section, the strategy of optimal control including installing immune patches or removing viruses is proposed to control the malware.
The parameters γ 1 = u 1 (t) and γ 2 = u 2 (t) are introduced. Thus, the optimal control system is given by: is a control variable, and it satisfies the Lebesgue measurable.
Thus, the objective function is given by: m and n are positive weighting factors. The Lagrangian function is given by: The Hamiltonian function can be given as: The indicative function is introduced as follows: According to the Pontryagin's maximum principle, the necessary conditions of the optimal control are the state equationẊ(t) = H λ (X, X τ , u, λ)(t), the optimal condition 0 = H u (X, X τ , u, λ)(t), and the adjoint equation − dλ dt = H x (X, X τ , u, λ)(t) + λ(t + τ) H X τ (X, X τ , u, λ)(t).

Simulation
In this section, the theoretical results are illustrated and completed by simulation, including three parts. The first part explains the relationship between the parameters and the reproductive numbers. The second part illustrates the trend in the quantity of state nodes and the bifurcation properties of the system. The last part describes the result of the optimal control.

Relationship between Parameters and Basic Reproductive Numbers
The relationship between basic reproductive numbers and other parameters are shown in Figure 1. The initial parameters are set as: A = 10, β 1 = 0.04, γ 1 = 0.2, c 1 = 0.3, c 2 = 0.2, d = 0.1, and ε = 0.4. As shown in Figure 1a, it can be observed that with the increase in d, the basic reproductive number R 1 decreases, and with the increase in β 1 , R 1 increases. When c 1 = 0, γ 1 = 0, or ε = 0, with the increase in β 1 , R 1 increases. Therefore, there is a linear relationship among β 1 , d, and R 1 . As shown in Figure 1b, the initial parameters are set as: A = 10, β 2 = 0.03, γ 2 = 0.4, c 2 = 0.2, and d = 0.1. Let β 2 and d vary. With the increase in d, the basic reproductive number R 2 decreases, and with the increase in β 2 , R 2 increases. There is a linear relationship among β 2 , d, and R 2 .

The Trend of the Quantity of State Nodes
In Sections 2 and 3, the three solutions of the system (2) are calculated by analyzing and demonstrating. The following part will illustrate the theoretical results. The influence of different Hopf bifurcation points (τ) are illustrated. Through analyzing three conditions, including E 0 ,Ē, and E * , we can obtain the variation in the state nodes and the stable situation.

Optimal Control Strategy Analysis
In this part, the control factor is added to this system, and the optimal control method is proposed to consider the cost of treating the infected nodes. The theoretical analysis for the optimal control is demonstrated in Section 5. In this part, there are different methods of control, including maximum control, minimum control, average control, and optimal control [43].
The parameters are set as A = 20, β 1 = 0.389, β 2 = 0.1, c 1 = 0.15, c 2 = 0.22, µ = 0.56, = 0.2, α = 0.115, d = 0.2, and τ = 2.2. The maximum control means full treatment; that is, the parameters' values are set as γ 1 = 1 and γ 2 = 1. The minimum control means no treatment; that is, the parameters' values are set as γ 1 = 0 and γ 2 = 0. The optimal control means γ 1 = u 1 (t) and γ 2 = u 2 (t); that is, u 1 (t) and u 2 (t) are the control curves of the patch input. Then, the average control is the average of the optimal control curve. When m = 30 and n = 20, the control curves are shown in Figure 6. Obviously, u 1 (t) is always greater than u 2 (t). As shown in Figure 7, we can obtain the curves of the cost for different methods. It is obvious that the cost of optimal control is the minimum. When τ = 3.2, the optimal control curve is shown in Figure 8. When τ = 5.2, the optimal control curve is shown in Figure 9. Thus, the larger the charging delay, the more oscillation appears in the optimal control curve.  When τ = 2.2, Figure 10 shows the ratio of different nodes for four control cases. For four cases, the proportion of S is the smallest and the proportion of R is the largest than the other state nodes of all the nodes. It is obvious that the ratios of I 1 and I 2 are lower for the full treatment than those for the no treatment. The ratio of I 1 decreases to 15.32%, and I 2 decreases to 3.11%. With the optimal control, the ratio of I 1 decreases to 18.67% and I 2 decreases to 7.31%. Then, the average control can also control the ratios of I 1 and I 2 .
As shown in Figure 11, the variation in nodes is shown. In addition to no treatment, the recovered nodes curve is the highest, which means that the number of recovered nodes is more than that of the other two control methods. The numbers I 1 and I 2 are less than that of no treatment. From Figure 11a, at the end of time, there is a noticeable wobble in the end of the optimal control curve.  In summary, according to the above analysis, the full treatment has the best effect to control the numbers of I 1 and I 2 when the control strategy of the system does not consider the cost. However, optimal control has the best effect to control the numbers of I 1 and I 2 when the control strategy of the system considers the cost.

Conclusions
In this paper, to investigate the malware propagation with charging delay and virus mutation in WRSNs, a new model is constructed and analyzed which improves the SIR model. The model involves a mutation which is caused by the infected variation. There is a charging delay in sensors. Unlike the SIR model, there are three equilibrium points including disease free equilibrium, boundary equilibrium, and endemic equilibrium. The analysis reveals the local stability of equilibrium points when the charging delay τ equals to zero. If R 1 < 1 and R 2 < 1, the malware does not affect the WRSNs and the malware will reduce to zero. In this charging delay model, the locally asymptotic stability is proven when τ is less than τ 0 . There are three different Hopf bifurcations.The system goes through a Hopf bifurcation, and the system goes into oscillation when τ is larger than τ 0 .
Meanwhile, the intent to study the SIIRL model of malware propagation in the WSN is used for controlling the spreading of mutation malware. We propose control strategies including full treatment control, average control, no treatment control, and optimal control. The cost relationships among the four control methods are Optimal > Average > Full > No. In addition, the malware proportion relationships among the four control methods are Full > Optimal > Average > No. Thus, considering the malware proportion and control cost, the best control method is optimal control. At the same time, the optimal control of the mutation model is also affected by the charging delay. The bigger the charging delay, the more oscillation appears in the optimal control curve. Finally, all of the theoretical analyses are proven by numerical simulation.
This paper discusses a mutation model with a charging delay. Although the dynamic propagation model is analyzed and the control strategy is proposed, there still exists an interesting work on how to establish rate of infection with actual situation. This work needs to be further studied in the future. However, with the development of networks, the heterogeneous network is one of our future research directions. Meanwhile, the stochastic model [50] and fractional order model [51] are also the trends in our future research.