Modeling and control of computer virus attack on a targeted network

Physica A 538 (2020) 122617

Contents lists available at ScienceDirect

Physica A journal homepage: www.elsevier.com/locate/physa

Modeling and control of computer virus attack on a targeted network Ranjit Kumar Upadhyay a , Prerna Singh b , a b

∗

Department of Mathematics & Computing, Indian Institute of Technology (Indian School of Mines), Dhanbad 826004, India School of Mathematics & Statistics, University of Sheffield, Sheffield S10 2TN, United Kingdom

article

info

Article history: Received 14 March 2019 Received in revised form 25 June 2019 Available online 8 October 2019 Keywords: Distributed attack Firewall Optimal control Sensitivity analysis Targeted network

a b s t r a c t In this paper, a model with two different frameworks of an attacking class and a targeted class is introduced to study the virus mobility of an attack in a targeted network. Existence and stability of equilibrium points are discussed in equivalence with the basic reproduction number R0 . It is observed that as R0 crosses unity, transcritical bifurcation arises in the system. The effect of firewall security is also shown graphically which has been taken as a media coverage factor in this work and it is discovered that firewall helps to bring down the virus propagation by reducing the infection peak. Optimal control concept is then introduced to propose one more measure for regulating the virus propagation. Numerical experimentations have been done to rationalize the analytical conclusions. In the end, a sensitivity analysis is executed that provides the information about relevance of parameters in deciding the virus spread in networks. © 2019 Published by Elsevier B.V.

1. Introduction Dynamics of computer viruses have always drawn resemblance with modeling of biological epidemiology [1]. Keeping in mind the understanding of infection spread in computer networks, some mathematical models have been discussed by several authors [2,3]. Cyber attack is now one of the biggest threat to economy and on personnel levels. Malignant codes not just simply influence a single PC, it may likewise get into whole of the system and can spread. It is able to send messages via email and obtain data or cause more harm by erasing the records. A type of malignant code is called the infection, which is a small program adjoining to various documents and will replicate itself in one PC and then spread to several other PCs. Infections can thus go from being little dangerous to significantly harmful for a network. All this had lead to an urgent call to deal with these computer viruses. The infections in computer networks are complex. However, unlike contagious diseases, viruses in a computer or computer networks can spread and cause destruction within a few seconds. Mathematical modeling has proved to be an important tool in the analysis of spread and control of virus in computer networks. Such models take into consideration the key factors that administer the virus spread, and prognosticate how the infection will profuse over time period. One of the serious threats of the time is the Distributed Denial of Service (DDoS) attack [4]. According to the 2010 report prepared by Arbor Networks, there is an alarming 102 percent increment of DDoS attacks in 2010 in comparison to 2009 [5]. In a classic DDoS attack, a huge number of attacking hosts are assembled to send meaningless packets to trouble a victim, or his Internet connection, and sometimes even both. In these few years, it has been observed that DDoS attack mechanisms are becoming more refined, practical, and more challenging to trace [6]. A Denial of Service attack that aims to stop a rightful ∗ Corresponding author. E-mail addresses: [email protected] (R.K. Upadhyay), [email protected] (P. Singh). https://doi.org/10.1016/j.physa.2019.122617 0378-4371/© 2019 Published by Elsevier B.V.

2

R.K. Upadhyay and P. Singh / Physica A 538 (2020) 122617

user to receive or use a regular service from internet can be carried out in several ways. A DDoS attack chiefly consists of four factors as was discussed in [7] and [8]. Initial one is the target host that takes in the hardship of attack. Second is the attacking agent, the program that carries out the attack on specific targeted resources. Another component is the control master program that regulates the attack. Last and fourth one is the actual attacker that enforces the attack [9,10]. Recent advancements in DDoS defense techniques have marked an end to the time when intruders could download some tool and carry out an attack against almost any website. In DDoS attacks nowadays, attackers make use of highly complicated mechanism to launch the attacks. Even with so much efforts towards lessening the sum of DDoS attack cases, they have increased to a great extent in frequency and also in size of the targeted networks and computers [11]. On one hand, carrying out a DDoS attack has been made relatively easier by the presence of various user-friendly attack means. However, there is yet lack of adequate solutions to guard against them in terms of interrupting a continuing attack and tracing backwards to the attacking sources. Thus, DDoS problems are likely to become more serious in the upcoming future. For example, they can be used in cyber battle to disrupt business, governmental agencies, industries, or military sites [6]. Previously, some efforts have been done to mathematically analyze and understand such type of attacks [12,13]. However, the affiliation between distributed methods of attack and targeted attacks still needs to be examined, keeping in mind the scene that targeted attacks are now becoming much more recurrent nowadays, and seriousness of selected targets is also expanding. In accordance with Symantec reports, overall balance of reported targeted attacks has heightened from 77 cases in 2010, to 82 in 2011, and then to 116 cases in 2012. Another quite famous targeted attack was called Stuxnet attack, that was endorsed in 2010 [2]. In this paper, an attempt has been made to obtain a mathematical model so as to study the dynamics of distributed attacks on fixed targeted resources for the success or failure of these attacks in leading to consequential damage of the resources. We separate targeted resources from the hosts used to fire the distributed attack. We propose an epidemic model that consists of a transmission of infection in one direction from particular hosts towards the targeted networks. On the same lines as in [2], we take into consideration both the chances of classifying new hosts for widening the spread of attack as well as causing damage to as many targets as possible. The targeted systems do not participate in spreading the attack, unlike the classical epidemiological models. We obtain a threshold on which the dynamics of the model depends, known as reproduction number. We also discuss the effect of firewall and optimal control exclusively on our model system of distributed attack and that has not been done so far. Section 2 presents the model formulation of the system. Analysis of threshold and equilibrium points have been obtained and discussed in Section 3. Section 4 comprises of numerical simulations, bifurcations analysis and firewall effect. In Section 5, theory of optimal control is introduced and its effect in reducing the virus prevalence is discussed. In Section 6, a sensitivity analysis is performed to ascertain the relative significance of every parameter in governing the dynamics of the system. Finally, Section 7 concludes the paper with results and discussions. 2. Model configuration A network model consisting of two separate classes of computer nodes, namely attacking class and a targeted class has been formulated. Attacking class launches the attack whereas targeted class is the victim of attack. Moreover, targeted class play no role in spread of infection. There is a further division of attacking and targeted systems into three classes. Attacking class is split into Susceptible(St ), Infectious(It ), and Recovered(Rt ) targeted nodes whereas targeted class is divided into Susceptible(Sa ), Infectious(Ia ), and External(Ea ) attacking nodes. The total targeted and attacking node populations at any instant of time are Nt and Na such that Nt = St + It + Rt , Na = Sa + Ia + Ea . Before moving towards the model formulation, we briefly discuss the underlying basic assumptions. The basic aim of distributed type of attack is to collect as many hosts as possible to launch the attack and carry out the attack at a large scale. The total number of targeted nodes and total number of attacking nodes remains constant overall. Recovered targeted nodes move to susceptible compartment at a constant rate. Crashing rate for all the attacking nodes is also assumed to be a constant. The class of external attacking nodes consists of those susceptible and infected attacking nodes which are removed from the internet. This external node compartment is also assumed to have nodes recruitment from recovered crashed nodes after certain point of time. On the same lines as in [10], we introduce a new mathematical model with the certain presumptions: (i) η > 0 is the rate at which each attacking node (susceptible/infectious) is detached from internet connection and join external node compartment. (ii) Each node from external attacking class is joined to the internet to become susceptible or infectious node(attacking) at the rate k1 > 0 or k2 > 0. (iii) We consider the crashing rates for each attacking node to be µ > 0. (iv) Each targeted node from susceptible class gets infected by an infectious attacking node at rate β1 > 0. (v) Each attacking node from susceptible class gets infected by an infectious attacking node at rate β2 > 0.

