Mathematical analysis of an influenza A epidemic model with discrete delay

Accepted Manuscript Mathematical analysis of an influenza A epidemic model with discrete delay P. Krishnapriya, M. Pitchaimani, Tarynn M. Witten PII: DOI: Reference:

S0377-0427(17)30201-7 http://dx.doi.org/10.1016/j.cam.2017.04.030 CAM 11110

To appear in:

Journal of Computational and Applied Mathematics

Received date: 12 January 2017 Please cite this article as: P. Krishnapriya, M. Pitchaimani, T.M. Witten, Mathematical analysis of an influenza A epidemic model with discrete delay, Journal of Computational and Applied Mathematics (2017), http://dx.doi.org/10.1016/j.cam.2017.04.030 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.

Manuscript Click here to view linked References

Mathematical analysis of an influenza A epidemic model with discrete delay P. Krishnapriya1 , M.Pitchaimani2 , and Tarynn M. Witten3 1,2

Ramanujan Institute for Advanced Study in Mathematics University of Madras, Chennai-600005, India, 3 Center for the Study of Biological Complexity Virginia Commonwealth University, Richmond, Virginia 23284-2030, USA [email protected]

Abstract Recently, a large number of mathematical models that are described by delay differential equations (DDEs) have appeared in the life sciences. In this paper, we present a delay differential model to describe influenza A (H1N1) dynamics. We begin by presenting the model with a brief discussion, followed by proving the positivity and boundedness of the model solution. We establish sufficient conditions for the global stability of the equilibria (the infection free equilibrium and the infected equilibrium), these conditions are obtained by means of the Lyapunov LaSalle invariance principle of the system. Also we have carried out bifurcation analysis along with an estimated length of delay to preserve the stability behaviour. In particular, we show the threshold dynamics in the sense that if (reproduction number) R0 < 1 the infectious population disappear so the disease dies out, while if R0 > 1 the infectious population persist. Sensitivity analysis of the influenza A (H1N1) model reveals which parameter values have a major impact on the model dynamics. Numerical simulations with application to H1N1 infection are given to verify the analytical results.

Keywords: H1N1 Model, Delays, Stability analysis, Sensitivity analysis. Mathematical Subject Classification (2010): 92B05, 93D05, 34D23, 90C31.

1

Introduction

Over the past one hundred years, mathematics has been used to understand and predict the spread of diseases, relating important public-health questions to basic transmission parameters. From prehistory to the present day, diseases have been a source of fear and superstition. A comprehensive picture of disease dynamics requires a variety of mathematical tools, from model creation to solving differential equations to statistical analysis. Although mathematics has been so far done quite well in dealing with epidemiology but there is no denying that there are certain factors which still lack proper mathematization.

1

Almost all mathematical models of diseases start from the same basic premise: that the population can be subdivided into a set of distinct classes, dependent upon their experience with respect to the disease [1, 2]. Most of them are described by nonlinear, differential equations (delay-difference, stochastic, etc.). The introduction of time delay is often used to model the latent period, that is, the time from the acquisition of infection to the time when the host becomes infectious [3]. Most of the authors assume that the latent period of diseases is negligible, i.e., once infected, each susceptible individual (S) instantaneously becomes infectious (I), and later recovers (R) with a permanent or temporary acquired immunity. These epidemic models are customarily called SIR (susceptible, infectious, recovered) models [4–7]. It is known that for some diseases, such as influenza and tuberculosis, on adequate contact with an infectious individual, a susceptible becomes exposed for a while; that is, infected but not yet infectious (incubating the disease). Thus it is realistic to introduce a latent compartment (also known as an incubation compartment). Time delays have also been incorporated into mathematical models to study H1N1 dynamics. Further within-host H1N1 models ODE and including time delays can be found in [8–11]. Human cases of swine influenza A (H1N1) virus infection have been recently identified in several countries. Cases of novel influenza (H1N1) were first identified in mid-April 2009 in California and soon thereafter in Texas and Mexico [12]. A simple ordinary differential equation model to study the epidemiological consequences of the drift mechanism for influenza A viruses was proposed in [13]. An extended SEIR model to describe the immigration of infected people and the stochastic variation of the infectious efficiency in [14]. The potential problems of the standard maximum likelihood estimate (MLE) approaches are demonstrated and possible remedies suggested on the epidemic data from the 2009 outbreak of H1N1 influenza on the campus of Washington State University [15]. Actually, Influenza is highly contagious and is easily transmitted through contact with droplets from the nose and throat of an infected person who is coughing and sneezing. The disease infects the nose, throat or lungs. Typically, an area can have epidemic conditions for a period of 46 weeks before it eases off. Onset of symptoms range from 18 to 72 h. The H1N1 influenza virus that caused the outbreak had never been observed before and was described as having components of human, avian and swine sources in its genome [16]. Swine serve as a mixing-vessel [17], since they are susceptible to infection with viruses from birds and other mammals, thereby providing an opportunity for genetic reassortment between influenza viruses during a mixed infection [18]. In general, delay-differential equations exhibit much more complicated dynamics than ordinary differential equations since a time delay could cause a stable equilibrium to become unstable and cause the populations to fluctuate. In 2012, Xueyong et.al. [19] developed a simple ode model for influenza A (H1N1) dynamics. They proposed on ODE model with vaccination. In their study shows that higher values of vaccination rate Φ significantly reduce the number of infected individuals, and lead to disease eradication. They defined the model as five compartments: but they were considering the four compartments: cells that are suspected, cells that are vaccination, cells that are exposed, and cells that are infectious. They described the dynamics of these populations by a system of four ordinary differential equations.

2

In this paper, we incorporate a discrete time delay to the model to describe the time necessary for exposed individual to become a infectious; the incubation period. The resulting model is a system of delay-differential equations. To understand the dynamics of the delay model, we study the transcendental characteristic equation of the linearized system at the positive infected steady state and obtain analytic conditions on the parameters under which the infected steady state is globally asymptotically stable. Numerical simulations are carried out to illustrate the obtained results. The paper is organized as follows. In Section 2 we first introduce the SEIR epidemic model with time delay system. Next, we analyze the model dynamics. Results on the existence of infection free and infection steady state and the threshold value are presented and local stability of infection free steady state and also global behavior of infection free and infected steady state by constructing suitable Lyapunov functional are analyzed in Section 3. In section 4 describe the estimation of length of delay preserve the stability. Section 5 addresses the sensitivity functions to evaluate the sensitivity of the parameters. In Section 6 we provide some numerical simulations of the system. Lastly, in Section 7, we close with discussion and conclusion.

2

Model Description

