Front. Phys. Frontiers in Physics Front. Phys. 2296-424X Frontiers Media S.A. 1098040 10.3389/fphy.2022.1098040 Physics Original Research Dynamic analysis and optimal control of worm propagation model with saturated incidence rate Wang et al. 10.3389/fphy.2022.1098040 Wang Ruiling 1 Xue Yakui 1 * Xue Kailin 2 1 School of Mathematics, North University of China, Taiyuan, China 2 Department of Mathematics, University of South Australia, Adelaide, Australia

Edited by: Xiaoming Zheng, Central Michigan University, United States

Reviewed by: Zizhen Zhang, Anhui University of Finance and Economics, China

Xue Ren, Harbin Engineering University, China

 *Correspondence: Yakui Xue, xyk5152@163.com

This article was submitted to Social Physics, a section of the journal Frontiers in Physics

 10 01 2023 2022 10 1098040 14 11 2022 09 12 2022 Copyright © 2023 Wang, Xue and Xue. 2023 Wang, Xue and Xue

This is an open-access article distributed under the terms of the Creative Commons Attribution License (CC BY). The use, distribution or reproduction in other forums is permitted, provided the original author(s) and the copyright owner(s) are credited and that the original publication in this journal is cited, in accordance with accepted academic practice. No use, distribution or reproduction is permitted which does not comply with these terms.

In order to prevent the propagation of computer worms effectively, based on the latent character of worms, the exposed compartments of computer and USB device are introduced respectively, and a series of computer worm propagation models with saturation incidence rate are proposed. The qualitative behavior of the proposed model is studied. Firstly, the threshold R 0 of the model is derived by using the next-generation matrix method, which completely characterized the stability of disease free equilibrium and endemic equilibrium. If R 0 < 1, the disease free equilibrium is asymptotically stable, implying that the worm dies out eventually and its attack remains under control; if R 0 > 1, the asymptotic stability of endemic equilibrium under certain conditions is proved, which means that the worm is always persistent and uncontrollable under such conditions. Secondly, the theoretical results are verified by numerical study, in which the relative importance of each parameter in worm prevalence is evaluated by sensitivity analysis. Finally, so as to minimize the number of computer and USB device carrying computer worms in short span of time, the worm propagation model is extended to incorporate three control strategies. The Pontryagin’s maximum principle is used to characterize the controls’ optimal levels. According to the control effect diagram, the combined strategy is effective in minimizing the transmission dynamics of worm virus in computer and USB devices populations respectively.

 computer worms saturation incidence rate stability analysis sensitivity analysis optimal control numerical simulation 11971278 National Natural Science Foundation of China10.13039/501100001809 1 Introduction

As a storage medium, USB device is favored by many office workers because of its large storage capacity and portability. It is often used to copy data between different computers. In fact, the transmission of data via USB device may be accompanied by the spread of worms [1]. Computer worm is a common computer virus, which has the characteristics of wide spread and serious harm [2]. On 13 January 2021, the news of incaseformat worm flooded wechat circle of friends, and infection cases were found in many provinces, cities, and industries, with a trend of large-scale outbreak. After execution, the worm will automatically copy to the windows directory of the system disk and create a registry to start itself [3]. Once the user restarts the computer, this will cause the worm host to execute from the windows directory, and then the worm process will traverse all disks except the system disk and delete the files, causing irreparable losses to the user. Therefore, it has become an urgent problem in the field of network security to study the propagation rule of computer worm and effectively control its propagation.

A computer worm is a group of programmed code that can reproduce itself and affect the normal use of the computer [4]. A biological virus is a unique tiny living organism, which can make use of the nutrients of the host cell to copy its own essential constituent material DNA or RNA and proteins [5–7]. Computer worms and biological viruses are two concepts from different fields, but some of their properties have striking similarities [8]. Mainly in the following four aspects: firstly, computer worms and biological viruses are infectious. Computer worms can spread from an infected computer to an uninfected computer through network or other remavable devices. Biological viruses can also be transmitted by direct or indirect contact between living organisms. Then, computer worms and biological viruses are latent and generally not easy to find. Computer worms can integrate their own program fragments into the system programs, so that their most important code can be saved from the anti-virus software. Biological viruses can also integrate their genetic material into the host’s DNA or RNA, thereby escaping attack by the host’s immune system. In addition, computer worms and biological viruses will mutate in the process of spreading, making their species more and more diverse. Finally, both computer worms and biological viruses are destructive. Computer worms can damage the computer system and interfere with the normal operation of the computer. Biological viruses may damage the cells or organs of the host organism, which may pose a threat to the organism’s life. Because of the high similarity between computer worms and biological viruses, we used the method of studying the transmission of infectious diseases to explore the dynamic behavior of computer worms.

In epidemiology, dynamic model has been an important method in analysing the spread and control of infectious diseases qualitatively and quantitatively [9–11]. In 1991, Kephart and White first introduced the dynamic model of infectious disease into the study of the propagation of computer virus. Since then, the dynamic model has also become a vital tool for the study of various computer viruses [12–16]. The research on computer worms have experienced a rapid growth, and a large number of models have been established, mathematically analyzed, which provides some useful and valid references for the characteristics of computer worms transmission [17]. In 2011, Song et al. [18] studied removable device as an independent carrier for the first time. They hypothesized that computer worms spread via both web-based scanning and removable devices, established SIR model of computer population and SI model of removable devices population, studied their dynamic behavior, and gave corresponding control strategy. In 2012, Zhu et al. [19] proposed a new dynamic model to describe the spread of computer virus by using the same method as in Ref. [18]. In addition, by qualitative analysis, it is concluded that controlling R 0 below one is an effective means to extinguish virus, which provides a good start point for understanding the transmission of computer virus through such interactions. In 2015, Ma et al. [20] proposed SIBV model of fixed nodes and SI model of mobile nodes considering the influence of benign worms and mobile devices in the network environment, discussed the influences of removable devices on the interaction dynamics between malicious worms and benign worms. Through numerical analysis and simulation, it is proved that the anti-worm technology can effectively suppress the spread of malicious worms. In 2018, Zhu et al. [21] introduced two control strategies of disconnecting computers from removable devices and reorganizing computers on the basis of Ref. [19], and studied the optimal control of virus transmission between computers and removable devices. Different from previous studies, this paper does not use a fixed cost weight index in the objective function, but uses a state-based cost weight index. In 2020, Kim et al. [1] proposed and analyzed a scheme to control virus transmission through vaccination, that is, installing effective anti-virus software. The results showed that the higher the vaccination rate, the lower the number of infected computers.

Inspired by the models above, we establish a new worm transmission model. We adopt the saturation incidence rate [22, 23], which can present the inhibition effect of susceptible devices and the crowding effect of infectious devices and can also ensure the boundness of contact rate by choosing suitable parameters [24]. In addition, considering the latent character of computer worms, the exposed compartments of computer and USB are introduced [25–28]. The rest of the article is organized as follows: in Section 2, qualitative analysis of the dynamical system is carried out, and sufficient conditions for the existence of equilibria and asymptotic stability are given; numerical findings and discussion are conducted in Section 3; in Section 4, optimal control strategy is proposed; and conclusion is presented in Section 5.

 2 Dynamical model

The main purpose of this paper is to explore the dynamic behavior of computer worms transmission in the process of using USB devices to transmit data, so as to formulate effective strategies to control the spread of worms. This model is an improved version of Kim et al. [1]. In this model, the computer acts as the host and the USB device acts as a vector. Among them, the total computer population is divided into the following four subclasses, which are S c (t), E c (t), I c (t), and R c (t), representing the numbers of susceptible computers, exposed computers, infectious computers and recovered computers at time t. And similarly for USB devices population, which is classified into three vector population subclasses, namely S u (t), E u (t), and I u (t), represent the numbers of susceptible USB devices, exposed USB devices and infectious USB devices at time t, respectively. In addition, the following model assumptions are presented:

(H1) All newly launched computers and removable USB devices are susceptible.

