A novel dynamic model for web malware spreading over scale-free networks

Accepted Manuscript A novel dynamic model for web malware spreading over scale-free networks Wanping Liu, Shouming Zhong

PII: DOI: Reference:

S0378-4371(18)30443-6 https://doi.org/10.1016/j.physa.2018.04.015 PHYSA 19447

To appear in:

Physica A

Received date : 7 January 2018 Revised date : 10 March 2018 Please cite this article as: W. Liu, S. Zhong, A novel dynamic model for web malware spreading over scale-free networks, Physica A (2018), https://doi.org/10.1016/j.physa.2018.04.015 This is a PDF file of an unedited manuscript that has been accepted for publication. As a service to our customers we are providing this early version of the manuscript. The manuscript will undergo copyediting, typesetting, and review of the resulting proof before it is published in its final form. Please note that during the production process errors may be discovered which could affect the content, and all legal disclaimers that apply to the journal pertain.

*Highlights (for review)

Highlights 

The impact of network topology on web malware propagation is addressed.



A new dynamic model is developed for web malware spreading on scale-free networks.



The global stability of the malware-free equilibrium is proved.



Malware spread can be availably controlled by properly adjusting network structure.

*Manuscript Click here to view linked References

A novel dynamic model for web malware spreading over scale-free networks 1

Wanping Liu1,2,∗ and Shouming Zhong1 School of Mathematical Sciences, University of Electronic Science and Technology of China, Chengdu 611731, China 2 College of Computer Science and Technology, Chongqing University of Technology, Chonqing 400054, China (Dated: March 10, 2018) The complex properties of real-world systems have been found in a wide range of physical and social networks. In this paper, a novel network-based model is developed to theoretically address the structure impact of scale-free networks on epidemic web malware spreading. We further split the compartments of the SDIRS model in terms of node degree, and apply the degree-based mean-field approach to formulate the model. The spreading threshold of the new model is explicitly calculated, and the global stability of the malware-free equilibrium is theoretically addressed provided the threshold below unity. Numerical simulations are given to illustrate the theoretical results. Moreover, the dynamics around the malware equilibrium and parameter analysis are also discussed by designing some numerical examples. The predictions are generally in agreement with numerical simulation results. Consequently, we suggest that web malware can be availably contained by properly adjusting the network structure so that the propagation threshold is less than unity. Keywords: Malware spreading, scale-free networks, propagation dynamics, stability analysis

∗

[email protected]

2 I.

INTRODUCTION

Advanced network technologies achieve the integration of scattered resources all over the world, so that people can promptly share and access valuable resources through networks. In recent years, the rapidly popularized Internet-based networks, such as social networks and e-commerce networks, have brought us much convenience. However, Internet resources, both hardware and software components, can be the targets of malicious attempts to gain unauthorized control to commit financial deception or digital extortion, cause interruptions, and steal valuable information. For instance, cyber criminals can profit much from bugs in browsers and website plugins, and it was reported that a total number of 54 zero-day vulnerabilities were newly discovered in 2015, a double increase from the year before [1]. One million web attacks against people were checked each day in 2015, because of nearly three quarters of all legitimate websites having unpatched flaws. In terms of data leakage, a record-setting number of 429 million exposed identities were reported in 2015, whereas a conservative number of half a billion unreported stolen records were estimated, since more corporations avoided revealing all information breaches [1]. As network security incidents are frequently reported, the situation of cyber security is becoming more serious. Especially, the ecosystem of malicious software is becoming more and more diverse, such as computer viruses which copy with users’ intervention, Internet worms which replicate automatically, botnets, ransomware and spyware. Malware spreading will be affected by a lot of factors, such as infection rate, user’s security awareness, network topology, as well as human behavioral responses [2–4]. Therefore, it is of great significance to study the propagation mechanisms of malware and to effectively prevent and control mass invasion attacks by the dynamic modeling approach [5–7]. Contagion-like propagating processes can be taken to depict not only the evolution of infectious diseases but also of rumors, culture, ideas, information, and other spreading phenomena in fields as diverse as social, physical, and biological sciences [8]. Due to their fundamental nature and simplicity, two compartmental models have attracted extensive attention in the field of epidemic modeling, the susceptible-infected-susceptible (SIS) model and the susceptible-infected-recovered (SIR) model. In both models, every infected individual is able to successfully propagate the disease to its neighbors with a certain (spreading) probability, and each infected individual possibly recovers with some other (recovering) probability. The only difference between them is essentially that whether the immunity is temporary or not. Specifically, for the SIS model, once recovered, infected individuals return to the susceptible state because of possibly catching the disease again (temporary immunity), whereas recovered individuals of the SIR case are assumed to acquire permanent immunity and become lost in the subsequent spreading processes. Both the schemes are sufficient to capture most of disease dynamics and many other epidemic-like processes, involving disease outbreaks, however, they actually remain simple. The study of malware dynamics aims to understand the mechanisms governing the spread of malware over networks, so as to develop mathematical models for clearly predicting the spreading laws of epidemic malware. In the early stages, the modeling of malware spreading was confined to homogeneous systems, where any couple of nodes are supposed to possess the same contact probability. Epidemic models for infectious diseases were originally borrowed to study mathematical modeling of malware propagation [9, 10]. A lot of compartmental models, ranging from the SIR/SIRS models to the SEIRS models, were established based on the homogeneously mixed assumption of the propagation network, and the change in the fraction of each compartment is the focus of study [11]. Many other variations of these models, which are more realistic and intricate, have also been developed. For instance, by newly splitting the traditional susceptible compartment into two different (weakly-protected and strongly-protected) sub-compartments, the authors [12] developed a novel compartmental WSIS model (extending the SIS model) to address the impact of heterogeneous immunities of computers on the propagation of epidemic malware. Considering the propagation mechanisms of malicious hyperlinks and the impact of human intervention on web malware spreading, in Ref. [13], the authors extended the SIRS model by adding another delitescent compartment, proposing an SDIRS model to study the spreading behavior of web malware over complete or regular networks. However, these models suffer from a common defect that the knowledge concerning the propagation network structure can not be fully applied. As a result, it is extremely difficult to understand the impact of network topology on malware prevalence by solely studying such compartment-based models. In fact, a multitude number of real-world networks are heterogeneously organized, including the world-wide-web and the Internet, which have been empirically found to be highly structured rather than simply homogeneous [14, 15]. These interesting findings motivate researchers to reexamine the previous models by considering nontrivial patterns among network nodes, such as power-law degree distributions. In Ref. [16], the authors addressed by the heterogeneous mean-field method (HMF) that the epidemic threshold approaches to zero in the thermodynamic limit on scale-free networks with characteristic exponents less than 3. This observation concerning the role of network topology totally changed our previous understanding of how to model and contain the onset of epidemics, concentrating the focus of attention not only on novel methodologies to model epidemic dynamics but also on the incorporation of real contact patterns in the dynamical settings. Since then, much attention was paid to spreading dynamics of epidemic malware over complex networks [16–19]. Moreover, many computational and theoretical frameworks have been further