In this paper, we propose the following SEIR epidemic model for (H1N1) influenza with discrete delay. From [19], we have the following set of equation as our basic model dS dt dV dt dE dt dI dt dR dt

= A − β(I(t) + ηE(t))S(t) − (Φ + µ)S(t), = ΦS(t) − (1 − σ)β(I(t) + ηE(t))V (t) − µV (t), = β(I(t) + ηE(t))S(t) + (1 − σ)β(I(t) + ηE(t))V (t) − k1 E(t) − µE(t), = k1 E(t) − dI(t) − δI(t) − µI(t), = δI(t) − µR(t).

(1)

For the purposes of our discussion, we have modified the above model (1) and also introduce the delay into the above system as follows, dS dt dE dt dI dt dR dt

= A − βI(t)S(t) − µ1 S(t), = βI(t)S(t) − k1 E(t − τ ) − µ2 E(t), = k1 E(t − τ ) − dI(t) − δI(t) − µ3 I(t), = δI(t) − µ4 R(t).

(2)

where τ represents the time necessary for exposed to become a infectious. Here, ‘A’ represents the birth or immigration rate, and the positive constant µ1 , µ2 , µ3 and µ4 represent 3

the death rate of susceptible, exposed, infectious and recovered. β is the infection rate and k1 represents the transfer rate between the exposed and infectious also d represents the disease death rate, δ represents the rate of recovery from the disease. Now, we that we have simplified our model (2), we observe that it has two steady ¯ E, ¯ I, ¯ R) ¯ and the infected steady state states: The infection free steady state E0 = (S, ∗ ∗ ∗ ∗ E1 = (S , E , I , R ). The basic reproduction number: The basic reproduction number, denoted R0 , is ‘the expected number of secondary cases produced, in a completely susceptible population, by a typical infective individual’ [20]. If R0 < 1, then on average an infected individual produces less than one new infected individual over the course of its infectious period, and the infection cannot grow. Conversely, if R0 > 1, then each infected individual produces, on average, more than one new infection, and the disease can invade the population. For the case of single infected compartment, R0 is simply the product of the infection rate and the mean duration of the infection. Now, we will calculate the basic reproduction number of the system (2). Let X = (S, E, I, R)T , then the system (2) can be written as dX dt

= F (X) − V (X),

where   βI(t)S(t)   0 , F (x) =    0 0 

We can obtain,

 k1 E(t) + µ2 E(t) −k1 E(t) + dI(t) + δI(t) + µ3 I(t) . V (x) =    −A + βI(t)S(t) + µ1 S(t) −δI(t) + µ4 R(t)

F

=



V

−1

 0 β S¯ , 0 0

V =



 k1 + µ2 0 . −k1 d + δ + µ3

yielding =

"

1 k1 +µ2 −k1 (k1 +µ2 )(d+δ+µ3 )

0 1 d+δ+µ3

#

.

F V −1 is the next generation matrix for model (2). It then follows that the spectral βk1 A . According to Theorem radius of matrix F V −1 is, ρ(F V −1 ) = µ1 (k1 + µ2 )(d + δ + µ3 ) 2 in [21], the reproduction number of model (2) is

4

R0 =

βk1 A . µ1 (k1 + µ2 )(d + δ + µ3 )

(3)

¯ 0, 0, 0) (where It is easy to see that if R0 < 1, then the infection free steady state E0 (S, A ¯ S = µ1 ) is the unique steady state, corresponding to the extinction of the infection free state. Now we have the following result for concerning the existence of equilibria. Theorem 2.1. If R0 > 1, then the system (2) has an unique equilibrium E1 (S ∗ , E ∗ , I ∗ , R∗ ) (i.e., S ∗ > 0, E ∗ > 0, I ∗ > 0, R∗ > 0) where S ∗ , E ∗ , I ∗ and R∗ are given in the proof. Proof. If R0 > 1, then the system (2) becomes as follows, A − βI ∗ S ∗ − µ1 S ∗ = 0,

βI ∗ S ∗ − k1 E ∗ − µ2 E ∗ = 0,

k1 E ∗ − dI ∗ − δI ∗ − µ3 I ∗ = 0, δI ∗ − µ4 R∗ = 0.

(4)

From the above equation (4), we easily find that S∗ = E∗ = I∗ = R∗ =

(k1 + µ2 )(d + δ + µ3 ) , βk1 (d + δ + µ3 ) (R0 − 1), k1 β 1 (R0 − 1) β δ (R0 − 1). µ4 β

Hence, the system (2) has a unique equilibrium E1 (S ∗ , E ∗ , I ∗ , R∗ ) if R0 > 1.



Summarizing the previous analysis, we have the following result. Theorem 2.2. Consider the system (2) with R0 defined in (3). If R0 < 1, then there is unique equilibrium, which is the infection-free steady state E0 ; while if R0 > 1, then there is unique equilibrium, which is the infected steady state E1 . Table I

Parameters description Parameters A µ1 d µ2 k1 µ3 β µ4 δ

Description Birth of immigrant rate Natural death rate Disease-induced death rate Natural death rate Transfer rates between the exposed and the infectious Natural death rate infection rate Natural death rate Rate of recovery from the disease

5

Values variable 5.48 ×10−5 /day 0.001/day 5.48 ×10−5 /day 0.2/day 5.48 ×10−5 /day variable 5.48 ×10−5 /day 0.14/day

Reference Assumed [30] [31] [30] [32] [30] Assumed [30] [33]

3

Positivity and boundedness

It is important to show that the positivity and the boundedness for the system (2) as they represent cell populations. Positivity implies that the cell population survives and boundedness may be interpreted as a natural restriction to growth as a consequence of limited resources. In this Section, we present some basic results, such as the positivity and the boundedness of solution of model (2). Now, we find out the positive solution of the newly discussed model (2). We denote by C([−τ, 0], R4+ ) is the Banach space of continuous function mapping the interval [−τ, 0] into R4+ , equipped with the sup-norm. By the standard theory of functional differential equation [22–24], we know that for any φ ∈ C([−τ, 0], R4+ ), there exists a unique solution  C+ = φ ∈ C([−τ, 0], R4+ ) : S(ξ) = φ1 (ξ), E(ξ) = φ2 (ξ), I(ξ) = φ3 (ξ), R(ξ) = φ4 (ξ)} .

(5)

of the delayed system (2), where φ = (φ1 , φ2 , φ3 , φ4 ) ∈ R4+ with φi (ξ) ≥ 0 : (ξ ∈ [−τ, 0], i = 1, ..., 4), and φ1 (0), φ2 (0), φ3 (0), φ4 (0) > 0.

3.1