(H2) At time t, the infectious force of infected USB devices to susceptible computers is given by 1 − σ ν β 1 S c I u 1 +   α 1 I u , while the infectibility of infected computers to susceptible USB devices is 1 − σ ν β 1 S u I c 1 +   α 2 I c .

(H3) Infected computers will gain temporary immunity after recovery due to anti-virus software installed, but after a period of time will join again susceptible class because of the variability of computer worms.

(H4) Computer dysfunction is not only related to natural death, but also to worm invasion.

The improved model as follow: d S c d t = Λ 1 − 1 − σ ν β 1 S c I u 1 + α 1 I u − σ ν + μ 1 S c + τ R c , d E c d t = 1 − σ ν β 1 S c I u 1 + α 1 I u − η + ϵ + μ 1 E c , d I c d t = η E c − γ + μ 1 + θ I c , d R c d t = σ ν S c + ϵ E c + γ I c − τ + μ 1 R c , d S u d t = Λ 2 − β 2 S u I c 1 + α 2 I c − μ 2 S u + ϕ E u + p I u , d E u d t = β 2 S u I c 1 + α 2 I c − ξ + ϕ + μ 2 E u , d I u d t = ξ E u − p + μ 2 I u .

For the computer population in system (Eq. 1), Λ1 is the recruitment rate of computer population; σ is the efficiency coefficient of the anti-virus software; ν is installation coverage rate coefficient of anti-virus software; β 1 is infection rate coefficient of susceptible computers by infectious USB devices; α 1 is the degree coefficient of protection measures taken for susceptible computers in case of worms outbreak; μ 1 is the natural elimination rate coefficient of computer population; θ is the elimination rate coefficient of computer population caused by worm invasion; τ, η, ϵ, γ are the transition rate coefficients among states in the computer population.

As for the USB population in system (Eq. 1), Λ2 is the recruitment rate of USB population; β 2 is infection rate coefficient of susceptible USB devices by infectious computers; α 2 is the degree coefficient of protection measures taken for susceptible USB devices in case of worms outbreak; μ 2 is the natural elimination rate coefficient of USB population; ϕ, p, ξ are the transition rate coefficients among states in the USB population.

Suppose N c and N u are the total number of computer population and USB population at time t, respectively, satisfying N c = S c + E c + I c + R c and N u = S u + E u + I u . Add the equations in system (Eq. 1) to get d N c d t = Λ 1 − μ 1 N c − θ I c ≤ Λ 1 − μ 1 N c and d N u d t = Λ 2 − μ 2 N u . The operable region for system (Eq. 1) is given below: Δ = S c , E c , I c , R c , S u , E u , I u ∈ R + 7 | 0 ≤ N c ≤ Λ 1 μ 1 , 0 ≤ N u ≤ Λ 2 μ 2 .

In this section, we firstly prove the existence of disease free equilibrium and endemic equilibrium of system (Eq. 1), and obtain the basic reproduction number R 0 through calculation which determines the propagation dynamics of system (Eq. 1). Secondly, we provide the local and global stability of disease free equilibrium and endemic equilibrium under certain conditions, which establishes the theoretical foundation for the control strategies of worm attack.

2.1 Existence of basic reproduction number and equilibria

Our aim is to study the characteristics of disease free equilibrium and endemic equilibrium. Let the right-hand side of system (Eq. 1) equals zero, and the disease free equilibrium can be obtained: P 0 = S c 0 , E c 0 , I c 0 , R c 0 , S u 0 , E u 0 , I u 0 = Λ 1 τ + μ 1 μ 1 σ ν + τ + μ 1 , 0,0 , Λ 1 σ ν μ 1 σ ν + τ + μ 1 , Λ 2 μ 2 , 0,0 .

Next, we calculate the basic reproduction number R 0, which is essentially the number of secondary infectious cases produced in a susceptible class by a single infectious node in its overall infectious duration. This is a decisive threshold in epidemiological systems that can be used to predict whether computers will be continuously attacked by worms [29]. According to the next-generation matrix method in literature [30], we have the following: F = 0 0 0 1 − σ ν β 1 S c 0 0 0 0 0 0 β 2 S u 0 0 0 0 0 0 0 , V = η + ϵ + μ 1 0 0 0 − η γ + μ 1 + θ 0 0 0 0 ξ + ϕ + μ 2 0 0 0 − ξ p + μ 2 .

Thus, we obtain the basic reproduction number: R 0 = 1 − σ ν β 1 β 2 η ξ S c 0 S u 0 η + ϵ + μ 1 γ + μ 1 + θ ξ + ϕ + μ 2 p + μ 2 = 1 − σ ν β 1 β 2 η ξ τ + μ 1 Λ 1 Λ 2 η + ϵ + μ 1 γ + μ 1 + θ ξ + ϕ + μ 2 p + μ 2 σ ν + τ + μ 1 μ 1 μ 2 .

Theorem 1 When R 0 > 1, system (Eq. 1) has a unique endemic equilibrium.

Proof: Make the right side of system (Eq. 1) equal to zero, and substitute the endemic equilibrium P * = ( S c * , E c * , I c * , R c * , S u * , E u * , I u * ) into the equation: E c ∗ = γ + μ 1 + θ η I c ∗ , R c ∗ = σ ν S c ∗ + ϵ E c ∗ + γ I c ∗ τ + μ 1 , E u * = p + μ 2 ξ I u * , S c ∗ = Λ 1 τ + μ 1 + τ ϵ E c ∗ + γ I c ∗ μ 1 σ ν + τ + μ 1 + 1 − σ ν β 1 I u * τ + μ 1 1 + α 1 I u * , S u * = Λ 2 + p I u * μ 2 + ξ + μ 2 ξ + ϕ + μ 2 ⋅ β 2 I c ∗ 1 + α 2 I c ∗ , I u * = ξ β 2 Λ 2 I c ∗ ξ + ϕ + μ 2 p + μ 2 1 + α 2 I c ∗ μ 2 + μ 2 β 2 p + ξ + μ 2 I c ∗ .

We get the expression G I c * = a I c * + b for I c * , where a = α 2 μ 2 ξ + ϕ + μ 2 p + μ 2 + μ 2 β 2 p + ξ + μ 2 + ξ β 2 Λ 2 α 1 μ 1 η + ϵ + μ 1 γ + μ 1 + θ σ ν + τ + μ 1 + ξ β 2 Λ 2 η β 1 μ 1 γ 1 − σ ν + μ 1 + θ 1 − σ ν τ + μ 1 + ξ β 2 Λ 2 β 1 μ 1 γ + μ 1 + θ 1 − σ ν τ + μ 1 + ϵ > 0 , b = μ 1 μ 2 η + ϵ + μ 1 γ + μ 1 + θ σ ν + τ + μ 1 ξ + ϕ + μ 2 p + μ 2 1 − R 0 2 . When R 0 > 1, b < 0, so I c * = − b a > 0 and G I c * = 0 , that is, when R 0 > 1, there is a unique positive equilibrium P*.

2.2 Stability of disease free equilibrium

In this subsection, we analyze the asymptotical stability of disease free equilibrium P 0.

Theorem 2 The disease free equilibrium P 0 of system (Eq. 1) is locally asymptotically stable if R 0 < 1 and is unstable if R 0 > 1.

Proof: The Jacobian matrix of system (Eq. 1) at P 0 is J 0 = a 11 0 0 τ 0 0 − 1 − σ ν β 1 S c 0 0 a 22 0 0 0 0 1 − σ ν β 1 S c 0 0 η a 33 0 0 0 0 σ ν ϵ γ a 44 0 0 0 0 0 − β 2 S u 0 0 − μ 2 ϕ p 0 0 β 2 S u 0 0 0 a 66 0 0 0 0 0 0 ξ a 77 , where a 11 = − σ ν + μ 1 , a 22 = − η + ϵ + μ 1 , a 33 = − γ + μ 1 + θ , a 44 = − τ + μ 1 , a 66 = − ξ + ϕ + μ 2 , a 77 = − p + μ 2 . Obviously, one of the eigenvalues of J |0| is λ 1 = −μ 2. In order to obtain the rest eigenvalues, the following simplified matrix is considered: J 1 0 = a 11 0 0 τ 0 − 1 − σ ν β 1 S c 0 0 a 22 0 0 0 1 − σ ν β 1 S c 0 0 η a 33 0 0 0 σ ν ϵ γ a 44 0 0 0 0 β 2 S u 0 0 a 66 0 0 0 0 0 ξ a 77 .