3 proposed, which indubitably made contagion modeling as an active research field and provided new phenomenological insights and useful approaches for studying real outbreaks. These epidemic models indicate that the structure of the propagation network plays a significant role in determining virus dynamics, i.e., heavily affecting whether viruses approach extinction. For more information about this topic, see Refs [20, 21]. Generally, network-based models can accommodate more knowledge of the network topology, thus the propagation network’s topological impact on virus spreading can be more deeply revealed by studying such models. In view of many real Internet-based networks exhibiting heterogeneous connectivity patterns, this paper is intended to formulate a new network-based model capturing the influence of network structure on web malware spreading. Our methodology allows for several new results. For that purpose, we incorporate a standard power-law distribution into the network structure. Each compartment of the SDIRS model is further split into a number of sub-compartments according to the maximum node degree. Then, a new dynamic model is proposed by the degree-based mean-field approach. Our model is able to describe the evolutions of the number of each sub-compartment with a specific node degree, largely extending the previous SDIRS model which can only roughly depict the evolutions of the four states. We compute the explicit propagating threshold for this case, and theoretically confirm the global stability of the malware-free equilibrium by establishing a theorem with a sufficient condition. Many numerical examples are also designed to illustrate our results. Finally, we also present results from extensive numerical simulations. II.

MODEL FORMULATION

Empirical research shows that a number of Internet-related networks, ranging from computer networks to social networks, are scale-free networks. In this section, we aim to apply the complex network approach to study the dynamics of web malware spreading over scale-free networks [22, 23]. Therefore, the structure of the propagating network is considered to admit a prescribed degree distribution. That is, the nodes asymptotically follow a power-law distribution P (k) ∼ k −r , k = 1, . . . , ∆, where P (k) stands for the probability that a randomly chosen node is of degree k, r is an indicator and ∆ is the maximum node degree of the complex network, indicating that P (k) = 0 for all k > ∆. In consideration of the essential differences between the spreading mechanisms of web malware and epidemic diseases, delitescent properties were introduced for the nodes receiving messages containing malicious links but not being clicked or activated yet [13]. Thus, for epidemic malware prevalence in the SDIRS scenarios, the whole network nodes are divided into four groups, susceptible nodes, delitescent nodes, intruded nodes, and recovered nodes. For convenience, several notations and quantities are introduced as follows. (1) Sk (t): the number of S-nodes with degree k at time t. (2) Dk (t): the number of D-nodes with degree k at time t. (3) Ik (t): the number of I-nodes with degree k at time t. (4) Rk (t): the number of R-nodes with degree k at time t. (5) Nk (t): the number of all network nodes with degree k at time t. That is, Nk (t) := Sk (t) + Dk (t) + Ik (t) + Rk (t). (6) Sek (t): the relative density of S-nodes with degree k at time t, i.e., Sek (t) = Sk (t)/Nk (t).

e k (t): the relative density of D-nodes with degree k at time t, i.e., D e k (t) = Dk (t)/Nk (t). (7) D

(8) Iek (t): the relative density of I-nodes with degree k at time t, i.e., Iek (t) = Ik (t)/Nk (t).

ek (t): the relative density of R-nodes with degree k at time t, i.e., R ek (t) = Rk (t)/Nk (t). (9) R

Next, we introduce several model parameters and impose some necessary assumptions, see also [13].

(A1) The parameter η represents the probability that a D-node turns into the recovered compartment per unit time. (A2) The parameter ζ depicts the probability of a R-node leaving for the susceptible compartment owning to immunity failure. (A3) The infection rate λ is the probability that an S-node receives malicious links sent by a neighboring I-node within a unit time. (A4) The parameter ε describes the probability of that a D-node leaves the delitescent compartment for the infected compartment.

4 (A5) The parameter γ describes the probability that an I-node system is fully scanned and immune, immediately turning to be recovered. Here, we only consider malware spreading over static networks, thus none of nodes will be newly-connected to or disconnected from the network, i.e., both the connection rate and the disconnection probability are zero. Besides, we also assume that the probability that a link has an I-node as one endpoint is independent on the degree of the other endpoint of the link and, thus, this probability is only a function of the components of Iek (t), k = 1, . . . , ∆. Let Θ(t) denote this probability, and direct probabilistic calculations yield Θ(t) =

1 X kP (k)Iek (t). k hki

Following the above assumptions and the compartment-based model in [13], we derive the following new degree-based model  dSek (t)      dt    e    dDk (t)  dt  e d I  k (t)     dt     ek (t)  d R  dt

ek (t), 1 ≤ k ≤ ∆, = −kλΘ(t)Sek (t) + ζ R

e k (t), 1 ≤ k ≤ ∆, = kλΘ(t)Sek (t) − (η + ε)D

(1)

e k (t) − γ Iek (t), 1 ≤ k ≤ ∆, = εD

e k (t) + γ Iek (t) − ζ R ek (t), 1 ≤ k ≤ ∆, = ηD

e k (0), Iek (0), R ek (0) ≥ 0, k = 1, . . . , ∆. The dynamical transfers among the compartments where the initial states Sek (0), D are depicted in Fig. 1.

FIG. 1. Schematic diagram for the state transitions of the nodes over scale-free networks.