R.K. Upadhyay and P. Singh / Physica A 538 (2020) 122617

3

Fig. 1. Simplified diagram of the proposed model.

(vi) As the immunity is not assumed to be permanent, recovered targeted node becomes susceptible again at some constant rate α > 0. (vii) Recruitment rate for external attacking class is taken to be ξ . (viii) The recovery rate of infected targeted nodes is taken to be γ . The transference between various nodes from different classes is revealed by the following sets of the nonlinear differential equations: Targeted Class dSt dt dIt dt dRt dt

= −β1 St Ia + α Rt , = β1 St Ia − γ It ,

(1)

= γ It − α Rt ,

with the initial conditions St (0) > 0, It (0) ≥ 0 and Rt (0) ≥ 0 and St (t) + It (t) + Rt (t) = 1. Attacking Class dSa dt dIa dt dEa dt

= −β2 Sa Ia + k1 Ea − ηSa − µSa , = β2 Sa Ia + k2 Ea − ηIa − µIa ,

(2)

= ηSa + ηIa − k1 Ea − k2 Ea − µEa + ξ ,

where Sa (t) + Ia (t) + Ea (t) = 1, with the initial conditions Sa (0) > 0, Ia (0) ≥ 0 and Ea (0) ≥ 0. Each parameter of targeted system (1) and attacking system (2) is positive and are illustrated in Table 1. Refer Fig. 1 for the model diagram. Systems (1) and (2) can be shortened to a system which we will further investigate is given by: dSt dt dIt dt dIa dt dEa dt

= −β1 St Ia + α (1 − St − It ), = β 1 S t Ia − γ It , (3)

= β2 (1 − Ia − Ea )Ia − µIa + k2 Ea − ηIa , = ξ + η(1 − Ea ) − (µ + k1 + k2 )Ea ,

4

R.K. Upadhyay and P. Singh / Physica A 538 (2020) 122617 Table 1 Characterization of parameters and variables. Variable/Parameter

Description

St It Rt Sa Ia Ea

Targeted susceptible nodes Targeted infectious nodes Targeted recovered nodes Attacking susceptible nodes Attacking infectious nodes Attacking external nodes Infection transmission rate for targeted class. Infection transmission rate for attacking class. Recovery rate of infected targeted nodes. Natural crashing rate for all the nodes of attacking system. Recovery rate for susceptible targeted nodes. Rate by which nodes from Sa or Ia class are detached from internet to join Ea class. Rate by which nodes from Ea class are attached to internet and join Sa class. Rate by which nodes from Ea class are attached to internet and join Ia class. Recruitment rate for external attacking nodes.

β1 β2 γ µ α η k1 k2

ξ

The operable region for this system (3) is given below:

τ = {(St , It , Ia , Ea ) ∈ R4 : St > 0, It ≥ 0, Ia ≥ 0, Ea ≥ 0, St + It ≤ 1, Ia + Ea ≤ 1}. We can easily observe that the region τ is positively invariant with respect to system (3), and the system has bounded solutions as well. We will study the above system for some dynamical properties in the upcoming sections.

3. Stability analysis By solving the equations

− ηξ

dSt dt

= 0,

dIt dt

= 0,

dIa dt

= 0,

dEa dt

= 0 of system (3), we get E0 = (− ηξ , 0, 0, 0) as one of the

equilibrium point, where = 1. Since all the infection terms are at zero level, we call it the infection free equilibrium of the system. We now calculate the Basic reproduction number R0 which is defined as the number of secondary infectious cases produced in a susceptible class by a single infectious node in its overall infectious duration [14]. This is a decisive threshold that can be used to predict whether the infection will exist in a network or it will diminish with time. Now using the approach of Next generation matrix [15], we formulate F and V as follows and extract the formulas of reproduction number for both the targeted and attacking systems respectively. For targeted class, F is derived from second equation of system (3) by choosing the infectious terms and V is derived using the remaining terms. The non-negative matrices F and V are given below:

] β1 St Ia , β2 (1 − Ia − Ea )Ia

