Next Article in Journal
Leveraging Elasticity to Uncover the Role of Rabinowitsch Suspension through a Wavelike Conduit: Consolidated Blood Suspension Application
Previous Article in Journal
Group Analysis of the Plane Steady Vortex Submodel of Ideal Gas with Varying Entropy
 
 
Font Type:
Arial Georgia Verdana
Font Size:
Aa Aa Aa
Line Spacing:
Column Width:
Background:
Article

Dynamical Behavior Analysis of a Time-Delay SIRS-L Model in Rechargeable Wireless Sensor Networks

School of Mechanical and Electric Engineering, Guangzhou University, Guangzhou 510006, China
*
Authors to whom correspondence should be addressed.
Mathematics 2021, 9(16), 2007; https://doi.org/10.3390/math9162007
Submission received: 4 June 2021 / Revised: 19 August 2021 / Accepted: 19 August 2021 / Published: 22 August 2021

Abstract

:
The traditional SIRS virus propagation model is used to analyze the malware propagation behavior of wireless rechargeable sensor networks (WRSNs) by adding a new concept: the low-energy status nodes. The SIRS-L model has been developed in this article. Furthermore, the influence of time delay during the charging behavior of the low-energy status nodes needs to be considered. Hopf bifurcation is studied by discussing the time delay that is chosen as the bifurcation parameter. Finally, the properties of the Hopf bifurcation are explored by applying the normal form theory and the center manifold theorem.

1. Introduction

With the development and application of communication technology and sensor technology, wireless sensor networks (WSNs) have come into being. WSNs are widely used globally, such as in the military industry, agricultural production, ecological monitoring, disaster warning, and infrastructure status monitoring [1]. Due to the self-organized network of WSNs, sensor nodes are vulnerable to various attacks in the process of information transmission [2,3,4,5,6,7,8,9,10]. Thus, it is of practical significance to prevent and deal with malicious attacks in WSNs.
The characteristics of WSNs can be listed as follows: adapting wireless communication, the topology being able to be changed dynamically, the security mechanism being imperfect, and the energy of nodes being limited, etc. Thus, the security of wireless sensor networks is facing great challenges. The attacks against sensor network nodes can be roughly divided into the following types: network worms, bots, and Trojans, among others [11,12,13,14]. Most malware attacks in WSNs mainly aim at paralyzing them or losing their power. In order to overcome the defect of short life cycle of WSNs, wireless rechargeable sensor networks (WRSNs) have been developed and have caught more attention. Consequently, a relatively special attack has appeared to against the nodes’ power, namely, the denial of charge (DOC) [15]. The prevention and control of malware spread in WRSNs have become a hot issue.
The propagation mechanism of malware is similar to that of viruses in biological groups. Many scholars have made great contributions to malware propagation dynamics. In the beginning, the spread of malware mainly focuses on the traditional Internet. For example, Zou et al. put forward the worm propagation models represented by SI (susceptible–infected), SIS (susceptible–infected–susceptible), and SIR (susceptible–infected–recovered) [16,17,18,19]. These studies lay a foundation for the later studies on the propagation behavior of malware in WSNs. However, it is worth noting that there are some differences between WSNs and the traditional Internet: (1) Because of their topological structure, the communication range of WSNs has physical limitations. (2) The limited energy characteristics lead to uneven communication capabilities. (3) The link-layer access conflict and avoidance mechanism of WSNs sharing wireless channel also introduces space-time correlation for the propagation dynamics of WSNs. (4) WSNs have a high degree of self-organization. Therefore, it is of great significance to study the malware propagation model in WSNs, which has been proposed in recent years. For example, Tanachaiwiwat began to study the propagation behavior of network worms in 2009 [20]. Khayam et al. used signal processing technology to study the propagation model of worms in 2006 [21]. Furthermore, the SEIR (susceptible–exposed–infected–recovered) model is adapted to study the spread of malware [22], while the SEIRS-V (susceptible–exposed–infected–recovered–susceptible and vaccinated) model explores the application of Hopf bifurcation during the malware spreading [23,24]. With the development, the research on WRSNs was pushed forward recently [25,26,27,28,29,30,31], considering the new concept of L (low-energy status) nodes. It is worth noting that the transition from the low-energy status nodes to the normal status nodes is called charging. At present, there are many charging methods for WRSNs, such as unmanned aerial vehicle (UAV) charging [32].
The related techniques and modeling basis of this paper mainly refer to the papers in Table 1:
In order to promote the study of virus dynamics in WRSNs, the modeling in this paper is mainly based on the SIRS (susceptible–infected–recovered–susceptible) infectious disease model. Combined with the charging delay in WRSNs, the bifurcation in the process of charging needs to be considered. Therefore, the general contributions of this paper are as follows:
  • Establishing the SIRS-L (susceptible–infected–recovered–susceptible–low-energy) model.
  • The equilibrium solutions of the SIRS-L model are obtained, and the basic reproductive number R0 is defined through the regeneration matrix [34].
  • Revealing of the stability of the SIRS-L model when the charging delay is ignored.
  • The variation of the solutions of the characteristic equation are discussed if the charging delay is considered through the theory in [35], and the occurrence conditions of Hopf bifurcation are figured out.
  • The properties of the Hopf bifurcation are explored by applying the normal form theory and the center manifold theorem [36].
The content of this paper is organized as follows: In Section 2, we establish the SIRS-L (susceptible–infected–recovered–susceptible–low-energy) model combined with the concept of low-energy nodes and the charging delay; in Section 3, the equilibrium points analysis and the corresponding Hopf bifurcation conditions are presented; the properties of the Hopf bifurcation are given in Section 4; the simulation analysis are revealed in Section 5; the conclusions are given in Section 6.

2. Modeling

On the basis of the classical SIRS model, we obtained the SIRS-L model as follows (1), with the parameters being explained in Table 2:
d S ( t ) d t = α 1 S ( t ) I ( t ) N ( t ) α 4 S ( t ) β 1 S ( t ) b S ( t ) + α 2 I ( t ) + α 5 R ( t ) + β 2 L S ( t τ )   d I ( t ) d t = α 1 S ( t ) I ( t ) N ( t ) α 2 I ( t ) b I ( t ) β 1 I ( t ) α 3 I ( t ) + β 2 L I ( t τ )   d R ( t ) d t = α 3 I ( t ) + α 4 S ( t ) b R ( t ) β 1 R ( t ) α 5 R ( t ) + β 2 L R ( t τ ) d L S ( t ) d t = β 1 S ( t ) b L S ( t ) β 2 L S ( t τ ) d L I ( t ) d t = β 1 I ( t ) b L I ( t ) β 2 L I ( t τ ) d L R ( t ) d t = β 1 R ( t ) b L R ( t ) β 2 L R ( t τ )
The following formula is obtained through the injection rate and the deactivation rate: d N ( t ) d t = b N ( t ) . Therefore, we can obtain N ( ) = b if the total number of nodes is stable. In order to make the analysis more clearer, the system (1) is simplified as follows (2) with the meaning x ( t ) = X ( t ) N ( ) , where x ( t ) = ( s ( t ) ,   i ( t ) ,   r ( t ) ,   l s ( t ) ,   l i ( t ) ,   l r ( t ) ) T and X ( t ) = ( S ( t ) ,   I ( t ) ,   R ( t ) ,   L S ( t ) ,   L I ( t ) ,   L R ( t ) ) T . Thus, x ( t ) represents the proportion of nodes in the system (1).
d s ( t ) d t = b α 1 s ( t ) i ( t ) α 4 s ( t ) β 1 s ( t ) b s ( t ) + α 2 i ( t ) + α 5 r ( t ) + β 2 l s ( t τ ) d i ( t ) d t = α 1 s ( t ) i ( t ) α 2 i ( t ) b i ( t ) β 1 i ( t ) α 3 i ( t ) + β 2 l i ( t τ ) d r ( t ) d t = α 3 i ( t ) + α 4 s ( t ) b r ( t ) β 1 r ( t ) α 5 r ( t ) + β 2 l r ( t τ ) d l s ( t ) d t = β 1 s ( t ) b l s ( t ) β 2 l s ( t τ ) d l i ( t ) d t = β 1 i ( t ) b l i ( t ) β 2 l i ( t τ ) d l r ( t ) d t = β 1 r ( t ) b l r ( t ) β 2 l r ( t τ )

3. Local Stability and Analysis of Hopf Bifurcation

