key: cord-0006186-8n85abve
authors: Imran, Mudassar; Malik, Tufail; Ansari, Ali R; Khan, Adnan
title: Mathematical analysis of swine influenza epidemic model with optimal control
date: 2016-01-21
journal: Jpn J Ind Appl Math
DOI: 10.1007/s13160-016-0210-3
sha: 5b3c68dbde7a6273c41ed5db9898399e32d4fbcb
doc_id: 6186
cord_uid: 8n85abve

A deterministic model is designed and used to analyze the transmission dynamics and the impact of antiviral drugs in controlling the spread of the 2009 swine influenza pandemic. In particular, the model considers the administration of the antiviral both as a preventive as well as a therapeutic agent. Rigorous analysis of the model reveals that its disease-free equilibrium is globally asymptotically stable under a condition involving the threshold quantity-reproduction number [Formula: see text] . The disease persists uniformly if [Formula: see text] and the model has a unique endemic equilibrium under certain condition. The model undergoes backward bifurcation if the antiviral drugs are completely efficient. Uncertainty and sensitivity analysis is presented to identify and study the impact of critical model parameters on the reproduction number. A time dependent optimal treatment strategy is designed using Pontryagin’s maximum principle to minimize the treatment cost and the infected population. Finally the reproduction number is estimated for the influenza outbreak and model provides a reasonable fit to the observed swine (H1N1) pandemic data in Manitoba, Canada, in 2009.

Influenza is a viral infection which is primarily transmitted through respiratory droplets produced by coughing and sneezing by an infected person [13, 41, 46] . Fever, cough, headache, muscle and joint pain, and reddened eyes are common symptoms of influenza. The infection, mostly prevalent in young children, elderly, pregnant women and patients with particular medical conditions such as chronic heart disease, has an initial incubation period of 1-4 days followed by a latency period of 2-10 days [33, 41, 46] .

For the past few centuries, influenza remains a serious threat to public health around the globe [2, 13, 41, 46] . It is reported that an influenza epidemic might have occurred in the sixteenth century [38] . During the past century, thousands of people lost their lives during three disastrous pandemics which include the Spanish flu, Asian flu, and Hong Kong flu during the year 1918, 1958, and 1968 respectively [2, 38, 46] . Most recently in 2014, India had reported 937 cases and 218 deaths from swine flu. The reported cases and deaths in 2015 had surpassed the previous numbers. The total number of confirmed cases crossed 33,000 with death of more than 2000 [43] .

In 2009, the world experienced the H1N1 influenza, also known as the Swine Influenza, an epidemic which led to over 16455 deaths globally [51] , 25 ,828 confirmed cases [36] and 428 deaths [37] in Canada alone. In Canada, this epidemic occurred in two waves, first wave occurred during the spring of 2009 whereas the second wave started around early October 2009 [5] and spanned about 3 months as it diminished by the end of January 2010 [31] . It is believed that genetic re-assortment involving the swine influenza virus lineages might be the main factor behind the emergence of H1N1 pandemic [22] . Various behavioral factors and chronic health conditions are considered to impose an increased risk of infection severity among H1N1 individuals. In particular, pregnant women, especially in their third trimester and infants are at increased risk of hospitalization and intensive care unit (ICU) admission [8, 10, 23, 50] . Moreover, individuals with pre-existing chronic conditions such as asthma, chronic kidney and heart diseases, and obesity are at far greater risk of death and ICU admission in comparison to healthy individuals [11, 28] . Additionally, geographical location may also impart a great risk to the severity of the infection. It is observed that people residing in remote, isolated communities and aboriginals, in the Canadian province of Manitoba, experienced severe illness due to the pandemic H1N1 infections [45] . Most recent outbreak of 2014, in India, had reported 937 cases and 218 deaths from swine flu. The reported cases and deaths in 2015 had surpassed the previous numbers. The total number of confirmed cases crossed 33,000 with death of more than 2000 [43] .