By elementary row operations, we get J 2 0 = a 11 0 0 τ 0 − 1 − σ ν β 1 S c 0 0 a 22 0 0 0 1 − σ ν β 1 S c 0 0 0 b 33 0 0 1 − σ ν η β 1 S c 0 0 0 0 b 44 0 b 46 0 0 0 0 b 55 1 − σ ν η β 1 β 2 S c 0 S u 0 0 0 0 0 0 b 66 ,

where b 33 = − η + ϵ + μ 1 γ + μ 1 + θ , b 44 = − μ 1 τ + σ ν + μ 1 η + ϵ + μ 1 2 γ + μ 1 + θ , b 55 = − η + ϵ + μ 1 γ + μ 1 + θ ξ + ϕ + μ 2 , b 66 = 1 − σ ν ξ η β 1 β 2 S c 0 S u 0 η + ϵ + μ 1 γ + μ 1 + θ ξ + ϕ + μ 2 − p + μ 2 , b 46 = 1 − σ ν η β 1 S c 0 γ σ ν + μ 1 η + ϵ + μ 1 + ϵ μ 1 − σ ν η + μ 1 1 − σ ν β 1 S c 0 η + ϵ + μ 1 γ + μ 1 + θ .

When ( 1 − σ ν ) ξ η β 1 β 2 S c 0 S u 0 ( η + ϵ + μ 1 ) ( γ + μ 1 + θ ) ( ξ + ϕ + μ 2 ) < ( p + μ 2 ) , that is, when R 0 < 1, the eigenvalues of J 2 | 0 | are all negative, and therefore the eigenvalues of J |0| are all negative. It shows that when R 0 < 1, the disease free equilibrium P 0 of system (Eq. 1) is locally asymptotically stable in Δ, whereas P 0 is unstable in Δ if R 0 > 1.

Theorem 3 The disease free equilibrium P 0 of system (Eq. 1) is globally asymptotically stable in region Δ if R 0 < 1 and is unstable if R 0 > 1.

Proof: The Castillo-Chavez method in Refs. [31] is used to prove it. First, let χ 1 = (S c , R c , S u ), χ 2 = (E c , I c , E u , I u ) and define P 0 = ( χ 1 0 , 0 ) , where χ 1 0 = Λ 1 τ + μ 1 μ 1 σ ν + τ + μ 1 , Λ 1 σ ν μ 1 σ ν + τ + μ 1 , Λ 2 μ 2 .

We decompose system (Eq. 1) into two subsystems, d χ 1 d t = G χ 1 , χ 2 , d χ 2 d t = H χ 1 , χ 2 ,

where G χ 1 , χ 2 = Λ 1 − 1 − σ ν β 1 S c I u 1 + α 1 I u − σ ν + μ 1 S c + τ R c σ ν S c + E c + γ I c − τ + μ 1 R c Λ 2 − β 2 S u I c 1 + α 2 I c − μ 2 S u + ϕ E u + p I u , H χ 1 , χ 2 = 1 − σ ν β 1 S c I u 1 + α 1 I u − η + + μ 1 E c η E c − γ + μ 1 + θ I c β 2 S u I c 1 + α 2 I c − ξ + ϕ + μ 2 E u ξ E u − p + μ 2 I u . Because S c = S c 0 , R c = R c 0 , S u = S u 0 , and G χ 1 , 0 = 0 , we have G χ 1 , 0 = Λ 1 − σ ν + μ 1 S c 0 + τ R c 0 σ ν S c 0 − τ + μ 1 R c 0 Λ 2 − μ 2 S u 0 = 0 . According to (Eq. 10), when t → ∞, χ 1 → χ 1 0 . Therefore, χ 1 = χ 1 0 is globally asymptotically stable.

Moreover, it is obvious that H χ 1 , χ 2 = B χ 2 − H ̄ χ 1 , χ 2 = a 22 0 0 1 − σ ν β 1 S c 0 η a 33 0 0 0 β 2 S u 0 a 66 0 0 0 ξ a 77 E c I c E u I u − 1 − σ ν β 1 S c 0 I u − 1 − σ ν β 1 S c I u 1 + α 1 I u 0 β 2 S u 0 I c − β 2 S u I c 1 + α 2 I c 0 .

The population limits of computer and USB device are S c 0 , S u 0 , that is, S c ≤ S c 0 , S u ≤ S u 0 , so 1 − σ ν β 1 S c 0 I u ≥ 1 − σ ν β 1 S c I u 1 + α 1 I u , β 2 S u 0 I c ≥ β 2 S u I c 1 + α 2 I c , which means H ̄ χ 1 , χ 2 ≥ 0 . Also from (Eq. 11), B is an M-matrix. Thus condition 1 and 2 are satisfied in section 4.1 of Refs. [32]. Therefore, Theorem 3 is proved.

2.3 Stability of endemic equilibrium

In this subsection, we analyze the asymptotical stability of endemic equilibrium P*.

Theorem 4 The endemic equilibrium P* of system Eq. 1 is locally asymptotically stable if max ϵ , γ < σ ν and 1 < R 0 < 1 + α 1 I u * 1 + α 2 I c * τ + μ 1 Λ 1 Λ 2 σ ν + τ + μ 1 μ 1 μ 2 S c * S u * .

Proof: The Jacobian matrix of system (Eq. 1) at P* is given by J | * | = c 11 0 0 τ 0 0 − 1 − σ ν β 1 S c ⋆ 1 + α 1 I u ∗ 2 1 − σ ν β 1 I u ∗ 1 + α 1 I u ∗ a 22 0 0 0 0 1 − σ ν β 1 S c ∗ 1 + α 1 I u ∗ 2 0 η a 33 0 0 0 0 σ ν ϵ γ a 44 0 0 0 0 0 − β 2 S u ∗ 1 + α 2 I c ∗ 2 0 c 55 ϕ p 0 0 β 2 S u ∗ 1 + α 2 I c ∗ 2 0 β 2 I c ∗ 1 + α 2 I c ∗ a 66 0 0 0 0 0 0 ξ a 77 , where c 11 = − σ ν + μ 1 − 1 − σ ν β 1 I u * 1 + α 1 I u * , c 55 = − μ 2 − β 2 I c * 1 + α 2 I c * .

After simplification we get J 1 * = c 11 0 0 τ 0 0 − 1 − σ ν β 1 S c ∗ 1 + α 1 I u ∗ 2 0 d 22 0 1 − σ ν β 1 I u ∗ τ 1 + α 1 I u ∗ 0 0 1 − σ ν β 1 S c ∗ σ ν + μ 1 1 + α 1 I u ∗ 2 0 0 d 33 1 − σ ν β 1 I u ∗ η τ 1 + α 1 I u ∗ 0 0 1 − σ ν β 1 S c ∗ η σ ν + μ 1 1 + α 1 I u ∗ 2 0 0 0 d 44 0 0 d 47 0 0 0 0 − μ 2 − ξ + μ 2 p 0 0 0 0 0 d 66 d 67 0 0 0 0 0 0 d 77 ,