Remark 1. Compared with the SDIRS model proposed in [13] which can only roughly describe the evolutions of four basic states, this new model (1) is high-dimensional, consisting of 4∆ equations, and thus is able to deeply depict the dynamics of the number of each group with a specific node degree. e k (t) + Iek (t) + R ek (t) = 1 holds for each k = 1, . . . , ∆, thus the first set of equations in system Note that Sek (t) + D (1) can be ignored and system (1) is reduced as follows  e k (t)  dD  e k (t) − Iek (t) − R ek (t)) − (η + ε)D e k (t), 1 ≤ k ≤ ∆,  = kλΘ(t)(1 − D   dt    e dIk (t) e k (t) − γ Iek (t), 1 ≤ k ≤ ∆, (2) = εD  dt     e   e k (t) + γ Iek (t) − ζ R ek (t), 1 ≤ k ≤ ∆,  dRk (t) = η D dt e k (0), Iek (0), R ek (0) ≥ 0, and D e k (0) + Iek (0) + R ek (0) ≤ 1, k = 1, . . . , ∆. where the initial states D

5 3∆ Denote R3∆ |xj ≥ 0, j = 1, 2, . . . , 3∆}, and the meaningful domain for system (2) is + = {(x1 , x2 , . . . , x3∆ ) ∈ R 3∆ Ω1 = {y = (y1 , y2 , . . . , y3∆ ) ∈ R+ |yj + yj+∆ + yj+2∆ ≤ 1, j = 1, . . . , ∆}, which can be confirmed as positively invariant for system (2). The following lemma is necessary to prove that Ω1 is a positive invariant for system (2), also see [17, 24].

Lemma II.1 Consider a system dx/dt = f(x) which is defined at least in a compact set C. Then, C is invariant if, for every point y on ∂C (the boundary of C), the vector f(y) is tangent to or pointing into C. Then we have the following lemma which concerns the set Ω1 . Lemma II.2 The set Ω1 is positively invariant for system (2). Proof. Note that ∂Ω1 (the boundary of Ω1 ) consists of the following 4∆ sets: Aj Bj Cj Dj

= {y ∈ Ω|yj = 0}, j = 1, 2, . . . , ∆, = {y ∈ Ω|yj+∆ = 0}, j = 1, 2, . . . , ∆, = {y ∈ Ω|yj+2∆ = 0}, j = 1, 2, . . . , ∆, = {y ∈ Ω|yj + yj+∆ + yj+2∆ = 1}, j = 1, 2, . . . , ∆,

which have j

φj = {0, . . . , 0, −1, 0, . . . , 0}, j = 1, 2, . . . , ∆, j+∆

ϕj = {0, . . . , 0, −1 , 0, . . . , 0}, j = 1, 2, . . . , ∆, j+2∆

ψj = {0, . . . , 0, −1 , 0, . . . , 0}, j = 1, 2, . . . , ∆, j

j+∆

j+2∆

ηj = {0, . . . , 0, 1, 0, . . . , 0, 1 , 0, . . . , 0, 1 , 0, . . . , 0}, j = 1, 2, . . . , ∆, as their outer normal vectors, respectively. Next, consider system (2), for j = 1, . . . , ∆, calculations yield P kP (k)yj+∆ dy ( |y∈Aj , φj ) = −jλ k (1 − yj+∆ − yj+2∆ ) ≤ 0, dt hki dy ( |y∈Bj , ϕj ) = −εyj ≤ 0, dt dy ( |y∈Cj , ψj ) = −(ηyj + γyj+∆ ) ≤ 0, dt dy ( |y∈Dj , ηj ) = −ζyj+2∆ ≤ 0. dt From the above inequalities, we know that dy/dt is tangent to or pointing into Ω1 for arbitrary y ∈ ∂Ω1 . Thus, the claimed result follows from Lemma II.1.  III.

PROPAGATION THRESHOLD AND EQUILIBRIA OF SYSTEM (2)

The propagation threshold of system (2), denoted by R0 , can be computed as (see Appendix A) s λε hk 2 i R0 = , hki γ(ε + η)

(3)

P P where hk 2 i = k k 2 P (k), hki = k kP (k). Remark 2. The propagation threshold R0 plays an important role in determining the dynamics of system (2). It can be observed that the model parameters λ, ε, γ and η obviously affect the value of R0 . Besides, the part hk 2 i/hki in the expression of R0 involves the degree distribution of the scale-free network, implying the topological effect on R0 . To study the dynamics of model (2), we begin with the following definitions of two equilibria.

6 Definition III.1 Let ∗ e∗ ∗ e∗ ∗ T e 1∗ , . . . , D e∆ e∆ E = (D , I1 , . . . , Ie∆ , R1 , . . . , R )

be an equilibrium of model (2). Then, we define ∗ (a) E is said to be a malware-free equilibrium if Ie1∗ = · · · = Ie∆ = 0. (b) E is said to be a malware equilibrium if there exists at least one j ∈ {1, . . . , ∆} such that Iej∗ > 0.

Obviously, system (2) always admits a malware-free equilibrium E0 = (0, . . . , 0, 0, . . . , 0, 0, . . . , 0)T , i.e., all infection components being zero. Certain calculations yield the following lemma. Lemma III.1 Consider system (2). The following assertions hold: 3∆

z }| { (1) There always exists a malware-free equilibrium E0 = (0, 0, · · · , 0)T . (2) There exists no malware equilibrium if R0 ≤ 1.

e 1∗ , . . . , D e ∗ , Ie1∗ , . . . , Ie∗ , R e1∗ , . . . , R e∗ )T if R0 > 1, i.e., (3) There exists a unique malware equilibrium E ∗ = (D ∆ ∆ ∆ e k∗ = γ Iek∗ , R ek∗ = γ (1 + η )Iek∗ , D ε ζ ε kλΘ∗ , Iek∗ = γ + ηγ/ε + kλΘ∗ (γ/ε + 1 + γ(1 + η/ε)/ζ)

where Θ∗ is the unique positive zero of the following function

f (x) =

λ X k 2 P (k) − 1. hki γ + ηγ/ε + kλx(γ/ε + 1 + γ(1 + η/ε)/ζ)

(4)

k

e 1, . . . , D e ∆ , Ie1 , . . . , Ie∆ , R e1 , . . . , R e∆ )T be an equilibrium of system (2). Then we have Proof. Let E = (D  e k − Iek − R ek ) − (ε + η)D e k = 0,  kλΘ(1 − D e k − γ Iek = 0, εD  e ek = 0. η Dk + γ Iek − ζ R

e k = 0, R ek = 0 for each k = 1, . . . , ∆, and thus Sek = 1, k = By letting Iek = 0, k = 1, . . . , ∆, then we can easily get D e k + Iek + R ek = 1. Therefore, system (2) always possesses a unique malware-free equilibrium 1, . . . , ∆ because of Sek + D 3∆

z }| { E0 = (0, 0, · · · , 0)T . In the sequel, let

∗ e∗ ∗ e∗ ∗ T e 1∗ , . . . , D e∆ e∆ E ∗ = (D , I1 , . . . , Ie∆ , R1 , . . . , R )

be a malware equilibrium of system (2). Then, it follows by system (2) that  e ∗ − Ie∗ − R e ∗ ) − (ε + η)D e ∗ = 0,  kλΘ∗ (1 − D k k k k ∗ ∗ e − γ Ie = 0, , εD k k  e∗ e∗ = 0. η Dk + γ Iek∗ − ζ R k where

Θ∗ =

1 X kP (k)Iek∗ . k hki

e∗ = γ (1 + η )Ie∗ , and e ∗ = γ Ie∗ , R From the above identities, we derive D k ε k k ζ ε k Iek∗ =

kλΘ∗

γ + ηγ/ε +

kλΘ∗ (γ/ε +

1 + γ(1 + η/ε)/ζ)

, k = 1, . . . , ∆.

7 Substituting the above equation into the expression of Θ∗ yields 1 X kP (k)Ik∗ Θ∗ = k hki k 2 λP (k)Θ∗ 1 X . = k γ + ηγ/ε + kλΘ∗ (γ/ε + 1 + γ(1 + η/ε)/ζ) hki

Thus Θ∗ is a positive root of the equation f (x) = 0 shown in (4). It is obvious that the function f (x) is monotonically decreasing in the interval [0, +∞), i.e., f ′ (x) < 0 for all x ≥ 0. Note that λ X k 2 P (k) f (0) = − 1 = R20 − 1, hki γ + ηγ/ε k

and f (+∞) = −1 < 0, thus f (x) = 0 has a unique positive root if R0 > 1. Hence, assertion (3) also holds. The proof is complete. ∗ Easily, we can get the inequality Ie1∗ < Ie2∗ < · · · < Ie∆ , which shows that, in the steady state, the infection density of states with higher node degree is higher than that of states with lower node degree. IV.

DYNAMIC ANALYSIS OF THE MODEL

It is of great importance to fully analyze the dynamics around the malware-free equilibrium of the model. For instance, when web malware has already intruded most of the network, the stability results will motivate us some proper measures to control or even eradicate malware spreading, such as adjusting related parameters so as to satisfy the condition. In this section, we aim to theoretically address the dynamics of system (2) with respect to Ω1 when the threshold R0 < 1. In the sequel, denote e 1 (t), . . . , D e ∆ (t), Ie1 (t), . . . , Ie∆ (t), R e1 (t), . . . , R e∆ (t))T y(t) = (D

and rewrite system (2) in matrix-vector notation as

dy(t) = Ay(t) + H(y(t)), dt where the initial state y(0) ∈ Ω1 , and

(5)



 −(η + ε)E∆ B 0 εE∆ −γE∆ 0 , A = ηE∆ γE∆ −ζE∆ 1 ijλP (j), 1 ≤ i, j ≤ ∆, B = (Bij )∆×∆ , Bij = hki

T   z 2∆ }| { e 1 (t) + Ie1 (t) + R e1 (t), . . . , ∆ D e ∆ (t) + Ie∆ (t) + R e∆ (t) , 0, . . . , 0 . H(y(t)) = −λΘ(t)D 

A.

Local stability of E0

The following lemma is necessary for proving the local stability of the malware-free equilibrium. Lemma IV.1 Denote a square matrix of size ∆ × ∆ by  (1 × 1)P (1) (1 × 2)P (2) . . . (1 × ∆)P (∆)  (2 × 1)P (1) (2 × 2)P (2) . . . (2 × ∆)P (∆)  M∆ =  .. .. ..  . . . (∆ × 1)P (1) (∆ × 2)P (2) . . . (∆ × ∆)P (∆)



  , 

8 P∆ where ∆ is the maximum node degree of the static scale-free network involved in the system (2), P (k) = k −r / k=1 k −r , k = 1, . . . , ∆. Then, we get that det(M∆ − xE∆ ) = (−x)∆−1 (hk 2 i − x), where E∆ represents the ∆ × ∆ identity matrix. The recursive relationship can be derived by certain algebraic calculations. Then, this Lemma IV.1 can be simply proved by induction on ∆ (also see Lemma 4.1 in Ref. [17]). Theorem IV.1 Consider system (2). The malware-free equilibrium E0 is locally asymptotically stable if R0 < 1, whereas it is a saddle point if R0 > 1. Proof. The jacobian matrix of system (2) evaluated at the malware-free equilibrium E0 is   −(η + ε)E∆ B 0 εE∆ −γE∆ 0 , J= ηE∆ γE∆ −ζE∆

Then, it follows by Lemma IV.1 and certain calculations that the characteristic equation of the above jacobian matrix is (x + η + ε)E∆ −B 0 −εE∆ (x + γ)E∆ 0 det (xE3∆ − J) = −ηE∆ −γE (x + ζ)E∆ ∆ (x + η + ε)E∆ −B = (x + ζ)∆ −εE∆ (x + γ)E∆ = (x + ζ)∆ det((x + η + ε)(x + γ)E∆ − εB) = (x + ζ)∆ (x + η + ε)∆−1 (x + γ)∆−1

2 ! k × (x + η + ε) (x + γ) − ελ hki = (x + ζ)∆ (x + η + ε)∆−1 (x + γ)∆−1

2 ! k 2 = 0. × x + (η + γ + ε)x + (η + ε)γ − ελ hki

The above equation has negative roots −ζ with multiplicity ∆, and negative roots −(η +ε) and −γ with multiplicity ∆ − 1, respectively. Now let

2 k 2 g(x) = x + (η + γ + ε)x + (η + ε)γ − ελ = 0. hki hk2 i Obviously, the above identity possesses two real roots. If R0 < 1, then (η + ε)γ − ελ hki > 0. Thus, it follows by the Hurwitz criterion that all roots of the characteristic equation have negative real parts, implying that E0 is locally hk2 i asymptotically stable. For the other case R0 > 1, then (η + ε)γ − ελ hki < 0. In this case, the characteristic equation exactly obtains one positive root, implying that E0 is a saddle point.  B.

Global stability of E0

Next, we will further prove the global stability of the malware-free equilibrium E0 . Consider an n × n matrix M, denote s(M) = max1≤i≤3∆ {Re λi }, where λi , i = 1, . . . , 3∆ are the eigenvalues of the matrix M, and Re(·) denotes the real part. To examine the global stability for E0 , we need the following lemma. Lemma IV.2 [24] Consider an n-dimensional autonomous system dx(t) = Bx(t) + H(x(t)), x(t) ∈ D, dt where D is a region containing the origin, H(x(t)) ∈ C 1 (D), limx→0 kH(x)k/kxk = 0. Assume that there are a positively invariant compact convex set C ⊂ D containing the origin, a positive number r, and a real eigenvector ω of BT such that

9 (C1) (x, ω) ≥ rkxk for all x ∈ C, (C2) (H(x), ω) ≤ 0 for all x ∈ C, (C3) the origin forms the largest positively invariant set included in the set {x ∈ C|(H(x), ω) = 0}. Let s(B) denote the maximum real part of all eigenvalues of B. Then we get (1) s(BT ) < 0 implies that the origin is globally asymptotically stable in C, (2) s(BT ) > 0 implies that there exists m > 0 such that x(0) ∈ C − {0} implies lim inf t→∞ kx(t)k ≥ m. Theorem IV.2 Consider system (2). Then, E0 is globally asymptotically stable with respect to Ω1 provided that R0 < 1. Proof. Let C = Ω1 . Note that the matrix AT is irreducible, and all of its non-diagonal entries are nonnegative, thus it follows by [24] that AT has a positive eigenvector ω = (ω1 , . . . , ω3∆ ) corresponding to its eigenvalue s(AT ). Let ω0 = min1≤j≤3∆ {ωj } > 0. Then, for all y ∈ Ω, we have (y, ω) ≥ ω0

3∆ X i=1

(H(y), ω) = −λΘ

yi ≥ ω0

∆ X i=1

3∆ X i=1

yi2

! 12

= ω0 kyk,

iωi (yi + yi+∆ + yi+2∆ ) ≤ 0.

Moreover, it follows by (H(y), ω) = 0 that yk+∆ = 0, k = 1, . . . , ∆, which together with system (2) leads to that the origin forms the largest positively invariant set in {y ∈ Ω|(H(y), ω) = 0}. In the above, we have proved that s(AT ) < 0 if R0 < 1, hence, the claimed result follows from assertion (1) of Lemma IV.2. V.

NUMERICAL SIMULATIONS

e e In this section, several numerical examples are designed to illustrate the dynamics of system (2). Let D(t), I(t), e R(t) denote the percentages of D-(I-, R-)nodes in all network nodes, respectively. We first show the formulae of e e and R(t) e e k (t), Iek (t) and R ek (t), respectively. Note that D e k (t) = Dk (t)/Nk (t), D(t), I(t) which can be expressed by D ek (t) = Rk (t)/Nk (t), then we get Iek (t) = Ik (t)/Nk (t), R P e X k Dk (t)Nk (t) e k (t)P (k), = D N k P P e X I (t) I (t)N (t) k k e = k k I(t) = k = Iek (t)P (k), N N k P P X e ek (t)P (k), e = k Rk (t) = k Rk (t)Nk (t) = R R(t) N N

e D(t) =

P

k

Dk (t) = N

k

P

where N denotes the whole number of network nodes, i.e., N = k Nk . In order to verify the correctness of Theorem IV.2 theoretically proved above, we first design the following experiment. Example V.1 Consider a scale-free network with ∆ = 50, r = 1.5 and the parameters of system (2) are specifically given by λ = 0.0005, ζ = 0.01, η = 0.1, ε = 0.2, γ = 0.02. Then, it follows by (3) that R0 = 0.5589 < 1, which together with Lemma III.1 implies that system (2) has the only malware-free equilibrium E0 (see Fig. 2). For the case shown in Example V.1, Theorem IV.2 implies the global stability of the malware-free equilibrium E0 . By fixing the initial vector of system (2), Fig. 2(a) shows the dynamics of system (2), where it can be seen that e e e e k (t) and D(t) e all components of system (2) and D(t), I(t), R(t) finally tend to zero. Specially, the solutions of D e e R(t)) e extremely rapidly converge to zero. Fig. 2(b) shows the evolutions of (D(t), I(t), by randomly giving 40 sets of

10 1

1 0.8

ek (t) R

0.6

0.6

e D(t)

0.4

e I(t)

0.5

e R(t)

0.2 0.4 0

0.3

Percentage of R−nodes

Iek (t)

0.8 0.7

0.5

e k (t) D

0.9

0

5

10

0.4

0.3

0.2

0.1

15

0.2

0 1

0.1 0.5 0

0

100

200

300

t

400

500

600

0 Percentage of I−nodes

(a)

0.2

0

0.4

0.6

0.8

1

Percentage of D−nodes

(b)

FIG. 2. Malware dynamics of system (2) with the only malware-free equilibrium. a) Evolutions of all solutions to system (2) e e and R(t) e with a set of specific randomly-given initial values. b) Phase diagram of D(t), I(t) with 40 different initial vectors (randomly given).