Positive Invariance

Let us put a system (2) in a vector form by setting X = (S, E, I, R)T ∈ R4 ,

(6)



   F1 (X) A − βI(t)S(t) − µ1 S(t), F2 (X)  βI(t)S(t) − k1 E(t − τ ) − µ2 E(t),     F (X) =  F3 (X) = k1 E(t − τ ) − dI(t) − δI(t) − µ3 I(t), F3 (X) δI(t) − µ4 R(t).

where F : C+ → R4 and F ∈ C ∞ (R4 ). Then equation (2) becomes X˙ = F (Xt )

(7)

where . ≡ d/dt and with Xt (θ) = X(t + θ), θ ∈ [−τ, 0] in [25]. It is easy to check in equation (7), that whenever we choose X(θ) ∈ C+ such that Xi = 0, then we obtain Fi (X)|Xi (t)=0 , Xt ∈ C+ ≥ 0, i = 1, 2, 3. Due to lemma [26], any solution of equation of (7) with Xt (θ) ∈ C+ , say, X(t) = X(t, X(0)), is such that X(t) ∈ R40+ for all t > 0. For boundedness of the solution, we define V (t) = S(t) + E(t) + I(t) + R(t) and p = min{µ1 , µ2 , d + µ3 , µ4 }. By non-negativity of the solution, it follows that dV (t) = A − p {S(t) + E(t) + I(t) + R(t)} dt ≤ A − pV (t). This implies that V (t) is bounded, and so are S(t), E(t), I(t) and R(t). This completes the proof. 6

The system (2) will be analyzed in biologically feasible region as follows. Consider the feasible region   A , E, I, R ≥ 0 Ω = (S, E, I, R) ∈ R4+ : S ≤ µ1 From the above, we can establish the following theorem. Theorem 3.1. The region Ω is positively invariant for the system (2) with initial conditions in R4+ . Thus, the epidemic model of the H1N1 will be considered in C([−τ, 0], R4+ ). Further we will find the sufficient conditions on the parameters for the stability of the infection free and the infected steady states.

3.2

Local stability analysis

In this section we study the local stability analysis of model (2). Theorem 3.2. The infection free steady state of model (2) is locally asymptotically stable at E0 when R0 < 1 and unstable when R0 > 1 in the case of τ ≥ 0. When R0 > 1, the system (2) has a infected steady state E1 = (S ∗ , E ∗ , I ∗ , R∗ ). Then the linearized system (2) at E1 yields −(µ1 + βI ∗ ) − λ 0 −βS ∗ 0 ∗ −λτ βI −(ke + µ2 ) − λ βS ∗ 0 −λτ 0 ke −(d + δ + µ3 ) − λ 0 0 0 δ −µ4 − λ

= 0,

Thus the characteristic polynomial of the above determinant as follows, λ3 + A1 λ2 + A2 λ + A3 + e−λτ (B1 λ2 + B2 λ + B3 ) = 0.

(8)

Theorem 3.3. The infected steady state of model (2) is locally asymptotically stable when R0 > 1 in the case of τ > 0. Proof. In the case of τ > 0, the above characteristic equation (8) an be rewritten as H(λ, τ ) = P (λ) + Q(λ)e−λτ = 0, where P (λ) = λ3 + A1 λ2 + A2 λ + A3 and Q(λ) = B1 λ2 + B2 λ + B3 . Also, A1 = µ2 + d + δ + µ3 + µ1 + R0 − 1,

A2 = µ2 (d + δ + µ3 + µ1 + R0 − 1) + (d + δ + µ3 )(µ1 + R0 − 1), A3 = µ2 (d + δ + µ3 )(µ1 + R0 − 1), B1 = k1 ,

B2 = k1 (d + δ + µ3 + µ1 + R0 − 1),

B3 = k1 (d + δ + µ3 )(µ1 + R0 − 1) + βA/µ1 . 7

If the equation (8) has a purely imaginary root λ = iω, with ω > 0, separating the real and imaginary parts, we obtain the system of transcendental equations as A1 ω 2 − A3 = (B3 − B1 ω 2 ) cos(ωτ ) + B2 ω sin(ωτ ),

(9)

ω 3 − A2 ω = B2 ω cos(ωτ ) − (B3 − B1 ω 2 ) sin(ωτ ).

(10)

Squaring and adding (9) and (10), we obtain ω 6 + C1 ω 4 + C2 ω 2 + C3 = 0.

(11)

where C1 = A21 − 2A2 − B12 ,

C2 = A22 − B22 − 2A1 A3 + 2B1 B3 , C3 = A23 − B32 .

Put ω 2 = ρ in (11), we get h(ρ) = ρ3 + C1 ρ2 + C2 ρ + C3 = 0.

(12)

From the hypothesis of R0 > 1, we deduce that C3 = A23 − B32 > 0, since for R0 > 1, C1 and C2 are also positive. Thus the equation (12) has no positive solution for R0 > 1. This completes the proof.  Now, we assume that C3 < 0, by using Descartes rule of signs, equation (12) has positive root ρ and thus equation (11) has a pair of purely imaginary roots iω ∗ . From equation (9) and (10), we obtain   (A1 ω ∗2 − A3 )(B3 − B1 ω ∗2 ) + B2 ω ∗2 (A2 − ω ∗2 ) 1 arccos τ∗ = ω∗ (B3 − B1 ω ∗2 )2 + B22 ω ∗2 2jπ where j = 0, 1, 2, ..., + ∗, ω   dRe(λ) Now we determine sign ∗ , where sign is the signum function and Re(λ) is dτ τ =τ the real part of λ. By using the mathematical calculation we can say that the infected steady state of model (2) remains stable for τ < τ ∗ and Hopf bifurcation occurs when τ = τ ∗. We first look for purely imaginary roots of λ = iω ∗ of equation (8) implies |P (iω ∗ )| = |Q(iω ∗ )| Differentiating (8) with respect to τ , we obtain n o (3λ2 + 2A1 λ + A2 ) + e−λτ (2B1 λ + B2 ) − τ e−λτ (B1 λ2 + B2 λ + B3 )

dλ = λe−λτ (B1 λ2 + B2 λ + B3 ) dτ

8

which implies,  −1 dλ = dτ = = Therefore,