The disease-free solution e 0 and endemic solution e + of the system (2) can be obtained as follows:
e 0 ( s 0 , i 0 , r 0 , l s 0 , l i 0 , l r 0 ) = ( b α 4 + β 1 + b β 2 ε α 4 α 5 b + β 1 + α 5 β 2 ε ,     0 ,     α 4 b ( α 4 + β 1 + b β 2 ε ) ( b + β 1 + α 5 β 2 ε ) α 4 α 5 ,    ε b α 4 + β 1 + b β 2 ε α 4 α 5 b + β 1 + α 5 β 2 ε ,     0 ,     ε α 4 b ( α 4 + β 1 + b β 2 ε ) ( b + β 1 + α 5 β 2 ε ) α 4 α 5 ) ,
e + ( s + , i + , r + , l s + , l i + , l r + ) = ( α 2 + b + β 1 + α 3 β 2 ε α 1           , ( b + β 1 + α 5 β 2 ε ) r + α 4 s + α 3           ,     ( α 2 α 1 s + ) α 4 s + + α 3 ( α 4 + β 1 + b β 2 ε ) s + α 3 b ( α 2 α 1 s + ) ( b + β 1 + α 5 β 2 ε ) + α 3 α 5   ,    ε ( α 2 + b + β 1 + α 3 β 2 ε ) α 1 ,     ε [ ( b + β 1 + α 5 β 2 ε ) r + α 4 s + ] α 3 ,     ε [ ( α 2 α 1 s + ) α 4 s + + α 3 ( α 4 + β 1 + b β 2 ε ) s + α 3 b ] ( α 2 α 1 s + ) ( b + β 1 + α 5 β 2 ε ) + α 3 α 5 ) ,
where ε = β 1 b + β 2 .
The basic reproductive number R0 is obtained through the next generation matrix method [34].
Set:
F = ( α 1 s ( t ) 0 0 0 ) ,
and
V = ( α 2 + b + β 1 + α 3 β 2 β 1 b + β 2 ) .
Thus,
R 0 = ρ ( F V 1 ) = α 1 s 0 ( b + β 2 ) ( α 2 + b + β 1 + α 3 ) ( b + β 2 ) β 1 β 2 .
It can be found that if R 0 > 1 , the model (2) has the unique endemic solution e + . By linearization technique, the characteristic matrix of system (2) about e + can be obtained as follows.
J ( e + ) = ( λ a 11 a 12 b 13 e λ τ b 14 e λ τ 0 0 a 21 λ a 22 0 0 b 25 e λ τ 0 a 31 a 32 λ a 33 b 33 e λ τ 0 0 b 36 e λ τ a 41 0 0 λ a 44 b 44 e λ τ 0 0 0 a 52 0 0 λ a 55 b 55 e λ τ 0 0 0 a 63 0 0 λ a 66 b 66 e λ τ ) .
where
a 11 = α 1 i + α 4 β 1 b , a 12 = α 1 s + + α 2 , a 13 = α 5 ,     b 14 = β 2 a 21 = α 1 i + ,     a 22 = α 1 s + α 2 b β 1 α 3 , b 25 = β 2 a 31 = α 4 ,     a 32 = α 3 ,     a 33 = b β 1 α 5 ,     b 36 = β 2 a 41 = β 1 ,     a 44 = b ,     b 44 = β 2 a 52 = β 1 ,     a 55 = b ,     b 55 = β 2 a 63 = β 1 ,     a 66 = b ,     b 66 = β 2 .
The characteristic equation of the model (2) at e + is shown as:
λ 6 + A 5 λ 5 + A 4 λ 4 + A 3 λ 3 + A 2 λ 2 + A 1 λ + A 0 + ( B 5 λ 5 + B 4 λ 4 + B 3 λ 3 + B 2 λ 2 + B 1 λ + B 0 ) e λ τ + ( C 4 λ 4 + C 3 λ 3 + C 2 λ 2 + C 1 λ + C 0 ) e 2 λ τ + ( D 3 λ 3 + D 2 λ 2 + D 1 λ + D 0 ) e 3 λ τ = 0 .
The detailed representations of X n ( X = A , B , C , D ; n = 1 , 2 , 3 ) are given in Appendix A.
Firstly, we study the local stability of system (2) at e + when τ = 0 .
If τ = 0 , Equation (8) can be rebuilt as:
λ 6 + E 5 λ 5 + E 4 λ 4 + E 3 λ 3 + E 2 λ 2 + E 1 λ + E 0 = 0 ,
where
E 5 = A 5 + B 5 ,             E 4 = A 4 + B 4 + C 4 ,             E 3 = A 3 + B 3 + C 3 + D 3 , E 2 = A 2 + B 2 + C 2 + D 2 ,             E 1 = A 1 + B 1 + C 1 + D 1 ,             E 0 = A 0 + B 0 + C 0 + D 0 .
The following formulas can be obtained:
{ 1   =   E 5 , 2   =   | E 5 1 E 3 E 4 | , 3   =   | E 5 1 0 E 3 E 4 E 5 E 1 E 2 E 3 | , 4   =   | E 5 1 0 0 E 3 E 4 E 5 1 E 1 E 2 E 3 E 4 0 E 0 E 1 E 2 | , 5   =   | E 5 1 0 0 0 E 3 E 4 E 5 1 0 E 1 E 2 E 3 E 4 E 5 0 E 0 E 1 E 2 E 3 0 0 0 E 0 E 1 | , 6   =   E 0 5 .
If hypothesis (H1) E 5 ,   E 4 ,   E 3 ,   E 2 ,   E 1 ,   E 0 > 0 ,   1 , 2 , 3 , 4 , 5 , 6 > 0 and τ = 0 , in Equation (10) holds, and all the characteristic solutions of Equation (9) are negative. Thus, e + is locally asymptotically stable under the theory of the Hurwitz criterion.
Then, we begin to consider the change of the characteristic solutions of Equation (9) when the charging time delay τ is introduced.
If τ > 0 , multiplying e λ τ on both sides of Equation (8), we can obtain the following equation:
B 5 λ 5 + B 4 λ 4 + B 3 λ 3 + B 2 λ 2 + B 1 λ + B 0
Referring to the processing method of the paper [35], let λ = (ω > 0) be the root of Equation (11). After separating the imaginary and real parts, we obtain
{ F 1 cos τ ω F 2 sin τ ω + F 3 = G 1 sin 2 τ ω + G 2 cos 2 τ ω F 4 sin τ ω F 5 cos τ ω + F 6 = G 1 cos 2 τ ω G 2 sin 2 τ ω ,
where
F 1 = ω 6 + A 4 ω 4 A 2 ω 2 + A 0 + C 4 ω 4 C 2 ω 2 + C 0 , F 2 = A 5 ω 5 A 3 ω 3 + A 1 ω + C 3 ω 3 C 1 ω , F 3 = B 4 ω 4 B 2 ω 2 + B 0 , F 4 = A 5 ω 5 A 3 ω 3 + A 1 ω C 3 ω 3 + C 1 ω , F 5 = ω 6 + A 4 ω 4 A 2 ω 2 + A 0 C 4 ω 4 + C 2 ω 2 C 0 , F 6 = B 5 ω 5 B 3 ω 3 + B 1 ω , G 1 = D 3 ω 3 D 1 ω , G 2 = D 2 ω 2 D 0 .
Squaring and adding the two Equations in Equation (12), we obtain
( F 1 cos τ ω F 2 sin τ ω + F 3 ) 2 + ( F 4 sin τ ω F 5 cos τ ω + F 6 ) 2 = G 1 2 + G 2 2
Because of sin τ ω = ± 1 cos 2 τ ω , the following two cases are discussed:
Case 1: If sin τ ω = 1 cos 2 τ ω , Equation (13) could be represented as follows:
( F 1 cos τ ω F 2 1 cos 2 τ ω + F 3 ) 2 + ( F 4 1 cos 2 τ ω F 5 cos τ ω + F 6 ) 2 = G 1 2 + G 2 2
Equation (14) can be calculated as follows:
H 4 cos 4 τ ω + H 3 cos 3 τ ω + H 2 cos 2 τ ω + H 1 cos τ ω + H 0 = 0 ,
where
H 4 = ( F 1 2 + F 5 2 F 2 2 F 4 2 ) 2 + 4 ( F 4 F 5 F 1 F 2 ) 2 , H 3 = ( F 1 F 3 + F 5 F 6 ) ( F 1 2 + F 5 2 F 2 2 F 4 2 ) + 8 ( F 4 F 5 F 1 F 2 ) ( F 4 F 6 F 2 F 3 ) , H 2 = 4 ( F 1 F 3 + F 5 F 6 ) 2 4 ( F 4 F 5 F 1 F 2 ) 2 + 4 ( F 4 F 6 F 2 F 3 ) 2           + 2 ( F 1 2 + F 5 2 F 2 2 F 4 2 ) ( F 2 2 + F 3 2 + F 4 2 + F 6 2 G 1 2 G 2 2 ) , H 1 = 4 ( F 2 2 + F 3 2 + F 4 2 + F 6 2 G 1 2 G 2 2 ) 2 ( F 1 F 3 + F 5 F 6 )                                                                     8 ( F 4 F 5 F 1 F 2 ) ( F 4 F 6 F 2 F 3 ) , H 0 = ( F 2 2 + F 3 2 + F 4 2 + F 6 2 G 1 2 G 2 2 ) 2 4 ( F 4 F 6 F 2 F 3 ) 2 .
Let cos τ ω = y , and the following equation can be obtained:
f ( y ) = y 4 + H 3 H 4 y 3 + H 2 H 4 y 2 + H 1 H 4 y + H 0 H 4 ,
and
f ˙ ( y ) = 4 y 3 + 3 H 3 H 4 y 2 + 2 H 2 H 4 y 1 + H 1 H 4 .
Set
y 3 + I 2 y 2 + I 1 y 1 + I 0 = 0 ,
where
I 2 = 3 H 3 4 H 4 ,             I 1 = H 2 2 H 4 ,             I 0 = H 1 4 H 4 .
The following solutions can be calculated:
y 1 = σ 1 I 1 3 I 2 2 9 σ 1 I 2 3 , y 2 = I 1 3 I 2 2 9 2 σ 1 I 2 3 σ 1 2 3 I 1 3 I 2 2 9 σ 1 + σ 1 2 i , y 3 = I 1 3 I 2 2 9 2 σ 1 I 2 3 σ 1 2 + 3 ( I 1 3 I 2 2 9 σ 1 + σ 1 ) 2 i ,
where
σ 1 = ( ( I 2 3 27 I 1 I 2 6 + I 0 2 ) 2 + ( I 1 3 I 2 2 9 ) 3 I 0 2 I 2 3 27 + I 1 I 2 6 ) 1 3 .
Then, the expression of cos τ ω can be obtained as f 1 ( ω ) = cos τ ω ; combined with Equation (13), we can obtain f 2 ( ω ) = sin τ ω . A function concerning ω can be obtained by
f 1 2 ( ω ) + f 2 2 ( ω ) = 1 .
Thus, (H2) Equation (20) has finite positive roots ω 1 i ( i = 1 , 2 , , n ). The corresponding time delay can be obtained as follows:
τ 1 i ( j ) = 1 ω 1 i a r c c o s f 1 ( ω 1 i ) + 2 j π ω 1 i ,   i = 1 , 2 , , n .     j = 0 , 1 , 2
Case 2: If sin τ ω = 1 cos 2 τ ω , Equation (13) could be represented as follows:
( F 1 cos τ ω F 2 1 cos 2 τ ω + F 3 ) 2 + ( F 4 1 cos 2 τ ω F 5 cos τ ω + F 6 ) 2 = G 1 2 + G 2 2 .
Similar to Case 1, a function concerning ω can be obtained by
f 3 2 ( ω ) + f 4 2 ( ω ) = 1 .
Thus, (H3) Equation (23) has finite positive roots ω 2 i ( i = 1 , 2 , , n ). The corresponding time delay can be obtained as follows:
τ 2 i ( j ) = 1 ω 2 i a r c c o s f 3 ( ω 2 i ) + 2 j π ω 2 i ,   i = 1 , 2 , , n .     j = 0 , 1 , 2
Let
τ 0 = m i n { τ 1 i ( 0 ) , τ 2 i ( 0 ) } , i = 1 , 2 , , n .
Then, if τ = τ 0 , Equation (13) has a pair of purely imaginary roots ± i ω 0 . Next, the conditions for bifurcation are obtained in the following analysis. From Equation (11), we obtain
[ d λ d τ ] 1 = J 1 ( λ ) + J 2 ( λ ) e λ τ + J 3 ( λ ) e λ τ + J 4 ( λ ) e 2 λ τ K 1 ( λ ) e λ τ + K 2 ( λ ) e λ τ + K 3 ( λ ) e 2 λ τ τ λ ,
with
J 1 ( λ ) = 5 B 5 λ 4 + 4 B 4 λ 3 + 3 B 3 λ 2 + 2 B 2 λ + B 1 , J 2 ( λ ) = 6 λ 5 + 5 A 5 λ 4 + 4 A 4 λ 3 + 3 A 3 λ 2 + 2 A 2 λ 1 + A 1 , J 3 ( λ ) = 4 C 4 λ 3 + 3 C 3 λ 2 + 2 C 2 λ + C 1 , J 4 ( λ ) = 3 D 3 λ 2 + 2 D 2 λ + D 1 , K 1 ( λ ) = λ 7 A 5 λ 6 A 4 λ 5 A 3 λ 4 A 2 λ 3 A 1 λ 2 A 0 λ , K 2 ( λ ) = C 4 λ 5 + C 3 λ 4 + C 2 λ 3 + C 1 λ 2 + C 0 λ , K 3 ( λ ) = 2 ( D 3 λ 4 + D 2 λ 3 + D 1 λ 2 + D 0 λ ) .
Thus,
[ d λ d τ ] λ = i ω 0 τ = τ 0 1 = M R + M I i N R + N I i ,
where
M R = 5 B 5 ω 0 4 3 B 3 ω 0 2 + B 1 + ( 5 A 5 ω 0 4 3 A 3 ω 0 2 + A 1 3 C 3 ω 0 2 + C 1 ) cos τ 0 ω 0 + ( 6 ω 0 5 + 4 A 4 ω 0 3 2 A 2 ω 0 4 C 4 ω 0 3 + 2 C 2 ω 0 ) sin τ 0 ω 0 + ( 3 D 3 ω 0 2 + D 1 ) cos 2 τ 0 ω 0 + 2 D 2 ω 0 sin 2 τ 0 ω 0 , M I = 4 B 4 ω 0 3 + 2 B 2 ω 0 + ( 6 ω 0 5 4 A 4 ω 0 3 + 2 A 2 ω 0 4 C 4 ω 0 3 + 2 C 2 ω 0 ) cos τ 0 ω 0 + ( 5 A 5 ω 0 4 3 A 3 ω 0 2 + A 1 + 3 C 3 ω 0 2 C 1 ) sin τ 0 ω 0 + 2 D 2 ω 0 cos 2 τ 0 ω 0 + ( 3 D 3 ω 0 2 D 1 ) sin 2 τ 0 ω 0 , N R = ( A 5 ω 0 6 A 3 ω 0 4 + A 1 ω 0 2 + 3 C 3 ω 0 4 C 1 ω 0 2 ) cos τ 0 ω 0 + ( ω 0 7 + A 4 ω 0 5 A 2 ω 0 3 + A 0 ω 0 + C 4 ω 0 5 C 2 ω 0 3 + C 0 ω 0 ) sin τ 0 ω 0 + 2 ( D 3 ω 0 4 D 1 ω 0 2 ) cos 2 τ 0 ω 0 + 2 ( D 2 ω 0 3 + D 0 ω 0 ) sin 2 τ 0 ω 0 , N I = ( ω 0 7 A 4 ω 0 5 + A 2 ω 0 3 A 0 ω 0 + C 4 ω 0 5 C 2 ω 0 3 + C 0 ω 0 ) cos τ 0 ω 0 + ( A 5 ω 0 6 A 3 ω 0 4 + A 1 ω 0 2 3 C 3 ω 0 4 + C 1 ω 0 2 ) sin τ 0 ω 0 + 2 ( D 2 ω 0 3 + D 0 ω 0 ) cos 2 τ 0 ω 0 + 2 ( D 3 ω 0 4 + D 1 ω 0 2 ) sin 2 τ 0 ω 0 .
If the hypothesis (H4) s i g n ( R e [ d λ d τ ] λ = i ω 0 1 ) = s i g n ( M R N R + M I N I ) 0 holds, the following Theorem 1 can be obtained by the Hopf bifurcation theorem [36].
Theorem 1.
If (H1)–(H4) hold, the endemic solution  e + of the model (2) is locally asymptotically stable for τ ϵ [ 0 , τ 0 ) , and the model (2) undergoes a Hopf bifurcation at the endemic solution e + if τ = τ 0 .

4. Properties of the Hopf bifurcation

The properties of the Hopf bifurcation are explored in this section by applying the normal form theory and the center manifold theorem [36]. Let τ = τ 0 + μ ,   μ ϵ R . The transformations are given as u 1 ( t ) = s ( t ) s + ,   u 2 ( t ) = i ( t ) i + ,   u 3 ( t ) = r ( t ) r + ,   u 4 ( t ) = l s ( t ) l s + ,   u 5 ( t ) = l i ( t ) l i + , u 6 ( t ) = l r ( t ) l r +   and   t ( t τ ) . Thus, the model (2) can be written as:
u ˙ ( t ) = L μ u t + F ( μ , u t ) ,
where u t = ( u 1 ( t )   u 2 ( t )   u 3 ( t )   u 4 ( t )   u 5 ( t ) ) T C ( [ 1 , 0 ] , R 6 ) and
L μ Φ = ( τ 0 + μ ) ( A + Φ ( 0 ) + B + Φ ( 1 ) ) , F ( μ , u t ) = ( τ 0 + μ ) ( α 1 Φ 1 ( 0 ) Φ 2 ( 0 ) α 1 Φ 1 ( 0 ) Φ 2 ( 0 ) 0 0 0 0 ) ,
where
A + = ( a 11 a 12 a 13 0 0 0 a 21 a 22 0 0 0 0 a 31 a 32 a 33 0 0 0 a 41 0 0 a 44 0 0 0 a 52 0 0 a 55 0 0 0 a 63 0 0 a 66 ) ,
and
B + = ( 0 0 0 b 14 0 0 0 0 0 0 b 25 0 0 0 0 0 0 b 36 0 0 0 b 44 0 0 0 0 0 0 b 55 0 0 0 0 0 0 b 66 )
According to the application of Riesz representation theorem, there exists a 6 × 6 matrix function η ( θ ,   μ ) : [ 1 , 0 ] R 6 × 6 such as:
L μ Φ 1 0 d η ( θ ,   μ ) Φ ( θ ) , Φ C .
The following equation is chosen:
η ( θ ,   μ ) = ( τ 0 + μ ) ( A + δ ( θ ) + B + δ ( θ + 1 ) ) .
with δ is the Dirac delta function.
For Φ C ( [ 1 , 0 ] , R 6 ) , the following equations are defined
A ( μ ) Φ = { d Φ ( θ ) d θ , 1 θ < 0 , 1 0 d η ( θ ,   μ ) Φ ( θ ) , θ = 0 ,
and
R ( μ ) Φ = { 0 , 1 θ < 0 , F ( μ ,   θ ) , θ = 0 .
The model (28) is equivalent to
u ˙ ( t ) = A ( μ ) u t + R ( μ ) u t .
A * is defined as
A * ( φ ) = { d φ ( s ) d s , 0 s < 1 , 1 0 d η T ( s ,   0 ) φ ( s ) , s = 0 ,
with a bilinear form
φ ( s ) , Φ ( θ ) = φ ¯ ( 0 ) Φ ( 0 ) θ = 1 0 ξ = 0 θ φ ¯ ( ξ θ ) d η ( θ ) Φ ( ξ ) d ξ ,
where η ( θ ) = η ( θ , 0 ) .
The eigenvector of A ( 0 ) is shown as h ( θ ) = ( 1 , h 2 , h 3 , h 4 , h 5 , h 6 ) T e i ω 0 θ   with   i ω 0 τ 0 , and the eigenvector of A * ( 0 ) is shown as h * ( s ) = D ( 1 , h 2 * , h 3 * , h 4 * , h 5 * , h 6 * ) e i ω 0 s   with i ω 0 τ 0 . Then the following results can be obtained:
{ h 2 = ( i ω 0 a 55 b 55 e i ω 0 τ 0 ) h 5 a 52 h 3 = ( i ω 0 a 66 b 66 e i ω 0 τ 0 ) h 6 a 63 h 4 = a 41 i ω 0 a 44 b 44 e i ω 0 τ 0 h 5 = a 21 a 52 ( i ω 0 a 22 ) ( i ω 0 a 55 b 55 e i ω 0 τ 0 ) b 25 a 52 e i ω 0 τ 0 h 6 = ( a 31 + a 32 h 2 ) a 63 ( i ω 0 a 33 ) ( i ω 0 a 66 b 66 e i ω 0 τ 0 ) b 36 a 63
and
{ h 2 * = ( a 55 + b 55 e i ω 0 τ 0 + i ω 0 ) h 5 * b 25 e i ω 0 τ 0 h 3 * = ( a 66 + b 66 e i ω 0 τ 0 + i ω 0 ) h 6 * b 36 e i ω 0 τ 0 h 4 * = b 14 e i ω 0 τ 0 a 44 + b 44 e i ω 0 τ 0 + i ω 0 h 5 * = ( a 12 + a 32 h 3 * ) b 25 e i ω 0 τ 0 ( a 22 + i ω 0 ) ( a 55 + b 55 e i ω 0 τ 0 + i ω 0 ) b 25 a 52 e i ω 0 τ 0 h 6 * = a 13 b 36 e i ω 0 τ 0 ( a 33 + i ω 0 ) ( a 66 + b 66 e i ω 0 τ 0 + i ω 0 ) b 36 a 63 e i ω 0 τ 0
We obtain Equation (38) from Equation (35):
h * ( s ) , h ( θ ) = D ¯ [ 1 + h 2 h ¯ 2 * + h 3 h ¯ 3 * + h 4 h ¯ 4 * + h 5 h ¯ 5 * + h 6 h ¯ 6 * + h ¯ 2 * b 25 h 5 + h ¯ 3 * b 36 h 6 + h ¯ 4 * b 44 h 4 + h ¯ 5 * b 55 h 5 + h ¯ 6 * b 66 h 6 ) ]
Equation (39) is chosen as:
D ¯ = [ 1 + h 2 h ¯ 2 * + h 3 h ¯ 3 * + h 4 h ¯ 4 * + h 5 h ¯ 5 * + h 6 h ¯ 6 * + τ 0 e i ω 0 τ 0 ( h 4 b 14 + h ¯ 2 * b 25 h 5 + h ¯ 3 * b 36 h 6    + h ¯ 4 * b 44 h 4 + h ¯ 5 * b 55 h 5 + h ¯ 6 * b 66 h 6 ) ] 1 ,
with h * , h = 1   and   h * , h ¯ = 0 .
Following the computation introduced in [36], we obtain
{ g 20 = 2 α 1 τ 0 D ¯ h 3 ( h ¯ 2 * 1 ) , g 11 = α 1 τ 0 D ¯ ( h 3 + h ¯ 3 ) ( h ¯ 2 * 1 ) , g 02 = 2 α 1 τ 0 D ¯ h ¯ 3 ( h ¯ 2 * 1 ) , g 21 = 2 α 1 τ 0 D ¯ h 3 ( h ¯ 2 * 1 ) ( W 11 ( 1 ) ( 0 ) h 3 + 1 2 W 20 ( 1 ) ( 0 ) h ¯ 3 + W 11 ( 3 ) ( 0 ) + 1 2 W 20 ( 3 ) ( 0 ) ) .
with
{ W 20 ( θ ) = i g 20 h ( 0 ) ω 0 τ 0 e i ω 0 τ 0 θ + i g ¯ 02 h ¯ ( 0 ) 3 ω 0 τ 0 e i ω 0 τ 0 θ + E 1 e 2 i ω 0 τ 0 θ W 11 ( θ ) = i g 11 h ( 0 ) ω 0 τ 0 e i ω 0 τ 0 θ + i g ¯ 11 h ¯ ( 0 ) ω 0 τ 0 e i ω 0 τ 0 θ + E 2 ,
and
E 1 = [ 2 i ω 0 a 11 a 12 a 13 b 14 e 2 i ω 0 τ 0 0 0 a 21 2 i ω 0 a 22 0 0 b 25 e 2 i ω 0 τ 0 0 a 31 a 32 2 i ω 0 a 33 0 0 b 36 e 2 i ω 0 τ 0 a 41 0 0 a 44 0 0 0 a 52 0 0 a 55 0 0 0 a 63 0 0 a 66 ] 1 [ E 2 ( 1 ) E 2 ( 2 ) 0 0 0 0 ]
E 2 = [ a 11 a 12 a 13 b 14 0 0 a 21 a 22 0 0 b 25 0 a 31 a 32 a 33 0 0 b 36 a 41 0 0 a 44 + b 44 0 0 0 a 52 0 0 a 55 + b 55 0 0 0 a 63 0 0 a 66 + b 66 ] 1 [ E 2 ( 1 ) E 2 ( 2 ) 0 0 0 0 ]
where
a 44 = 2 i ω 0 a 44 b 44 e 2 i ω 0 τ 0 ,             a 55 = 2 i ω 0 a 55 b 55 e 2 i ω 0 τ 0 , a 66 = 2 i ω 0 a 66 b 66 e 2 i ω 0 τ 0 , E 1 ( 1 ) = α 1 h 2 ,             E 1 ( 2 ) = α 1 h 2 ,             E 2 ( 1 ) = α 1 ( h 2 + h ¯ 2 ) ,             E 2 ( 2 ) = α 1 ( h 2 + h ¯ 2 )  
Thus, the following values can be obtained:
{ C 1 ( 0 ) = i 2 τ 0 ω 0 ( g 11 g 20 2 | g 11 | 2 | g 02 | 2 3 ) + g 21 2 ,      μ 2 = R e { C 1 ( 0 ) } R e { λ ( τ 0 ) } ,      β 2 * = 2 R e { C 1 ( 0 ) } ,      T 2 = I m { C 1 ( 0 ) } + μ 2 I m { λ ( τ 0 ) } τ 0 ω 0 .
Theorem 2 is given as:
Theorem 2.
For the model (2), if μ 2 > 0 ( μ 2 < 0 ), the Hopf bifurcation is supercritical (subcritical). If β 2 * < 0 ( β 2 * > 0 ), the bifurcating periodic solutions are stable (unstable). If T 2 > 0 ( T 2 < 0 ), the bifurcating periodic solutions increase (decrease).

5. Simulation

5.1. Parameter Dependence of R 0

In this section, prevention methods are analyzed in system (2).
The basic reproduction number R 0 is generally used to describe the infectivity of virus. Therefore, it is also applicable to the spread of malware in WRSNs. On the assumption that the system is stable, if R 0 1 , malware will be cleared eventually. If R 0 > 1 , malware will always exist.
The size of R 0 is affected by different parameters. What is worth concerning about is the influence of the parameter ( α 2 , α 3 , α 4   or   α 5 ) on the malware spread under different diffusion rates α 1 . Here, we first assume that only one parameter ( α 2 , α 3 , α 4   or   α 5 ) and α 1 are changed while the other parameters are unchanged. We observe the spread of malware in the system.
(1)
The parameters are as follows in Figure 1a: b = 0.05 ,   α 3 = 0.15 ,   α 4 = 0.3 ,   α 5 = 0.7 ,   α 1 [ 0 ,   1 ] , and   α 2 [ 0 ,   1 ] . Figure 1a shows that with the increase of the diffusion rate α 1 , the malware is more easier prevailing. However, the self-healing rate of the infected nodes α 2 has a inhibitory effect on the spread of malware. Figure 1a also shows that the charging behavior ( β 2 = 0.2 ) makes the spread of malware easier.
(2)
The parameters are as follows in Figure 1b: b = 0.05 ,   α 2 = 0.02 ,   α 4 = 0.3 ,   α 5 = 0.7 ,   α 1 [ 0 ,   1 ] , and   α 3 [ 0 ,   1 ] . Figure 1b shows that with the increase of the recovery rate of the infected nodes α 3 , the spread of malware can be effectively suppressed, which will provide us with the reference value in the control of malware. Similarly, Figure 1b also shows that without the charging behavior ( β 2 = 0 ), the control of malware will become easier.
(3)
The parameter settings are as follows in Figure 1c: b = 0.05 ,   α 2 = 0.02 ,   α 3 = 0.15 ,   α 5 = 0.7 ,   α 1 [ 0 ,   1 ] , and   α 4 [ 0 ,   1 ] , and in Figure 1d: b = 0.05 ,   α 2 = 0.02 ,   α 3 = 0.15 ,   α 4 = 0.3 ,   α 1 [ 0 ,   1 ] , and   α 5 [ 0 ,   1 ] . It is shown that the increase of the immune prevention rate α 4 can suppress the spread of malware but the increase of the immune failure rate α 5 helps the malware to spread.
(4)
Figure 1a–d together reflect that the charging behavior ( β 2 = 0.2 ) will encourage the spread of malware, which will provide us with the data reference for the control of the malware spread in the SIRS-L model.
Next, the influence on the malware spread under different low-energy-node conversion rate β 1 and charging success rate β 2 are also discussed.
(5)
The parameter settings are as follows in Figure 1e: b = 0.05 ,   α 1 = 0.87 , α 2 = 0.02 ,   α 3 = 0.15 ,   α 4 = 0.3 ,   α 5 = 0.7 ,   β 1 [ 0 ,   1 ] , and   β 2 [ 0 ,   1 ] . Figure 1e reflects the influence of the low-energy node conversion rate β 1 and the charging success rate β 2 on malware propagation. It is shown that if the low-energy node conversion rate β 1 is less than 0.2, the increase of the charging success rate β 2 can easily lead to the prevalence of malware. On the other side, if the low-energy node conversion rate β 1 is larger than 0.4, we can appropriately reduce the charging success rate β 2 to suppress the spread of malware.

5.2. Analysis and Display of Equilibrium Solutions

In this section, the distribution of nodes’ states under different parameters of R 0 and τ are discussed to verify the stability of the system and to discuss the bifurcation.
(1)
The parameter settings are as follows in Figure 2: b = 0.05 ,   α 1 = 0.87 , α 2 = 0.02 ,   α 3 = 0.15 ,   α 4 = 0.3 ,   α 5 = 0.7 ,   β 1 = 0.2 , β 2 = 0 ,   a n d   R 0 = 0.3149 ,   τ = 0 with the initial value of nodes’ scale: (s(0), i(0), r(0), ls(0), li(0), lr(0))= (0.9, 0.1, 0, 0, 0, 0). It is shown that (s(∞), i(∞), r(∞), ls(∞), li(∞), lr(∞)) = (0.1520, 0, 0.0480, 0.6079, 0, 0.1920). It is shown that if the malware appears in the model (2), it will gradually disappear if R 0 < 1 . It also can be shown that if β 2 = 0 , the total proportion of the low-energy nodes (ls, li, lr) is larger than that of the general nodes (s, i, r).
(2)
The parameter settings are as follows in Figure 3: b = 0.05 ,   α 1 = 0.87 , α 2 = 0.02 ,   α 3 = 0.15 ,   α 4 = 0.3 ,   α 5 = 0.7 ,   β 1 = 0.02 , β 2 = 0 ,   a n d   R 0 = 1.8633 , τ = 0 , with the initial value of nodes’ scale: (s(0), i(0), r(0), ls(0), li(0), lr(0))= (0.9, 0.1, 0, 0, 0, 0). It is shown that (s(∞), i(∞), r(∞), ls(∞), li(∞), lr(∞)) = (0.2756, 0.2770, 0.1616, 0.1103, 0.1109, 0.0646). It shows that the malware will prevail if R 0 > 1 .
(3)
The parameter settings are as follows in Figure 4: b = 0.05 ,   α 1 = 0.87 , α 2 = 0.02 ,   α 3 = 0.15 ,   α 4 = 0.3 ,   α 5 = 0.7 ,   β 1 = 0.2 , β 2 = 0.2 ,   a n d   R 0 = 0.5743 , τ = 0 , with the initial value of nodes’ scale: (s(0), i(0), r(0), ls(0), li(0), lr(0))= (0.9, 0.1, 0, 0, 0, 0). It is shown that (s(∞), i(∞), r(∞), ls(∞), li(∞), lr(∞)) = (0.4027, 0, 0.1529, 0.3221, 0, 0.1223). It shows that if the malware appears in the model (2), it will gradually disappear if R 0 < 1 . It also can be shown that if β 2 = 0.2 , the total proportion of the low-energy nodes (ls, li, lr) is smaller than the proportion of the low-energy nodes in Figure 2 if the charging operation is not performed β 2 = 0 .
(4)
The parameter settings are as follows in Figure 5: b = 0.05 ,   α 1 = 0.87 , α 2 = 0.02 ,   α 3 = 0.15 ,   α 4 = 0.3 ,   α 5 = 0.7 ,   β 1 = 0.2 , β 2 = 0.2 ,   a n d   R 0 = 1.3473 , τ = 0 , with the initial value of nodes’ scale: (s(0), i(0), r(0), ls(0), li(0), lr(0)) = (0.9, 0.1, 0, 0, 0, 0). It is shown that (s(∞), i(∞), r(∞), ls(∞), li(∞), lr(∞)) = (0.2988, 0.1204, 0.1364, 0.2391, 0.0963, 0.1091). It shows that the malware will prevail if R 0 > 1 . Similarly, it also can be shown that if β 2 = 0.2 , the total proportion of the low-energy nodes (ls, li, lr) is smaller than the proportion of the low-energy nodes in Figure 2 if the charging operation is not performed β 2 = 0 . The total proportion of the low-energy nodes is similar to the corresponding proportion in Figure 4, which represents the proportion of low-energy nodes is only related to β 1 and β 2 .
The above results (1)–(4) can be summarized as follows: the system is locally stable on the premise that the parameters meet the hypothesis (H1). In addition, if R 0   1 , the malicious software will persist in the system, while if R 0 < 1 , the malicious software will be finally cleared. Next, we verify the bifurcation of the system.
(5)
The parameter settings are as follows in Figure 6 and Figure 7: b = 0.05 ,   α 1 = 0.9059 , α 2 = 0.02 ,   α 3 = 0.15 ,   α 4 = 0.3 ,   α 5 = 0.7 ,   β 1 = 0.2 , β 2 = 0.2 ,   a n d   R 0 = 1.4029 , with the initial value of nodes’ scale: (s(0), i(0), r(0), ls(0), li(0), lr(0)) = (0.9, 0.1, 0, 0, 0, 0). It can be calculated that τ 0 = 11.5370 and ω 0 = 0.2254 . In Figure 6, it can be seen that the model (2) is asymptotically stable if τ = 10.25 < τ 0 = 11.5370 . What’s more, in Figure 7, the model (2) undergoes a Hopf bifurcation if τ = 11.74 > τ 0 = 11.5370 . According to Equation (41) and Theorem 2, the following parameters can be obtained: λ ( τ 0 ) = 0.0292 + 0.0328 i ,   C 1 ( 0 ) = 1.4520 + 0.7622 i ,   μ 2 = 49.7260 > 0 ,   β 2 * = 2.9040 < 0 ,   and   T 2 = 0.9203 < 0 . The result can be concluded that the Hopf bifurcation is supercritical, the bifurcating periodic solutions are stable and the period of the periodic solutions decreases.

6. Conclusions

This article puts forward a novel SISR-L model that considers the low-energy status nodes. Combined with the reality, the charging delay τ is introduced in the model to study the bifurcation. We have proved that the system is locally stable under certain conditions. In addition, the influence of different parameters on the spread of malware is displayed in the simulation. It is worth noting that in the SISR-L system, the charging behavior will make it more difficult to control the spread of malware. On the analysis of bifurcation, we reveal that the SISR-L model is locally asymptotically stable if the charging delay is less than τ 0 , and the model undergoes a Hopf bifurcation at the endemic solution e + if the charging delay is larger than τ 0 . Finally, the properties of Hopf bifurcation were explored through applying the normal form method and the center manifold theorem. All the analyses can provide theoretical reference for the control of malware propagation in the SISR-L model.

Author Contributions

Conceptualization, G.L., J.L. and Z.P.; methodology, G.L. and J.L.; software, G.L. and J.L.; validation, G.L. and J.L.; formal analysis, G.L. and J.L.; investigation, G.L., J.L. and Z.L.; writing—original draft preparation, J.L.; writing—review and editing, G.L., J.L. and Z.L. All authors have read and agreed to the published version of the manuscript.

Funding

The authors acknowledge funding received from the following science foundations: the National Natural Science Foundation of China (61403089, 51975136, 52075109); the 2020 Department of Education of Guangdong Province Innovative and Strong School Project (Natural Sciences)—Young Innovators Project (Natural Sciences) under grant 2020KQNCX054; the National Key Research and Development Program of China (2018YFB2000501); the Science and Technology Innovative Research Team Program in Higher Educational Universities of Guangdong Province (2017KCXTD025); the Innovative Academic Team Project of Guangzhou Education System (1201610013); the Special Research Projects in the Key Fields of Guangdong Higher Educational Universities (2019KZDZX1009); the Science and Technology Research Project of Guangdong Province (2017A010102014, 2016A010102022); and the Science and Technology Research Project of Guangzhou (201707010293). They are all appreciated for supporting this work.

Institutional Review Board Statement

Not applicable.

Informed Consent Statement

Not applicable.

Data Availability Statement

Not applicable.

Conflicts of Interest

No conflict of interest.

Appendix A

A 5 = a 11 a 22 a 33 a 44 a 55 a 66 , A 4 = a 66 ( a 11 + a 22 + a 33 + a 44 + a 55 ) + a 11 a 22 a 12 a 21 a 13 a 31 + a 44 ( a 11 + a 22 + a 33 ) + a 55 ( a 11 + a 22 + a 33 + a 44 ) + a 33 ( a 11 + a 22 ) , A 3 = a 13 ( a 22 a 31 a 21 a 32 ) ( a 33 + a 44 + a 55 + a 66 ) ( a 11 a 22 a 12 a 21 ) + a 44 ( a 13 a 31 a 33 ( a 11 + a 22 ) ) + a 55 ( a 13 a 31 a 44 ( a 11 + a 22 + a 33 ) a 33 ( a 11 + a 22 ) ) + a 66 ( a 13 a 31 a 44 ( a 11 + a 22 + a 33 ) a 55 ( a 11 + a 22 + a 33 + a 44 ) a 33 ( a 11 + a 22 ) ) , A 2 = a 44 ( a 13 a 21 a 32 a 13 a 22 a 31 ) + ( a 33 ( a 44 + a 55 + a 66 ) + ( a 44 + a 55 ) a 66 + a 44 a 55 ) ( a 11 a 22 a 12 a 21 ) a 66 ( a 44 ( a 13 a 31 a 33 ( a 11 + a 22 ) ) + a 55 ( a 13 a 31 a 44 ( a 11 + a 22 + a 33 ) a 33 ( a 11 + a 22 ) ) + ( a 55 a 44 a 13 ) ( a 21 a 32 a 22 a 31 ) ) a 55 a 44 · ( a 13 a 31 a 33 ( a 11 + a 22 ) ) , A 1 = a 66 ( a 44 ( a 33 ( a 11 a 22 a 12 a 21 ) + a 13 ( a 21 a 32 a 22 a 31 ) ) + a 55 ( a 33 ( a 11 a 22 a 12 a 21 ) + a 44 ( a 11 a 22 a 12 a 21 a 13 · a 31 + a 33 ( a 11 + a 22 ) ) + a 13 ( a 21 a 32 a 22 a 31 ) ) ) a 44 a 55 ( a 33 ( a 11 a 22 a 12 a 21 ) + a 13 ( a 21 a 32 a 22 a 31 ) ) , A 0 = a 44 a 55 a 66 ( a 33 ( a 11 a 22 a 12 a 21 ) + a 13 ( a 21 a 32 a 22 a 31 ) ) , B 5 = b 44 b 55 b 66 , B 4 = b 44 ( a 11 + a 22 + a 33 + a 55 + a 66 ) + b 55 ( a 11 + a 22 + a 33 + a 44 + a 66 ) + b 66 ( a 11 + a 22 + a 33 + a 44 + a 55 ) a 41 b 14 a 52 b 25 a 63 b 36 , B 3 = ( b 66 + b 44 + b 44 b 55 ) ( a 11 a 22 a 12 a 21 a 13 a 31 ) ( b 66 a 55 + b 44 a 66 + b 55 a 66 ) ( a 11 + a 22 + a 33 + a 44 ) ( b 66 a 44 + b 55 a 44 + b 44 a 55 ) ( a 11 + a 22 + a 33 ) ( b 66 a 33 + b 44 a 33 b 44 b 55 a 33 ( a 11 + a 22 ) + ( a 22 + a 55 ) ( a 41 b 14 + a 63 b 36 ) + ( a 11 + a 44 ) ( a 52 b 25 + a 63 b 36 ) + ( a 33 + a 66 ) ( a 41 b 14 + a 52 b 25 ) , B 2 = b 55 a 33 ( a 11 a 22 a 12 a 21 ) + ( a 44 b 55 + a 55 + a 66 b 55 + a 44 b 66 + a 55 b 44 ) ( a 11 a 22 a 12 a 21 a 13 a 31 ) + ( a 44 a 33 b 55 + a 55 a 33 + a 66 a 33 b 55 + a 33 b 66 + a 55 a 33 b 44 + a 66 a 33 b 44 ) ( a 11 + a 22 ) + ( a 66 a 44 b 55 + a 55 a 44 + a 55 a 66 b 44 ) ( a 11 + a 22 + a 33 ) + b 55 a 13 ( a 21 a 32 a 22 a 31 ) a 33 ( a 22 a 41 b 14 + a 11 a 52 b 25 ) a 44 ( a 11 ( a 52 b 25 + a 63 b 36 ) + a 33 a 52 b 25 + a 22 a 63 b 36 ) a 66 ( a 33 ( a 41 b 14 + a 52 b 25 ) + a 22 a 41 b 14 + a 11 a 52 b 25 + a 41 a 55 b 14 + a 44 a 52 b 25 ) a 55 ( a 41 b 14 + a 63 b 36 ) + a 33 a 41 b 14 + b 36 ( a 11 a 63 + a 44 a 63 ) ) ) + ( b 44 a 13 a 13 a 66 ) ( a 21 a 32 a 22 a 31 ) + ( a 44 b 66 a 33 a 63 b 36 + b 44 a 33 ) ( a 11 a 22 a 12 a 21 ) + a 13 a 31 a 52 b 25 , B 1 = a 66 ( a 33 ( a 22 a 41 b 14 + a 11 a 52 b 25 ) + a 55 ( a 22 a 41 b 14 + a 33 a 41 b 14 ) + a 44 ( a 11 a 52 b 25 + a 33 a 52 b 25 ) a 13 a 31 a 52 b 25 ) b 55 · ( a 44 ( a 33 ( a 11 a 22 a 12 a 21 ) ) + a 66 ( a 33 ( a 11 a 22 a 12 a 21 ) ) ) b 66 ( a 44 ( a 33 ( a 11 a 22 a 12 a 21 ) + a 13 ( a 21 a 32 a 22 a 31 ) ) + a 55 ( a 33 ( a 11 a 22 a 12 a 21 ) ) ) b 44 ( a 55 ( a 33 ( a 11 a 22 a 12 a 21 ) + a 13 ( a 21 a 32 a 22 a 31 ) ) ) ) + a 44 ( a 11 a 33 a 52 b 25 a 13 a 31 · a 52 b 25 + a 11 a 22 a 63 b 36 a 12 a 21 a 63 b 36 ) + a 55 ( a 44 ( a 11 a 63 b 36 + a 22 a 63 b 36 ) + a 22 a 33 a 41 b 14 + a 11 a 22 a 63 b 36 a 12 a 21 · a 63 b 36 ) + ( b 55 a 44 a 13 + a 55 a 66 a 33 b 44 b 66 a 13 a 44 + ( 1 a 66 ) a 55 a 13 b 44 ) ( a 21 a 32 a 22 a 31 ) + ( a 44 a 44 b 55 a 66 + b 66 a 44 · a 44 + a 66 a 55 a 55 b 44 ) ( a 11 a 22 a 12 a 21 a 13 a 31 + a 33 ( a 11 + a 22 ) ) , B 0 = a 66 ( a 44 ( ( a 11 a 33 a 13 a 31 ) a 52 b 25 ) + a 22 a 33 a 41 a 55 b 14 ) a 44 a 55 ( a 11 a 22 a 12 a 21 ) a 63 b 36 + a 44 a 55 b 66 ( a 33 ( a 11 a 22 a 12 · a 21 ) + a 13 ( a 21 a 32 a 22 a 31 ) ) + a 44 a 66 b 55 ( a 33 ( a 11 a 22 a 12 a 21 ) + a 13 ( a 21 a 32 a 22 a 31 ) ) + a 55 a 66 b 44 ( a 33 ( a 11 a 22 a 12 · a 21 ) + a 13 ( a 21 a 32 a 22 a 31 ) ) , C 4 = b 44 b 55 + b 66 ( b 44 + b 55 ) , C 3 = b 55 ( a 41 b 14 + a 63 b 36 b 44 ( a 11 + a 22 + a 33 + a 66 ) ) + b 44 ( a 52 b 25 + a 63 b 36 ) + b 66 ( a 41 b 14 + a 52 b 25 b 44 ( a 11 + a 22 + a 33 + a 55 ) b 55 ( a 11 + a 22 + a 33 + a 44 ) ) , C 2 = a 41 a 52 b 14 b 25 b 44 ( a 11 ( a 52 b 25 + a 63 b 36 ) + ( a 33 a 52 b 25 + a 22 a 63 b 36 + a 52 a 66 b 25 + a 55 a 63 b 36 ) b 66 ( a 33 ( a 41 b 14 + a 52 · b 25 ) b 44 ( a 11 a 22 a 12 a 21 a 13 a 31 + a 55 ( a 11 + a 22 + a 33 ) + a 33 ( a 11 + a 22 ) ) b 55 ( a 11 a 22 a 12 a 21 a 13 a 31 + a 44 ( a 11 + a 22 + a 33 ) + a 33 ( a 11 + a 22 ) ) + a 22 a 41 b 14 + a 11 a 52 b 25 + a 41 a 55 b 14 + a 44 a 52 b 25 ) b 55 ( a 22 ( a 41 b 14 + a 63 b 36 ) b 44 · ( a 11 a 22 a 12 a 21 a 13 a 31 + a 66 ( a 11 + a 22 + a 33 ) + a 33 ( a 11 + a 22 ) ) + a 33 a 41 b 14 + a 11 a 63 b 36 + a 41 a 66 b 14 + a 44 a 63 b 36 ) + a 41 a 63 b 14 b 36 + a 52 a 63 b 25 b 36 , C 1 = b 44 ( a 66 ( a 11 a 52 b 25 + a 33 a 52 b 25 ) + a 55 ( a 11 a 63 b 36 + a 22 a 63 b 36 ) + a 11 a 33 a 52 b 25 a 13 a 31 a 52 b 25 + a 11 a 22 a 63 b 36 a 12 · a 21 a 63 b 36 ) + b 55 ( a 66 ( a 22 a 41 b 14 + a 33 a 41 b 14 ) + a 44 ( a 11 a 63 b 36 + a 22 a 63 b 36 ) b 44 ( a 33 ( a 11 a 22 a 12 a 21 ) + a 66 ( a 11 · a 22   a 12 a 21 a 13 a 31 + a 33 ( a 11 + a 22 ) ) + a 13 a 21 a 32 a 13 a 22 a 31 ) + a 22 a 33 a 41 b 14 + a 11 a 22 a 63 b 36 a 12 a 21 a 63 b 36 ) + b 66 · ( a 33 ( a 22 a 41 b 14 + a 11 a 52 b 25 ) + a 55 ( a 22 a 41 b 14 + a 33 a 41 b 14 ) + a 44 ( a 11 a 52 b 25 + a 33 a 52 b 25 ) b 44 ( a 33 ( a 11 a 22 a 12 a 21 ) + a 55 ( a 11 a 22 a 12 a 21 a 13 a 31 + a 33 ( a 11 + a 22 ) ) + a 13 a 21 a 32 a 13 a 22 a 31 ) b 55 ( a 33 ( a 11 a 22 a 12 a 21 ) + a 44 ( a 11 a 22 a 12 · a 21 a 13 a 31 + a 33 ( a 11 + a 22 ) ) + a 13 a 21 a 32 a 13 a 22 a 31 ) a 13 a 31 a 52 b 25 ) a 33 a 41 a 52 b 14 b 25 a 22 a 41 a 63 b 14 b 36 a 11 a 52 · a 63 b 25 b 36 a 41 a 52 a 66 b 14 b 25 a 41 a 55 a 63 b 14 b 36 a 44 a 52 a 63 b 25 b 36 , C 0 = b 66 ( a 44 ( a 11 a 33 a 52 b 25 a 13 a 31 a 52 b 25 ) a 44 b 55 ( a 33 ( a 11 a 22 a 12 a 21 ) + a 13 a 21 a 32 a 13 a 22 a 31 ) a 55 b 44 ( a 33 ( a 11 a 22   a 12 a 21 ) + a 13 a 21 a 32 a 13 a 22 a 31 ) + a 22 a 33 a 41 a 55 b 14 ) b 55 ( a 44 ( a 11 a 22 a 63 b 36 a 12 a 21 a 63 b 36 ) a 66 b 44 ( a 33 ( a 11 a 22 a 12 · a 21 ) + a 13 a 21 a 32 a 13 a 22 a 31 ) + a 22 a 33 a 41 a 66 b 14 ) b 44 ( a 66 ( a 11 a 33 a 52 b 25 a 13 a 31 a 52 b 25 ) + a 55 ( a 11 a 22 a 63 b 36 a 12 · a 21 a 63 b 36 ) ) + a 11 a 44 a 52 a 63 b 25 b 36 + a 22 a 41 a 55 a 63 b 14 b 36 + a 33 a 41 a 52 a 66 b 14 b 25 , D 3 = b 44 b 55 b 66 ,   D 2 = b 66 ( b 55 ( a 41 b 14 b 44 ( a 11 + a 22 + a 33 ) ) + a 52 b 25 b 44 ) a 63 b 36 b 44 b 55 , D 1 = b 55 ( b 44 ( a 11 a 63 b 36 + a 22 a 63 b 36 ) a 41 a 63 b 14 b 36 ) + b 66 ( b 44 ( a 11 a 52 b 25 + a 33 a 52 b 25 ) + b 55 ( a 22 a 41 b 14 b 44 ( a 11 a 22 a 12 a 21 a 13 a 31 + a 33 ( a 11 + a 22 ) ) + a 33 a 41 b 14 ) a 41 a 52 b 14 b 25 ) a 52 a 63 b 25 b 36 b 44 , D 0 = b 55 ( b 44 ( a 11 a 22 a 63 b 36 a 12 a 21 a 63 b 36 ) a 22 a 41 a 63 b 14 b 36 ) + b 66 ( b 55 ( b 44 ( a 33 ( a 11 a 22 a 12 a 21 ) + a 13 a 21 a 32 a 13 · a 22 a 31 ) a 22 a 33 a 41 b 14 ) b 44 ( a 11 a 33 a 52 b 25 a 13 a 31 a 52 b 25 ) + a 33 a 41 a 52 b 14 b 25 ) + a 11 a 52 a 63 b 25 b 36 b 44 a 41 a 52 a 63 · b 14 b 25 b 36 ,

References

  1. Rashid, B.; Rehmani, M.H. Applications of wireless sensor networks for urban areas: A survey. J. Netw. Comput. 2016, 60, 192–219. [Google Scholar] [CrossRef]
  2. Oliviero, F.; Romano, S.P. A reputation-based metric for secure routing in wireless mesh networks. In Proceedings of the IEEE GLOBECOM 2008–2008 IEEE Global Telecommunications Conference, New Orleans, LA, USA, 30 November–4 December 2008; pp. 1–5. [Google Scholar]
  3. Rehman, A.U.; Rehman, S.U.; Raheem, H. Sinkhole attacks in wireless sensor networks: A survey. Wirel. Pers. Commun. 2019, 106, 2291–2313. [Google Scholar] [CrossRef]
  4. Alajmi, N. Wireless sensor networks attacks and solutions. arXiv 2014, arXiv:1407.6290. [Google Scholar]
  5. Ngai, E.C.; Liu, J.; Lyu, M.R. On the intruder detection for sinkhole attack in wireless sensor networks. In Proceedings of the 2006 IEEE International Conference on Communications, Istanbul, Turkey, 11–15 June 2006; Volume 8, pp. 3383–3389. [Google Scholar]
  6. Hu, Y.C.; Perrig, A.; Johnson, D.B. Packet leashes: A defense against wormhole attacks in wireless networks. In Proceedings of the IEEE INFOCOM 2003. Twenty-second Annual Joint Conference of the IEEE Computer and Communications Societies (IEEE Cat. No. 03CH37428), San Francisco, CA, USA, 30 March–3 April 2003; Volume 3, pp. 1976–1986. [Google Scholar]
  7. Yu, B.; Xiao, B. Detecting selective forwarding attacks in wireless sensor networks. In Proceedings of the 20th IEEE International Parallel & Distributed Processing Symposium, Rhodes Island, Greece, 25–29 April 2006; p. 8. [Google Scholar]
  8. Salehi, S.A.; Razzaque, M.A.; Naraei, P.; Farrokhtala, A. Detection of sinkhole attack in wireless sensor networks. In Proceedings of the 2013 IEEE International Conference on Space Science and Communication (IconSpace), Melaka, Malaysia, 1–3 July 2013; pp. 361–365. [Google Scholar]
  9. Sundararajan, R.K.; Arumugam, U. Intrusion detection algorithm for mitigating sinkhole attack on LEACH protocol in wireless sensor networks. J. Sens. 2015, 2015. [Google Scholar] [CrossRef] [Green Version]
  10. Chen, C.; Song, M.; Hsieh, G. Intrusion detection of sinkhole attacks in large-scale wireless sensor networks. In Proceedings of the 2010 IEEE International Conference on Wireless Communications, Networking and Information Security, Beijing, China, 25–27 June 2010; pp. 711–716. [Google Scholar]
  11. Liang, C.J.M.; Musăloiu-e, R.; Terzis, A. Typhoon: A Reliable Data Dissemination Protocol for Wireless Sensor Networks. In European Conference on Wireless Sensor Networks; Springer: Berlin/Heidelberg, Germany, 2008; pp. 268–285. [Google Scholar]
  12. Song, Y.; Jiang, G.P. Model and Dynamic Behavior of Malware Propagation over Wireless Sensor Networks. In International Conference on Complex Sciences; Springer: Berlin/Heidelberg, Germany, 2009; pp. 487–502. [Google Scholar]
  13. Yetgin, H.; Cheung, K.T.K.; El-Hajjar, M.; Hanzo, L.H. A Survey of Network Lifetime Maximization Techniques in Wireless Sensor Networks. IEEE Commun. Surv. Tutor. 2017, 19, 828–854. [Google Scholar] [CrossRef] [Green Version]
  14. Rasheed, A.; Mahapatra, R.N. The Three-Tier Security Scheme in Wireless Sensor Networks with Mobile Sinks. IEEE Trans. Parallel Distrib. Syst. 2010, 23, 958–965. [Google Scholar] [CrossRef]
  15. Lin, C.; Shang, Z.; Du, W.; Ren, J.K.; Wang, L.; Wu, G.W. CoDoC: A Novel Attack for Wireless Rechargeable Sensor Networks through Denial of Charge. In Proceedings of the IEEE INFOCOM, Paris, France, 29 April–2 May 2019. [Google Scholar]
  16. Desmedt, Y.; Frankel, Y. Shared generation of authenticators and signatures. In Annual International Cryptology Conference; Springer: Berlin/Heidelberg, Germany, 1991; pp. 457–469. [Google Scholar]
  17. Weidong, C.; Dengguo, F. A group of threshold group-signature schemes with privilege subsets. In Progress on Cryptography; Springer: Boston, MA, USA, 2004; pp. 81–88. [Google Scholar]
  18. Yi, S.; Dengguo, F. The Design and analysis of a new group of (tj, t, n) threshold group signature scheme. China Crypto. 2000. Available online: https://en.cnki.com.cn/Article_en/CJFDTotal-ZKYB200102002.htm (accessed on 20 July 2021).
  19. Wang, X.; Dong, Y. Threshold group signature scheme with privilege subjects based on ECC. In Proceedings of the 2010 International Conference on Communications and Intelligence Information Security, Nanning, China, 13–14 October 2010; pp. 84–87. [Google Scholar]
  20. Tanachaiwiwat, S.; Helmy, A. Encounter-based worms: Analysis and defense. Ad Hoc Netw. 2009, 7, 1414–1430. [Google Scholar] [CrossRef] [Green Version]
  21. Khayam, S.A.; Radha, H. Using signal processing techniques to model worm propagation over wireless sensor networks. IEEE Signal Process. Mag. 2006, 23, 164–169. [Google Scholar] [CrossRef]
  22. Batista, F.K.; Martin del Rey, A.; Queiruga-Dios, A. A new individual-based model to simulate malware propagation in wireless sensor networks. Mathematics 2020, 8, 410. [Google Scholar] [CrossRef] [Green Version]
  23. Zhang, Z.; Si, F. Dynamics of a delayed SEIRS-V model on the transmission of worms in a wireless sensor network. Adv. Differ. Equ. 2014, 2014, 1–15. [Google Scholar] [CrossRef] [Green Version]
  24. Zhang, Z.; Wang, Y. Bifurcation Analysis for an SEIRS-V Model with Delays on the Transmission of Worms in a Wireless Sensor Network. Math. Probl. Eng. 2017, 2017. [Google Scholar] [CrossRef] [Green Version]
  25. Liu, G.; Peng, B.; Zhong, X. A Novel Epidemic Model for Wireless Rechargeable Sensor Network Security. Sensors 2021, 21, 123. [Google Scholar] [CrossRef]
  26. Liu, G.; Peng, B.; Zhong, X. Epidemic Analysis of Wireless Rechargeable Sensor Networks Based on an Attack–Defense Game Model. Sensors 2021, 21, 594. [Google Scholar] [CrossRef]
  27. Liu, G.; Peng, B.; Zhong, X.; Cheng, L.; Li, Z. Attack-Defense Game between Malicious Programs and Energy-Harvesting Wireless Sensor Networks Based on Epidemic Modeling. Complexity 2020, 2020, 1–19. [Google Scholar] [CrossRef]
  28. Liu, G.; Peng, B.; Zhong, X.; Lan, X. Differential Games of Rechargeable Wireless Sensor Networks against Malicious Programs Based on SILRD Propagation Model. Complexity 2020, 2020, 13. [Google Scholar] [CrossRef]
  29. Liu, G.; Li, J.; Liang, Z.; Peng, Z. Analysis of Time-Delay Epidemic Model in Rechargeable Wireless Sensor Networks. Mathematics 2021, 9, 978. [Google Scholar] [CrossRef]
  30. Liu, G.; Peng, Z.; Liang, Z.; Li, J.; Cheng, L. Dynamics Analysis of a Wireless Rechargeable Sensor Network for Virus Mutation Spreading. Entropy 2021, 23, 572. [Google Scholar] [CrossRef]
  31. Liu, G.; Huang, Z.; Wu, X.; Liang, Z.; Hong, F.; Su, X. Modelling and Analysis of the Epidemic Model under Pulse Charging in Wireless Rechargeable Sensor Networks. Entropy 2021, 23, 927. [Google Scholar] [CrossRef]
  32. Liu, G.; Shu, C.; Liang, Z.; Peng, B.; Cheng, L. A Modified Sparrow Search Algorithm with Application in 3d Route Planning for UAV. Sensors 2021, 21, 1224. [Google Scholar] [CrossRef]
  33. Zhu, L.; Guan, G. Dynamical Analysis of a Rumor Spreading Model with Self-Discrimination and Time Delay in Complex Networks. Phys. A Stat. Mech. Appl. 2019, 533, 121953. [Google Scholar] [CrossRef]
  34. Diekmann, O.; Heesterbeek, H.; Britton, T. Mathematical Tools for Understanding Infectious Disease Dynamics; Princeton University Press: Princeton, NJ, USA, 2012; Volume 7. [Google Scholar]
  35. Cooke, K.L.; Van Den Driessche, P. On zeroes of some transcendental equations. Funkc. Ekvacioj 1986, 29, 77–90. [Google Scholar]
  36. Hassard, B.D.; Kazarinoff, N.D.; Wan, Y.H.; Wan, Y.W. Theory and Applications of Hopf Bifurcation; CUP Archive: Cambridge, UK, 1981; Volume 41. [Google Scholar]
Figure 1. Mutual influence of parameters. (a) The relationships between α 1 ,   α 2 , and R 0 . (b) The relationships between α 1 ,   α 3 , and R 0 . (c) The relationships between α 1 , α 4 , and R 0 . (d) The relationships between α 1 , α 5 , and R 0 . (e) The relationships between β 1 , β 2 , and R 0 .
Figure 1. Mutual influence of parameters. (a) The relationships between α 1 ,   α 2 , and R 0 . (b) The relationships between α 1 ,   α 3 , and R 0 . (c) The relationships between α 1 , α 4 , and R 0 . (d) The relationships between α 1 , α 5 , and R 0 . (e) The relationships between β 1 , β 2 , and R 0 .
Mathematics 09 02007 g001
Figure 2. The distribution of nodes if R 0 = 0.3149 and β 2 = 0 .
Figure 2. The distribution of nodes if R 0 = 0.3149 and β 2 = 0 .
Mathematics 09 02007 g002
Figure 3. The distribution of nodes if R 0 = 1.8633 and β 2 = 0 .
Figure 3. The distribution of nodes if R 0 = 1.8633 and β 2 = 0 .
Mathematics 09 02007 g003
Figure 4. The distribution of nodes if R 0 = 0.5743 and β 2 = 0.2 .
Figure 4. The distribution of nodes if R 0 = 0.5743 and β 2 = 0.2 .
Mathematics 09 02007 g004
Figure 5. The distribution of nodes if R 0 = 1.3473 and β 2 = 0.2 .
Figure 5. The distribution of nodes if R 0 = 1.3473 and β 2 = 0.2 .
Mathematics 09 02007 g005
Figure 6. The distribution of nodes if R 0 = 1.4029 and τ = 10.25 < τ 0 = 11.5370 .
Figure 6. The distribution of nodes if R 0 = 1.4029 and τ = 10.25 < τ 0 = 11.5370 .
Mathematics 09 02007 g006
Figure 7. The distribution of nodes if R 0 = 1.4029 and τ = 11.74 > τ 0 = 11.5370 .
Figure 7. The distribution of nodes if R 0 = 1.4029 and τ = 11.74 > τ 0 = 11.5370 .
Mathematics 09 02007 g007
Table 1. Related research on the modeling and analysis of this paper.
Table 1. Related research on the modeling and analysis of this paper.
AuthorsModelCharacteristicsReference Content
Zhu et al. [33]SIRS (susceptible–infected–recovered–susceptible)The authors put forward the time delay of the immune validity; the SIRS model is applied to WSNs analysis.SIRS model is taken as the premise of modeling; the SIRS-L model considering the low-energy state nodes is proposed in this paper.
Zhang et al. [23]SEIRS-V (susceptible–exposed–infected-recovered–susceptible and vaccinated)Time delay is applied to SEIRS-V model.It provides the analysis reference of Hopf bifurcation and the corresponding mathematical processing method.
Liu et al. [25]SIS-L
(susceptible–infected–susceptible–low-energy status)
It first proposes the low-energy state nodes and combines them into the research of WSRNs.It provides a theoretical basis for the low-energy status modeling.
Liu et al. [26]SIAS-L
(susceptible–infected–anti-malware–susceptible–low-energy status)
The status of anti-malware is proposed and the optimal strategy is considered.It provides a theoretical basis for the low-energy status modeling.
Liu et al. [29]SIS-L
(susceptible–infected–susceptible–low-energy status)
Time delay is considered for the first time in the model with the low-energy state nodes. However, the bifurcation is not discussed.It provides a theoretical reference for time-delay analysis and the feasibility in modeling.
Table 2. Interpretation of the model parameters.
Table 2. Interpretation of the model parameters.
S(t) The Susceptible Nodes
I(t) The infected nodes
R(t) The recovered nodes
LS(t) The low-energy status susceptible nodes
LI(t) The low-energy status infected nodes
LR(t) The low-energy status recovered nodes
N(t) Total number of nodes
Injection rate of new sensor nodes
α 1 Diffusion rate of malware
α 2 Self-healing rate of the infected nodes
α 3 Recovery rate of the infected nodes
α 4 Immune rate of the susceptible nodes
α 5 Immune failure rate of the recovered nodes; the recovered nodes will be re-exposed to the malware and may be infected again
β 1 Low-energy node conversion rate, which is to describe the process of the general nodes dropping to the low-energy nodes
β 2 Charging success rate
bNode deactivation rate
Publisher’s Note: MDPI stays neutral with regard to jurisdictional claims in published maps and institutional affiliations.

Share and Cite

MDPI and ACS Style

Liu, G.; Li, J.; Liang, Z.; Peng, Z. Dynamical Behavior Analysis of a Time-Delay SIRS-L Model in Rechargeable Wireless Sensor Networks. Mathematics 2021, 9, 2007. https://doi.org/10.3390/math9162007

AMA Style

Liu G, Li J, Liang Z, Peng Z. Dynamical Behavior Analysis of a Time-Delay SIRS-L Model in Rechargeable Wireless Sensor Networks. Mathematics. 2021; 9(16):2007. https://doi.org/10.3390/math9162007

Chicago/Turabian Style

Liu, Guiyun, Junqiang Li, Zhongwei Liang, and Zhimin Peng. 2021. "Dynamical Behavior Analysis of a Time-Delay SIRS-L Model in Rechargeable Wireless Sensor Networks" Mathematics 9, no. 16: 2007. https://doi.org/10.3390/math9162007

APA Style

Liu, G., Li, J., Liang, Z., & Peng, Z. (2021). Dynamical Behavior Analysis of a Time-Delay SIRS-L Model in Rechargeable Wireless Sensor Networks. Mathematics, 9(16), 2007. https://doi.org/10.3390/math9162007

Note that from the first issue of 2016, this journal uses article numbers instead of page numbers. See further details here.

Article Metrics

Back to TopTop