initial values for system (2). It can be observed in Fig. 2(b) that all trajectories starting from different initial points eventually approach to the green zero point, agreeing with the result of Theorem IV.2. In the previous section, we just theoretically address the global stability for the malware-free equilibrium E0 of model (2) with parameters implying R0 < 1. Due to the complexity, the other parameter cases such that R0 > 1 are not theoretically considered yet. Thus, the following experiment is designed to show the dynamics of system (2) with a malware equilibrium. Example V.2 Consider a scale-free network with ∆ = 50, r = 1.5 and the parameters of system (2) are specifically given by λ = 0.005, ζ = 0.01, η = 0.01, ε = 0.03, γ = 0.004. Then, it follows by (3) that R0 = 4.1920, which is obviously greater than unity. Thus, in this case, Lemma III.1 implies that system (2) has a malware equilibrium E ∗ . By randomly fixing a set of initial values for system (2), Fig. 3(a) shows evolutions of all components of system (2), where it can be seen that all of them finally tend to their corresponding constants. Fig. 3(b) shows the evolutions of e e R(t)) e (D(t), I(t), by randomly fixing 20 sets of initial values for system (2). It can be observed in Fig. 3(b) that all trajectories starting from different initial points finally approach to the same red point. Based on a multitude number of numerical simulations such as Example V.2 and Fig. 3, we can present a conjecture that the malware equilibrium E ∗ of system (2) is globally asymptotically stable with respect to Ω1 − {0} provided R0 > 1. Next, we aim to examine the impact of model parameters as well as the network topology on malware spreading. Deep understanding and analysis of the propagation thresholds of compartmental models may provide some insights to develop feasible measures to effectively inhibit or even control malware prevalence over networks. The model parameters λ, ε, γ and η in (3) are viewed as variables, respectively. Then, obviously, R0 is strictly decreasing with respect to γ and η, whereas R0 is strictly increasing with respect to λ. For the remaining parameter ε, it follows by some calculations that ∂R0 λη ηR0 1 hk 2 i = = > 0. 2 ∂ε 2R0 hki γ(ε + η) 2ε(ε + η)

Thus, R0 is monotonically increasing with respect to ε. From the above analysis, we can infer that reducing the parameter ε and the invasion rate λ is helpful to inhibit malware spreading, whereas raising the parameters γ and η is conducive to control malware prevalence. Next, an experiment is designed to illustrate how the propagation threshold R0 evolves over scale-free networks with different topologies. Except for ∆ and r, the other parameters in (3) are specifically fixed. Thus, without loss of generality, λ, ε, γ and η can be chosen such that λε/[γ(ε + η)] = 0.1. Fig. 4 shows the values of R0 with varying ∆ ∈ [40, 100] and r ∈ [1.2, 2.6]. It is easily found that R0 is greater over a scale-free network with larger maximum node degree ∆ and smaller exponent r.

11 1

0.5 e k (t) D

e D(t)

Iek (t)

0.8

e I(t)

ek (t) R

0.7 0.6

0.45 0.4 Percentage of R−nodes

0.9

e R(t)

0.5 0.4

0.3 0.25 0.2 0.15 0.1

0.3

0.05

0.2

0 0

0.1 0

0.35

0.5 0

100

200

300

t

(a)

400

500

600

0

0.2

0.4

0.6

0.8

Percentage of I−nodes

Percentage of D−nodes

(b)

FIG. 3. Malware dynamics of system (2) having a malware equilibrium. a) Solutions of system (2) with a set of specific initial e e and R(t) e values (randomly given). b) Phase diagram of D(t), I(t) with 20 different sets of initial vectors (randomly given).