3λ2 + 2A1 λ + A2 2B1 λ + B2 τ + − , −λτ 2 2 λe (B1 λ + B2 λ + B3 ) λ(B1 λ + B2 λ + B3 ) λ 2B1 λ + B2 τ 3λ2 + 2A1 λ + A2 + − , 3 2 2 −λ(λ + A1 λ + A2 λ + A3 ) λ(B1 λ + B2 λ + B3 ) λ 2λ3 + A1 λ2 − A3 B 1 λ2 − B 3 τ + − −λ2 (λ3 + A1 λ2 + A2 λ + A3 ) λ2 (B1 λ2 + B2 λ + B3 ) λ

   2λ3 + A1 λ2 − A3 B 1 λ2 − B 3 τ Ξ = sign Re + − −λ2 (λ3 + A1 λ2 + A2 λ + A3 ) λ2 (B1 λ2 + B2 λ + B3 ) λ λ=iω ∗    (A3 + A1 ω ∗2 ) + i2ω ∗3 1 B1 ω ∗2 + B3 = + sign Re ω ∗2 (A1 ω ∗2 − A3 ) + i(ω ∗3 − A2 ω ∗ ) (B3 − B1 ω ∗2 ) + i(B2 ω ∗ )   (A3 + A1 ω ∗2 )(A1 ω ∗2 − A3 ) + 2ω ∗3 (ω ∗3 − A2 ω ∗ ) (B1 ω ∗2 + B3 )(B3 − B1 ω ∗2 ) 1 sign + = ω ∗2 (A1 ω ∗2 − A3 )2 + (ω ∗3 − A2 ω ∗ )2 (B3 − B1 ω ∗2 )2 + (B2 ω ∗ )2   1 (A3 + A1 ω ∗2 )(A1 ω ∗2 − A3 ) + 2ω ∗3 (ω ∗3 − A2 ω ∗ ) + (B1 ω ∗2 + B3 )(B3 − B1 ω ∗2 ) = sign ω ∗2 (B3 − B1 ω ∗2 )2 + (B2 ω ∗ )2  ∗6  1 2ω + (A21 − 2A2 − B12 )ω ∗4 + (A23 − B32 ) = sign ω ∗2 (B3 − B1 ω ∗2 )2 + (B2 ω ∗ )2 and this determines a set of possible eigenvalues of ω ∗ . Our aim is to determine the direction of motion of λ as τ is varied. i.e., we determine (   )   dλ −1 d(Re(λ)) = sign Re Ξ = sign dτ dτ λ=iω ∗ ∗ λ=iω

A21 − 2A2

− and − are both positive by virtue of equation (11), we have  As dRe(λ) > 0. Thus, the solution curve of the characteristic equation (11) dτ ω=ω ∗ ,τ =τ ∗ crosses the imaginary axis. This shows that a Hopf bifurcation occurs at 0 < τ = τ ∗ . By continuity, the infected steady state is locally asymptotically stable when τ < τ ∗ . By summarizing the above analysis, we arrive at the following theorem. B12

A23

B32

Theorem 3.4. Suppose R0 > 1, the following result can be obtained. 1. The infected equilibrium E1 is stable when τ ∈ [0, τ ∗ ) and unstable when τ > τ ∗ . τ is the Hopf bifurcation value, which means that periodic solutions will bifurcate from this infected equilibrium as τ passes through the critical value τ ∗ .

3.3

Global stability analysis

We shall consider the global stability of the infection free and the infected steady state of model (2) by the Lyapunov direct method. We define a function g : R>0 → R≥0 as g(u) = u − 1 − ln u. We note that g(u) ≥ 0 for any u > 0 and has the global minimum 0 at u = 1. 9

Theorem 3.5. The infected free steady state of model (2) is globally asymptotically stable when R0 < 1. Proof. We consider a Lyapunov functional U1 (t) as follows:   S ¯ U1 = SF + E + I. S¯

(13)

Here U1 is well-defined, continuous and positive definite for all (S, E, I, R) > 0 and θ ∈ [0, τ ]. Also, the global minimum U1 (= 0) occurs at the infection free steady state E0 . Thus every solution tends to the infection free steady state E0 . Further functions U1 along the trajectories of the system (2) satisfies   dU1 S¯ dS dE dI = 1− + + dt S dt dt dt   ¯ S (A − βIS − µ1 S) + βIS − k1 E(t − τ ) − µ2 E = 1− S +k1 E(t − τ ) − (d + δ + µ3 )I (14) Using the infection free steady state condition of model (2) A = µS¯ in (14), then the equation (14) becomes dU1 dt

µ1 ¯ 2 + βI S¯ − µ2 E − (d + δ + µ3 )I, (S − S) S
µ1 ¯ 2 + β S¯ − (d + δ + µ3 ) I − µ2 E, = − (S − S) S   k1 µ1 2 ¯ R0 − 1 (d + δ + µ3 )I − µ2 E, = − (S − S) + S k1 + µ2 ≤ 0.

= −

(15)

dU1 dU1 ¯ E(t) = is negative if R0 < 1. Also note that = 0 if and only if S(t) = S, dt dt I(t) = R(t) = 0. Hence by LaSalle invariance principle [25], the infection free equilibrium point E0 is globally asymptotically stable on Ω.  Thus

Theorem 3.6. If infected steady state E1 exists, then it is globally asymptotically stable for any τ < τ ∗ , when R0 > 1. Proof. Define a Lyapunov functional as follows:    Z t  E(θ) I ∗ ∗ + k1 E dθ U2 = S(t) + E(t) + I g I∗ E∗ t−τ

(16)

It is easy to see that U2 > 0, moreover U2 = 0 if and only if (S, E, I, R) takes the steady state value (S ∗ , E ∗ , I ∗ , R∗ ) and S(t − θ) = S ∗ , E(t − θ) = E ∗ , I(t − θ) = I ∗ and R(t − θ) = R∗ for all θ ∈ [0, τ ].     dU2 dS dE I ∗ dI E(t − τ ) ∗ E(t) − E(t − τ ) = + + 1− + k1 E + ln dt dt dt I dt E∗ E

10



 I∗ = A − βIS − µ1 S + βIS − k1 E(t − τ ) − µ2 E + 1 − (k1 E(t − τ ) I   E(t − τ ) ∗ −(d + δ + µ3 ))I + k1 E − k1 E(t − τ ) + k1 E ln , E

(17)