where d 22 = − a 22 c 11 , d 33 = a 22 a 33 c 11 , d 44 = μ 1 τ + σ ν + μ 1 a 22 2 a 33 c 11 2 + μ 1 1 − σ ν β 1 I u ∗ 1 + α 1 I u ∗ a 22 2 a 33 c 11 2 + η μ 1 + θ − μ 1 a 33 τ 1 − σ ν β 1 I u ∗ 1 + α 1 I u ∗ a 22 c 11 2 , d 47 = σ ν η θ + ϵ − σ ν μ 1 a 33 − η μ 1 γ − σ ν 1 − σ ν β 1 S c ∗ 1 + α 1 I u ∗ 2 a 22 c 11 2 , d 66 = − β 2 I c ∗ 1 + α 2 I c ∗ ξ + μ 2 + μ 2 a 66 a 22 a 33 c 11 d 44 , d 67 = β 2 I c ∗ p 1 + α 2 I c ∗ a 22 a 33 c 11 d 44 + μ 1 μ 2 η τ + σ ν + μ 1 × 1 − σ ν β 1 S c ∗ 1 + α 1 I u ∗ 2 β 2 S u ∗ 1 + α 2 I c ∗ 2 a 11 a 22 2 a 33 c 11 2 + μ 1 2 μ 2 η τ a 33 + μ 1 μ 2 η τ ϵ a 33 + μ 1 μ 2 η a 11 a 22 a 33 − μ 1 μ 2 η 2 τ μ 1 +   θ − μ 1 μ 2 η 2 τ γ ⋅ 1 − σ ν β 1 I u ∗ 1 + α 1 I u ∗ β 2 S u ∗ 1 + α 2 I c ∗ 2 1 − σ ν β 1 S c ∗ 1 + α 1 I u ∗ 2 a 22 c 11 2 , d 77 = − μ 2 p + ξ + μ 2 β 2 I c ∗ 1 + α 2 I c ∗ a 22 a 33 c 11 d 44 + 1 − σ ν β 1 I u ∗ 1 + α 1 I u ∗ μ 2 a 22 a 33 a 66 a 77 d 44 + ξ μ 2 η τ σ ν η θ + ϵ − σ ν μ 1 a 33 1 − σ ν β 1 I u ∗ 1 + α 1 I u ∗ β 2 S u ∗ 1 + α 2 I c ∗ 2 1 − σ ν β 1 S c ∗ 1 + α 1 I u ∗ 2 a 22 c 11 2 + ξ μ 2 η τ − η μ 1 γ − σ ν 1 − σ ν β 1 I u ∗ 1 + α 1 I u ∗ β 2 S u ∗ 1 + α 2 I c ∗ 2 1 − σ ν β 1 S c ∗ 1 + α 1 I u ∗ 2 a 22 c 11 2 + a 11 a 44 d 44 R 0 2 1 − σ ν β 1 β 2 η ξ Λ 1 Λ 2 μ 1 σ ν + τ + μ 1 1 − σ ν + τ + μ 1 μ 1 μ 2 S c ∗ S u ∗ R 0 2 τ + μ 1 Λ 1 Λ 2 1 + α 1 I u ∗ 2 1 + α 2 I c ∗ 2 . + a 11 a 44 d 44 R 0 2 1 − σ ν β 1 β 2 η ξ Λ 1 Λ 2 μ 1 σ ν + τ + μ 1 1 − σ ν + τ + μ 1 μ 1 μ 2 S c ∗ S u ∗ R 0 2 τ + μ 1 Λ 1 Λ 2 1 + α 1 I u ∗ 2 1 + α 2 I c ∗ 2 .

It is obvious that a 22, a 33, a 66, c 11, c 55 < 0, and then, d 22, d 33, d 44, d 66 < 0 can be obtained according to the above expressions. Besides, when max ϵ , γ < σ ν and 1 < R 0 < 1 + α 1 I u ∗ 1 + α 2 I c ∗ τ + μ 1 Λ 1 Λ 2 σ ν + τ + μ 1 μ 1 μ 2 S c ∗ S u ∗ , we get d 77 < 0. Just as we have shown that the real parts of all eigenvalues of J ∗ are negative, so the endemic equilibrium J ∗ is locally asymptotically stable under condition max ϵ , γ < σ ν and 1 < R 0 < 1 + α 1 I u ∗ 1 + α 2 I c ∗ τ + μ 1 Λ 1 Λ 2 σ ν + τ + μ 1 μ 1 μ 2 S c ∗ S u ∗ .

Next, the geometrical approach in Refs. [33] is used to determine the condition that system (Eq. 1) achieves global asymptotic stability at the endemic equilibrium.

Theorem 5 If max ϕ S c ∗ E c ∗ − β 2 I c ∗ 1 + α 2 I c ∗ , β 2 I c ∗ 1 + α 2 I c ∗ ⋅ E c ∗ S c ∗ < 1 − σ ν β 1 I u ∗ 1 + α 1 I u ∗ < β 2 I c ∗ 1 + α 2 I c ∗ ⋅ S u ∗ E u ∗ and R 0 > 1, then system (Eq. 1) is globally asymptotically stable at endemic equilibrium P∗ and unstable otherwise.

Proof: We prove the global stability of endemic equilibrium P∗ under the above conditions by using geometrical approach [32, 34, 35], which is the generalization of Lyapunov theory in essence, by using the third additive compound matrix. Consider the subsystems of system (Eq. 1) d S c d t = Λ 1 − 1 − σ ν β 1 S c I u 1 + α 1 I u − σ ν + μ 1 S c + τ R c , d E c d t = 1 − σ ν β 1 S c I u 1 + α 1 I u − η + ϵ + μ 1 E c , d S u d t = Λ 2 − β 2 S u I c 1 + α 2 I c − μ 2 S u + ϕ E u + p I u , d E u d t = β 2 S u I c 1 + α 2 I c − ξ + ϕ + μ 2 E u .

Let J be the Jacobian matrix of system (Eq. 14) given by J = − σ ν + μ 1 − 1 − σ ν β 1 I u ∗ 1 + α 1 I u ∗ 0 0 0 1 − σ ν β 1 I u ∗ 1 + α 1 I u ∗ − η + ϵ + μ 1 0 0 0 0 − μ 2 − β 2 I c ∗ 1 + α 2 I c ∗ ϕ 0 0 β 2 I c ∗ 1 + α 2 I c ∗ − ξ + ϕ + μ 2 .

Based on [35], the third additive compound matrix of J is denoted by J 3 = A 11 ϕ 0 0 β 2 I c ∗ 1 + α 2 I c ∗ A 22 0 0 0 0 A 33 0 0 0 1 − σ ν β 1 I u ∗ 1 + α 1 I u ∗ A 44 ,

where A 11 = − σ ν + 2 μ 1 + μ 2 + η + ϵ − 1 − σ ν β 1 I u ∗ 1 + α 1 I u ∗ − β 2 I c ∗ 1 + α 2 I c ∗ , A 22 = − σ ν + 2 μ 1 + μ 2 + η + ϵ + ξ + ϕ − 1 − σ ν β 1 I u ∗ 1 + α 1 I u ∗ , A 33 = − σ ν + μ 1 + 2 μ 2 + ξ + ϕ − 1 − σ ν β 1 I u ∗ 1 + α 1 I u ∗ − β 2 I c ∗ 1 + α 2 I c ∗ , A 44 = − μ 1 + 2 μ 2 + η + ϵ + ξ + ϕ − β 2 I c ∗ 1 + α 2 I c ∗ .

Now let us choose a function P χ = P S c , E c , S u , E u so that P χ = diag S c , E c , S u , E u , and accordingly, P − 1 χ = diag 1 S c , 1 E c , 1 S u , 1 E u . The time derivative of P χ yields that P f χ = diag S ̇ c , E ̇ c , S ̇ u , E ̇ u .

