Mathematical analysis of swine influenza epidemic model with optimal control.

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 R c . The disease persists uniformly if R c > 1 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. © The JJIAM Publishing Committee and Springer Japan 2016.

Introduction

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–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, Bolle 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 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 , exhibits disease free equilibrium as well as the phenomenon of backward bifurcation. We will perform uncertainty and sensitivity analysis of 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.

Model formulation

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 , high risk susceptible individuals , susceptible individuals who were given antiviral drugs/vaccines , latent individuals , symptomatic individuals at early stage , symptomatic individuals at later stage , hospitalized individuals , treated individuals and recovered individuals , so thatThe high risk susceptible population includes pregnant women, children, health-care workers and providers (including all front-line workers), the elderly and other immune-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 bywhere is the effective contact rate, are the modification parameters accounting for relative infectiousness of individuals in the and T classes respectively in comparison to those in the class. High risk susceptible individuals acquire infection at a rate , where 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 (), and individuals in all epidemiological classes suffer natural death at a rate .

Fig. 1

Schematic diagram of the model (1)

Table 1

Description and nominal values of the model parameters

	Description	Value	Ref
π	Birth rate	1,119,583/80*365	Assumed
1/μ	Average human lifespan	80*365	Assumed
p	Fraction of susceptible individuals at high risk for contracting infection	0.4	Assumed
β	Effective contact rate for transmitting H1N1 influenza	0.9 (initial)	To be estimated
σ_L	Antiviral coverage rate for low risk susceptible individuals	0.3	[32]
σ_H	Antiviral coverage rate for high risk susceptible individuals	0.5	[32]
α	Rate at which latent individuals become infectious	1 / 1.9	[32]
τ_1	Treatment rate for individuals in the early stage of infection	1 / 5	[32]
τ_2	Treatment rate for individuals in the later stage of infection	1 / 3	[32]
ϕ_I_2	Recovery rate for symptomatic infectious individuals in the	1/5	[32]
	Later stage
ϕ_H	Recovery rate for hospitalized individuals	1/5	[32]
ϕ_T	Recovery rate treated individuals	1/3	[32]
η_1	Modification parameter (see text)	0.1	[32]
η_2	Modification parameter (see text)	1/2	[32]
η_3	Modification parameter (see text)	1.2	[32]
η_4	Modification parameter (see text)	1	[32]
θ_H	Modification parameter for infection rate of high risk	1.2	[32]
	Susceptible individuals
1-θ_P	Drug efficacy against infection	0.5	[32]
ψ	Hospitalization rate of individuals in I_2 class	0.5	[32]
γ	Progression rate from I_1 to I_2 classes	0.06	[32]
δ	Disease-induced death rate of individuals in I_2 class	1/100	[32]
θ_1δ	Disease-induced death rate for hospitalized individuals	1/100	[32]

Susceptible individuals who received prophylaxis can become infected at a reduced rate , where 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 . Individuals in the class receive antiviral treatment at a rate . These individuals progress to the later infectious class () at a rate . Similarly, individuals in the class are treated (at a rate ), hospitalized (at a rate ), recover (at a rate ) and suffer disease induced death (at a rate ). Hospitalized individuals, recover (at a rate ) and suffer disease-induced death (at a reduced rate , where accounts for the assumption that hospitalized individuals, in the H class, are less likely to die than uncapitalized infectious individuals in the class).

Finally, treated individuals recover (at a rate ). It is assumed that recovery confers permanent immunity against re-infection with H1N1. The model (1) is an extension of the model presented by Sharomi et al. [39] by:The basic qualitative properties of the model (1) will now be analyzed.

administering antivirals to susceptible individuals.

stratifying the susceptible population based on the risk of infection.

Basic properties of the model

Lemma 1

The closed setis positively-invariant and attracting with respect to the model (1).

Thus, in the region , the model is well-posed epidemiologically and mathematically [20]. Hence, it is sufficient to study the qualitative dynamics of the model (1) in .

Existence and stability of equilibria

Local stability of disease-free equilibrium (DFE)

The model (1) has a DFE, given bywith,Following [44], the linear stability of 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 , is given bywhere, is the spectral radius (dominant eigenvalue in magnitude) of the next generation matrix and Hence, using Theorem 2 of [44], the following result is established.

The epidemiological significance of the control reproduction number, , 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 to a value less than unity. Biologically-speaking, Lemma 2 implies that the H1N1 pandemic can be eliminated from the population (when ) if the initial sizes of the sub-populations in various compartments of the model are in the basin of attraction of the DFE . 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.

Lemma 2

The DFE, of the model (1), is locally asymptotically stable (LAS) if , and unstable if .

Global stability of DFE

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

Theorem 1

The DFE, , of the model (1) is GAS in if and ,

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 .

Theorem 2

If , then the disease is uniformly persistent: there exists an such thatfor all solutions of (1) with .

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

It is convenient to define the following quantitieswhere, . The quantities and represent the control reproduction numbers associated with the prophylactic () or therapeutic () use of antivirals in the community only.