Using the infected steady state condition A = βI ∗ S ∗ + µS ∗ and d + δ + µ3 = k1 E ∗ in (17), then equation (17) becomes    E(t − τ ) I I∗ ∗ ∗ ∗ ∗ ∗ − ∗ = −µ1 (S − S ) + BI S − k1 E(t − τ ) − µ2 E + k1 E I 1 − I E∗I ∗ I   E(t − τ ) +k1 E − k1 E(t − τ ) + k1 E ∗ ln , E After some little algebra we obtain   E(t − τ ) E∗ E(t − τ )I ∗ ∗ ∗ E(t − τ ) , − 1 − ln + ln = −µ1 (S − S ) − k1 E − k 1 E∗ E∗ E I ≤ 0. (18)    ∗ E(t − τ ) E(t − τ ) E since ln = ln + ln . It follows that U2 (t) is bounded and ∗ E E E dU2 (t) = 0 if and only if E = non-increasing and thus limt→∞ U2 (t) exists. Note that dt dU2 (t) E ∗ , E(t − τ ) = E ∗ . Therefore, the maximal compact invariant set in = 0 is the dt singleton E1 . By the LaSalle invariance principle for delay systems in [25], the infected steady state of model (2) is globally attracting. Hence the infected steady state of model (2) is globally asymptotically stable when R0 > 1, in the case of τ < τ ∗ . 

4

Estimation of the Length of Delay to Preserve Stability

Following lines of Erbe et al. [27] and using the Nyquist criterion [28], it can be shown that the sufficient conditions for the local asymptotic stability of E1 (S ∗ , E ∗ , I ∗ , R∗ ) are given by, ImH(iω0 ) > 0,

(19)

ReH(iω0 ) = 0,

(20)

where H(ν) = ν 3 + A1 ν 2 + A2 ν + A3 + e−ντ (B1 ν 2 + B2 ν + B3 ) and ω0 is the smallest positive root of (20). Inequality (19) and (20) can alternatively be written as A2 ω0 − ω03 > −B2 ω0 cos(ω0 τ ) + B3 sin(ω0 τ ) − B1 ω02 sin(ω0 τ ),

(21)

A3 − A1 ω02 = B1 ω02 cos(ω0 τ ) − B3 cos(ω0 τ ) − B2 ω0 sin(ω0 τ ).

(22)

Now if (21) and (22) are satisfied simultaneously, they are sufficient conditions to guarantee stability. These are now used to get an estimate to the length of the time delay. The aim is 11

to find an upper bound ω+ to ω0 independent of τ from (22) and then to estimate τ so that (21) holds true for all values of ω such that 0 ≤ ω ≤ ω+ , and hence, in particular at ω = ω0 . Equation (22) can be rewritten as A1 ω02 = A3 − B1 ω02 cos(ω0 τ ) + B3 cos(ω0 τ ) + B2 ω0 sin(ω0 τ ).

(23)

Maximizing the right hand side of (23) subject to, | sin(ω0 τ )| ≤ 1,

| cos(ω0 τ )| ≤ 1,

(24)

we obtain |A1 |ω02 ≤ |A3 | + |B3 | + |B1 |ω02 + |B2 |ω0 .

(25)

Hence if, ω+ =

1 2(|A1 | − |B1 |)

  q 2 |B2 | + B2 + 4(|A1 | − |B1 |)(|A3 | + |B3 |) ,

(26)

then clearly from (25) we have ω0 ≤ ω+ . From (21), we obtain

ω02 < A2 + B2 cos(ω0 τ ) + B1 ω0 cos(ω0 τ ) −

B3 sin(ω0 τ ) . ω0

(27)

Since E1 (S ∗ , E ∗ , I ∗ , R∗ ) is locally asymptotically stable for τ = 0, the inequality (27) will continue to hold for sufficiently small τ > 0. Using (23) and (27) can be rearranged as  
A1 B3 2 sin(ω0 τ ) B3 − B1 ω0 − A1 B2 (cos(ω0 τ ) − 1) + (B2 − A1 B1 )ω0 + ω0 < A1 A2 − A3 − B3 + B1 ω02 + A1 B2 .

(28)

Using the bound ω τ 
0 B3 − B1 ω02 − A1 B2 (cos(ω0 τ ) − 1) = −(B3 − B1 ω02 − A1 B2 )2sin2 2 1 2 2 ≤ |(−B3 + B1 ω02 + A1 B2 )|ω+ τ , 2  
A1 B3 2 (B2 − A1 B1 )ω0 + sin(ω0 τ ) ≤ |(B2 − A1 B1 )|ω+ + |A1 ||B3 | τ, ω0

(29)

we obtain from (27)

L1 τ 2 + L2 τ < L3 ,

(30)

where, 1 2 |(−B3 + B1 ω02 + A1 B2 )|ω+ , 2 2 = |(B2 − A1 B1 )|ω+ + |A1 ||B3 |,

L1 = L2

2 + A1 B2 . L3 = A1 A2 − A3 − B3 + B1 ω+

12

(31)

Hence if, τ+ =

1 2L1



L2 +

q

L22



+ 4L1 L3 ,

(32)

then for 0 ≤ τ ≤ τ+ , the Nyquist criterion holds true and τ+ estimates the maximum length of the delay preserving the stability.

5

Parameter sensitivity analysis

The sensitivity and identifiability analysis were done by estimating at each data point the derivative of the clinical score Q = [S, E, I, R]T with respect to the vector parameter q = [µ1 , µ2 , µ3 , µ4 , k1 , d, β, δ, A]T . The sensitivities are dimensionless, as they are scaled with respect to the variables and parameter values. In the following the sensitivity functions with respect to an arbitrary parameter p, for the system (2) are denoted by, Qi,q =

∂Qi (t) , ∂q

i = 1, ...4

(33)

These sensitivities allow to determine which parameter is least important to the model output and these least sensitive parameters can be fixed and not used for calibration. There are different approaches to find the sensitivity functions of DDEs [29]. However, for simplicity we will use the so called direct approach to find sensitivity functions of system (2). The corresponding sensitivity system (2), with respect to the parameter ‘A’ is as follows,   dS = −βI ∗ SA (t, A) − βIA (t, A)S ∗ − µ1 SA (t, A), dt t,A   dE = βI ∗ SA (t, A) + βIA (t, A)S ∗ − k1 EA (t − τ, A) − µ2 EA (t, A), dt t,A   dI = k1 EA (t − τ, A) − dIA (t, A) − δIA (t, A) − µ3 IA (t, A), dt t,A   dR = δIA (t, A) − µ4 RA (t, A). (34) dt t,A The corresponding sensitivity system (2), with respect to the parameter ‘δ’ is as follows,   dS = −βI ∗ Sδ (t, δ) − βIδ (t, δ)S ∗ − µ1 Sδ (t, δ), dt t,δ   dE = βI ∗ Sδ (t, δ) + βIδ (t, δ)S ∗ − k1 Eδ (t − τ, δ) − µ2 Eδ (t, δ), dt t,δ   dI = k1 Eδ (t − τ, δ) − dIδ (t, δ) − I(t) − µ3 Iδ (t, δ), dt t,δ   dR = I(t) − µ4 Rδ (t, δ). (35) dt t,δ