Therefore, P f P − 1 = diag S ̇ c S c , E ̇ c E c , S ̇ u S u , E ̇ u E u , and P J 3 P − 1 = A 11 ϕ S c ∗ E c ∗ 0 0 β 2 I c ∗ 1 + α 2 I c ∗ ⋅ E c ∗ S c ∗ A 22 0 0 0 0 A 33 0 0 0 1 − σ ν β 1 I u ∗ 1 + α 1 I u ∗ ⋅ E u ∗ S u ∗ A 44 . So that B = P f P − 1 + P J 3 P − 1 , B = B 11 B 12 0 0 B 21 B 22 0 0 0 0 B 33 0 0 0 B 43 B 44 , where B 11 = S ̇ c S c − σ ν + 2 μ 1 + μ 2 + η + ϵ − 1 − σ ν β 1 I u ∗ 1 + α 1 I u ∗ − β 2 I c ∗ 1 + α 2 I c ∗ , B 12 = ϕ S c ∗ E c ∗ , B 22 = E ̇ c E c − σ ν + 2 μ 1 + μ 2 + η + ϵ + ξ + ϕ − 1 − σ ν β 1 I u ∗ 1 + α 1 I u ∗ , B 21 = β 2 I c ∗ 1 + α 2 I c ∗ ⋅ E c ∗ S c ∗ , B 33 = S ̇ u S u − σ ν + μ 1 + 2 μ 2 + ξ + ϕ − 1 − σ ν β 1 I u ∗ 1 + α 1 I u ∗ − β 2 I c ∗ 1 + α 2 I c ∗ , B 44 = E ̇ u E u − μ 1 + 2 μ 2 + η + ϵ + ξ + ϕ − β 2 I c ∗ 1 + α 2 I c ∗ , B 43 = 1 − σ ν β 1 I u ∗ 1 + α 1 I u ∗ ⋅ E u ∗ S u ∗ .

Consequently, if 1 − σ ν β 1 I u ∗ 1 + α 1 I u ∗ + β 2 I c ∗ 1 + α 2 I c ∗ > ϕ S c ∗ E c ∗ , then h 1 t = B 11 + ∑ j = 2 4 B 1 j , h 1 t = S ̇ c S c − σ ν + 2 μ 1 + μ 2 + η + ϵ − 1 − σ ν β 1 I u ∗ 1 + α 1 I u ∗ − β 2 I c ∗ 1 + α 2 I c ∗ + ϕ S c ∗ E c ∗ = S ̇ c S c − σ ν + 2 μ 1 + μ 2 + η + ϵ − 1 − σ ν β 1 I u ∗ 1 + α 1 I u ∗ + β 2 I c ∗ 1 + α 2 I c ∗ − ϕ S c ∗ E c ∗ ≤ S ̇ c S c − σ ν + 2 μ 1 + μ 2 + η + ϵ . and if 1 − σ ν β 1 I u ∗ 1 + α 1 I u ∗ > β 2 I c ∗ 1 + α 2 I c ∗ ⋅ E c ∗ S c ∗ , then h 2 t = B 22 + ∑ j = 1 j ≠ 2 4 B 2 j , h 2 t = E ̇ c E c − σ ν + 2 μ 1 + μ 2 + η + ϵ + ξ + ϕ − 1 − σ ν β 1 I u ∗ 1 + α 1 I u ∗ + β 2 I c ∗ 1 + α 2 I c ∗ ⋅ E c ∗ S c ∗ = E ̇ c E c − σ ν + 2 μ 1 + μ 2 + η + ϵ + ξ + ϕ − 1 − σ ν β 1 I u ∗ 1 + α 1 I u ∗ − β 2 I c ∗ 1 + α 2 I c ∗ ⋅ E c ∗ S c ∗ ≤ E ̇ c E c − σ ν + 2 μ 1 + μ 2 + η + ϵ + ξ + ϕ .

Similarly, h 3 t = B 33 + ∑ j = 1 j ≠ 3 4 B 3 j , h 3 t = S ̇ u S u − σ ν + μ 1 + 2 μ 2 + ξ + ϕ − 1 − σ ν β 1 I u ∗ 1 + α 1 I u ∗ − β 2 I c ∗ 1 + α 2 I c ∗ , ≤ S ̇ u S u − σ ν + μ 1 + 2 μ 2 + ξ + ϕ . and if β 2 I c ∗ 1 + α 2 I c ∗ > 1 − σ ν β 1 I u ∗ 1 + α 1 I u ∗ ⋅ E u ∗ S u ∗ , then h 4 t = B 44 + ∑ j = 1 3 B 4 j , h 4 t = E ̇ u E u − μ 1 + 2 μ 2 + η + ϵ + ξ + ϕ − β 2 I c ∗ 1 + α 2 I c ∗ + 1 − σ ν β 1 I u ∗ 1 + α 1 I u ∗ ⋅ E u ∗ S u ∗ = E ̇ u E u − μ 1 + 2 μ 2 + η + ϵ + ξ + ϕ − β 2 I c ∗ 1 + α 2 I c ∗ − 1 − σ ν β 1 I u ∗ 1 + α 1 I u ∗ ⋅ E u ∗ S u ∗ ≤ E ̇ u E u − μ 1 + 2 μ 2 + η + ϵ + ξ + ϕ .

If h 1 , h 2 , h 3 , h 4 ∈ R 4 and we define the Lozinski measure μ B as μ B = h i , i = 1,2,3,4 then integration of μ B by taking limits as t → ∞ the following equations follow: lim t → ∞ sup 1 t ∫ 0 t h 1 t d t ≤ 1 t log S c t S c 0 − σ ν + 2 μ 1 + μ 2 + η + ϵ < − σ ν + 2 μ 1 + μ 2 + η + ϵ , lim t → ∞ sup 1 t ∫ 0 t h 2 t d t ≤ 1 t log E c t E c 0 − σ ν + 2 μ 1 + μ 2 + η + ϵ + ξ + ϕ < − σ ν + 2 μ 1 + μ 2 + η + ϵ + ξ + ϕ , lim t → ∞ sup 1 t ∫ 0 t h 3 t d t ≤ 1 t log S u t S u 0 − σ ν + μ 1 + 2 μ 2 + ξ + ϕ < − σ ν + μ 1 + 2 μ 2 + ξ + ϕ , lim t → ∞ sup 1 t ∫ 0 t h 4 t d t ≤ 1 t log E u t E u 0 − μ 1 + 2 μ 2 + η + ϵ + ξ + ϕ < − μ 1 + 2 μ 2 + η + ϵ + ξ + ϕ .

Thus, combining the above four inequalities, we get the following inequality: q ̄ = lim t → ∞ sup 1 t ∫ 0 t μ B d t < 0 .

This shows that the subsystem of system (Eq. 1) containing four non-linear differential equations is globally asymptotically stable around the equilibrium S c ∗ , E c ∗ , S u ∗ , E u ∗ . The solutions to the rest of the equations of system (Eq. 1) are I c → I c ∗ , R c → R c ∗ , I u → I u ∗ when t → ∞. Therefore, the endemic equilibrium P∗ is globally asymptotically stable.

3 Numerical simulation

In this section, we refer to the data of Refs. [36–38] to analyze the performance and dynamic behavior of the studied model through numerical simulation, and demonstrate the feasibility of the theoretical results obtained.

3.1 Stability of equilibria

The stability of disease free equilibrium and endemic equilibrium of system (Eq. 1) is verified by numerical simulation. Figure 1 and Figure 2 show the stability simulation diagram of disease free equilibrium and endemic equilibrium. Initially the number of susceptible computers, exposed computers, infected computers, recovery computers, susceptible USB devices, exposed USB devices and infected USB devices are S c (0) = 50,000, E c (0) = 20,000, I c (0) = 10,000, R c (0) = 10,000, S u (0) = 5,000, E u (0) = 3,000, and I u (0) = 1,000, separately. We take the parameters Λ1 = 0.75, Λ2 = 0.1, σ = 0.6, ν = 0.1, β 1 = 0.035, β 2 = 0.035, α 1 = 0.8, α 2 = 0.3, μ 1 = 0.1, μ 2 = 0.1, τ = 0.1, η = 0.45, ϵ = 0.25, ϕ = 0.05, p = 0.003, ξ = 0.005, θ = 0.001, γ = 0.05, at which time, R 0 ≈ 0.0039 < 1. It can be seen from Figure 1 that the number of exposed computers E c , infected computers I c , exposed USB devices E u and infected USB devices I u gradually approach zero, which shows the correctness of Theorem 2 and Theorem 3. When β 1 = 0.053, μ 1 = 0.01, ξ = 0.05, and other parameters remain unchanged, R 0 ≈ 1.0670 > 1. Compared with Figure 1, the number of exposed computers E c , infected computers I c , exposed USB devices E u and infected USB devices I u in Figure 2 are not zero and reach a stable value, which is also consistent with Theorems 4 and 5.