FIG. 4. Values of R0 as a function of the maximum node degree ∆ and the scale-free network index r. The other parameters in (3) are specifically fixed such that λε/[γ(ε + η)] = 0.1.

Remark 2. The propagation threshold R0 involves the power law degree distribution of scale-free networks characterized by ∆ and r. For the sake of inhibiting malware prevalence over scale-free networks, it follows by the above numerical simulations that we can reduce the threshold value not only by properly adjusting the parameter values, but also by changing the network topology, such as reducing the maximum node degree ∆ and increasing the index r. Example V.3 Consider a scale-free network with ∆ = 100, r = 2.5 and consider system (2) with parameters λ = 0.005, ε = 0.2, η = 0.3, ζ = 0.01. e of system (2) with several different values of γ. Different colors represent Figure 5 shows the evolutions of I(t) distinct values of the parameter γ, while different types of lines show the evolutions with several sets of initial values. e eventually converge to corresponding equilibrium points at higher speeds as As shown in Figure 5, the values of I(t) the parameter γ grows. This implies that the number of I-nodes will become stable more quickly through increasing the parameter γ, such as timely and comprehensive security scan for the systems soon after the outbreak of malware. It can be also observed in Figure 5 that the final density of intruded nodes possesses significant decrease when the value of parameter γ increases from 0.001 to 0.007, while it has inconspicuous decrease when γ values from 0.007 to