13

The corresponding sensitivity system (2), with respect to the parameter ‘k1 ’ is as follows,   dS = −βI ∗ Sk1 (t, k1 ) − βIk1 (t, k1 )S ∗ − µ1 Sk1 (t, k1 ), dt t,k1   dE = βI ∗ Sk1 (t, k1 ) + βIk1 (t, k1 )S ∗ − E(t − τ ) − µ2 Ek1 (t, k1 ), dt t,k1   dI = E(t − τ ) − dIk1 (t, k1 ) − δIk1 (t, k1 ) − µ3 Ik1 (t, k1 ), dt t,k1   dR = δIk1 (t, k1 ) − µ4 Rk1 (t, k1 ). (36) dt t,k1 The corresponding sensitivity system (2), with respect to the parameter ‘µ1 ’ is as follows,   dS = −βI ∗ Sµ1 (t, µ1 ) − βIµ1 (t, µ1 )S ∗ − S(t), dt t,µ1   dE = βI ∗ Sµ1 (t, µ1 ) + βIµ1 (t, µ1 )S ∗ − k1 Eµ1 (t − τ, µ1 ) − µ2 Eµ1 (t, µ1 ), dt t,µ1   dI = k1 Eµ1 (t − τ, µ1 ) − dIµ1 (t, µ1 ) − δIµ1 (t, µ1 ) − µ3 Iµ1 (t, µ1 ), dt t,µ1   dR = δIµ1 (t, µ1 ) − µ4 Rµ1 (t, µ1 ). (37) dt t,µ1 The corresponding sensitivity of system (2), with respect to the parameter ‘µ2 ’ is as follows,   dS = −βI ∗ Sµ2 (t, µ2 ) − βIµ2 (t, µ2 )S ∗ − µ1 Sµ2 (t, µ2 ), dt t,µ2   dE = βI ∗ Sµ2 (t, µ2 ) + βIµ2 (t, µ2 )S ∗ − k1 Eµ2 (t − τ, µ2 ) − E(t), dt t,µ2   dI = k1 Eµ2 (t − τ, µ2 ) − dIµ2 (t, µ2 ) − δIµ2 (t, µ2 ) − µ3 Iµ2 (t, µ2 ), dt t,µ2   dR = δIµ2 (t, µ2 ) − µ4 Rµ2 (t, µ2 ). (38) dt t,µ2 The corresponding sensitivity of system (2), with respect to the parameter ‘µ3 ’ is as follows,   dS = −βI ∗ Sµ3 (t, µ3 ) − βIµ3 (t, µ3 )S ∗ − µ1 Sµ3 (t, µ3 ), dt t,µ3   dE = βI ∗ Sµ3 (t, µ3 ) + βIµ3 (t, µ3 )S ∗ − k1 Eµ3 (t − τ, µ3 ) − µ2 Eµ3 (t, µ3 ), dt t,µ3   dI = k1 Eµ3 (t − τ, µ3 ) − dIµ3 (t, µ3 ) − δIµ3 (t, µ3 ) − I(t), dt t,µ3   dR = δIµ3 (t, µ3 ) − µ4 Rµ3 (t, µ3 ). (39) dt t,µ3

14

The corresponding sensitivity of system (2), with respect to the parameter ‘µ4 ’ is as follows,   dS = −βI ∗ Sµ4 (t, µ4 ) − βIµ4 (t, µ4 )S ∗ − µ1 Sµ4 (t, µ4 ), dt t,µ4   dE = βI ∗ Sµ4 (t, µ4 ) + βIµ4 (t, µ4 )S ∗ − k1 Eµ4 (t − τ, µ4 ) − µ2 Eµ4 (t, µ4 ), dt t,µ4   dI = k1 Eµ4 (t − τ, µ4 ) − dIµ4 (t, µ4 ) − δIµ4 (t, µ4 ) − µ3 Iµ4 (t, µ4 ), dt t,µ4   dR = δIµ4 (t, µ4 ) − R(t). (40) dt t,µ4 The corresponding sensitivity system (2), with respect to the parameter ‘d’ is as follows,   dS = −βI ∗ Sd (t, d) − βId (t, d)S ∗ − µ1 Sd (t, d), dt t,d   dE = βI ∗ Sd (t, d) + βId (t, d)S ∗ − k1 Ed (t − τ, d) − µ2 Ed (t, d), dt t,d   dI = k1 Ed (t − τ, d) − I(t) − δId (t, d) − µ3 Id (t, d), dt t,d   dR = δId (t, d) − µ4 Rd (t, d). (41) dt t,d The corresponding sensitivity system (2), with respect to the parameter ‘β’ is as follows,   dS = −I ∗ S(t) − I(t)S ∗ − µ1 Sβ (t, β), dt t,β   dE = I ∗ S(t) + I(t)S ∗ − k1 Eβ (t − τ, β) − µ2 Eβ (t, β), dt t,β   dI = k1 Eβ (t − τ, β) − dIβ (t, β) − δIβ (t, β) − µ3 Iβ (t, β), dt t,β   dR = δIβ (t, β) − µ4 Rβ (t, β). (42) dt t,β The semi-relative sensitivity solutions (depicted in Figures.1-9) are calculated by simply multiplying the unmodified sensitivity solutions by a chosen parameter which provides information concerning the amount the state will change when that parameter is doubled. It is best to calculate this type of sensitivity solution to obtain a more thorough understanding of the dynamics. Using the parameter values (from Table I) and sensitivity results are obtained in this section. These shown in the following Figures. 1-9.

15

Fig. 1-2 show the semi-relative sensitivity analysis for the system (34) and (35) for Table I values respectively. From Fig.1, the parameter ‘A’ is very sensitive in the susceptible and infectious system (34). From Fig.2, the parameter ‘δ’ is very sensitive in the susceptible and infectious system (35). But, those parameters are negatively proportion with increasing the initial function and it is very sensitive in the early time intervals and the sensitivity decreases by time to be insensitive in the infected steady state in exposed and recovered levels.

Fig. 3-4 show the semi-relative sensitivity analysis for the system (36) and (37) for Table I values respectively. From Fig.3, the parameter ‘k1 ’ is very sensitive in the susceptible, exposed and recovered system (36). But, it is negatively proportion with increasing the initial function and it is very sensitive in the early time intervals and the sensitivity decreases by time to be insensitive in the infected steady state at infectious levels. From Fig.4, the parameter

16

‘µ1 ’ is very sensitive in the early time intervals at susceptible, exposed, infectious and recovered system (37).