Time series plot of system (1) with Λ1 = 0.75, Λ2 = 0.1, σ = 0.6, ν = 0.1, β 1 = 0.035, β 2 = 0.035, α 1 = 0.8, α 2 = 0.3, μ 1 = 0.1, μ 2 = 0.1, τ = 0.1, η = 0.45, ϵ = 0.25, ϕ = 0.05, p = 0.003, ξ = 0.005, θ = 0.001, γ = 0.05. (A) the tendency of the worm propagation in a short period, (B) the tendency of the worm propagation in a later period.

Time series plot of system (1) with Λ1 = 0.75, Λ2 = 0.1, σ = 0.6, ν = 0.1, β 1 = 0.053, β 2 = 0.035, α 1 = 0.8, α 2 = 0.3, μ 1 = 0.01, μ 2 = 0.1, τ = 0.1, η = 0.45, ϵ = 0.25, ϕ = 0.05, p = 0.003, ξ = 0.05, θ = 0.001, γ = 0.05. (A) the tendency of the worm propagation in a short period, (B) the tendency of the worm propagation in a later period.

 3.2 Performance comparison of the SEIRSEI model with the SIRSI model

To evaluate the performance of the SEIRSEI model, numerical methods are used to compare it with the SIRSI model [1]. We simulate SEIRSEI and SIRSI models by Matlab. Both models share the same parameters as in Section 3.1. Figure 3 shows the number of infected computers and infected USB devices for both models. The results show that, compared with SIRSI model, the propagation speed of worms and the number of infected computers and USB devices in SEIRSEI model are significantly decreased. This implies that the SEIRSEI model takes much smaller time to combat the worm than the SIRSI model. Obviously, the results shown in Figure 3 verify that the proposed model is more effective in controlling computer worms attack than the previous model.

Performance comparison of the SEIRSEI model with the SIRSI model. (A) graph of infectious number of computers, (B) graph of infectious number of USB devices.

 3.3 Sensitivity analysis

Sensitivity analysis is used to determine which parameters have a significant effect on decreasing the spread of disease. Forward sensitivity analysis is considered an important component of epidemic modeling, although it is tedious to calculate for complex biological models. The ecologist and epidemiologist paid a lot of attention to the sensitivity study of R 0.

Definition. The normalized forward sensitivity index of R 0 that depends differentiability on a parameter ω is defined as S ω = ω R 0 ∂ R 0 ∂ ω

The following three methods are normally used to calculate the sensitivity indices, (i) by direct differentiation, (ii) by a Latin hypercube sampling method, (iii) by linearizing system (Eq. 1), and then solving the obtain set of linear algebraic equations. We will apply the direct differentiation method as it gives analytical expressions for the indices. The indices not only shows us the influence of various aspects associated with the spreading of infectious disease but also gives us important information regarding the comparative change between R 0 and different parameters. Consequently, it helps in developing the control strategies.

Table 1 demonstrates that the parameters Λ1, Λ2, β 1, β 2, τ, η and ξ have a positive influence on the reproduction number R 0, which describe that the growth or decay of these parameters say by 10% will increase or decrease the reproduction number by 10%, 10%, 10%, 10%, 1.1%, 4.3% and 9.6% respectively. But on the other hand, the index for parameters σ, ν, μ 1, μ 2, ϵ, ϕ, p, θ and γ illustrates that increasing their values by 10% will decrease the values of reproduction number R 0 by 2.9%, 2.9%, 16.7%, 26.1%, 3.1%, 3.2%, 0.29%, 0.06%, and 3.3% respectively. α 1 and α 2 have no impact on the reproduction number R 0. Figures 4A–H are plotted to demonstrate the variations in the basic reproductive number R 0 with respect to the different model parameters.

Sensitivity indices of the basic reproduction number R 0 against mentioned parameters.

Parameter S.Index Value Parameter S.Index Value
Λ1 S Λ 1 1 Λ2 S Λ 2 1
σ S σ −0.2946 ν S ν −0.2946
β 1 S β 1 1 β 2 S β 2 1
α 1 S α 1 0 α 2 S α 2 0
μ 1 S μ 1 −1.6719 μ 2 S μ 2 −2.6160
τ S τ 0.1154 η S η 0.4375
ϵ S ϵ −0.3125 ϕ S ϕ −0.3226
p S p −0.0291 ξ S ξ 0.9677
θ S θ −0.0066 γ S γ −0.3311

Sensitivity analysis of different parameters. (A) R 0 versus sensitive parameters p and ϕ, (B) R 0 versus sensitive parameters γ and ϵ, (C) R 0 versus sensitive parameters σ and η, (D) R 0 versus sensitive parameters ν and η, (E) R 0 versus sensitive parameters τ and μ 1, (F) R 0 versus sensitive parameters θ and μ 1, (G) R 0 versus sensitive parameters τ and η, (H) R 0 versus sensitive parameters τ and γ.

 4 Optimal control strategy formulation

In order to explore how to effectively control the propagation of computer worms, system (Eq. 1) is improved in this section. Three control variables u i t i = 1,2,3 are introduced into system (Eq. 1), which denote the effectiveness of the corresponding control strategy to inhibit worm propagation. u 1 t 0 ≤ u 1 t ≤ 1 represents the effectiveness of improving the installation coverage rate of anti-virus software, u 2 t 0 ≤ u 2 t ≤ 1 and u 3 t 0 ≤ u 3 t ≤ 1 represent the effectiveness of disinfecting computers and USB devices, respectively. Establish the optimal control system: d S c d t = Λ 1 − 1 − u 1 β 1 S c I u 1 + α 1 I u − u 1 + μ 1 S c + τ R c , d E c d t = 1 − u 1 β 1 S c I u 1 + α 1 I u − η + u 2 + μ 1 E c , d I c d t = η E c − u 2 + μ 1 + θ I c , d R c d t = u 1 S c + u 2 E c + u 2 I c − τ + μ 1 R c , d S u d t = Λ 2 − β 2 S u I c 1 + α 2 I c − μ 2 S u + u 3 E u + u 3 I u , d E u d t = β 2 S u I c 1 + α 2 I c − ξ + u 3 + μ 2 E u , d I u d t = ξ E u − u 3 + μ 2 I u .

Defining the objective function J u 1 , u 2 , u 3 = ∫ 0 T A 1 E c + A 2 I c + A 3 E u + A 4 I u + 1 2 B 1 u 1 2 + B 2 u 2 2 + B 3 u 3 2 d t , where A 1, A 2, A 3 and A 4 are positive weight constants of E c , I c , E u and I u respectively, and T is the final time to implement the control strategies.

The control set is defined as U = u 1 t , u 2 t , u 3 t u i t i = 1,2,3 , which is a Lebesgue measurable control quantity satisfying u i t ∈ 0,1 and t ∈ 0 , T . The purpose of this paper is to find the optimal control u 1 ∗ , u 2 ∗ , u 3 ∗ so that J u 1 ∗ t , u 2 ∗ t , u 3 ∗ t = min J u 1 t , u 2 t , u 3 t u 1 t , u 2 t , u 3 t ∈ U .

Based on the boundedness of the right end of system (Eq. 21) and the convexity of the integrand of (Eq. 22), the existence of the optimal control u 1 ∗ , u 2 ∗ , u 3 ∗ can be obtained [39, 40].

To find the optimal solution of the optimal control problem (Eqs 21, 22), the Hamilton function is defined as follows H = A 1 E c + A 2 I c + A 3 E u + A 4 I u + 1 2 B 1 u 1 2 + B 2 u 2 2 + B 3 u 3 2 + λ S c d S c d t + λ E c d E c d t + λ I c d I c d t + λ R c d R c d t + λ S u d S u d t + λ E u d E u d t + λ I u d I u d t , where λ S c , λ E c , λ I c , λ R c , λ S u , λ E u , λ I u are the adjoint variables. According to Pontryagin’s maximum principle [41], the following theorem is obtained.