12 0.02. This indicates that increasing the recovering rate will contribute less to further keep the number of I-nodes at a lower level when the recovering rate γ is large enough.

0.6

0.5 γ=0.001 γ=0.004

0.4

e I(t)

γ=0.007 γ=0.01

0.3

γ=0.02 0.2

0.1

0

0

500

1000

1500

2000

2500

3000

t

e in system (2) with varying γ and other parameters given in Example V.3. Colors represent different FIG. 5. Evolutions of I(t) values of γ; linetypes show the evolutions with distinct starting initial vectors.

Example V.4 Consider a scale-free network with ∆ = 50, r = 2.1 and consider system (2) with parameters λ = 0.01, ε = 0.1, γ = 0.01, ζ = 0.01. e in system (2) with several different values of η will finally tend to It is shown in Figure 6 that the evolutions of I(t) the equilibrium, respectively. Furthermore, the I-component of the equilibrium reduces more slowly as the parameter grows. More clearly, it decreases more quickly as η increases in the interval (0, 0.2] than (0.2, 0.9]. 0.8 η=0 0.7

η=0.1 η=0.2

0.6

η=0.3 0.5

e I(t)

η=0.9

0.4 0.3 0.2 0.1 0

0

200

400

600

800

1000

t

e in system (2) with varying η and other parameters given in Example V.4. Colors represent different FIG. 6. Evolutions of I(t) values of η; linetypes represent evolutions starting from distinct initials.

Example V.5 Consider a scale-free network with ∆ = 100, r = 2.5 and consider system (2) with parameters ε = 0.2, γ = 0.01, η = 0.1, ζ = 0.01. e evolves It is well known that greater intrusion rate will benefit web malware propagation. Figure 7 shows how I(t) e under several different cases of λ, where I(t) undergoes a nearly common jump for all cases, and then decrease and finally remain constant at a higher number of I-nodes for larger value of λ.

13

0.7 λ=0.001 λ=0.005

0.6

λ=0.01 λ=0.05

0.5

λ=0.1

e I(t)

0.4

0.3

0.2

0.1

0

0

100

200

300

400

500

600

700

800

900

t

e in system (2) with varying λ and other parameters given in Example V.5. FIG. 7. Evolutions of I(t)

The parameter ε describes the probability with which a user will click on the malicious hyperlink after reading the malware message. Obviously, greater value of ε is helpful to spread web malware. The next example is designed to show how the dynamics of model (2) evolve with several varying ε-values. Example V.6 Consider a scale-free network with ∆ = 60, r = 2.5 and consider system (2) with parameters λ = 0.01, γ = 0.02, η = 0.2, ζ = 0.2. Figure 8 quantitatively shows how the parameter ε contributes web malware spreading. It can be seen that the percentage of I-nodes will finally keep at a higher number as ε increases. This illustrates the positive impact of ε on web malware spreading.

0.6 ε=0.1 ε=0.2

0.5

ε=0.4 ε=0.6 0.4

e I(t)

ε=0.8

0.3

0.2

0.1

0

0

200

400

600

800

1000

t

e in system (2) with several cases of ε and other parameters specifically shown in Example V.6. FIG. 8. Evolutions of I(t)

Note that the parameter ζ is not included in the expression of the propagation threshold R0 , i.e., the value of R0 will not be affected by ζ. This indicates that the dynamics of model (2) is generally determined by the other parameters ε, η, γ, λ and the topology of the spreading network. However, it is clearly noted that ζ is incorporated in the expressions of the malware equilibrium E ∗ when R0 > 1. Thus, the value of ζ will affect the final states at which the percentages of each kind of nodes will keep. Example V.7 Consider a scale-free network with ∆ = 100, r = 2.5 and consider system (2) with parameters λ = 0.01, ε = 0.3, γ = 0.01, η = 0.2.

14 It can be seen in Figure 9 that the percentages of I-nodes will increase dramatically in the early time and then decrease until achieve a steady state. Ie∗ values a higher number as ζ increases, exhibiting significant change when ζ ∈ [0.001, 0.1], while slight variety for ζ ∈ [0.1, 0.4]. 0.45 ζ=0.001 0.4

ζ=0.01 ζ=0.1

0.35

ζ=0.2 0.3 ζ=0.4

e I(t)

0.25

0.2

0.15 0.1 0.05 0

0

100

200

300

400

500

600

700

800

t

e in system (2) with different ζ and other parameters specifically fixed in Example V.7. FIG. 9. Evolutions of I(t)

The following example is designed to illustrate the influence of network structure on malware spreading. Example V.8 Consider system (2) with parameters λ = 0.005, ε = 0.2, γ = 0.02, η = 0.1, ζ = 0.2.

e of system (2) with several different values of r and ∆, respectively. It is Figure 10 shows the evolutions of I(t) e eventually converge to corresponding equilibrium points at slower observed from Figure 10(a) that the values of I(t) rates as the parameter r increases. It is also indicated in Figure 10(a) that a scale-free network with lower power-law e with several exponent would benefit the spread of malware. By fixing r = 2.5, Figure 10(b) gives the evolutions of I(t) different ∆-values, showing that the final percentage of I-nodes will be higher over networks with greater ∆. Therefore, except for changing related parameters, these numerical simulations indicate that it is conducive to inhibit malware spreading over scale-free networks by properly adjusting the network structure, such as increasing the exponent r and applicably deleting the edges around the maximum-degree nodes (without destroying the connectivity) so as to make the greatest degree ∆ smaller. 0.5

0.7 r=2

0.45

∆=50 0.6

r=2.2 0.4

r=2.4

0.35

r=2.6

∆=100 ∆=200

0.5

∆=300

r=2.8

∆=400

0.3

e I(t)

e I(t)

0.4

0.25

0.3

0.2 0.15

0.2

0.1 0.1 0.05 0

0

100

200

300

400

500

600

700

800

900

0

0

100

200

300

400

500

t

t

(a)

(b)

600

700

800

900

FIG. 10. Malware dynamics over different scale-free networks. a) The exponent r is varying, and the maximum node degree is fixed as ∆ = 150. b) ∆ is varying, and the exponent r is fixed as r = 2.5.