Fig. 5-6 show the semi-relative sensitivity analysis for the system (38) and (39) for Table I values respectively. From Fig.5, the parameter ‘µ2 ’ is very sensitive in the exposed and recovered system (38). But, it is negatively proportion with increasing the initial function and it is very sensitive in the early time intervals and the sensitivity decreases by time to be insensitive in the infected steady state at infectious levels. From Fig.6, the parameter ‘µ3 ’ is very sensitive in the susceptible and infectious system (39). But, it is negatively proportion with increasing the initial function and it is very sensitive in the early time intervals and the sensitivity decreases by time to be insensitive in the infected steady state at exposed and recovered levels.

17

Fig. 7-8 show the semi-relative sensitivity analysis for the system (40) and (41) for Table I values respectively. From Fig.7, the parameter ‘µ4 ’ is very sensitive susceptible and infectious system (40). But, it is negatively proportion with increasing the initial function and it is very sensitive in the early time intervals and the sensitivity decreases by time to be insensitive in the infected steady state at exposed and recovered levels. From Fig.8, the parameter ‘d’ is very sensitive to small perturbation level in infected steady state at all level of susceptible, exposed, infectious and recovered system (41).

Fig. 9 show the semi-relative sensitivity analysis for the system (42) for Table I values respectively. From Fig.9, the parameter ‘β’ is sensitive only in recovered level. Since ‘β’ oscillate with increasing the initial function and it is very sensitive in the early time intervals and the sensitivity decreases by time to be insensitive in susceptible, exposed and infectious level. Also, from the above plots, we observe that a small change in any parameter can produce significant changes in the levels of system (34-42). Here the parameter ‘β’ plays a vital role in the model dynamics. Since it is very insensitive except in the recovered level. Therefore the small perturbation in the level of infection rate ‘β’, there is no change in that parameter. The size of the perturbation should be selected with care when determining the effect of the parameter on the outputs.

6

Application with Numerical simulation

In this section, we carry out some numerical simulations to display the qualitative behaviour of H1N1 model (2)(see Figs. 11 (a) - 12 (d)) fit the model to given data. The numerical simulations confirmed the theoretical results obtained in Section 3. We have seen in previous sections that the basic reproduction number R0 plays a decisive rule in determining the H1N1 dynamics. It is known that infected free and infected steady state is locally asymptotically stable if all the roots of the corresponding characteristic equa18

tions (8) has negative real parts (roots are located on the open left side of the complex plane i.e., Reλi < 0, i = 1, 2, ...). The precise values of the time delays intervals for some parameter based on information, we assumed that τ is ≤ 6 days, and estimated both of these by given data

6.1

Application

We linearized the system (2) at infected steady state E1 (7838.524944, 6426.052489, 9111.426889, 6377.998822), when R0 = 1.164005684 > 1 and using A = 0.5 β = 0.000018, τ = 6, we obtain, dS dt dE dt dI dt dR dt

= −0.000018I(t)S ∗ − 0.000018S(t)I ∗ − 0.0000548S(t), = 0.000018I(t)S ∗ + 0.000018S(t)I ∗ − 0.2E(t − τ ) − 0.0000548E(t), = 0.2E(t − τ ) − 0.001I(t) − 0.14I(t) − 0.0000548I(t), = 0.14I(t) − 0.0000548R(t).

(43)

Then the characteristic equation of (43), L(λ, τ ) = λ3 + 0.505170084λ2 + 0.0231505093λ + 0.000001268155228 +(0.2λ2 + 0.0610230568λ + 0.1360151651)e−6λ = 0.

(44)

The eigenvalues corresponding to variational matrix of infected equilibrium are −0.2657299313−0.6317097714I, −0.2657299313+0.6317097714I, −0.3885483113+1.318614981I, −0.4330897050 − 2.345248641I, −0.3885483113 − 1.318614981I

Finally, we can see that the Theorem 3.3 is satisfied, when R0 > 1. Hence the infected steady state of model (43) is locally asymptotically stable. This, together with Theorems 3.6 implies that the global stability of equilibria for system (43) is completely determined by the reproduction number R0 . From Table I, now we conclude that the model (43), takes the given data , all the characteristic roots of the characteristic equation for the model (44) has negative real parts, when R0 > 1 . Therefore our model (2) is always stable. The plots of the characteristic equations as shown in Fig. 10 (a)- 10 (b).

19

Fig. 10 (a) shows that the plot of the characteristic equation of (43) with given data for delay case (R0 > 1). Fig. 10 (b) shows that the plot of characteristic equation of (43) with given data for non-delay case (R0 > 1). To see the behaviour of system (2) and to verify stability results obtained for the given data. The solution trajectories are plotted in Fig. 11 (a) - 12(d).

20

21

Fig. 11 (a) - 11 (d): Solution trajectories of the system (43) with respect to the value of β = 0.000018 and using the initial conditions (S(0) = 1200, E(0) = 70, I(0) = 110, R(0) = 40) for the delay case (R0 > 1). Fig. 12 (a) - 12 (d): Solution trajectories of the system (43) with respect to the value of β = 0.000018 and using the initial conditions (S(0) = 1200, E(0) = 70, I(0) = 110, R(0) = 40) for the non-delay case (R0 > 1). Now, we compared our model for both delay and non-delay case at the same infection rate, which can be shown in Fig. 11 (a) - 12 (d). The delay model seems to fit the data, it could be better than using without delay model. Thus, our numerical results are reasonable representatives of our model. Moreover, small perturbations of parameters will give small perturbation of the matrix entries used in finding eigenvalues for determining the stability of two steady state points. The stability results in our numerical calculation are biologically reasonable and represent a qualitative outcome that is possible for the parameter values.

7

Discussion and Conclusion

In this paper, we have incorporated time delay into our H1N1 models. We showed that the threshold value R0 plays an important role in determining the stability of the steady states of the delay model of H1N1 model dynamics. The global stability of the infected free and the infected steady states have been established by using suitable Lyapunov functionals and LaSalle invariant principle. More specifically, we have proven that, if the threshold value R0 is greater than unity the infected steady state is locally asymptotically stable and globally attracting. For using sensitivity analysis, we can easily see that the infection rate parameter β plays a vital role in the model dynamics. The sensitivity functions are useful to evaluate which parameters of a significant uncertainty effect in the model. The allowable time delay for activation of the immune system and the estimation of the length of delay to preserve stability may be the two important parameters β and A that may help decide the mode of action for controlling the disease.