Endemic equilibrium and backward bifurcation

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. Letrepresents an arbitrary endemic equilibrium of the model (1). Further, letbe 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 getwhere,The large number of new coefficients () 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 , of (5), is always positive, and the signs of the remaining parameters are dependent upon the value of . 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.

Theorem 3

If , then the model (1) exhibits a unique endemic equilibrium. However in case of no endemic equilibrium exists.

Proof

If , it is easy to verify that all the coefficients of (5) are positive. Therefore by Descartes’s sign rule there is no positive root and hence no endemic equilibrium.

However if , we have . Considering the other parameters, there is exactly one sign change among consecutive parameters. According to Descartes sign rule, we have exactly one positive root. Thus we have a unique endemic equilibrium. This completes the argument.

Now the question arises what happens when . In general, we cannot make any comment about this scenario. However the numerical simulations (using actual parameter estimates) suggest that there is no endemic equilibrium for and there exists a unique endemic equilibrium when .

The most important case arises when we assume that vaccination was perfect and that infection rate of high risk susceptibles is unity . These assumptions reduce (5) to a quadratic equation given belowwhere,

Lemma 3

If (along with perfect vaccination, and infection rate of high risk susceptibles, relative to that of low risk susceptibles, unity , then Backward Bifurcation is observed.

Clearly . Now if , we have . Thus, when and , we have two positive roots and hence we have Backward Bifurcation. The actual parameter estimates clearly show this phenomenon (Fig. 2).

Fig. 2

Backward bifurcation

Uncertainty and sensitivity analysis

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 , 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 as the model output. The parameter estimates and the assigned distributions are given in the Appendix.

The estimate of from uncertainty analysis is 1.32 with 95  CI (1.13, 1.81). The probability that 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 are (the ones having PRCC value above 0.5 or below 0.5) , , and . 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).

Fig. 3

Uncertainty analysis

Fig. 4

Sensitivity analysis

Optimal control

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 and 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 and be functions of time. The control set isThe goal is to minimize the cost function defined asThis 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 and are balancing cost factors due to scales and importance of the cost part of the objective functional, and are the maximum attainable values of the control, is the positive constants to keep balance in the population size in each class. We attempt to find an optimal control pair such thatThe 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 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 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 which minimizes J over . Also, there exists adjoint system of ’s such thatwhere

The transversality condition givesThe vaccination control is characterized as

Proof

It is easy to verify that the integrand of J is convex with respect to and . Also the solutions of the (1) are bounded as 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 obtainevaluated at the optimal control, which results in the stated adjoint system . The optimality condition isTherefore on the set , we obtainConsidering the bounds on and , 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

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 , 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 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 . 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 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 the infected people. 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.

Fig. 5

Simulations show the optimal treatment control and

Fig. 6

Infected population under different treatment control strategies

Fig. 7

Simulation of accumulated cost of different control strategies

Fig. 8

Optimal control strategies for different , a optimal control , b optimal control

Estimating the basic reproduction number for the 2009 H1N1 influenza outbreak in Manitoba, Canada

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 algorithm 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 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 .

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 .

Table 2

Estimates of the initial conditions for model (1)

Initial condition	Value
S_L(0)	1,119,583∗0.6
S_H(0)	1,119,583∗0.4
P(0)	0
L(0)	0
I_1(0)	1
I_2(0)	0
H(0)	0
T(0)	0
R(0)	0

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 of April 2009 to the of January 2010. The data is displayed in Fig. 9.

Fig. 9

The daily confirmed cases of H1N1 influenza in Manitoba, Canada

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 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 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.

Fig. 10

The best-fit trajectory of model (1) along with the daily confirmed cases of H1N1 influenza in Manitoba, Canada. Dashed lines represent post-antiviral and therapeutic interventions for 1st and 2nd wave

For the pre-control wave, we obtain ;

For the post-control wave, we obtain ;

For the pre-control wave, we obtain ;

For the post-control wave, we obtain .

Summary

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 mutually-exclusive epidemiological compartments, shows the following:

The disease-free equilibrium of the model is shown to be globally-asymptotically stable when .

The disease persists in the community if .

The model has a unique endemic equilibrium when the associated reproduction threshold .

Backward bifurcation is observed when all susceptible individuals are equally likely to acquire infection () and the vaccination is perfect .

Uncertainty analysis suggests that with mean value of and 95  confidence interval (1.13, 1.81), the disease will get endemic under the given conditions and data.

Sensitivity analysis recommends that we must estimate the parameter values of with great precision as they greatly influence the magnitude of which in turn controls the spread of the disease.

So far all the analysis is dependent on 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.

Table 3

Mean values of the model parameters with their assigned distributions

Parameter	Mean	Distribution	Parameter	Mean	Distribution
μ	3e-5	N	η_1	0.1	G
β	0.5	G	η_2	1.2	G
σ_L	0.08	N	η_3	0.2	G
σ_H	0.14	N	η_4	0.2	G
α	0.35	G	θ_H	1.2	U
τ_1	0.5	N	θ_P	0.3	N
τ_2	0.6	N	ψ_1	0.5	N
ϕ_I_2	0.2	N	γ_1	0.06	G
ϕ_H	0.2	N	δ	0.01	N
ϕ_T	0.25	N	θ_1δ	0.01	N

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