[ F =

γ It

[ V =

]

µIa − k2 Ea + ηIa

,

Now F and V are constructed by obtaining the derivatives with respect to Ia and It respectively, with their values obtained at infection-free equilibrium E0 , and are given as:

[ F =

V =

0 0

] β1 , β2

[ γ

µ+η

0

[ FV

−1

=

]

0

0

0

β1

,

]

(η+µ)

β2 η+µ

Now the effective reproduction number for whole system is given by the spectral radius of FV −1 i.e ρ (FV −1 ), R0 = ρ (FV −1 ) =

β2 . η+µ

(4)

R.K. Upadhyay and P. Singh / Physica A 538 (2020) 122617

5

Jacobian matrix is given by:

−(β1 Ia + α ) β 1 Ia ⎢ J =⎣ 0

⎡

0

−β1 St β1 St β2 (1 − 2Ia − Ea ) − µ − η

−α −γ 0 0

⎤

0 0

⎥

⎦ k2 − β2 Ia −(µ + k1 + k2 + η)

0

.

Jacobian matrix at E0 = (1, 0, 0, 0) is obtained as follows:

⎡ −α ⎢ 0 J0 = ⎣ 0

−α −γ 0 0

0

⎤

0 0 k2

−β1 β1 β2 − µ − η

⎥ ⎦.

−(µ + k1 + k2 + η)

0

Now, eigenvalues of the Jacobian matrix(upper triangular) J0 , are given as follows:

λ1 λ2 λ3 λ4

= = = =

−α, −γ , β 2 − (µ + η ) , −(µ + k1 + k2 + η),

(5)

Theorem 3.1. The infection-free equilibrium E0 (1, 0, 0, 0) is locally asymptotically stable for R0 < 1, and goes unstable for R0 > 1. β

2 Proof. For R0 = µ+η < 1, the eigenvalue λ3 < 0, hence we can see from (5) that all the eigenvalues of Jacobian matrix J0 are negative. So, equilibrium state E0 is locally asymptotically stable for R0 < 1. On the other hand, for R0 > 1, λ3 > 0, i.e., positive. Hence IFE, E0 is unstable for R0 > 1. □

β

2 = 1, we have β2 = µ + η. Then taking β2 as the bifurcation parameter at R0 = 1, we get λ3 = 0. At R0 = 1, i.e., µ+η Hence J0 must have a zero eigenvalue for R0 = 1. In this case, E0 (1, 0, 0, 0) is a non hyperbolic equilibrium point.

3.1. Endemic equilibrium point To determine another equilibrium point, we start with obtaining the values of Ia∗ and Ea∗ by solving the last two equations of system (3). The obtained values are shown below: Ia∗ =

−(η + µ)2 − (η + µ)(k1 + k2 ) + β2 (µ − ξ + k1 + k2 ) + 2β2 (η + µ + k1 + k2 )

√ b

,

where b = (η + µ)2 (η + µ + k1 + k2 )2 − 2(η + µ + k1 + k2 )((η + µ)(µ − ξ + k1 ) + (−η + µ − 2ξ )k2 )β2 + (µ − ξ + k1 + k2 )2 β22 , and Ea∗ =

η+ξ . η + µ + k1 + k2

A fundamental condition for Ia∗ to exist at a non-negative level is Then, solving for It∗ from system (3), we have, It∗ = −

√

b + β2 (µ − ξ + k1 + k2 ) > (µ + η)2 + (µ + η)(k1 + k2 ),

αβ1 ((η + µ)2 + (η + µ)(k1 + k2 ) − (µ − ξ + k1 + k2 )β2 −

√ b)

2αγ (η + µ + k1 + k2 )β2 + (α + γ )β1 (−(η + µ)2 − (η + µ)(k1 + k2 ) + (µ − ξ + k1 + k2 )β2 + ∗

∗

√ . b)

∗

Now, with these values of It and Ia , we obtain St which is given by: St∗ =

2αγ (η + µ + k1 + k2 )β2 2αγ (η + µ + k1 + k2 )β2 + (α + γ )β1 (−(η + µ)2 − (η + µ)(k1 + k2 ) + β2 (µ − ξ + k1 + k2 ) +

√ . b)

And a sufficient condition for St and It to exist at a positive level is (η + µ) + (η + µ)(k1 + k2 ) ≤ β2 (µ − ξ + k1 + k2 ) + Hence, the endemic equilibrium E ∗ (St∗ , It∗ , Ia∗ , Ea∗ ) exists only if the above quoted conditions are satisfied. ∗

∗

2

√ b.

Theorem 3.2. The endemic equilibrium point E ∗ (St∗ , It∗ , Ia∗ , Ea∗ ) is locally asymptotically stable when effective reproduction number, R0 > 1. Proof. To prove this theorem, we apply the approach established on Central Manifold theory, taking β2 as the bifurcation parameter [16]. At R0 = 1, value of the bifurcation parameter is given by

β2 = µ + η = βc .

6

R.K. Upadhyay and P. Singh / Physica A 538 (2020) 122617

The Jacobian matrix J, at IFE E0 and at β2 = βc is given by

⎡ −α ⎢ 0 J0 = ⎣ 0

−α −γ

−β1 β1

0 0

0 0

0

0 0 k2

⎤

−(µ + k1 + k2 + η)

⎥ ⎦.

Now, the eigenvalues of J0 are −α , −γ , 0 and −(µ + k1 + k2 + η). Clearly, J0 has a zero now, corresponding ( eigenvalue )

to which, the right eigenvector is given by W = (w1 , w2 , w3 , w4 )T , where w1 = −w3

w3 is arbitrary. Let w3 = 1, then )⎤ ( ⎡ − βα1 γα + 1 ⎢ ⎥ ⎢ ⎥ β1 ⎢ ⎥ γ W =⎢ ⎥, ⎢ ⎥ 1 ⎣ ⎦

β1 α

α γ

+ 1 , w2 =

w3 β1 , γ

w4 = 0 and

0 such that det(J0 − 0.I).W = [0]4×4 . Furthermore, the left eigenvector corresponding to zero eigenvalue, V = (v1 , v2 , v3 , v4 ) must meet the equalities V .J0 = 0, V .W = 1. Thereafter,

( )⎤ ⎡ − βα1 γα + 1 ⎥ ⎢ ⎥ ⎢ β1 ⎥ ⎢ γ V .W = 1 ⇒ (v1 , v2 , v3 , v4 ) ⎢ ⎥ = 1, ⎥ ⎢ 1 ⎦ ⎣ 0

(

⇒ v1 − and, V .J0 = 0

β1 α

(

α +1 γ

))

⎡ −α ⎢ 0 ⇒ (v1 , v2 , v3 , v4 ) ⎣ 0 0

+ v2 .

β1 + v3 .1 + v4 .0 = 1. γ

−α −γ

−β1 β1

0 0

0 0

0 0 k2

⎤ ⎥ ⎦ = 0.

−(µ + k1 + k2 + η)

Upon solving the above systems for V , we obtain v1 = 0, v2 = 0, v3 = 1, and v4 =

k2

µ+η+k1 +k2

(

. So, V = 0, 0, 1,

k2 k1 +k2 +η+µ

)

is the left eigenvector corresponding to zero eigenvalue. Denote x1 = St , x2 = It , x3 = Ia , and x4 = Ea and then we define, a=

4 ∑

vk wi wj

k,i,j=1

b=

4 ∑ k,i=1

vk wi

∂ 2 fk (0, 0) , ∂ xi ∂ xj

∂ 2 fk (0, 0) , ∂ xi ∂β2

where fk = kth component of f . Non-zero derivatives are given by

∂ 2 f3 = − 2β2 , ∂ x23

∂ 2 f3 = −β2 , ∂ x2 ∂ x4

∂ 2 f3 ∂ 2 f1 = − β2 , = −β1 , ∂ x4 ∂ x3 ∂ x1 ∂ x3 2 2 ∂ f2 ∂ f3 =β1 , = −β1 , ∂ x1 ∂ x3 ∂ x3 ∂ x1 ∂ 2 f3 =1. ∂ x3 ∂β2

(6)

(7)

R.K. Upadhyay and P. Singh / Physica A 538 (2020) 122617

7

Inserting all second order derivatives calculated at infection free equilibrium E0 in (6) and (7), and for β2 = βc , we get a = −β2 v3 (v3 w32 + w4 ), b = v3 w 3 . Lastly, substituting the above calculated values of V and W , we get a = −β2 ,

(8)

b = 1.

(9)

As we can see that a < 0 and b > 0 at β2 = βc , therefore applying theorem 4.1 from [16], we ascertain that transcritical bifurcation appears at R0 = 1. At β2 = βc or R0 = 1, a non hyperbolic equilibrium point appears and when R0 crosses 1, another equilibrium(endemic), turns up which is unstable for R0 < 1, and grow into locally stable when R0 > 1. □ 3.2. Global stability of endemic equilibrium We use Li and Muldowney’s geometric approach [17] to establish the global stability of endemic equilibrium. Since the dimension of matrix n > 2, this method will be easier to approach to prove stability. Theorem 3.3. For R0 > 1, the endemic equilibrium is globally asymptotically stable for the region inside τ . Proof. If R0 > 1, endemic equilibrium is locally asymptotically stable and unstable for R0 < 1 from Theorem 3.2. Also, from Theorem 3.1, infection-free equilibrium E0 is unstable for R0 > 1, which suggests that system (3) is uniformly persistent in τ [18]. This implies that there must exist some constant c > 0 so that for any initial point (St (0), It (0), Ia (0), Ea (0)) ∈ τ , each solution of system (3) in τ satisfies min[lim inf St (t), lim inf It (t), lim inf Ia (t), lim inf Ea (t)] ≥ c . t →∞

t →∞

t →∞

t →∞

By the theory of Li and Muldowney [17], if x → f (x) ∈ Rn is a continuously differentiable function in an open subset E of Rn and x′ = f (x), then (i) E is a simply connected set. (ii) There must exist some compact absorbing set D in E. (iii) The system x′ = f (x) has some unique equilibrium in E. As the endemic equilibrium is locally asymptotically stable, it will also be globally stable provided (i) and (ii) holds. Now the Jacobian associated with endemic equilibrium point E ∗ (St∗ , It∗ , Ia∗ , Ea∗ ) is given by

⎡−β I ∗ − α 1 a ⎢ β1 Ia∗ J =⎢ ⎣ 0

0

−β1 St∗ β1 St∗ ∗ β2 (1 − Ia − Ea∗ ) − µ − η

0

0

−α −γ

0

0

⎤

⎥ 0 ⎥. −β2 Ia∗ − k2 ⎦ −k1 − k2 − η − µ

The second compound matrix of Jacobian is given by

⎡

J

[2]

a11 ⎢0 ⎢ ⎢0 =⎢ ⎢0 ⎣0 0

β1 St∗

0

a22 0 β1 Ia∗ 0 0

−β2 Ia − k2

−β1 St∗ −α

0 0

a33 0 β1 Ia∗ 0

0 a44 0 0

−α −β2 Ia∗ − k2

∗

a55 0

0 0

⎤

⎥ ⎥ −β1 St∗ ⎥ ⎥, 0 ⎥ ∗ ⎦ β1 S t a66

where a11 = −β1 Ia∗ − α − γ , a22 = −β1 Ia∗ − α + β2 (1 − 2Ia∗ − Ea∗ ) − µ − η, a33 = −β1 Ia − α − η − k1 − k2 − µ, ∗

a44 = β2 (1 − 2Ia∗ − Ea∗ ) − µ − η − γ ,

a55 = −γ − η − k1 − k2 − µ, a66 = β2 (1 − 2Ia∗ − Ea∗ ) − 2µ − 2η − k1 − k2 . To obtain matrix A according to the approach, we shall define a 6 × 6 diagonal matrix P as

( P = diag

1,

It

,

It

Ia Ia

,...,

It Ia

)

.

Let Pf be the matrix obtained by taking directional derivative of each diagonal element of P.

8

R.K. Upadhyay and P. Singh / Physica A 538 (2020) 122617

Then, Pf P −1 = diag

(

0,

It′ It

−

Ia′ Ia

It′

,...,

It

−

Ia′

)

.

Ia

Hence matrix A is given as A = Pf P −1 + PJ [2] P −1 = Pf P −1 + J [2] , such that

[

]

A11 A= A21

A12 , A22

⎡ ⎤ 0

where A11 = −β1 Ia∗ − α − γ , A12 = β1 St∗

[

]

[

0

−β1 St∗

]

0

0 , A21

⎢0⎥ ⎢ ⎥ = ⎢0⎥ and ⎣0⎦ 0

⎡ A22

b11

∗

−β2 Ia − k2

−α

0

b22

0

0

b33

−α −β2 Ia∗ − k2

β 1 Ia

0

b44

0

0

0

⎢ 0 ⎢ ⎢ = ⎢β1 Ia∗ ⎢ ⎣ 0

∗

0

0

−β1 St∗ ⎥ ⎥ ⎥ 0 ⎥, ⎥ β1 St∗ ⎦ b55

where b11 = β1 Ia∗ − α + β2 (1 − Ia∗ − Ea∗ ) − µ − η + I′

⎤

I′

It′ It

−

Ia′ , Ia

b22 = β1 Ia∗ − α − η − µ − k1 − k2 +

I′

I′

It′ It

− I′

Ia′ , Ia Ia′ . Ia

b33 = β2 (1 − Ia∗ − Ea∗ ) −µ−η−γ + It − Ia , b44 = −µ−η−γ − k1 − k − 2 + It − Ia , b55 = β2 (1 − Ia∗ − Ea∗ ) − 2µ− 2η− k1 − k2 + It − t a t a t Let µ be the Lozinskii measure with respect to L1 norm. Then from assumptions (ii) and (iii) and using the theory from [19], we have the result that if µ(Pf P −1 + PJ [2] P −1 ) < 0, then follows the global stability of endemic equilibrium. We also have µ(A) ≤ sup{l1 , l2 }, where l1 = µ1 (A11 ) + |A12 |

(10)

l2 = |A21 | + µ1 (A22 )

(11)

where |A12 |, |A21 | are matrix norms with respect to L1 vector norm and µ1 (A11 ), µ1 (A22 ) have been obtained by taking the maximal of sums retrieved by adding the absolute value of non-diagonals in one column with the diagonal element in that column. So, we have

µ1 (A11 ) = −β1 Ia∗ − α − γ

(12)

|A12 | = β1 St∗

(13)

|A21 | = 0

(14)

And µ1 (A22 ) is calculated as

µ1 (A22 ) = max{2β1 Ia∗ − α + β2 (1 − Ia∗ − Ea∗ ) − µ − η + ′

α + β2 (1 − Ia∗ − Ea∗ ) − µ − η − γ +

It It

It′ It

−

Ia′ Ia

, β2 Ia∗ + 2β1 Ia∗ − α − µ − η − k1 +

′

−

Ia Ia

′

, α + β2 Ia∗ + k2 − µ − η − γ − k1 − k2 − 2 + ′

2β1 St∗ + β2 (1 − Ia∗ − Ea∗ ) − 2µ − 2η − k1 − k2 +

µ1 (A22 ) = 2β1 Ia∗ + β2 (1 − Ia∗ − Ea∗ ) − α − µ − η +

It It

It′ It

It It

It′ It

−

Ia′ Ia

,

′

−

Ia Ia

,

′

−

Ia

−

Ia′

Ia

Ia

}

,

(15)

if 2β1 Ia∗ + β2 Ia∗ − α − k1 > α − γ + β2 (1 − Ia∗ − Ea∗ ) and β2 Ia∗ − k1 > β2 (1 − Ia∗ − Ea∗ ). From system (3), we have It′ It Ia′ Ia

=

β1 St Ia It

− γ,

= β2 (1 − Ia − Ea ) − µ +

(16) k2 Ea Ia

− η.

Substituting Eqs. (12)–(15) in (10) and (11) respectively and using (16) and (17), we get l1 = β1 St∗ − β1 Ia∗ − α − γ ≤

It′ It

− ϵ,

(17)

R.K. Upadhyay and P. Singh / Physica A 538 (2020) 122617

9

Fig. 2. Time Series plot for the behavior of nodes, for the case (a) R0 = 0.5 < 1 at β2 = 0.01, ξ = 0 and η = 0 and (b) R0 = 2.8125 > 1, at β2 = 0.9. Rest of the parameters take the same value as given in (18).

l2 = β2 Ia∗ + 2β1 Ia∗ − α − µ − η − k1 +

It′ It

−

Ia′ Ia

It′

≤

It

− ϵ,

I′

for some ϵ > 0, when β2 Ia∗ + 2β1 Ia∗ < α + µ + η + k1 + Ia . Then we must have, µ(A) ≤ sup{l1 , l2 } according to the method. a

Now, a quantity q¯2 is defined as follows q¯2 = lim sup sup t →∞

x0 ∈τ

t

∫

1 t

µ(A(x(s, x0 )))ds. 0

Now, 1 t

t

∫

µ(A)ds ≤ 0

=

1

∫

t

(

It′

t 0 It logIt (t)

⇒ q¯2 < 0.

t

− ϵ )dt −ϵ

Hence the global stability of endemic equilibrium is verified.

□

4. Numerical simulations Now the authentication of our mathematical findings will be shown using numerical simulations. Existence, stability of both the steady states, and different dynamical properties of the model such as occurrence of transcritical bifurcation has been verified for the system (3) using MATLAB. The given set of parameters is been assumed for the experimentations: m = 0.2, β1 = 0.8, α = 0.3, η = 0.3, γ = 0.4, k1 = 0.2, k2 = 0.4, ξ = 4.

(18)

In Fig. 2, we have plot the time series by varying the values of different parameters to analyze the stability of both the equilibrium points and verifying all the theoretical conditions corresponding to each point. The system (3) is simulated taking initial node population sizes as St0 = 0.7, It0 = 0.1, Ia0 = 0.6, and Ea0 = 0.3. At β2 = 0.01, the value of effective reproduction number R0 is 0.5 (i.e, R0 < 1). We can see from Fig. 2(a) that the nodes approaches (1, 0, 0, 0) as the time increases. Whereas, for β2 = 0.9, we have R0 = 2.8125 > 1 and it is clear from Fig. 2(b) that the nodes move closer to unique endemic equilibrium E ∗ (0.3167, 0.2929, 0.4625, 4.6739), thus establishing the stability of endemic steady state for R0 > 1. 4.1. Bifurcation analysis Now we present in Fig. 3 the bifurcation diagram for the model that exhibits a forward transcritical bifurcation at R0 = 1. Remaining parameters are same as mentioned in Eq. (18). We observed that when reproduction number is less than 1, the system remains in infection free condition, but as it goes beyond 1, the endemic environment takes place i.e., the figure shows bifurcation at R0 = 1 along with the appearance of a positive(endemic) equilibrium as soon as R0 crosses 1. Thereby, Fig. 3 validates the accuracy of our analytical findings.

10

R.K. Upadhyay and P. Singh / Physica A 538 (2020) 122617

Fig. 3. Curve showing the transcritical bifurcation for system (3) at η = 0 and ξ = 0. Parameters are same as mentioned in (18).

Fig. 4. Effect of firewall security term m on the quantity of infectious targeted nodes, for (a) β2 = 0.01 and R0 < 1, (b) β2 = 0.9 and R0 > 1. Other parameters are same as given in (18).

4.2. Firewall effect Compelled by the need for rapid employment of defensive measures against DDoS and other type of targeted attacks, the notion of ‘‘Internet firewall’’ was suggested to secure the internet [6]. The main reason to employ firewall is to attain superior network protection and topmost user ease. Firewalls are typically employed at the edge of networks or at the entry point of some private network. Further, the internet traffic(incoming or outgoing) is investigated by the firewalls. In simpler way, firewall works as a system that manipulates data flow between two networks. Then employed firewall subsequently permits the following transmissions from one network to other [20]. We acknowledge firewall security as the media coverage factor in this model [10,21,22] and [23]. We employ firewall security coefficient in the transmission rate term of the model system (3) as β1 (t) = β1 e−mIa . Fig. 4 shows the effect of firewall security term m on the model dynamics. As the value of m is increased, the infection peak can be seen to be going down. So, we can say that firewall helps to control the virus proliferation to a symbolic level. In Fig. 4(a), for β2 = 0.01 and R0 < 1, the targeted infectious nodes tends to zero as m is increased, thus confirming the negative impact of m on the virus spread. On the other hand, for R0 > 1, the nodes seem to be approaching a finite value with time and in this case, we can see from Fig. 4(b) that higher the value of m, lower is the infection peak and thus firewall reduces the infection level up to some extent. In Table 2, simulated data is given for the case of a successful attack. It is observed that after a certain amount of time(t = 9 here), the fraction of nodes approaches a fixed value. However, in the unsuccessful attack scenario i.e., for parameters such that R0 < 1, we see from Table 3 that the fraction of infected targeted nodes decreases as the time increases. The starting point has been taken same for both the Tables.

R.K. Upadhyay and P. Singh / Physica A 538 (2020) 122617

11

Table 2 Distriution of node population for targeted and attacking class in a successful attack case for R0 > 1. t

St

It

Ia

Ea

% of infected targeted nodes

1 5 7 9 10 20 40 60 80

0.5013 0.3156 0.319 0.3141 0.3151 0.3167 0.3166 0.3167 0.3167

0.2777 0.3163 0.3013 0.2948 0.2936 0.2929 0.2929 0.2929 0.2928

0.4965 0.4623 0.4625 0.4623 0.4627 0.4625 0.4625 0.4626 0.4628

2.9308 4.6299 4.6669 4.6728 4.6735 4.6739 4.6739 4.6739 4.6739

>20% >30% >30% >20% >20% >20% >20% >20% >20%

Table 3 Distribution of node population for targeted and attacking class in an unsuccessful attack case for R0 < 1. t

St

It

Ia

Ea

% of infected targeted nodes

5 20 30 50 70 90 150 250 300

0.2212 0.2692 0.3001 0.3627 0.4252 0.4866 0.6537 0.8450 0.9010

0.3605 0.3107 0.2974 0.2706 0.2438 0.2175 0.1464 0.0653 0.0417

0.7218 0.5708 0.4899 0.3684 0.2827 0.2202 0.1099 0.0378 0.0226

0.0135 0 0 0 0 0 0 0 0

>30% >30% >20% >20% >20% >20% <20% <10% <10%

5. Optimal control analysis Now we examine the effect of given control parameters in the dynamics of virus spread in networks [24]. We modify the transmission rate given in (3) and develop a control problem. The aim is to find an optimal solution and minimize the infection transmission. We also aim for minimizing the cost of enforcing the control measures. In system (3), we change the transmission rate by including the term of control variable u(t), where u(t) indicates the use of antivirus measures to be installed in the targeted system. Thereby, in the targeted population, the related infection is reduced by a factor of (1 − u(t)). The purpose is to reduce the number of infected nodes from targeted class, and also the cost of employing the control using the least possible control variable u [25]. With control measures, the dynamics of system (3) is given by: dSt dt dIt dt dIa dt dEa

= −β1 St Ia (1 − u(t)) + α (1 − St − It ), = β1 St Ia (1 − u(t)) − γ It , (19)

= β2 (1 − Ia − Ea )Ia + k2 Ea − µIa − ηIa ,

= ξ + η(1 − Ea ) − (µ + k1 + k2 )Ea , dt Now, the objective functional with respect to system (19) is defined as follows: T

∫

(

J(u) =

A1 It + 0

1 2

B1 u2

)

dt .

(20)

where A1 and B1 are weight constants of infected targeted nodes and antivirus software implementations respectively. The control set U consists of Lebesgue measurable control measures u(t) on [0, 1] such that 0 ≤ t ≤ T , 0 < u(t) ≤ 1. We aim to find the required control function u∗ such that J(u∗ ) = min {J(u)} u∈U

We show that control problem exists for the system (19) at t = 0. Then using the concept used in [25] and [26], Lagrangian for the system (19) is given by 1

B1 u2 . 2 Now we proceed to define the Hamiltonian H for this control problem as follows: L = A 1 It +

H = L(It , u) + λ1

dSt dt

+ λ2

dIt dt

+ λ3

dIa dt

+ λ4

dEa dt

.

(21)

(22)

12

R.K. Upadhyay and P. Singh / Physica A 538 (2020) 122617

Fig. 5. Figures showing the behavior of (a) susceptible and (b) infected targeted nodes, with and without control. Parameters are same as mentioned in (18).

We use the results from [25–27] to prove that an optimal control u∗ ∈ U exists such that J(u∗ ) = minu∈U J(u). Now, we obtain the optimal solution using theory from [24]. as follows: Let (x, u) be the optimal solution, and λ = (λ1 , λ2 , . . . , λi , i = 1, 2, 3, 4) be non-trivial vector function that satisfies:

λ′ = −

∂ H(t , x, u, λ) , ∂x

(23)

Now we use (23) to obtain adjoint variables λi , i = 1, 2, 3, 4 that satisfies:

λ′1 = λ1 α + β1 Ia (1 − u) − λ2 β1 Ia (1 − u), λ′2 = −A1 + αλ1 + γ λ2 , λ′3 = λ1 β1 St Ia (1 − u) − λ2 β1 St (1 − u) − λ3 β2 (1 − 2Ia − Ea ) + (µ + η)λ3 , λ′4 = −λ3 (k2 − β2 Ia ) + λ4 (µ + k1 + k2 + η),

(24)

and transversality conditions

λi (tb ) = 0, i = 1, 2, 3, 4. Then solving the equation

(25)

dH du

= 0, we have the optimal solution: } } (λ2 − λ1 )β1 St∗ Ia∗ ∗ ,1 ,0 . u = max min β1