References [1] Junyuan Yanga, Maia Martcheva, Lin Wanga, Global threshold dynamics of an SIVS model with waning vaccine-induced immunity and nonlinear incidence, Math. Biosci., 268 : 18, (2015). [2] P. Krishnapriya, M. Pitchaimani, Optimal control of mixed immunotherapy and chemotherapy of tumours with discrete delay , Int. J. Dynam. Control, DOI 10.1007/s40435-015-0221-y, (2015). [3] B.D. Hassard, N.D. Kazarinoff, Y.H. Wan, Theory and Applications of Hopf Bifurcation, Cambridge University, Cambridge, (1981). [4] M.E. Alexander, C. Bowman, S.M. Moghadas, R. Summors, A.B. Gumel, B.M. Sahai, A vaccination model for transmission dynamics of influenza, SIAM. Appl.Dyn. Syst. 3 : 503-524, (2004).

22

[5] A. Korobeinikov, G.C. Wake, Lyapunov functions and global stability for SIR,SIRS & SIS epidemiological models, Appl. Math. Lett. 15 : 955-960, (2002). [6] Z. Ma & J. Li, Dynamical modeling and analysis of epidemics World scientific, (2009). [7] T. Zhang, Z. Teng, Global behaviour and permenance of SIRS epidemic models with time delay, Nonlinear Anal. Real World Appliations, 9 : 1409-1424, (2008). [8] Rukmini Kumar, Gilles Clermont, Yoram Vodovotz, Carson C. Chow, The dynamics of acute inflammation, J. Theo. Biol., 230 : 145-155, (2004). [9] Himanshu Manchandaa, Nora Seidel, Andi Krumbholz, Andreas Sauerbrei, Michaela Schmidtke, Reinhard Guthkea, Within-host influenza dynamics: A small-scale mathematical modeling approach, BioSyst. 118 : 51-59, (2014). [10] Xuhui Tan, Lingling Yuan, Jingjing Zhou, Yinan Zheng, Fen Yangd, Modeling the initial transmission dynamics of influenza A H1N1 in Guangdong Province, China, Int. J. Infec. Diseases, 17 : e479-e484, (2013). [11] Shailendra K. Guptaa, Shishir K. Gupta, Suchi Smita, Mugdha Srivastava, Xin Lai, Ulf Schmitz, Qamar Rahman, Olaf Wolkenhauer, Julio Vera, Computational analysis and modeling the effectiveness of Zanamivir targeting neuraminidase protein in pandemic H1N1 strains, Infect. Genetics and Evolution, 11 : 1072-1082, (2011). [12] Morb Mortal Wkly, Swine influenza A (H1N1) infection in two children-Southern California, Rep., 58 : 400-402, (2009). [13] Renato Casagrandi, Luca Bolzoni, Simon A. Levin, Viggo Andreasen, The SIRC model and influenza A, Math. Biosci. 200 : 152-169, (2006). [14] Masaya M. Saitoa, Seiya Imoto, Rui Yamaguchi, Hiroki Sato, Haruka Nakada, Masahiro Kami, Satoru Miyano, Tomoyuki Higuchi, Extension and verifcation of the SEIR model on the 2009 influenza A (H1N1) pandemic in Japan, Math. Biosci., 246 : 47-54, (2013). [15] Elissa J. Schwartza, Boseung Choib, Grzegorz A. Rempalac, Estimating epidemic parameters: Application to H1N1 pandemic data, Math. Biosci., Article in press, 1-6, (2015). [16] T.T. Wang , P. Palese, Unraveling the mystery of swine influenza virus, Cell 137 : 983-985, (2009). [17] R.J. Webby, S.L. Swenson, S.L. Krauss, P.J. Gerrish, S.M. Goyal, R.G. Webster, Evolution of swine H3N2 influenza viruses in the United States, J. Virol., 74 : 82438251, (2000). [18] A. Altmuller, M. Kunerl, K. Muller, V.S. Hinshaw, W.M. Fitch, C. Scholtissek, Genetic relatedness of the nucleoprotein (NP) of recent swine, turkey, and human infuenza A virus (H1N1) isolates, Virus Res., 22 : 79-87, (1992). [19] Xueyong Zhou, Zhen Guo, Analysis of an influenza A (H1N1) epidemic model with vaccination, Arab J. Math., 1 : 267-282, (2012).

23

[20] O. Diekmann, J.A.P. Heesterbeek, J.A.P. Metz, On the definition and the computation of teh basic reproduction ration R0 in models for infectious diseases in heterogeneous populations, J. Math. Biol., 28 : 365, (1990). [21] P. Van den Driessche, J. Watmough, Reproduction numbers and sub-threshold endemic equilibria for compartmental models of disease transmission, Math. Biosci. 180 : 29-48, (2002). [22] Y. Kuang, Delay differential equations with applications in population dynamics, Math. Sci. Eng., Academic Press, Boston, (1993). [23] N. MacDonald, Biological Delay Systems: Linear Stability Theory, Cambridge University, Cambridge, (1989). [24] J. Hale, Theory of Functional differential equations, Springer, New York, (1997). [25] J.K. Hale, G. Verduyn Lunel, Introduction to Functional Differential Equations, Springer, New York, (1993). [26] X. Yang, L. Chen, J. Chen , Permanence and positive periodic solution for the singlespecies non autonomous delay diffusive model, Comput. Math. Appl., 32(4): 109 116, (1996). [27] L. H. Erbe, H. I. Freedman, V. Sree Hari Rao, ‘Three-species food-chain models with mutual interference and time delays’, Math. Biosci., 80(1) : 57-80, (1986). [28] H. I. Freedman, V. S. H. Rao, ‘The trade-off between mutual interference and time lags in predator-prey systems’, Bull. Math. Biol., 45(6) : 991-1004, (1983). [29] F.A. Rihan, Sensitivity analysis of dynamical systems with time lags, J. Comput. Appl. Math., 151 : 445-463, (2003). [30] Z. Mukandavire, W. Garira, Sex-structured HIV/AIDS model to analyse the effects of condom use with application to Zimbabwe, J. Math. Biol., 54 : 669-699, (2007). [31] S.M. Tracht, S.Y. DelValle, J.M. Hyman, Mathematical modeling of the effectiveness of facemasks in reducing the spread of novel influenza A (H1N1), PLoS One, 5(2) : e9018, (2010). [32] A.R. Tuite, A.L. Greer, M. Whelan et. al., Estimated epidemiologic parameters and morbidity associated with pandemic H1N1 influenza, CMAJ. 182(2) : 131-136, (2009). [33] S. Leekha, N.L. Zitterkopf, M.J. Espy, T.F. Smith, R.L. Thompson et. al., Duration of influenza A virus shedding in hospitalized patients and implications for infection control, Infect. Control Hosp. Epidemiol., 28 : 1071-1076, (2007).

24