15 VI.

CONCLUSION

In this study, we mainly examine how the topology of propagation networks affects malware spreading processes. By incorporating the standard power-law distribution into the propagating network, a scale-free network-based model is newly developed. The propagating threshold of the model is explicitly computed, see (3), which includes hk 2 i/hki, implying the effect of node degree distribution. Theoretically, the dynamic properties of model (1) are carefully addressed, proving the global stability of the malware-free equilibrium if the threshold is less than unity. Because of the complexity of the system, theoretical dynamics of the malware equilibrium are extremely hard to address. But, we design a lot of numerical simulations which intuitively show how the system evolves around the malware equilibrium, motivating a conjecture of the global stability of the malware equilibrium. Our numerical simulations suggest that malware prevalence over networks, to a certain degree, can be effectively prevented and controlled by properly modifying the network topology such that the characteristics of the propagating networks satisfy the requirements. However, this model can be mainly applied to depict malware propagation over some real networks (e.g., webpage networks or social networks) that are scale-free or approximately follow the pow-law distribution. Our next work will further modify and improve the model by establishing a node-level system which involves in more accurate impact of the network structure on malware propagation. Moreover, we will also consider simulated data or even real data for such network topology and then consider the functioning of these models on that data.

ACKNOWLEDGMENTS

This research was supported by National Natural Science Foundation of China (Grant No. 61603065), Project of Humanities and Social Sciences of Ministry of Education of China (Grant No. 15YJC790061), Project of National Bureau of Statistics of China (Grant No. 2016LZ08), fund of Chongqing Science and Technology Committee (Grant Nos. cstc2016jcyjA0076, cstc2017jcyjAX0277), and China Postdoctoral Science Foundation (Grant No. 2017M612946).

Appendix A: Calculations for the propagation threshold of model (1)

In Ref. [25], the authors developed a standard method to address the problem of computing the threshold of compartmental models. For convenience, it is presented below. Consider an n−dimensional dynamical model for virus spreading, where the first m variables correspond to all infected compartments which are numbered as compartment 1 through m, and the left n − m compartments which are numbered as compartment m + 1 through n correspond to uninfected nodes. Denote a variable vector x = (x1 , x2 , . . . , xn ), where xi denote the number (or fraction) of nodes in the i-th compartment. Let Fi (x) be the rate of appearance of new infections into compartment i, Vi+ (x) be the rate of transfer of nodes into compartment i by all other means, and Vi− (x) be the rate of transfer of nodes out of compartment i. Then the considered model can be rewritten as the following system dxi = fi (x) = Fi (x) − Vi (x), i = 1, 2, . . . , n, dt where Vi (x) = Vi− (x) − Vi+ (x). Let F(x) = (F1 , . . . , Fn )T , f(x) = (f1 , . . . , fn )T , and V(x) = (V1 , . . . , Vn )T . For the functions in the above system, five assumptions (A1)–(A5) are described below. (A1) If x ≥ 0, then Fi , Vi− , Vi+ ≥ 0 for i = 1, 2, . . . , n. (A2) If xi = 0, then Vi− = 0. In particular, if x ∈ Xs := {x ≥ 0|xi = 0, i = 1, . . . , m}, which is defined as the set of all infection-free states, then Vi− = 0 for i = 1, . . . , m. (A3) Fi (x) = 0 if i > m. (A4) If x ∈ Xs , then Fi (x) = 0 and Vi+ (x) = 0 for i = 1, . . . , m. (A5) If F(x) = (F1 , . . . , Fn )T is set to zero, then all eigenvalues of Df(x0 ) have negative real parts, where Df(x0 ) is the derivative [∂fi /∂xi ](i.e., Jacobian matrix) evaluated at x0 which is a (locally asymptotically) stable equilibrium.