Theorem 6 Let u 1 ∗ t , u 2 ∗ t , u 3 ∗ t be optimal controls and S c ∗ , E c ∗ , I c ∗ , R c ∗ , S u ∗ , E u ∗ , I u ∗ be optimal state solutions of the optimal control problem (Eqs 21, 22), and then there exists adjoint variables λ S c , λ E c , λ I c , λ R c , λ S u , λ E u , λ I u , satisfying, d λ i d t = − ∂ H ∂ i , where i = S c , E c , I c , R c , S u , E u , I u . In addition, the corresponding optimal control variables u 1 ∗ t , u 2 ∗ t , u 3 ∗ t are u 1 ∗ = max 0 , min β 1 S c I u 1 + α 1 I u λ E c − λ S c − S c λ R c − λ S c B 1 , 1 , u 2 ∗ = max 0 , min − E c λ R c − λ E c − I c λ R c − λ I c B 2 , 1 , u 3 ∗ = max 0 , min E u λ E u − λ S u − I u λ I u − λ S u B 3 , 1 .

Proof: To prove the theorem, the differential of (Eq. 23), namely the Hamiltonian function H, with respect to the state variables S c , E c , I c , R c , S u , E u , I u are obtained, and the corresponding adjoint system is as follows: d λ S c d t       =       − ∂ H ∂ S c = − 1 − u 1 β 1 I u 1 + α 1 I u λ E c − λ S c − u 1 λ R c − λ S c + μ 1 λ S c , d λ E c d t       =       − ∂ H ∂ E c = − A 1 − η λ I c − λ E c − u 2 λ R c − λ E c + μ 1 λ E c , d λ I c d t       =       − ∂ H ∂ I c = − A 2 − β 2 S u 1 + α 2 I c 2 λ E u − λ S u − u 2 λ R c − λ I c + μ 1 + θ λ I c , d λ R c d t       =       − ∂ H ∂ R c = τ λ R c − λ S c + μ 1 λ R c , d λ S u d t       =       − ∂ H ∂ S u = − β 2 I c 1 + α 2 I c λ E u − λ S u + μ 2 λ S u , d λ E u d t       =       − ∂ H ∂ E u = − A 3 + u 3 λ E u − λ S u − ξ λ I u − λ E u + μ 2 λ E u , d λ I u d t       =       − ∂ H ∂ I u = − A 4 − 1 − u 1 β 1 S c 1 + α 1 I u 2 λ E c − λ S c + u 3 λ I u − λ E u + μ 2 λ I u , and with transversality conditions λ S c T = 0 , λ E c T = 0 , λ I c T = 0 , λ R c T = 0 , λ S u T = 0 , λ E u T = 0 , λ I u T = 0 . Also, the optimal controls u 1 ∗ , u 2 ∗ and u 3 ∗ satisfy ∂ H ∂ u i ∗ = 0 , i = 1, 2, 3. Therefore, the optimal controls are characterized by u 1 ∗ = max 0 , min β 1 S c I u 1 + α 1 I u λ E c − λ S c − S c λ R c − λ S c B 1 , 1 , u 2 ∗ = max 0 , min − E c λ R c − λ E c − I c λ R c − λ I c B 2 , 1 , u 3 ∗ = max 0 , min E u λ E u − λ S u − I u λ I u − λ S u B 3 , 1 .

For sufficiently small final time T, the uniqueness of the optimal control of the system has been achieved owing to the boundedness of the state variables and adjoint variables together with the Lipschitz property of systems (Eq. 21) and (Eq. 26).

In order to show the rationality and effectiveness of the control strategy mentioned in Section 4, numerical simulation is carried out. Figure 5 shows a substantial decrease in the population of exposed computers, infected computers, exposed USB devices and infected USB devices respectively relative to the case of no control, when the combined strategy is implemented, including the installation coverage rate of anti-virus software and the disinfection of computer and USB device. Parameter Λ1 = 0.75, Λ2 = 0.1, σ = 0.6, ν = 0.1, β 1 = 0.035, β 2 = 0.035, α 1 = 0.8, α 2 = 0.3, μ 1 = 0.1, μ 2 = 0.1, τ = 0.1, η = 0.45, ϵ = 0.25, ϕ = 0.05, p = 0.003, ξ = 0.005, θ = 0.001, γ = 0.05. We draw the conclusion that the spread of worm virus can be effectively controlled by strengthening users’ computer security education, improving their understanding of anti-virus software, and calling on users to disinfect computers and USB devices regularly.

Plots of different classes showing the difference between with control and without control, respectively: (A) exposed computers, (B) infected computers, (C) exposed USB devices, (D) infected USB devices.

 5 Conclusion and discussion

This paper mainly considers the dynamic behavior of computer worms in the process of using USB to transmit data. The novel ideas of our paper are as follows: (a) the exposed compartments of computer and USB are introduced because of the latent nature of worms; (b) infected computers will gain temporary immunity after recovery due to anti-virus software installed, but after a period of time will join again susceptible class because of the variability of computer worms; (c) the elimination rate of computers is not only related to natural mortality, but also to worm attacks; (d) the saturation incidence is used thanks to the inhibition effect of susceptible equipments and the crowding effect of infectious equipments.

Firstly, through theoretical analysis, the existence of disease free equilibrium and endemic equilibrium is studied, and the sufficient conditions for their asymptotic stability are given. The global stability of disease free equilibrium is carried out using the Castillo-Chavez method, whereas, the global stability of the unique endemic equilibrium is also investigated via using geometrical approach. It can be seen that when R 0 < 1, worm transmission is effectively controlled and verified by numerical simulation. Therefore, one effective means to extinguish worm is to keep R 0 below 1. Furthermore, the sensitivity analysis is complemented by numerical simulation to explore the influence of some parameter changes on the basic reproduction number R 0. It is seen that Λ1, Λ2, β 1, β 2, τ, η and ξ are positively correlated with R 0, while σ, ν, μ 1, μ 2, ϵ, ϕ, p, θ and γ are negatively correlated with R 0. Finally, we establish a optimal control system by taking the installation coverage rate of anti-virus software and the disinfection of computer and USB device as control variables. The numerical simulation results show that compared with the situation without control measures, the total number of computers and USB devices in the exposed state and infected state is significantly reduced after control measures are taken, that is to say, adding control could effectively reduce the spread of computer worms.

In this paper, the mathematical model is used to describe the propagation process of computer worms, which is helpful to analyze the propagation characteristics of worms, and then propose more effective defense measures and control strategies. In addition, the research of this subject is beneficial to predict the propagation trend of worms, test the effect of various defensive measures, and provide theoretical guidance for the defense and control of computer worms.

 Data availability statement

The original contributions presented in the study are included in the article/Supplementary Material, further inquiries can be directed to the corresponding author.

 Author contributions

RW, YX, and KX conceived the study. RW developed the mathematical model, computations, simulation coding, manuscript writing and data interpretation. YX contributed to supervision and validation. KX performed literature survey and computations. All authors contributed to the article and approved the submitted version.

 Funding

This work was supported by the National Natural Science Foundation of China (grant number 11971278); the Shanxi Scholarship Council of China (grant number 2022-150) and Natural Science Foundation of Shanxi Province (grant number 202203021211086).

 Conflict of interest

The authors declare that the research was conducted in the absence of any commercial or financial relationships that could be construed as a potential conflict of interest.

 Publisher’s note