Like seasonal influenza, H1N1 influenza is also transmitted through coughing and sneezing by infected individuals. It may also be transmitted through contact with contaminated objects. Main symptoms of H1N1 infection include fever, cough, sore throat, body aches, headache, chills and fatigue [7] . Various preventive measures were taken to minimize H1N1 infection cases. These measures included (a) social exclusion such as closing schools during the epidemic and banning public gatherings etc., (b) mass vaccination mainly to high risk groups due to limited supply and (c) the use of antiviral drugs such as oseltamivir (Tamiflu) and zanamivir (Relenza) [6] . Although the epidemic breakout occurred in 2009, H1N1 influenza still poses significant challenges to public health worldwide [9, 14, 18, [47] [48] [49] .

In view of the serious consequences posed by the H1N1 epidemic on the public health, various mathematical models have been proposed and analyzed in order to understand the transmission dynamics of the H1N1 influenza [3, 4, 17, 19, 21, 22, 34, 39] . In particular, Boëlle et al. [3] estimated the basic reproduction number for the H1N1 influenza outbreak in Mexico. Gojovic et al. [19] proposed mitigation strategies to overcome H1N1 influenza. Sharomi et al. [39] presented an H1N1 influenza model that accounts for the role of an imperfect vaccine and antiviral drugs administered to infected individuals. In the past, the basic reproduction number R c has also been calculated by various researchers to control the disease by considering most-sensitive parameters but these strategies were not time-dependent. Lenhart proposed HIV models [15, 27] , using which time-dependent optimal controls were proposed. In this paper, we will propose an optimal treatment strategy to prevent the H1N1 influenza epidemic.

The purpose of the current study is to complement the aforementioned studies by formulating a new model in which the susceptible population is stratified according to the risk of infection. Our model also accounts for the impact of singular use of the antiviral drugs, administered as a preventive agents only i.e., given to susceptible individuals or, in addition, as a therapeutic agent i.e., given to individuals with symptoms of the pandemic H1N1 infection in the early stage of illness. In this study we will present a rigorous mathematical analysis of our model and show that our model, under suitable conditions involving the basic reproduction number R c , exhibits disease free equilibrium as well as the phenomenon of backward bifurcation. We will perform uncertainty and sensitivity analysis of R c in order to identify crucial model parameters. Additionally we will also propose optimal strategies that will eliminate the infection with minimal application of resources.

The paper is organized as follows. The model is formulated in Sect. 2 and rigorously analyzed in Sect. 3. We perform uncertainty and sensitivity analysis in Sect. 4 and using ideas from optimal control theory, we present an optimal treatment strategy in Sects. 5 and 6 marks the end of our paper with a discussion on our results.

The total human population at time t, denoted by N (t), is sub-divided into nine mutually-exclusive sub-populations comprising of low risk susceptible individuals (S L (t)), high risk susceptible individuals (S H (t)), susceptible individuals who were given antiviral drugs/vaccines (P(t)), latent individuals (L(t)), symptomatic individuals at early stage (I 1 (t)), symptomatic individuals at later stage (I 2 (t)), hospitalized individuals (H (t)), treated individuals (T (t)) and recovered individuals (R(t)), so that

The high risk susceptible population includes pregnant women, children, health-care workers and providers (including all front-line workers), the elderly and other immune- (1) compromised individuals. The rest of the susceptible population is considered to be at low risk of acquiring H1N1 infection. The model to be considered in this study is given by the following deterministic system of non-linear differential equations (a schematic diagram of the model is given in Fig. 1 , and the associated variables and parameters are described in Table 1) . 

In (1), π represents the recruitment rate into the population (all recruited individuals are assumed to be susceptible) and p is the fraction of recruited individuals who are at high risk of acquiring infection. Low risk susceptible individuals acquire infection at a rate λ, given by

where β is the effective contact rate, η i (i = 1, . . . , 4) are the modification parameters accounting for relative infectiousness of individuals in the L , I 2 , H 1 and T classes respectively in comparison to those in the I 1 class. High risk susceptible individuals acquire infection at a rate θ H λ, where θ H > 1 accounts for the assumption that high risk susceptible individuals are more likely to get infected in comparison to low risk susceptible individuals. Low (high) risk susceptible individuals receive antiviral at a rate σ L (σ H ), and individuals in all epidemiological classes suffer natural death at a rate μ. Susceptible individuals who received prophylaxis (P) can become infected at a reduced rate θ P λ, where 1 − θ P (0 < θ P < 1) is the efficacy of the antiviral in preventing infection. Individuals in the latent class become infectious at a rate α and move to the symptomatic early stage of infection I 1 . Individuals in the I 1 class receive antiviral treatment at a rate τ 1 . These individuals progress to the later infectious class (I 2 ) at a rate γ . Similarly, individuals in the I 2 class are treated (at a rate τ 2 ), hospitalized (at a rate ψ), recover (at a rate φ I 2 ) and suffer disease induced death (at a rate δ). Hospitalized individuals, recover (at a rate φ H ) and suffer disease-induced death (at a reduced rate θ 1 δ, where 0 < θ 1 < 1 accounts for the assumption that hospitalized individuals, in the H class, are less likely to die than uncapitalized infectious individuals in the I 2 class).

Finally, treated individuals recover (at a rate φ T ). It is assumed that recovery confers permanent immunity against re-infection with H1N1. The model (1) The basic qualitative properties of the model (1) will now be analyzed. Thus, in the region D, the model is well-posed epidemiologically and mathematically [20] . Hence, it is sufficient to study the qualitative dynamics of the model (1) in D.

3 Existence and stability of equilibria

The model (1) has a DFE, given by

Following [44] , the linear stability of E 0 can be established using the next generation operator method on system (1). The calculation is presented in Appendix. It follows that the control reproduction number, denoted by R c , is given by

where, ρ is the spectral radius (dominant eigenvalue in magnitude) of the next generation matrix F V −1 and K 1 = α + μ,

Hence, using Theorem 2 of [44] , the following result is established. The epidemiological significance of the control reproduction number, R c , which represents the average number of new cases generated by a primary infectious individual in a population where some susceptible individuals receive antiviral prophylaxis, is that the H1N1 pandemic can be effectively controlled if the use of antiviral can bring the threshold quantity (R c ) to a value less than unity. Biologically-speaking, Lemma 2 implies that the H1N1 pandemic can be eliminated from the population (when R c < 1) if the initial sizes of the sub-populations in various compartments of the model are in the basin of attraction of the DFE (E 0 ). To ensure that disease elimination is independent of the initial sizes of the sub-populations of the model, it is necessary to show that the DFE is globally asymptotically stable (GAS). This is considered below, for a special case.

The GAS property of E 0 is presented here. We proved the following result: following:

The proof is presented in the appendix. Theorem 1 shows that the H1N1 pandemic can be eliminated from the community if the use of antivirals can lead to R c ≤ R * .

Proof is presented in the appendix. The epidemiological implication of the result is that disease will prevail in the community if R c > 1.

It is convenient to define the following quantities

where, ω = Ω σ L =σ H =0 . The quantities R P and R T represent the control reproduction numbers associated with the prophylactic (R P ) or therapeutic (R T ) use of antivirals in the community only.

In order to find endemic equilibria of the basic model (1) (that is, equilibria where atleast one of the infected components of the model (1) are non-zero), the following steps are taken. Let

represents an arbitrary endemic equilibrium of the model (1). Further, let

be the associated force of infection at steady-state. Solving the equations of the model at steady-state yields:

Using (8) and the expression for λ * * we get

where,

The large number of new coefficients (A i ) are defined in the Appendix.

The positive endemic equilibria of the model (1) are obtained by solving the cubic (5) for λ * * and substituting the results (positive values of λ * * ) into (4). The coefficient a 0 , of (5), is always positive, and the signs of the remaining parameters are dependent upon the value of R c . Since we have a cubic equation to solve, it might not be possible to capture the dynamics of the system in general. Therefore, we will try to make natural and reasonable assumptions to study the system. Now we will discuss the solutions of (5) under different conditions. 

The most important case arises when we assume that vaccination was perfect (θ P = 0) and that infection rate of high risk susceptibles is unity (θ H = 1). These assumptions reduce (5) to a quadratic equation given below where, 

This section attempts to discuss the possible means of uncertainty in the values of our model parameters and to quantify the results. In a deterministic model, uncertainty can be generated by only input (initial conditions, parameters) variation and model parameters are among the most important part of the input data. The input data, particularly parameter estimates, can be uncertain due to the lack of sophisticated measuring methods, error in collecting and interpreting data and natural variations. Therefore, in order to study this uncertainty we perform global uncertainty and sensitivity analysis on our model output. Uncertainty analysis gives a qualitative measure of parameters which are critical in the model output. This analysis also presents the degree of confidence on the model parameters. The Latin hypercube sampling will be used to perform the uncertainty analysis on the model parameters presented in Table 1 . On the other hand, sensitivity analysis presents a qualitative measure of important parameters and quantifies the impact of each parameter on the model output, keeping the other parameters constant. Sensitivity analysis comprises of calculating the Partial Rank Correlation Coefficient (PRCC). The rank correlation coefficient is a standard method of quantifying the degree of monotonicity between an input parameter and the model output.

So far the dynamics of the model are determined by the threshold quantity R c , which determines the prevalence of the disease and therefore we consider it as our model output. Both uncertainty and sensitivity analyses are performed on the model parameters as the input data and R c as the model output. The parameter estimates and the assigned distributions are given in the Appendix.

The estimate of R c from uncertainty analysis is 1.32 with 95 % CI (1.13, 1.81). The probability that R c > 1 is 99 %. This suggests that H1N1 will get endemic under the assumed conditions. However, the time taken to reach that state could be different.

From the sensitivity analysis, it is apparent that most significant parameters to R c are (the ones having PRCC value above 0.5 or below −0.5) β, τ 1 , φ T and η 4 . This result implies that these parameters need to be estimated with utmost precision and accuracy in order to capture the transmission dynamics of H1N1 (Figs. 3, 4) .

Optimal control theory has found use in making decisions involving epidemic and biological models. The desired results and performance of the control functions depend on the different situations. Lenhart's HIV models [15, 27] used optimal control to design the treatment strategies. Lenhart [25] provides a very good example of deciding how to divide the efforts between two treatment strategies (case holding and case finding) of the two strain TB model. Yan [52] used an optimal isolation strategy to fight the SARS epidemic. In [24] Joshi formulated two control functions as coefficients of the ODE system representing treatment effects in a two drug regime in an HIV immunology model. The goal was to maximize the concentration of T cells while minimizing the toxic effects of the drug. The analytic and numerical results illustrated the level of two drugs to be used over the chosen time interval. The required balancing effect between two competing goals was well predicted by optimal control theory.

While formulating an optimal control problem, deciding how and where to introduce the control (through vaccination, drug treatment etc) in the system of differential equation is very important. The formulation of the optimal control problem must be a reasonable and practical representation of the situation to be considered. The form of the optimal control depends heavily on the system being analyzed and the objective functional to be optimized. We will consider a quadratic dependence on the control in the objective functional, this will be briefly justified when we discuss our formulation of the problem.

The existence and uniqueness of the optimal control in an optimality system can be established using the Lipchitz properties of the differential equations [16, 35] . For a detailed example, see the work of Lenhart [15, 16] . After establishing the existence and the uniqueness results, we can confidently continue to numerically solving the optimality system to get the desired optimal control.

Control theory for models defined as a system of differential equations was formulated by Pontryagin [35] . Since his time, the theory and its application have grown considerably in number. Pontryagin's maximum principle is a technique to optimize the performance criterion, (mostly cost and efficiency functional) which depends on different control parameters. The control parameters introduced are mostly functions of time appearing as coefficients in the model. The technique involves reducing a problem in which an optimal function over an entire time domain is to be determined to a standard optimization problem. This is achieved by appending an adjoint system of differential equations (with terminal boundary conditions) to the original model (state system) of differential equations (with initial conditions). The adjoint functions behave very similar to Lagrange multipliers (appending constraints to the function of several variables to be maximized or minimized). The adjoint variables maximize or minimize the state variables with respect to the desired objective functional. The details of the necessary conditions for the adjoint and optimal controls are presented here [16, 26, 35] . For the application of these results see [15] .

We will attempt to optimize the treatment rates τ 1 and τ 2 of our model (1). We aim to design an optimal treatment strategy that minimizes an objective functional taking into account both the cost and the number of infectious individuals. Let the treatment rates τ 1 and τ 2 be functions of time. The control set V is

The goal is to minimize the cost function defined as

This objective functional (performance criterion) involves the numbers of infected individuals as well as the cost of treatment. Since the total cost includes not only the consumption for every individual but also the cost of organization, management, and cooperation, a non-linear cost function seems to be the right choice. In this paper, a quadratic function is implemented for measuring the control cost by reference to literature in epidemics control [15, 25, 27, 52] . The coefficients W 1 and W 2 are balancing cost factors due to scales and importance of the cost part of the objective functional, b 1 and b 2 are the maximum attainable values of the control, X i is the positive constants to keep balance in the population size in each class. We attempt to find an optimal control pair (τ * 1 , τ * 2 ) such that

The existence of an optimal solution to the optimal control problem can be established by Theorem III.4.1 and its corresponding Corollary in [16] . The boundedness of solutions to the system (3.1) for the finite time interval is needed to establish these conditions and D satisfies it. Pontryagin's maximum principle [35] provides the necessary conditions to be satisfied by the optimal vaccination v(t). As stated earlier this reduces the problem to one of minimizing point wise a Hamiltonian, H , with respect to v

where k i represents the right hand side of the modified model's ith differential equation. Using the Pontryagin's maximum principle [35] and the optimal control existence result from [16] , we have Theorem 4 There exists a unique optimal pair τ * 1 (t), τ * 2 (t) which minimizes J over V. Also, there exists adjoint system of ξ i 's such that

The transversality condition gives

The vaccination control is characterized as

Proof It is easy to verify that the integrand of J is convex with respect to τ 1 and τ 2 .

Also the solutions of the (1) are bounded as N (t) ≤ μ for all time. Also it is clear the state system (1) has the Lipschitz property with respect to the state variables. With these properties and using Corollary 4.1 of [16] , we have the existence of the optimal control. Since we have the existence of the optimal vaccination control. Using the Pontryagin's maximum principle, we obtain

evaluated at the optimal control, which results in the stated adjoint system . The optimality condition is

Therefore on the set {t :

Considering the bounds on τ 1 and τ 2 , we have the characterizations of the optimal control as in (14) and (15) . Clearly the state and the adjoint functions are bounded. Also it is easily verifiable that state system and adjoint system has Lipschitz structure with respect to the corresponding variables, from here we obtain the uniqueness of the optimal control for sufficiently small time T [16, 35] . The uniqueness of the optimal control pair follows from the uniqueness of the optimality system, which consists of (3.1) and (3.5), with characterizations (3.6). Note here the uniqueness of the optimal control is only for a certain length of time. This restriction on the length of the time interval arises due to the opposite time orientations of (3.1), and (3.5); the state system has initial conditions and the adjoint system has final conditions. For example see [15] , this restriction is very common in control problems (see [27] ).

The optimality system consisting of 9 state system of differential equations and an adjoint system differential equations with S L (0) = S L0 , S H (0) = S H 0 , P(0) = P 0 , L(0) = L 0 , I 1 (0) = I 10 ,

Next, we discuss the numerical solutions of the optimality system and the corresponding optimal control pair. The optimal treatment strategy is obtained by solving the optimality system. An iterative scheme is used for solving the optimality system. We start to solve the state equations with a guess for the control pair (τ 1 , τ 2 ), using a forward Runge-Kutta method. The adjoint functions have final time conditions. Because of this transversality conditions on the adjoint functions , the adjoint equations are solved by a backward Runge-Kutta method using the current iteration solution of the state equations. Then, the controls are updated by using a convex combination of the previous control and the value from the characterizations. This process is repeated and iteration is stopped if the values of unknowns at the previous iteration are very close to the ones at the present iteration.

The parameter values used have R c > 1 when the model without control is considered. Thus the disease is not expected to die out without intervention strategies. Figure 5 represents the control strategy to be employed for the optimal results. This control strategy minimizes both the cost and the infected population (I 1 + I 2 ). It is well understood that in order to eradicate an epidemic, a large fraction of infected population needs to be treated. Therefore, an upper bound of 0.7 was chosen for treatment control. The optimal controls remain at higher values initially and then steadily decreasing to 0. In fact, at the beginning of simulated time, the optimal control is staying at higher values in order to treat as many infected as possible to prevent the susceptible individuals from getting infected and epidemic to break out. The steadily decreasing of the v is determined by the balance between the cost of the infected individuals and the cost of the controls. Figure 6 shows the infected (I 1 + I 2 )population for the optimal control and constant control. It is easy to see that the optimal control is much more effective for reducing the number of infected individuals and decreasing the time-span of the epidemic. As normally expected, in the early phase of the epidemic breakouts, keeping the control at their upper bound will directly lead to the decreasing of the number of Figure 7 shows the cost associated with the optimal and constant control strategy. It is clear the cost of optimal strategy is much less than the cost of constant strategy. Moreover Fig. 8 represents control strategies for different values of the effective contact rate β. It can be observed that increased effective contact rate requires us to implement control strategy for a longer duration of time.

Our purpose in this section is to estimate the transmission ability of the H1N1 virus during the two waves of the 2009 H1N1 epidemic in Manitoba, Canada, using epidemic data in the form of the daily confirmed cases of influenza. We will make use of the previous work of Cintron-Arias et al. [12] , who have developed a standard algo-rithm involving the use of a generalized least squares (GLS) scheme to fit epidemic data to a proposed deterministic ODE model. The GLS scheme is executed in the context of a statistical model and using the assumption that the epidemic data contains random noise. We then use statistical asymptotic theory to obtain the limiting probability distribution of the unknown model parameters. A detailed description of the methodology can be found in [12] . A simple but effective measure of the transmission ability of an infectious disease is given by the basic reproduction number, defined as the total number of secondary infections produced by introducing a single infective in a completely susceptible population. We will use the methods developed by Cintron-Arias et al. [12] to estimate the parameters of model (1) and thereby estimate the effective reproduction number R c for the H1N1 epidemic.

Estimates for several of the model parameters used in model (1) can be obtained from existing studies on H1N1 influenza. Table 1 lists these parameters along with reasonable estimates of their values.

The effective contact rate for transmitting H1N1 influenza β, which is a measure of the rate at which contact between an infective and a susceptible individual occurs and the probability that such contact will lead to an infection, is extremely difficult to determine directly. Consequently, we adopt an indirect approach, similar to previous studies such as [12] , by first finding the value of the parameter β for which the model (1) has the best agreement with the epidemic data, and then using the resultant parameter values to estimate R c .

As mentioned before, for the purpose of simulating model (1) we require knowledge of the initial conditions. It is possible to consider the initial conditions as model parameters along with the effective contact rate β and estimate values for all parameters.

Such an approach produces slightly unreliable results. This is explained by the fact that the available epidemic data is restricted to the daily confirmed cases of H1N1, while the optimization schemes that we will employ produce estimates for ten variables. There are thus too many degrees of freedom and the 'best-fit' may result in unrealistic estimates for the initial conditions. We will therefore use reasonable estimates for the initial conditions given in Table 2 and restrict ourselves to optimizing only the effective contact rate β.

The epidemic data to be used in this study is in the form of daily confirmed cases of H1N1 influenza in Manitoba, Canada, reported over approximately eight months, from the 24th of April 2009 to the 3rd of January 2010. The data is displayed in Fig. 9 .

The best-fit trajectory obtained by applying the GLS methodology in [12] to model (1) is displayed in the figure below. We have assumed that the first wave began on the 24th of April 2009 and control measures involving antiviral and therapeutic agents was started 48 days later. Similarly, we assume that the second wave began on the 3rd of October and control measures for the second wave began 50 days later. Furthermore, we assume that the antiviral and therapeutic control interventions that were implemented during the first wave also continue to be implemented before, during and after the second wave. Figure 10 shows data fitting of daily confirmed cases of H1N1 influenza in Manitoba, Canada. 

A deterministic model is designed and rigorously analyzed to assess the transmission dynamics and impact of antiviral drugs in curtailing the spread of disease of the 2009 swine influenza pandemic. The analysis of the model, which consists of nine mutuallyexclusive epidemiological compartments, shows the following:

(i) The disease-free equilibrium of the model is shown to be globally-asymptotically stable when R c ≤ Ω θ H . (iv) Backward bifurcation is observed when all susceptible individuals are equally likely to acquire infection (θ H = 1) and the vaccination is perfect (θ P = 0). (v) Uncertainty analysis suggests that with mean value of R c = 1.32 and 95 % confidence interval (1.13, 1.81), the disease will get endemic under the given conditions and data. (vi) Sensitivity analysis recommends that we must estimate the parameter values of β, τ 1 , φ T and η 4 with great precision as they greatly influence the magnitude of R c , which in turn controls the spread of the disease. (vii) So far all the analysis is dependent on R c for the spread of disease. Therefore, we have designed an optimal treatment strategy (using Pontryagin's maximum principle) to prevent the epidemic. It is observed that the optimal strategy is much more efficient in combating infection with minimal cost and resources.

Since N (t) ≥ 0, it follows using a standard comparison theorem [29] that

Therefore, N (t) ≤ π/μ if N (0) ≤ π/μ. This proves the positive invariance of D.

To prove that D is attracting, from (12), it is clear that d N dt < 0, whenever N (t) > π/μ. Thus, either the solution enters D in finite time, or N (t) approaches π/μ, and the variables denoting the infected classes approach zero. Hence, D is attracting and all solutions in R 9 + eventually enter D.

The matrices F (for the new infection terms) and V (of the transition terms) are given,

Reproduction number R c is given as

where, ρ denotes the spectral radius.

Proof Consider the model (1) with θ H = 1 and θ P = 0. Further, consider the Lyapunov function

where,

The Lyapunov derivative is given by (where a dot represents differentiation with respect to t)

so that, F = g 1 λ[S L (t)+ θ H S H (t)+ θ P P]− K 1 K 2 K 3 K 4 K 5 (η 1 L + I 1 +η 2 I 2 + η 3 H +η 4 T ),

It can be shown that By setting x(t) = (L(t), I 1 (t), I 2 (t), H (t), T (t)) T , equations for the infected components of (1) can be written as

where Y (x) = S L +θ H S H +θ P P N Ω F − V and it is clear that Y (E 0 ) = F − V . Also it is easy to check that Y (E 0 ) is irreducible. We will apply the Lemma A.4 in [1] to show that M is a uniform weak repeller. Since E 0 is a steady state solution, we can consider it to be a periodic orbit of period T = 1. P(t, x), the fundamental matrix of the solutions for (7) is e tY . Since the spectral radius of Y (E 0 ) = R c − 1 > 0, the spectral radius of e Y (E 0 ) > 1. So condition 2 of Lemma A.4 is satisfied. Taking x = E 0 , we get P(T, E 0 ) = e Y (E 0 ) which is a primitive matrix, because Y (E 0 ) is irreducible, as mentioned in Theorem A.12(i) [40] . This satisfies the condition 1 of Lemma A.4. Thus, M is a uniform weak repeller and disease is weakly persistent. By definition, M = ∂D. M is trivially closed and bounded relative to D and hence, compact. Therefore by Theorem 1.3 of [42] , we have that M is a uniform strong repeller and disease is uniformly persistent.

The coefficients are defined as follows A 1 = πθ H θ P A 2 = π(1 − p)[θ P (σ H + μ) + θ H μ] + πθ H p(θ P (σ L + μ) + μ) + πθ P A 1 A 3 = πμ(1 − p)(σ H + μ) + πθ H pμ(σ L + μ) + πθ P B A 4 = π(1 − p)[θ P (σ H + μ) + μθ H ] + π p[μ + θ P (σ L + μ)] + π A 1 A 5 = πμ(1 − p)(σ H + μ) + π pμ(σ L + μ) + π A 2 A 6 = π(1 − p)θ H θ P + π pθ P A 7 = σ L θ H (1 − p) + σ H p

The estimated parameters are presented in Table 3 . 

N , G and U represents the normal, gamma and uniform distribution respectively

Persistence and global stability in a selection-mutation sizestructured model

Optimal isolation control strategies and cost-effectiveness analysis of a two-strain avian influenza model

A preliminary estimation of the reproduction ratio for new influenza A(H1N1) from the outbreak in Mexico

Modeling influenza epidemics and pandemics: insights into the future of swine flu (H1N1)

Canada enters 2nd wave of H1N1

Centers for Disease Control and Prevention. Three reports of oseltamivir resistant novel influenza A (H1N1) viruses

Outbreak of swine-origin influenza A (H1N1) virus infection-Mexico

Centers for Disease Control and Prevention. Pregnant women and novel influenza A (H1N1): considerations for clinicians

Information on people at high risk of developing flu-related complications

The esitimation of the effective reproductive number from disease outbreak data

What happened in China during the 1918 influenza pandemic?

Optimizing chemotherapy in an HIV model. Electron

Deterministic and Stochastic Optimal Control

Pandemic Potential of a Strain of Influenza A (H1N1): Early Findings

GenBank sequences from 2009 H1N1 influenza outbreak

Modelling mitigation strategies for pandemic (H1N1)

The mathematics of infectious diseases

Early epidemiological assessment of the virulence of emerging infectious diseases: a case study of an influenza pandemic, modelling mitigation strategies for pandemic (H1N1)

H1N1 2009 influenza virus infection during pregnancy in the USA

Optimal control of an immunology model

Optimal control of treatments in a two-strain tuberculosis model

Dynamic Optimisation

Optimal control of the chemotherapy of HIV

Critically ill patients with 2009 influenza A (H1N1) infection in Canada

Stability Analysis of Nonlinear Systems

The Stability of Dynamical Systems. Regional Conference Series in Applied Mathematics

Confirmed Cases of H1N1 Flu in Manitoba

Assessing transmission control measures, antivirals and vaccine in curtailing pandemic influenza: scenarios for the US, UK, and The Netherlands

Epidemiological evidence of an early wave of the 1918 influenza pandemic in New York City

The first influenza pandemic in the new millennium: lessons learned hitherto for current control efforts and overall pandemic preparedness

The Mathematical Theory of Optimal Processes

Numerical study of a diffusive epidemic model of influenza with variable transmission coefficient

The Theory of the Chemostat

The 1918 influenza virus: a killer comes into view

Persistence under relaxed point-dissipativity

Reproduction numbers and sub-threshold endemic equilibria for compartmental models of disease transmission

Winnipeg Regional Health Authority Report. Outbreak of novel H1N1 influenza A virus in the Winnipeg health region

World Health Organization. Influenza A (H1N1)-update 49. Global Alert and Response (GAR)

World Health Organization. Statement by Director-General

Human infection with new influenza A (H1N1) virus: clinical observations from Mexico and other affected countries

World Health Organization. Pandemic (H1N1) 2009-update 81

Optimal quarantine and isolation strategies in epidemics control

Proof Adding the equations of the model (1) gives