{

{

(26)

2

d L ∗ ∗ Also, du 2 > 0 implies that the optimal control problem is minimum at u . Substituting u in the control system (19), we obtain:

( { { } }) (λ2 − λ1 )β1 St∗ Ia∗ = −β1 St Ia 1 − max min ,1 ,0 + α (1 − St − It ), dt β1 ( { { } }) dIt∗ (λ2 − λ1 )β1 St∗ Ia∗ = β1 St Ia 1 − max min ,1 ,0 − γ It , dt β1

dSt∗

dIa∗ dt dEa∗

= β2 (1 − Ia − Ea )Ia + k2 Ea − µIa − ηIa ,

= ξ + η(1 − Ea ) − (µ + k1 + k2 )Ea , dt with H ∗ at (t , St∗ , It∗ , Ia∗ , Ea∗ , u∗ , λ1 , λ2 , λ3 , λ4 ) so that ∗

∗

H = A1 It +

1 2

(

( B1

{

max min

{

(λ2 − λ1 )β1 St∗ Ia∗

β1

(27)

}

,1 ,0

}))2

.

(28)

Now, we numerically obtain the optimal control solution by solving the systems (27) and (28). We adopt a numerical method, described in [28] and [29] to solve the optimal system. Then, we display the graphs plotted using MATLAB with A1 = 100 and B1 = 2. Remaining parameters are same as given in (19). In Fig. 5, the dotted line displays the behavior of nodes with control and solid line depicts the system without control. We can see from Fig. 5(a) that the level of susceptible population with control is significantly greater than the one without control and

R.K. Upadhyay and P. Singh / Physica A 538 (2020) 122617

13

Fig. 6. Optimal solution curves for systems (27) and (28). Table 4 R The sensitivity indices, γpi 0 =

∂ R0 ∂ pi

×

pi R0

, of reproduction number R0 for parameters pi .

Parameter (pi )

Value

Sensitivity indices

β2 η µ β1

0.9 0.3 0.02 0.9 0.2 0.4 0.4 4 0.2

+1.00 −0.9375 −0.0625

k1 k2

γ ξ

m

0 0 0 0 0 0

from Fig. 5(b), it is observed that the population of infected nodes with control declines more rapidly than the nodes without control. Next is the Fig. 6 which presents the optimal curves for the systems (27) and (28). 6. Sensitivity analysis Now, to find out the relative importance of each parameter in virus dispersion prevalence, we conduct the sensitivity analysis of each parameter using the method adopted in [30,31]. The initial infection transmission is directly related to reproduction number (R0 ), and pervasiveness of infection is strongly linked with the endemic equilibrium. The proportion of infected targeted nodes is most important being the focus of our work and is directly related to the nodes getting crashed because of infection and hence needs specific attention. Thus, we conduct the sensitivity analysis of R0 , for the parametric values mentioned in Table 4. Sensitivity indices let us know regarding the significance of each parameter in infection dispersion and widespread prevalence. 6.1. Sensitivity indices of R0 R

We use the formula γp 0 =

∂ R0 ∂p

× Rp , to determine the sensitivity of R0 with respect to the parameters given in Table 4. 0

Parameters taken are also mentioned in the table. From Table 4 and Fig. 7, we see that R0 increases on increasing the transmission rate β2 and decreases with η and µ. R R As γβ20 = +1, increasing or decreasing β2 for say 50%, leads to same change in proportion of R0 . γη 0 = −0.9375 signifies that 10% increment in η will lead to 9.375% decrement in value of R0 and vice versa. 6.2. Sensitivity indices of endemic steady state E ∗ After R0 , we determine the sensitivity investigation of state variables at the endemic steady state. The obtained calculations have been shown in Table 5. Fig. 8 also shows the effect of β1 on infected targeted nodes and support our I∗

findings. γηt = −0.1832058 indicates that increasing (or decreasing) the value of η by 10%, decreases (or increases) It∗ by 1.832058%.

14

R.K. Upadhyay and P. Singh / Physica A 538 (2020) 122617

Fig. 7. Variation of R0 with parameters (a) β2 , (b) η and (c) µ. Parameters are same as mentioned in (18). Table 5

E∗

The sensitivity indices, γpii =

∂ Ei∗ ∂ pi

×

pi Ei∗

, of state variables for the parameters pi given in

Table 4. I∗

Parameter (pi )

γpit

β1 β2 η µ

0.3166077 −0.2349634 −0.1832058 −0.116727 0.00252831 0.29211198 0.3904434 −0.70705223 −0.01081884

k1 k2

α γ ξ

I∗

E∗

γpia

γpia

0

0 0

−0.834893 −0.057865 −0.0036864 0.432478 0.922736 0 0 −0.0341711

−0.768957 −0.0217391 −0.217391 −0.434783 0 0 0.930237

Now, we have in Fig. 9(a) the contour graph for change in R0 with infection transmission rate and natural crashing rate of all attacking system. We find that low contact rate with susceptible attacking nodes and high natural crashing rate of infected nodes will bring down R0 less than one i.e will hinder infection spread. Contour graph from Fig. 9(b) displays the change in reproduction number R0 with infection transmission rate and the rate at which nodes from Sa and Ia class gets disconnected from the internet to join external attacking nodes. It can be inferred that higher value of η and low transmission rate will keep R0 less than 1. 7. Conclusions and discussions In this work, a network model with distributed attack on targeted resources has been proposed. Cautious mathematical inspections have been done to justify the analytical findings. A certain threshold has been formulated on which the whole dynamics of system depends. Stabilities of equilibrium points and bifurcation have been discussed in terms of the threshold, called reproduction number. After that, effect of using firewall security is also shown graphically. A sensitivity analysis is also executed to find the relative impact of various parameters on virus propagation. Here, we conclusively present some of the crucial findings of the model:

R.K. Upadhyay and P. Singh / Physica A 538 (2020) 122617

15

Fig. 8. Variation in infected targeted nodes with infection transmission rate β1 .

Fig. 9. Contour graphs displaying the variation in (a) R0 , with β2 and µ, (b) R0 , with β2 and η. Other parameters are same as given in (18).

• The model system (3) has an infection-free equilibrium which is locally stable when R0 ≤ 1, and becomes unstable for R0 > 1. • As R0 crosses unity, system approaches another equilibrium which is the endemic equilibrium point. Using the results from central manifold theory, we have shown that endemic equilibrium is unstable for R0 > 1 and a forward transcritical bifurcation occurs at R0 = 1. • The coefficient of firewall helps in mitigating the virus propagation in a computer network to a certain extent. • It is observed that with proper control variables, the density of susceptible targeted nodes increases and that of infected targeted nodes goes down. • Sensitivity analysis has been carried out for R0 and endemic state variables, that tells the relative importance of each parameter in virus spread. • For further study of models of this kind, the infection transmission rate can be taken to be time-dependent to study the effect of time dependent infectiousness on the infection duration. For instance, transmission rate is negligible up to certain point of time and after that attains some value. The study made in this work may be applicable for controlling DDoS attacks. Many intrusion-detectors nowadays are equipped with technology to prevent DDoS attacks by using connection authentication methods and stopping specific requests from reaching operation servers. A backup internet connection shall be maintained with different IP addresses in case the primary circuit is flooded by malicious requests. DDoS prevention products can also be used that are specifically designed to detect and obstruct such attacks. References [1] N.T. Bailey, The Mathematical Theory of Infectious Diseases and its Applications, Charles Griffin & Company Ltd, 5a Crendon Street, High Wycombe, Bucks HP13 6LE, 1975.

16

R.K. Upadhyay and P. Singh / Physica A 538 (2020) 122617

[2] K. Haldar, B.K. Mishra, A mathematical model for a distributed attack on targeted resources in a computer network, Commun. Nonlinear Sci. Numer. Simul. 19 (9) (2014) 3149–3160. [3] J.R. Piqueira, A.A. De Vasconcelos, C.E. Gabriel, V.O. Araujo, Dynamic models for computer viruses, Comput. Secur. 27 (7–8) (2008) 355–359. [4] C.E.R. Team, C. E. R. Team, Results of the distributed-systems intruder tools workshop, http://www.cert.org/reports/dsitworkshop-final.html. [5] E. Alomari, S. Manickam, B. Gupta, S. Karuppayah, R. Alfaris, Botnet-based distributed denial of service (DDoS) attacks on web servers: classification and art, arXiv preprint arXiv:1208.0403. [6] R.K. Chang, Defending against flooding-based distributed denial-of-service attacks: a tutorial, IEEE Commun. Mag. 40 (10) (2002) 42–51. [7] F. Lau, S.H. Rubin, M.H. Smith, L. Trajkovic, Distributed denial of service attacks, in: Systems, Man, and Cybernetics, 2000 IEEE International Conference on, Vol. 3, IEEE, 2000, pp. 2275–2280. [8] S. Bellovin, Distributed denial of service attacks, Talk, NANOG, San Jose. [9] C. C. Center, CERTR advisory CA-2000-01 denial-of-service developments (2000). [10] S. Kumari, P. Singh, R.K. Upadhyay, Virus dynamics of a distributed attack on a targeted network: Effect of firewall and optimal control, Commun. Nonlinear Sci. Numer. Simul. (2019). [11] S.T. Zargar, J. Joshi, D. Tipper, A survey of defense mechanisms against distributed denial of service (DDoS) flooding attacks, IEEE Commun. Surv. Tutor. 15 (4) (2013) 2046–2069. [12] E. Gelenbe, M. Gellman, G. Loukas, Defending networks against denial-of-service attacks, in: Unmanned/Unattended Sensors and Sensor Networks, Vol. 5611, International Society for Optics and Photonics, 2004, pp. 233–244. [13] F. Xing, W. Wang, Understanding dynamic denial of service attacks in mobile ad hoc networks, in: Military Communications Conference, 2006. MILCOM 2006. IEEE, IEEE, 2006, pp. 1–7. [14] O. Diekmann, J.A.P. Heesterbeek, J.A. Metz, On the definition and the computation of the basic reproduction ratio R0 in models for infectious diseases in heterogeneous populations, J. Math. Biol. 28 (4) (1990) 365–382. [15] P. Van den Driessche, J. Watmough, Reproduction numbers and sub-threshold endemic equilibria for compartmental models of disease transmission, Math. Biosci. 180 (1–2) (2002) 29–48. [16] C. Castillo-Chavez, B. Song, Dynamical models of tuberculosis and their applications, Math. Biosci. Eng. 1 (2) (2004) 361–404. [17] M.Y. Li, J.S. Muldowney, A geometric approach to global-stability problems, SIAM J. Math. Anal. 27 (4) (1996) 1070–1083. [18] S. Sun, Global dynamics of a seir model with a varying total population size and vaccination, Int. J. Math. Anal. 6 (40) (2012) 1985–1995. [19] G. Butler, P. Waltman, Persistence in dynamical systems, J. Differential Equations 63 (2) (1986) 255–263. [20] K. Salah, K. Elbadawi, R. Boutaba, Performance modeling and analysis of network firewalls, IEEE Trans. Netw. Serv. Manage. 9 (1) (2012) 12–21. [21] G.P. Sahu, J. Dhar, Dynamics of an SEQIHRS epidemic model with media coverage, quarantine and isolation in a community with pre-existing immunity, J. Math. Anal. Appl. 421 (2) (2015) 1651–1672. [22] J. Cui, Y. Sun, H. Zhu, The impact of media on the control of infectious diseases, J. Dynam. Differential Equations 20 (1) (2008) 31–53. [23] A. Bhargava, D.K. Soni, P. Jain, J. Dhar, Dynamics of attack of malicious codes on the targeted network: Effect of firewall, in: Recent Trends in Information Technology (ICRTIT), 2016 International Conference on, IEEE, 2016, pp. 1–6. [24] L.S. Pontryagin, Mathematical Theory of Optimal Processes, Routledge, 2018. [25] A.A. Lashari, G. Zaman, Optimal control of a vector borne disease with horizontal transmission, Nonlinear Anal. RWA 13 (1) (2012) 203–212. [26] E. Bonyah, M. Khan, K. Okosun, J. Gómez-Aguilar, Modelling the effects of heavy alcohol consumption on the transmission dynamics of gonorrhea with optimal control, Math. Biosci. (2018). [27] D.L. Lukes, Differential Equations: Classical to Controlled, Elsevier, 1982. [28] S. Lenhart, J.T. Workman, Optimal Control Applied to Biological Models, Crc Press, 2007. [29] H.S. Rodrigues, M.T.T. Monteiro, D.F. Torres, Optimal control and numerical software: an overview, arXiv preprint arXiv:1401.7279. [30] R.K. Upadhyay, P. Roy, Deciphering dynamics of recent epidemic spread and outbreak in west africa: the case of ebola virus, Int. J. Bifurcation Chaos 26 (09) (2016) 1630024. [31] M. Samsuzzoha, M. Singh, D. Lucy, Uncertainty and sensitivity analysis of the basic reproduction number of a vaccinated epidemic model of influenza, Appl. Math. Model. 37 (3) (2013) 903–915.