All claims expressed in this article are solely those of the authors and do not necessarily represent those of their affiliated organizations, or those of the publisher, the editors and the reviewers. Any product that may be evaluated in this article, or claim that may be made by its manufacturer, is not guaranteed or endorsed by the publisher.

 References Kim KS Ibrahim MM Jung IH Kim S . Mathematical analysis of the effectiveness of control strategies to prevent the autorun virus transmission propagation. Appl Maths Comput (2020) 371:124955. 10.1016/j.amc.2019.124955 Deng Y Pei Y Li C . Parameter estimation of a susceptible–infected–recovered–dead computer worm model. SIMULATION (2022) 98:209–20. 10.1177/00375497211009576 Piqueira JRC Cabrera MA Batistela CM . Malware propagation in clustered computer networks. Physica A: Stat Mech its Appl (2021) 573:125958. 10.1016/j.physa.2021.125958 Xiao X Fu P Dou C Li Q Hu G Xia S . Design and analysis of seiqr worm propagation model in mobile internet. Commun nonlinear Sci Numer simulation (2017) 43:341–50. 10.1016/j.cnsns.2016.07.012 Zheng W Wang T Liu C Yan Q Zhan S Li G Liver transcriptomics reveals microrna features of the host response in a mouse model of dengue virus infection. Comput Biol Med (2022) 106057:106057. 10.1016/j.compbiomed.2022.106057 Carrillo M Góngora PA Rosenblueth DA . An overview of existing modeling tools making use of model checking in the analysis of biochemical networks. Front Plant Sci (2012) 3:155. 10.3389/fpls.2012.00155 Yang G Gómez Tejeda Zañudo J Albert R . Target control in logical models using the domain of influence of nodes. Front Physiol (2018) 454:454. 10.3389/fphys.2018.00454 Ren J Xu Y . A compartmental model for computer virus propagation with kill signals. Physica A: Stat Mech its Appl (2017) 486:446–54. 10.1016/j.physa.2017.05.038 Razzaq OA Khan NA Faizan M Ara A Ullah S . Behavioral response of population on transmissibility and saturation incidence of deadly pandemic through fractional order dynamical system. Results Phys (2021) 26:104438. 10.1016/j.rinp.2021.104438 Paul AK Kuddus MA . Mathematical analysis of a Covid-19 model with double dose vaccination in Bangladesh. Results Phys (2022) 35:105392. 10.1016/j.rinp.2022.105392 Tiwari V Deyal N Bisht NS . Mathematical modeling based study and prediction of Covid-19 epidemic dissemination under the impact of lockdown in India. Front Phys (2020) 8:586899. 10.3389/fphy.2020.586899 Toutonji O Yoo S-M . Passive benign worm propagation modeling with dynamic quarantine defense. KSII Trans Internet Inf Syst (Tiis) (2009) 3:96–107. 10.3837/tiis.2009.01.005 Wang F Zhang Y Wang C Ma J Moon S . Stability analysis of a seiqv epidemic model for rapid spreading worms. Comput Security (2010) 29:410–8. 10.1016/j.cose.2009.10.002 Yao Y Xie X-w. Guo H Yu G Gao F-X Tong X-j. Hopf bifurcation in an internet worm propagation model with time delay in quarantine. Math Comput Model (2013) 57:2635–46. 10.1016/j.mcm.2011.06.044 Zhou P Gu X Nepal S Zhou J . Modeling social worm propagation for advanced persistent threats. Comput Security (2021) 108:102321. 10.1016/j.cose.2021.102321 Yang L-X Yang X . The impact of nonlinear infection rate on the spread of computer virus. Nonlinear Dyn (2015) 82:85–95. 10.1007/s11071-015-2140-z Zhang Z Zou J Upadhyay RK ur Rahman G . An epidemic model with multiple delays for the propagation of worms in wireless sensor networks. Results Phys (2020) 19:103424. 10.1016/j.rinp.2020.103424 Song L-P Jin Z Sun G-Q Zhang J Han X . Influence of removable devices on computer worms: Dynamic analysis and control strategies. Comput Maths Appl (2011) 61:1823–9. 10.1016/j.camwa.2011.02.010 Zhu Q Yang X Ren J . Modeling and analysis of the spread of computer virus. Commun Nonlinear Sci Numer Simulation (2012) 17:5117–24. 10.1016/j.cnsns.2012.05.030 Ma J Chen Z Wu W Zheng R Liu J . Influences of removable devices on the anti-threat model: Dynamic analysis and control strategies. Information (2015) 6:536–49. 10.3390/info6030536 Zhu Q Loke SW Zhang Y . State-based switching for optimal control of computer virus propagation with external device blocking. Security and Communication Networks (2018). 10.1155/2018/4982523 Upadhyay RK Kumari S . Detecting malicious chaotic signals in wireless sensor network. Physica A: Stat Mech its Appl (2018) 492:1129–52. 10.1016/j.physa.2017.11.043 Gao Q Zhuang J . Stability analysis and control strategies for worm attack in mobile networks via a veiqs propagation model. Appl Maths Comput (2020) 368:124584. 10.1016/j.amc.2019.124584 Hu R Gao Q Wang B . Dynamics and control of worm epidemic based on mobile networks by seiqr-type model with saturated incidence rate. Discrete Dyn Nat Soc (2021) 2021:1–22. 10.1155/2021/6637263 Madhusudanan V Srinivas M Nwokoye CH Murthy BN Sridhar S . Hopf-bifurcation analysis of delayed computer virus model with holling type iii incidence function and treatment. Scientific Afr (2022) 15:e01125. 10.1016/j.sciaf.2022.e01125 Mishra BK Pandey SK . Dynamic model of worm propagation in computer network. Appl Math Model (2014) 38:2173–9. 10.1016/j.apm.2013.10.046 Mishra BK Pandey SK . Dynamic model of worms with vertical transmission in computer network. Appl Math Comput (2011) 217:8438–46. 10.1016/j.amc.2011.03.041 Guillén JH Del Rey AM Encinas LH . Study of the stability of a seirs model for computer worm propagation. Physica A: Stat Mech its Appl (2017) 479:411–21. 10.1016/j.physa.2017.03.023 Upadhyay RK Singh P . Modeling and control of computer virus attack on a targeted network. Physica A: Stat Mech its Appl (2020) 538:122617. 10.1016/j.physa.2019.122617 Van den Driessche P Watmough J . Reproduction numbers and sub-threshold endemic equilibria for compartmental models of disease transmission. Math biosciences (2002) 180:29–48. 10.1016/S0025-5564(02)00108-6 Castillo-Chavez C Blower S Van den Driessche P Kirschner D Yakubu A-A . Mathematical approaches for emerging and reemerging infectious diseases: Models, methods, and theory, 126. Springer Science & Business Media (2002). Khan A Zarin R Inc M Zaman G Almohsen B . Stability analysis of leishmania epidemic model with harmonic mean type incidence rate. The Eur Phys J Plus (2020) 135:528–0. 10.1140/epjp/s13360-020-00535-0 Li MY Muldowney JS . A geometric approach to global-stability problems. SIAM J Math Anal (1996) 27:1070–83. 10.1137/S0036141094266449 Khan K Zarin R Khan A Yusuf A Al-Shomrani M Ullah A . Stability analysis of five-grade leishmania epidemic model with harmonic mean-type incidence rate. Adv Difference Equations (2021) 2021:86–27. 10.1186/s13662-021-03249-4 Khan A Zarin R Ahmed I Yusuf A Humphries UW . Numerical and theoretical analysis of rabies model under the harmonic mean type incidence rate. Results Phys (2021) 29:104652. 10.1016/j.rinp.2021.104652 Yang F Zhang Z Zeb A . Hopf bifurcation of a veiqs worm propagation model in mobile networks with two delays. Alexandria Eng J (2021) 60:5105–14. 10.1016/j.aej.2021.03.055 Haldar K Mishra BK . A mathematical model for a distributed attack on targeted resources in a computer network. Commun Nonlinear Sci Numer Simulation (2014) 19:3149–60. 10.1016/j.cnsns.2014.01.028 Khanh NH Huy NB . Stability analysis of a computer virus propagation model with antidote in vulnerable system. Acta Mathematica Scientia (2016) 36:49–61. 10.1016/S0252-9602(15)30077-1 Alexander ME Moghadas SM Röst G Wu J . A delay differential model for pandemic influenza with antiviral treatment. Bull Math Biol (2008) 70:382–97. 10.1007/s11538-007-9257-2 Lukes DL . Differential equations: Classical to controlled (1982). Pontryagin LS . Mathematical theory of optimal processes. New York, NY: CRC Press (1987).