16 Then, van den Driessche and Watmough proved a useful lemma (see Lemma 1 in the reference [25]). That is, if the above assumptions (A1)–(A5) are satisfied, then the derivatives DF(x0 ) and DV(x0 ) are partitioned as     V 0 F 0 , DV(x0 ) = , DF(x0 ) = J3 J4 0 0 where F and V are the m × m matrices defined by     ∂Vi ∂Fi (x0 ) and V = (x0 ) with 1 ≤ i, j ≤ m. F = ∂xj ∂xj Further, F is non-negative, V is a non-singular matrix and all eigenvalues of J4 have positive real part. Then, the threshold can be defined as ρ(F V −1 ), where ρ(A) denotes the spectral radius of a matrix A (refer to the reference [25]). Next, in order to formulate the threshold of the compartmental model, we denote x(t) = (x1 (t), x2 (t), . . . , x3∆ (t))T e 1 (t), . . . , D e ∆ (t), Ie1 (t), . . . , Ie∆ (t), R e1 (t), . . . , R e∆ (t))T . = (D

Then the scale-free network-based SDIRS model (2) can be rewritten as follows          

where Θ(t) =

1 hki

P

k

dxk (t) = kλΘ(t)(1 − xk (t) − x∆+k (t) − x2∆+k (t)) dt − (η + ε)xk (t), 1 ≤ k ≤ ∆, dx (t) ∆+k  = εxk (t) − γx∆+k (t), 1 ≤ k ≤ ∆,    dt      dx2∆+k (t) = ηxk (t) + γx∆+k (t) − ζx2∆+k (t), 1 ≤ k ≤ ∆, dt

kP (k)x∆+k (t). Note that the above system can be also rewritten as the following vector form dx = f(x) = F(x) − V(x), dt

− T + T where F = (F1 , F2 , . . . , F3∆ )T , V(x) = V− (x) − V+ (x), V− = (V1− , V2− , . . . , V3∆ ) , V+ = (V1+ , V2+ , . . . , V3∆ ) , and



λΘ(1 − x1 − x∆+1 − x2∆+1 ) 2λΘ(1 − x2 − x∆+2 − x2∆+2 ) .. .

      ∆λΘ(1 − x∆ − x2∆ − x3∆ )  εx1   εx2  F(x) =  ..  .   εx∆   0   0   ..  . 0





(η + ε)x1 (η + ε)x2 .. .

            (η + ε)x∆     γx∆+1    −  γx∆+2  , V (x) =  ..   .       γx∆+∆   ζx   2∆+1   ζx 2∆+2     ..   . ζx2∆+∆





 0   0     ..     .   0     0    +  0  , V (x) =    ..   .     0     ηx + γx 1 ∆+1     ηx + γx 2 ∆+2     ..  . ηx∆ + γx∆+∆



           .           

It is easy to verify that the functions satisfy assumptions (A1)–(A5). (A1) If x ≥ 0, then Fi , Vi− , Vi+ ≥ 0 for i = 1, . . . , 3∆. (A2) If xi = 0, then Vi− = 0 for any i ∈ {1, . . . , 3∆}. In particular, if x ∈ Xs , then Vi− = 0 for i = 1, . . . , 2∆.

17 (A3) Fi = 0 if i > 2∆. (A4) If x ∈ Xs , then Fi = 0 and Vi+ = 0 for i = 1, . . . , 2∆. (A5) If F(x) is set to zero, then all eigenvalues of Df(x0 ) have negative real parts. For the assumption (A5), if F(x) is set to zero, then we have   −(η + ε)E∆ 0 0 0 −γE∆ 0 , Df(x0 ) =  ηE∆ γE∆ −ζE∆

which obviously implies that the matrix Df(x0 ) has three different negative eigenvalues −(η + ε), −γ, and −ζ, each with multiplicity of ∆. Thus, the above (A5) holds.     (η + ε)E∆ 0 0 B ,V = , Then, it follows by the above result (see also Lemma 1 in Ref. [25]) that F = 0 γE∆ εE∆ 0 λ where B = (Bij )∆×∆ , Bij = hki ijP (j), 1 ≤ i, j ≤ ∆. Direct calculations yield that 

 F V −1 = 

0 ε E∆ η+ε

 λ M∆  γhki . 0

It follows by Lemma IV.1 and certain calculations that the characteristic equation of the matrix F V −1 is λ M xE − ∆ ∆
γhki det xE2∆ − F V −1 = ε − xE∆ η + ε E∆   ελ M∆ = det x2 E∆ − γ(η + ε) hki ∆    γ(η + ε) hki 2 ελ det M∆ − x E∆ = − γ(η + ε) hki ελ ∆  ∆−1  γ(η + ε) hki 2 ελ − x = −  γ(η + ε) hki  ελ

 γ(η + ε) hki 2 × k2 − x ελ   γ(η + ε) hki 2 2  ελ = 0. x2(∆−1) x − k = γ(η + ε) hki ελ Thus, the propagation threshold of model (2) is the spectral radius of the matrix F V −1 , i.e.,

R0 = ρ(F V −1 ) =

s

hk 2 i λε . hki γ(ε + η)

References

[1] 2016 Internet Security Threat Report (ISTR), Volume 21, April 2016. https://www.symantec.com/security-center/ threat-report (Date of access: 24th March, 2017) [2] S. Funk, E. Gilad, C. Watkins, and V. A. A. Jansen, The Spread of Awareness and Its Impact on Epidemic Outbreaks, Proc. Natl. Acad. Sci. U.S.A. 106, 6872 (2009). [3] M. Sandro, P. Nicola, A. Alex, G. Sergio, M. Yamir, and V. Alessandro, Modeling Human Mobility Responses to the Large-Scale Spreading of Infectious Diseases, Sci. Rep. 1, 62 (2011). [4] S. Funk, M. Salath´e, and V. A. A. Jansen, Modelling the Influence of Human Behaviour on the Spread of Infectious Diseases: A Review, J. R. Soc. Interface 7, 1247 (2010).

18 [5] K. Shi, Y. Tang, X. Liu, S. Zhong, Non-fragile sampled-data robust synchronization of uncertain delayed chaotic Lurie systems with randomly occurring controller gain fluctuation, ISA Transactions 66, 185–199 (2017) [6] B. Yan, X. Zhou, J. Cheng, F. Lang, Finite-Time Boundedness of Markov Jump System with Piecewise-Constant Transition Probabilities via Dynamic Output Feedback Control, Mathematical Problems in Engineering, 2015, Article ID 327947 (2015) [7] X. Liu, K. Zhang, Existence, uniqueness and stability results for functional differential equations on time scales, Dynamic Systems and Applications, 25(4), 501–530 (2016) [8] R. Pastor-Satorras, C. Castellano, P. Van Mieghem, and A. Vespignani, Epidemic Processes in Complex Networks, Rev. Mod. Phys. 87, 925 (2015). [9] W. Liu, C. Liu, X. Liu, A discrete dynamic model for computer worm propagation. In: Proceedings of the mathematics & statistics, vol. 150, pp. 119-131 (2015) [10] C. Gan, M. Yang, Z. Zhang, W. Liu, Global dynamics and optimal control of a viral infection model with generic nonlinear infection rate, Discrete Dynamics in Nature and Society, vol. 2017, Article ID 7571017 (2017) [11] M. E. J. Newman, Threshold Effects for Two Pathogens Spreading on a Network, Phys. Rev. Lett. 95, 108701 (2005). [12] W. Liu, C. Liu, X. Liu, S. Cui, and X. Huang, Modeling the spread of malware with the influence of heterogeneous immunization. Appl. Math. Model. 40(4), 3141-3152 (2016). [13] W. Liu, and S. Zhong, Web malware spread modelling and optimal control strategies, Sci. Rep., 7, Article number: 42308 (2017). [14] W. Willinger, D. Alderson, J. Doyle, Mathematics and the Internet: A Source of Enormous confusion and Great Potential, Notices of the American Mathematical Society, 56(5), 586-599 (2009) [15] A. Clauset, C.R. Shalizi, M.E.J. Newman, Power-law distributions in empirical data, SIAM Review, 51, 661-703 (2009) [16] R. Pastor–Satorras, and A. Vespignani, Epidemic Spreading in Scale-Free Networks, Phys. Rev. Lett., 86(14), 3200–3203 (2001) [17] W. Liu, C. Liu, Z. Yang, X. Liu, Y. Zhang, and Z. Wei, Modeling the propagation of mobile malware on complex networks. Commun. Nonlinear Sci. 37, 249-264 (2016). [18] G.F. Arruda, E. Cozzo, T.P. Peixoto, F.A. Rodrigues, and Y. Moreno, Disease Localization in Multilayer Networks, Phys. Rev. X 7, 011014 (2017). [19] A.V. Goltsev, S.N. Dorogovtsev, J.G. Oliveira, and J.F.F. Mendes, Localization and Spreading of Diseases in Complex Networks, Phys. Rev. Lett. 109, 128702 (2012). [20] L.-X. Yang, X. Yang, and Y. Wu, The impact of patch forwarding on the prevalence of computer virus: A theoretical assessment approach, Appl. Math. Model. 43, 110-125 (2017). [21] F. Radicchi, and G. Bianconi, Redundant Interdependencies Boost the Robustness of Multiplex Networks, Phys. Rev. X 7, 011013 (2017). [22] R. Albert, H. Jeong, A.-L. Barab´ asi, Diameter of the world wide web, Nature 401, 130-131 (1999). [23] B.A. Huberman, L.A. Adamic, Growth dynamics of the World-Wide Web, Nature 401, 131 (1999). [24] A. Lajmanovich, and J.A. Yorke, A deterministic model for gonorrhea in a nonhomogenous population, Math. Biosci. 28(3-4) 221-236 (1976). [25] P. van den Driessche, and J. Watmough, Reproduction numbers and sub-threshold endemic equilibria for compartmental models of disease transmission, Math. Biosci., 180(1–2) (2002) 29–48