\documentclass[reqno]{amsart}
\usepackage{hyperref}
\usepackage{graphicx}

\AtBeginDocument{{\noindent\small
\emph{Electronic Journal of Differential Equations},
Vol. 2011 (2011), No. 155, pp. 1--21.\newline
ISSN: 1072-6691. URL: http://ejde.math.txstate.edu or http://ejde.math.unt.edu
\newline ftp ejde.math.txstate.edu}
\thanks{\copyright 2011 Texas State University - San Marcos.}
\vspace{9mm}}

\begin{document}
\title[\hfilneg EJDE-2011/155\hfil H1N1 influenza pandemic]
{Deterministic model for the role of antivirals in controlling 
the spread of the H1N1 influenza pandemic}

\author[M. Imran, T. M. Malik, S. M. Garba \hfil EJDE-2011/155\hfilneg]
{Mudassar Imran, Mohammad T. Malik, Salisu M. Garba}  % in alphabetical order

\address{Mudassar Imran \newline
Department of Mathematics, Lahore University of Management Sciences, 
Lahore, Pakistan}
\email{mudassar.imran@gmail.com}

\address{Mohammad Tufail Malik \newline
Department of Mathematics, University of Manitoba, Winnipeg,
Manitoba, R3T 2N2, Canada}
\email{malik@cc.umanitoba.ca
\url{http://home.cc.umanitoba.ca/~malik/}}

\address{Salisu M. Garba \newline
Department of Mathematics and Applied Mathematics, University of Pretoria, 
Pretoria 0002, South Africa}
\email{Salisu.Garba@up.ac.za}

\thanks{Submitted July 6, 2011. Published November 18, 2011.}
\subjclass[2000]{92D30}
\keywords{H1N1 influenza; deterministic model; antiviral; Lyapunov function}

\begin{abstract}
 A deterministic model is designed and used to theoretically
 assess the impact of antiviral drugs in controlling the spread of
 the 2009 swine influenza pandemic. In particular, the model
 considers the administration of the antivirals 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 certain conditions involving
 having the associated reproduction number less than unity.
 Furthermore, the model has a unique endemic equilibrium if the
 reproduction threshold exceeds unity. The model provides a
 reasonable fit to the observed H1N1 pandemic data for the Canadian
 province of Manitoba. Numerical simulations of the model suggest
 that the singular use of antivirals as preventive agents only
 makes a limited population-level impact in reducing the burden of
 the disease in the population (except if the effectiveness level
 of this ``prevention-only'' strategy is high). On the other hand,
 the combined use of the antivirals (both as preventive and
 therapeutic agents) resulted in a dramatic reduction in disease
 burden. Based on the parameter values used in these simulations,
 even a moderately-effective combined treatment-prevention
 antiviral strategy will be sufficient to eliminate the H1N1
 pandemic from the province.
\end{abstract}

\maketitle
\numberwithin{equation}{section}
\newtheorem{theorem}{Theorem}[section]
\newtheorem{lemma}[theorem]{Lemma}
\newtheorem{remark}[theorem]{Remark}
\newtheorem{definition}[theorem]{Definition}
\allowdisplaybreaks

\section{Introduction}

Since its emergence in the Spring of 2009, the H1N1 Influenza  A
pandemic (also known as the Swine Influenza) continues to pose
significant challenges to public health around the world
\cite{cdc5,el,gen,who1,who2,who3}. For instance, the H1N1 pandemic
has accounted for 33,494 cases \cite{cancases} (8,669 hospitalized
cases including 1,472 admitted to ICU) and 429 deaths in Canada
(over 16,000 people have died globally \cite{who9}). The H1N1
pandemic  is believed to have resulted from a genetic reassortment
involving swine influenza virus lineages \cite{sty5}. Several
chronic conditions and behavioral and other factors have been
associated with increased risk of disease severity among
H1N1-infected individuals. Infants and pregnant women (especially
in the third trimester) are at increased risk of hospitalization
and intensive care unit (ICU) admission
\cite{cdc4,jamieson,uscdc1,whosalah1}. Furthermore, some studies
have shown that people with pre-existing chronic conditions (such
as asthma and other chronic lung diseases, chronic kidney and
heart diseases, obesity and conditions associated with immune
suppression) are prone to increased risk of death and ICU
admission \cite{kumar, uscdc2}. In the Canadian province of
Manitoba, Aboriginals and people residing in remote and isolated
communities are at increased risk of severe illness due to the
pandemic H1N1 infection \cite{wrha}.


Like in the case of seasonal flu, the H1N1 pandemic  is believed
to be spreading mainly through coughs and sneezes of people who
are infected with the pandemic  and touching of contaminated
objects.  H1N1 infection has been reported to cause a wide range
of flu-like symptoms, including fever, cough, sore throat, body
aches, headache, chills and fatigue. Additionally, many have
reported nausea, vomiting and/or diarrhea \cite{cdc1}.  The second
wave of H1N1  started around early October 2009 in Canada
\cite{cbc} and  data \cite{manitoba2}  suggests that it has
diminished by the end of January 2010 (after  peaking in November
2009).

 The pandemic H1N1 is being controlled using various measures
including social exclusion (e.g., closure of schools, banning
large public gatherings, etc.), mass vaccination (the limited
supply of the H1N1 vaccine, in the early stage of the second wave,
has forced its prioritization to high-risk groups, such as
children under the age of six, pregnant women, people with
weakened immune systems, etc.) and the use of antiviral drugs
(notably, {\it oseltamivir} (Tamiflu) and {\it zanamivir}
(Relenza)). It should be mentioned that a rare occurrence of
antiviral resistance (specifically to the use of Tamilflu) has
been reported \cite{resist}.

Some mathematical modeling studies have been carried out, aimed at
gaining insight into the transmission dynamics and control of the
H1N1 pandemic (see, for instance,
\cite{blower,marija,sty5,sty1,sty2,sty3,sty4,shar}). In
particular, Sharomi {\it et al.} \cite{shar} presented a model for
the spread of H1N1 that incorporates an imperfect vaccine and
antiviral drugs administered to infected individuals with disease
symptoms (the model in \cite{shar} stratifies the infected
population in terms of their risk of developing severe illness).
The purpose of the current study is to complement the
aforementioned studies by designing a new  model for assessing the
impact of singular use of the antiviral drugs, administered as a
preventive agents only (i.e., given to susceptible individuals)
or, as a therapeutic agent (i.e., given to individuals with
symptoms of the pandemic H1N1 infection in the early stage of
illness), in curtailing or mitigating the burden of the H1N1
pandemic in a population. An additional special feature of the
model to be designed is that it stratifies the susceptible
population according to  risk of infection. The model is used  to
theoretically compare the potential impact of the targeted
administration of the available antivirals (as a preventive agent
alone, or their combined use as both preventive and therapeutic
agent combined) in combatting the spread of the H1N1 pandemic in
the Canadian province of Manitoba.

The paper is organized as follows. The model is formulated  in
Section 2 and  rigorously analysed in Section 3. Numerical
simulations are reported  in Section 4.

\section{Model Formulation}\label{det_model}

\begin{figure}[tb]
\begin{center}
\includegraphics[width=0.9\textwidth]{fig1} % schematic
\end{center}
 \caption{Schematic diagram of the model \eqref{eqn1})}
\label{schematic}
\end{figure}

The total human population at time $t$, denoted by $N(t)$,  is
sub-divided into ten mutually-exclusive sub-populations of
low-risk susceptible individuals $(S_L(t))$, high-risk susceptible
individuals $(S_H(t))$, susceptible individuals who were given
antiviral drugs $(P(t))$, latent individuals $(L(t))$, infectious
individuals who show no disease symptoms $(A(t))$, symptomatic
individuals at early stage $(I_{1}(t))$, symptomatic individuals
at later stage $(I_{2}(t))$, hospitalized individuals $(H_1(t))$,
hospitalized individuals at ICU $(H_2(t))$, treated individuals
$(T(t))$ and recovered individuals $(R(t))$, so that
$$
N(t)=S_L(t)+S_H(t)+P(t)+L(t)+A(t)+I_{1}(t)+I_{2}(t)+H_1(t)
+H_2(t)+T(t)+R(t).
$$
The high-risk susceptible population includes pregnant women,
children, health-care workers and providers (including all
front-line workers), the elderly and other immuno-compromised
individuals. The rest of the susceptible population is considered
to be of 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 Figure \ref{schematic}, and the
associated variables and parameters are described in Table \ref{table1}):
\begin{equation} \label{eqn1}
\begin{gathered}
\frac{dS_L}{dt}=\Pi(1-p)-\lambda{S}_L-\sigma_LS_L-\mu{S_L},\\
\frac{dS_H}{dt}=\Pi{p}-\theta_H\lambda{S}_H-\sigma_HS_H-\mu{S_H},\\
\frac{dP}{dt}=\sigma_LS_L+\sigma_HS_H-\theta_P\lambda{P}-\mu{P},\\
\frac{dL}{dt}=\lambda(S_L+\theta_HS_H+\theta_PP)-(\alpha_1+\mu){L},\\
\frac{dA}{dt}=\alpha_1(1-f)L-(\phi_1+\mu-\delta)A,\\
\frac{dI_{1}}{dt}=\alpha_1f{L}-(\tau_1+\gamma_1+\mu)I_{1},\\
\frac{dI_{2}}{dt}=\gamma_1I_1-(\tau_2+\psi_1+\phi_2+\mu+\delta)I_{2},\\
\frac{dH_1}{dt}=\psi_1I_2-(\psi_2+\phi_3+\mu+\theta_1\delta)H_1,\\
\frac{dH_2}{dt}=\psi_2{H}_1-(\phi_4+\mu+\theta_2\delta)H_2,\\
\frac{dT}{dt}=\tau_1{I}_1+\tau_2{I}_2-(\phi_5+\mu)T,\\
\frac{dR}{dt}=\phi_1A+\phi_2I_2+\phi_3H_1+\phi_4H_2+\phi_5T-\mu{R}.
\end{gathered}
\end{equation} 
Here, the parameter $\Pi$ 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 $\lambda$, given by
\[
\lambda=\frac{\beta(\eta_1L+\eta_{2}A+I_{1}+\eta_{3}{I}_{2}
+\eta_{4}{H}_{1}+\eta_{5}{H}_{2}+\eta_{6}T)}{N},
\]
where $\beta$ is the effective contact rate, $\eta_i$
$(i=1,\dots, 6)$ are the modification parameters accounting for relative
infectiousness of individuals in the $L,A,I_2,H_1,H_2$ and $T$
classes in comparison with those in the $I_1$ class. high-risk
susceptible individuals acquire infection at a rate
$\theta_H\lambda$, where $\theta_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 antivirals at a rate
$\sigma_L$ ($\sigma_H$), and individuals in all epidemiological
classes suffer natural death at a rate $\mu$ (it is assumed, in
this study, that individuals in the latent class can transit
infection).

Susceptible individuals who received prophylaxis ($P$)  can become
infected  at a reduced rate $\theta_P\lambda$, where $1-\theta_P$
(with $0<\theta_P<1$) is the efficacy of the antiviral drugs in
preventing infection.  Individuals in the latent class become
infectious at a rate $\alpha_1$. A fraction, $f$, of these
individuals  display clinical symptoms of the disease and the
remaining fraction, $1-f$, does not show disease symptoms (and are
moved to the class $A$). Infectious individuals that show no
disease symptoms recover at a rate $\phi_1$ and die due to the
disease at a rate $\delta$. Individuals in the $I_1$ class receive
antiviral treatment at a rate $\tau_1$. These individuals progress
to the later infectious class ($I_2$) at a rate $\gamma_1$.
Similarly, individuals in the $I_2$ class are treated (at a rate
$\tau_2$), hospitalized (at a rate $\psi_1$), recover (at a rate
$\phi_2$) and suffer disease induced death (at a rate $\delta$).
Hospitalized individuals (not currently in ICU) are admitted into
ICU (at a rate $\psi_2$). Hospitalized individuals recover (at a
rate $\phi_3$) and suffer disease-induced death (at a reduced rate
$\theta_1\delta$, where $0<\theta_1<1$ accounts for the assumption
that hospitalized individuals, in the $H_1$ class, are less likely
to die than unhospitalized infectious individuals in the $I_2$
class).

Individuals in ICU recover (at a rate $\phi_4$) and die due  to
the H1N1  pandemic (at an increased rate $\theta_2\delta$, where
$\theta_2>1$ accounts for the assumption that those in ICU are
more likely to die than those in the $I_2$ class). Finally,
treated individuals recover (at a rate $\phi_5$). It is assumed
that recovery confers permanent immunity against re-infection with
H1N1.

The model \eqref{eqn1} is an extension of the model presented by
 Sharomi {\it et al}. \cite{shar} by:
\begin{enumerate}
\item[(i)] administering antivirals to susceptible individuals
 (only individuals with clinical symptoms of the disease are
  given antivirals in \cite{shar});
\item[(ii)] stratifying the susceptible population based on their
risk of acquiring infection.
\end{enumerate}


\begin{table}[t]\label{par_desc}
\caption{Description and nominal values of the model parameters}
\label{table1}
\scriptsize \centering
\begin{tabular}{|c|l|c|l|}
\hline
 & \textbf{Description} & \textbf{Value} & \textbf{Ref}\\
\hline
$\Pi$ & birth rate & 1119583/80*365 & \\
$1/\mu $ & average human lifespan & 80*365 &\\
$\beta$& probability for transmitting swine flu & 0.9 & assumed\\
$\sigma_L$ & antiviral coverage rate for low-risk susceptible individuals & 1/25 & assumed\\
$\sigma_H$ & antiviral coverage rate for high-risk susceptible individuals & 1/25& assumed\\
$\alpha_1$ & rate at which latent individuals become infectious & 0.35 & assumed\\
$f$ & fraction of latent individuals that progress to the& 0.2&\cite{Nuno2007}\\
& symptomatic class & &\\
$\tau_1$ & treatment rate for individuals in the early stage of infection & 0.5& assumed\\
$\tau_2$ & treatment rate for individuals in the later stage of infection & 0.5& assumed\\
$\phi_1$ & recovery rate for asymptomatic infectious individuals & 1/5&\cite{Nuno2007}\\
$\phi_2$ & recovery rate for symptomatic infectious individuals in the & 1/5&\cite{Nuno2007}\\
& later stage & &\\
$\phi_3$ & recovery rate for hospitalized individuals & 1/5&\cite{Nuno2007}\\
$\phi_4$ & recovery rate for individuals in ICU & 1/7& assumed\\
$\phi_5$ & recovery rate treated individuals & 1/7& assumed\\
$\eta_1$ & modification parameter (see text) & 0.1&\cite{Nuno2007}\\
$\eta_2$ & modification parameter (see text) & 1/2&\cite{Nuno2007}\\
$\eta_3$ & modification parameter (see text) & 1.2&\cite{Nuno2007}\\
$\eta_4$ & modification parameter (see text) & 1&\cite{Nuno2007}\\
$\eta_5$ & modification parameter (see text) & 0.01&\cite{Nuno2007}\\
$\eta_6$ & modification parameter (see text) & 0.2&\cite{Nuno2007}\\
$\theta_H$ & modification parameter for infection rate of high-risk & $\ge 1$ & assumed \\
& susceptible individuals & &\\
$1-\theta_P$ & drug efficacy in preventing infection & [0,1] &assumed\\
$\psi_1$ & hospitalization rate of individuals in $I_2$ class & 0.5&\cite{Nuno2007}\\
$\psi_2$ & rate of ICU admission of hospitalized individuals & 0.8&\cite{Nuno2007}\\
$\gamma_1$& progression rate from $I_1$ to $I_2$ classes & 0.06&\cite{Nuno2007}\\
$\delta$ & disease-induced death rate of individuals in $I_2$ class & 1/100&\cite{Nuno2007}\\
$\theta_1\delta$ & disease-induced death rate for hospitalized individuals & 1/100&\cite{Nuno2007}\\
$\theta_2\delta$ & disease-induced death rate for individuals in ICU & 1/100&\cite{Nuno2007}\\
\hline
\end{tabular}
\end{table}

The basic qualitative properties of the model \eqref{eqn1} will
now be analyzed.

\subsection{Basic properties of the model}

\begin{theorem}\label{thm1}
The variables of the model \eqref{eqn1} are non-negative for all time
$t>0$. In other words, the solutions of the model \eqref{eqn1} with
positive initial data will remain positive for all $t>0$.
\end{theorem}

\begin{proof}
Let,
\begin{align*}
t_1&=\sup\big\{t>0:S_L>0,\, S_H>0,\, P>0,\, L>0,\, A>0,\\
&\quad  I_1>0,\, I_2>0,\, H_1>0,\, H_2>0,\, T>0,\, R>0\big\}.
\end{align*}
Thus, $t_1>0$. The first equation of the model \eqref{eqn1} can be
rewritten as
\[
\frac{d}{dt}\big\{S_L(t)\exp\big[(\sigma+\mu)t
+\int^t_0\lambda(\psi)d\psi\big]\big\}
=\Pi (1-p)\exp\big[(\sigma+\mu)t+\int^t_0\lambda(\psi)d\psi\big],
\]
so that
\begin{align*}
&S_L(t_1){\exp}\big[(\sigma+\mu){t_1}+\int^{t_1}_0\lambda(\psi)
 d\psi\big]-S_L(0)\\
&=\int^{t_1}_0\big[\Pi(1-p){\exp}\Big((\sigma+\mu){y}
 +\int^y_0\lambda d\psi\Big)\big]dy.
\end{align*}
Therefore,
\begin{align*}
S_L(t_1)&\geq{S_L(0)}{\exp}\big[-(\sigma+\mu){t_1}
-\int^{t_1}_0\lambda d\psi\big]
+{\exp}\big[-(\sigma+\mu){t_1}
 -\int^{t_1}_0\lambda d\psi\big]\\
&\quad\times \int^{t_1}_0\big\{\Pi(1-p){\exp}\big[(\sigma+\mu){y}
+\int^y_0\lambda d\psi\big]\big\}dy
>0.
\end{align*}
Similarly, it can be shown that $ S_H(t)>0$, $P(t)>0$, $L(t)>0$,
$A(t)>0$, $I_1(t)>0$, $I_2(t)>0$, $H_1(t)>0$, $H_2(t)>0$,
$T(t)>0$, $R(t)>0$ for all $t>0$.
\end{proof}

Theorem \ref{thm1} can also be proven by applying a result from
Appendix A of \cite{thieme}.
Adding all the equations in the system \eqref{eqn1} gives
\begin{equation} %2
\frac{dN}{dt}=\Pi-\mu{N}-\delta(I_2+\theta_1H_1+\theta_2H_2),
\end{equation}
so that,
\begin{equation} %3
\frac{dN}{dt}\leq\Pi-\mu{N}.
\end{equation}
Since $N(t)\geq0$, it follows, using Gronwall inequality,
that
\[
N(t)\leq{N(0)e^{-\mu t}+\frac{\Pi}{\mu}(1-e^{-\mu t})}.
\]
Hence,
\begin{equation}
N(t)\leq{\Pi}/{\mu}\quad \text{if }  N(0)\leq{\Pi}/{\mu}.
\end{equation}
This result is summarized below.

\begin{lemma} \label{lem1}
The following biologically-feasible region of model \eqref{eqn1}
is positively-invariant:
\begin{align*}
\mathcal{D}=\Bigl\{&(S_L,S_H,P,L,A,I_1,I_2,H_1,H_2,T,R)\in\mathbb{R}_+^{11}:\\
&S_L+S_H+P+L+A+I_1+I_2+H_1+H_2+T+R\leq\frac{\Pi}{\mu}\Big\}.
\end{align*}
\end{lemma}

Thus, in the region $\mathcal{D}$, the model is well-posed epidemiologically
and mathematically \cite{hethcote2000}. Hence, it is sufficient to
study the qualitative dynamics of the model \eqref{eqn1} in $\mathcal{D}$.

\section{Existence and Stability of Equilibria}

\subsection{Local stability of disease-free equilibrium (DFE)}

The model \eqref{eqn1} has a DFE, given by
\begin{align*}
\mathcal{E}_0&=(S_L^*,S_H^*,P^*,L^*,A^*,I_1^*,I_2^*,H_1^*,H_2^*,T^*,R^*)\\
&=\Big(S_L^*,S_H^*,\frac{\sigma_LS_L^*+\sigma_HS_H^*}{\mu},
 0,0,0,0,0,0,0,0\Big),
\end{align*}
with
$$
S_L^*=\frac{\Pi(1-p)}{\sigma_L+\mu}, \quad
S_H^*=\frac{\Pi{p}}{\sigma_H+\mu}.
$$
Following \cite{vanden2002}, the linear stability of $\mathcal{E}_0$
can be established using the next generation operator method on
system \eqref{eqn1}. The matrices $F$ (for the new infection terms)
and $V$ (of the transition terms) are
 given, respectively, by
$$
F=\begin{bmatrix}
  \beta\eta_1\Omega & \beta\eta_2\Omega & \beta\Omega
  &  \beta\eta_3\Omega & \beta\eta_4\Omega & \beta\eta_5\Omega
  & \beta\eta_6\Omega\\
  0 & 0 & 0 & 0 & 0 & 0 & 0\\
  0 & 0 & 0 & 0 & 0 & 0 & 0\\
  0 & 0 & 0 & 0 & 0 & 0 & 0\\
  0 & 0 & 0 & 0 & 0 & 0 & 0\\
  0 & 0 & 0 & 0 & 0 & 0 & 0\\
  0 & 0 & 0 & 0 & 0 & 0 & 0\\
\end{bmatrix},
$$
$$
V=\begin{bmatrix}
  K_1 & 0 & 0 & 0 & 0 & 0 & 0\\
  -\alpha_1(1-f) & K_2 & 0 & 0 & 0 & 0 & 0\\
  -\alpha_1{f} & 0 & K_3 & 0 & 0 & 0 & 0\\
  0 & 0 & -\gamma_1 & K_4 & 0 & 0 & 0\\
  0 & 0 & 0 & -\psi_1 & K_5 & 0 & 0\\
  0 & 0 & 0  & 0 & -\psi_2 & K_6 & 0\\
  0 & 0 & -\tau_1 & -\tau_2 & 0 & 0 & K_7\\
\end{bmatrix},
$$
where, $\Omega=(S_L^*+\theta_HS_H^*+\theta_PP^*)/N^*$,
$K_1=\alpha_1+\mu$, $K_2=\phi_1+\mu+\delta$,
$K_3=\tau_1+\gamma_1+\mu$, $K_4=\tau_2+\psi_1+\phi_2+\mu+\delta$,
$K_5=\psi_2+\phi_3+\mu+\theta_1\delta$,
$K_6=\phi_4+\mu+\theta_2\delta$,
$K_7=\phi_5+\mu$  and $N^*=\frac{\Pi}{\mu}$.

It follows then that the \emph{control reproduction number},
denoted by $\mathcal{R}_c $, is given by
\begin{equation} \label{eqn2}
\begin{split}
\mathcal{R}_c &=\rho(FV^{-1})\\
&=\frac{\beta\Omega}{K_1}\big\{\eta_1+\big\{\frac{\alpha_1
\eta_2(1-f)}{K_2}+\frac{\alpha_1 f}{K_3K_4}\big[K_4+\eta_3\gamma_1
+\frac{\eta_6}{K_7}(\tau_1K_4+\tau_2\gamma_1)+Q\big]\big\}\big\},
\end{split}
\end{equation}
where,  $Q=(\psi_1\gamma_1(K_6\eta_4+\psi_2\eta_5))/(K_5K_6)$ and
$\rho$ is the spectral radius (dominant eigenvalue in
magnitude) of the next generation matrix $FV^{-1}$. Hence, using
\cite[Theorem 2]{vanden2002}, the following result is established.

\begin{lemma} \label{DFE-LAS}
The DFE, $\mathcal{E}_0$, of the model
\eqref{eqn1}, is locally asymptotically stable (LAS) if $\mathcal{R}_c <1$, and
unstable if $\mathcal{R}_c >1$.
\end{lemma}

The epidemiological significance of the \emph{control reproduction
number}, $\mathcal{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 $(\mathcal{R}_c )$ to a value less
than unity. Biologically-speaking, Lemma \ref{DFE-LAS} implies
that the H1N1 pandemic  can be eliminated from the population
(when $\mathcal{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 $(\mathcal{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.


It is convenient to define the following quantities
\begin{equation} \label{eqn4}
\begin{split}
\mathcal{R}_{P}
&=\mathcal{R}_c \big|_{\tau_1=\tau_2=0}
 =\frac{\beta\Omega}{K_1}\big\{\eta_1
 +\big\{\frac{\alpha_1\eta_2(1-f)}{K_2}+\frac{\alpha_1f}{K_3K_4}
 [K_4+\eta_3\gamma_1+Q]\big\}\big\},\\
\mathcal{R}_T
&=\mathcal{R}_c \big|_{\sigma_L=\sigma_H=0}\\
&=\frac{\beta\omega}{N^*K_1}\big\{\eta_1
 +\big\{\frac{\alpha_1\eta_2(1-f)}{K_2}+\frac{\alpha_1f}{K_3K_4}
 [K_4+\eta_3\gamma_1+\frac{\eta_6}{K_7}(\tau_1K_4+\tau_2\gamma_1)
 +Q]\big\}\big\}.
\end{split}
\end{equation}
The quantities, $\mathcal R_P$ and $\mathcal R_T$, represent the
control reproduction numbers associated with the singular
prophylactic ($\mathcal{R}_P$) or therapeutic ($\mathcal{R}_T$)
use of antivirals in the community, respectively.


\begin{figure}[tb]
\begin{center}
\includegraphics[width=0.6\textwidth]{fig2} %dfe.eps
\end{center}
\caption{Simulations of the model \eqref{eqn1} showing the total number of
infected individuals as a function of time, using various initial conditions. 
Parameter values used are as in Table \ref{table1} with $p=0.5$, 
$\theta_H=1$, $\theta_P=0$, $f=0.5$, $\phi_4=1/7$, $\phi_5=1/7$, $c=1$, 
$\beta=0.4$, $\eta_1=0.1$ (so that, $\mathcal{R}_c =0.004<1$)}
\label{dfe}
\end{figure}

\subsection{Global stability of DFE}

The GAS property of the DFE ($\mathcal{E}_0$) of the model is
considered for the special case where all the susceptible
individuals (i.e., there is no stratification of the susceptible
individuals based on risk of infection) are equally likely to
acquire infection (i.e., $\theta_H=1$) and the susceptible
individuals who received antivirals are fully protected against
infection (i.e., $\theta_P=0$). Define
$\tilde{\mathcal{R}_c} = \mathcal{R}_c \big|_{\theta_{H=1}, \theta_{P=0}}$.
We claim the following result.

\begin{theorem}\label{thm_dfe}
The DFE, $\mathcal{E}_0$, of the model \eqref{eqn1} with $\theta_H=1$
and $\theta_P=0$, is GAS in $\mathcal{D}$ if $\tilde{\mathcal{R}}_c \leq1$.
\end{theorem}

\begin{proof}
Consider the model \eqref{eqn1} with $\theta_H=1$ and $\theta_P=0$.
Further, consider the Lyapunov function
$$
\mathcal{F}=g_1L+g_2A+g_3I_1+g_4I_2+g_5H_1+g_6H_2+g_7T,
$$
where
\begin{gather*}
\begin{aligned}
g_1&=\eta_1K_2K_3K_4K_5K_6K_7+\eta_2\alpha_1K_3K_4K_5K_6K_7(1-f)
 +\alpha_1fK_2K_4K_5K_6K_7\\
&\quad +\eta_3\gamma_1\alpha_1fK_2K_5K_6K_7
 +\eta_4\gamma_1\psi_1\alpha_1fK_2K_6K_7
 +\eta_5\psi_1\psi_2\gamma_1\alpha_1fK_2K_7\\
&\quad +\eta_6\alpha_1\tau_1fK_2K_4K_5K_6+\eta_6\alpha_1\tau_2\gamma_1fK_2K_5K_6,
\end{aligned}\\
g_2=\eta_2K_1K_3K_4K_5K_6K_7,\\
\begin{aligned}
g_3&=K_1K_2(K_4K_5K_6K_7+\eta_3\gamma_1K_5K_6K_7+\eta_4\psi_1\gamma_1K_6K_7
 +\eta_5\psi_1\psi_2\gamma_1K_7\\
&\quad +\eta_6\tau_1K_4K_5K_6+\eta_6\tau_2\gamma_1K_5K_6),
\end{aligned}\\
g_4=K_1K_2K_3(\eta_3K_5K_6K_7+\eta_4\psi_1K_6K_7+\eta_5\psi_1\psi_2K_7+\eta_6\tau_2K_5K_6),\\
g_5=K_1K_2K_3K_4K_7(\eta_4K_6+\eta_5\psi_2),\\
g_6=\eta_5K_1K_2K_3K_4K_5K_7,\\
g_7=\eta_6K_1K_2K_3K_4K_5K_6.
\end{gather*}
The Lyapunov derivative is given by (where a dot represents
differentiation with respect to $t$)
\begin{align*}
\mathcal{\dot{F}}
&= g_1\dot{L}+g_2\dot{A}+g_3\dot{I}_1+g_4\dot{I}_2+g_5\dot{H}_1
 +g_6\dot{H}_2+g_7\dot{T}\\
&= g_1[\lambda(S_L+S_H+\theta_PP)-(\alpha_1+\mu){L}]
 +g_2[\alpha_1(1-f)L-(\phi_1+\mu)A]\\
&\quad+g_3[\alpha_1f{L}-(\tau_1+\gamma_1+\mu)I_{1}]
 +g_4[\gamma_1I_1-(\tau_2+\psi_1+\phi_2+\mu+\delta)I_{2}]\\
&\quad+g_5[\psi_1I_2-(\psi_2+\phi_3+\mu+\theta_1\delta)H_1]
 +g_6[\psi_2{H}_1-(\phi_4+\mu+\theta_2\delta)H_2]\\
&\quad+g_7[\tau_1{I}_1+\tau_2{I}_2-(\phi_5+\mu)T],
\end{align*}
so that
\begin{align*}
\mathcal{\dot{F}}
&= g_1\lambda[S_L(t)+S_H(t)]-K_1K_2K_3K_4K_5K_6K_7\big(\eta_1L+\eta_{2}A
 +I_{1}+\eta_{3}{I}_{2}\\
 &\quad +\eta_{4}{H}_{1}+\eta_{5}{H}_{2}\big)\\
&= g_1\lambda[S_L(t)+S_H(t)]-K_1K_2K_3K_4K_5K_6K_7
 \frac{\lambda N}{\beta}\\
&\leq {g}_1\lambda{N}-K_1K_2K_3K_4K_5K_6K_7
 \frac{\lambda{N}}{\beta},\,
\end{align*}
since $S_L(t)+S_H(t)\leq{N(t)}$  in $\mathcal{D}$.

It can be shown that $g_1=\frac{\beta}K_1K_2K_3K_4K_5K_6K_7$.
Hence,
\begin{align*}
\dot{\mathcal{F}}
&\leq \frac{\tilde{\mathcal{R}_c}}{\beta}K_1K_2K_3K_4K_5K_6K_7\lambda{N}
 -K_1K_2K_3K_4K_5K_6K_7\frac{\lambda{N}}{\beta}\\
&= K_1K_2K_3K_4K_5K_6K_7\frac{\lambda{N}}{\beta}
(\tilde{\mathcal{R}_c}-1).
\end{align*}
Thus, $\dot{\mathcal{F}}\leq0$ if $\tilde{\mathcal{R}_c}\leq1$  with
$\dot{\mathcal{F}}=0$ if and only if $L=A=I_1=I_2=H_1=H_2=T=0$.
Further, the largest compact invariant set in
$\{(S_L,S_H,P,L,A,I_1,I_2,H_1, H_2, T,\\ R)\in \mathcal{D}:
\dot{\mathcal{F}}=0\}$ 
is the singleton $\{\mathcal{E}_0\}$. It
follows from the LaSalle Invariance Principle 
(\cite[Chapter 2, Theorem 6.4]{lasalle1976}) 
that every solution to the equations in
\eqref{eqn1} with $\theta_H=1$ and $\theta_P=0$ and  with initial
conditions in $\mathcal{D}$ converge to
DFE, $\mathcal{E}_0$, as $t\to \infty$. That is,
$[L(t),A(t),I_1(t),I_2(t),H_1(t),H_2(t),T(t)]\to (0,0,0,0,0,0,0)$
as $t\to \infty$. Substituting $L=A=I_1=I_2=H_1=H_2=T=0$ into
the first three equations of the model \eqref{eqn1} gives
$S_L(t)\to {S_L^*},\,S_H(t)\to {S_H^*}$ and
$P(t)\to {P^*}$ as $t\to \infty$. Thus,
$[S_L(t),S_H(t),P(t),L(t),A(t),I_1(t),I_2(t),H_1(t),H_2(t),
T(t),R(t)]\to (S_L^*,S_H^*,P^*,0,0,0,0,0,0,0,0)$
as $t\to \infty$ for $\tilde{\mathcal{R}_c}\leq1$, so that
the DFE, $\mathcal{E}_0$, is GAS in
$\mathcal{D}$ if $\tilde{\mathcal{R}_c}\leq1$.
\end{proof}

Theorem \ref{thm_dfe} shows that the H1N1 pandemic can be eliminated
from the comunity if the use of antivirals can lead to
$\tilde{\mathcal{R}_c}\leq1$ for the case when $\theta_H=1$ and
$\theta_p=0$. This result is illustrated in Figure \ref{dfe}.


\subsection{Existence and stability of an endemic equilibrium
point (EEP)}

In order to find endemic equilibria of the basic model \eqref{eqn1}
(that is, equilibria where the infected components of the model
\eqref{eqn1} are non-zero), the following steps are taken. Let
$$
E_1=(S_L^{**},S_H^{**},P^{**},L^{**},A^{**},I_1^{**},I_2^{**},
H_1^{**},H_2^{**},T^{**},R^{**})
$$
represents an arbitrary endemic equilibrium of the model \eqref{eqn1}.
Further, let
\[
\lambda^{**}=\frac{\beta(\eta_1L^{**}+\eta_{2}A^{**}
+I_{1}^{**}+\eta_{3}{I}_{2}^{**}+\eta_{4}{H}_{1}^{**}
+\eta_{5}{H}_{2}^{**}+\eta_{6}T^{**})}{N^{**}},
\]
be the associated force of infection at steady-state.
Solving the equations of the model at steady-state gives
\begin{gather} 
S_L^{**} = \frac{\Pi(1-p)}{\lambda^{**}+\sigma_L+\mu},\quad
S_H^{**}=\frac{\Pi{p}}{\lambda^{**}+\sigma_H+\mu}, \nonumber \\
P^{**}=\frac{\Pi\{[\sigma_L(1-p)+\sigma_H{p}](\lambda^{**}+\mu)
 +\sigma_L\sigma_H\}}{(\lambda^{**}+\sigma_H+\mu)(\lambda^{**}
 +\sigma_L+\mu)\mu}, \nonumber \\
L^{**}= \frac{\lambda^{**}\Pi[\sigma_H(1-p)+\sigma_L
 +\lambda^{**}+\mu]}{K_1(\lambda^{**}+\sigma_H+\mu)
  (\lambda^{**}+\sigma_L+\mu)}, \nonumber \\
A^{**}= \frac{\lambda^{**}\Pi\alpha_1(1-f)[\sigma_H(1-p)
 +\sigma_L+\lambda^{**}+\mu]}{K_1K_2(\lambda^{**}
  +\sigma_H+\mu)(\lambda^{**}+\sigma_L+\mu)}, \nonumber \\
I_1^{**}= \frac{\lambda^{**}\Pi\alpha_1{f}[\sigma_H(1-p)
 +\sigma_L+\lambda^{**}+\mu]}{K_1K_3(\lambda^{**}
  +\sigma_H+\mu)(\lambda^{**}+\sigma_L+\mu)}, \nonumber \\
I_2^{**}= \frac{\lambda^{**}\Pi\alpha_1{f}\gamma_1[\sigma_H(1-p)
        +\sigma_L+\lambda^{**}+\mu]}{K_1K_3K_4(\lambda^{**}
        +\sigma_H+\mu)(\lambda^{**}+\sigma_L+\mu)}, \nonumber \\
H_1^{**}= \frac{\lambda^{**}\Pi\alpha_1{f}\gamma_1\psi_1
 [\sigma_H(1-p)+\sigma_L+\lambda^{**}+\mu]}{K_1K_3K_4K_5(\lambda^{**}
 +\sigma_H+\mu)(\lambda^{**}+\sigma_L+\mu)}, \nonumber  \\
H_2^{**}= \frac{\lambda^{**}\Pi\alpha_1{f}\gamma_1\psi_1\psi_2
  [\sigma_H(1-p)+\sigma_L+\lambda^{**}+\mu]}{K_1K_3K_4K_5K_6
  (\lambda^{**}+\sigma_H+\mu)(\lambda^{**}+\sigma_L+\mu)}, \nonumber \\
T^{**}=\frac{\lambda^{**}\Pi\alpha_1{f}(K_4\tau_1
 +\gamma_1\tau_2)[\sigma_H(1-p)+\sigma_L+\lambda^{**}
  +\mu]}{K_1K_3K_4K_7(\lambda^{**}+\sigma_H+\mu)(\lambda^{**}
   +\sigma_L+\mu)}, \nonumber \\
R^{**}=\frac{\lambda^{**}\Pi\alpha_1[\sigma_H(1-p)+\sigma_L
 +\lambda^{**}+\mu]X}{\mu{K_1}K_2K_3K_4K_5K_6K_7(\lambda^{**}
  +\sigma_H+\mu)(\lambda^{**}+\sigma_L+\mu)},
  \label{eq2}
\end{gather}
where
\begin{align*}
X&=K_5K_6 [ K_7K_3K_4\phi_1{\it (1-f)}
 +K_2\phi_5f(K_4\tau_1+\gamma_1\tau_2)]\\
&\quad -fk_2\gamma_1K_7 ( \psi_2\phi_4\psi_
{{1}}+K_6\phi_3\psi_1+K_6\phi_2K_5 ).
\end{align*}
It can be shown that the non-zero
equilibria of the model satisfy the following quadratic (in terms of
$\lambda^{**}$):
 \begin{equation}\label{eq3}
 a_0(\lambda^{**})^2+b_0\lambda^{**}+c_0=0,
\end{equation}
where,
\begin{align*}
a_0&= \mu\,K_2K_6K_5K_4K_3K_7+\alpha_1{\it(1-f)}
 K_7K_3K_4K_5K_6\mu+\alpha_1K_7\psi_1\gamma_1fk_2K_6\phi_3\\
&\quad+\psi_1\gamma_1\alpha_1fk_7K_6K_2\mu
 +\alpha_1K_5K_6K_2 \phi_5K_4\tau_1f
 +\alpha_1fk_5K_6K_2\mu\,\gamma_1\tau_2\\
&\quad +\alpha_1fk_7K_4K_5K_6K_2\mu
 +\alpha_1K_7K_3K_4K_5K_6\phi_1{\it(1-f)}
 +\psi_2\psi_1\gamma_1\alpha_1fk_7K_2\mu\\
&\quad +\alpha_1K_5K_6K_2  \phi_5f\gamma_1\tau_2
 +\alpha_1K_7\psi_1\gamma_1fk_2\phi_4\psi_2
 +\gamma_1\alpha_1fk_7K_5K_6K_2\mu \\
&\quad +\alpha_1K_7 \gamma_1fk_2K_6\phi_2K_5
 +\alpha_1fk_5K_6K_2\mu\,K_4\tau_1,
\end{align*}

\begin{align*}
b_0&= K_6K_7K_5K_4K_3K_1K_2\mu - \mu\beta \eta_1K_2K_3K_4K_5K_6K_7\\
&\quad - \mu\beta  \eta_2\alpha_1(1-f)K_3K_4K_5K_6K_7
 -\mu\beta \alpha_1fK_2K_4K_5K_6K_7\\
&\quad -\mu\beta \eta_3\gamma_1\alpha_1fK_2K_5K_6K_7
 -\mu\beta \eta_4\psi_1 \gamma_1\alpha_1fK_2K_6K_7
 -\mu\beta \eta_5\psi_2\psi_1\gamma_1\alpha_1fK_2K_7\\
&\quad -\mu\beta \eta_6\alpha_1fK_2K_5K_6\tau_2\gamma_1
  +\mu K_2K_3K_4K_5K_7K_6\sigma_lp
  -\mu\beta \eta6\alpha_1fK_2K_5K_6\tau_1K_4\\
&\quad +K_2K_1K_3K_4K_5K_7K_6\sigma_l+K_6K_7K_5K_4K_3K_2\mu^2\\
&\quad+K_2K_1K_3K_4K_5K_7K_6\sigma_hp-K_2K_1K_3K_4K_5K_7K_6\sigma_lp\\
&\quad -\alpha_1(1-f)\mu K_3K_4K_5K_7K_6\sigma_hp
 -\mu K_2K_3K_4K_5K_7K_6\sigma_hp\\
&\quad  +\mu K_2K_3K_4K_5K_7K_6\sigma_h
 +\alpha_1(1-f)\mu K_3K_4K_5K_7K_6\sigma_lp\\
&\quad +\alpha_1(1-f)\mu K_3K_4K_5K_7K_6\sigma_h
  +\alpha_1f\mu K_2K_4K_5K_7K_6\sigma_lp\\
&\quad +\psi_2\psi_1\gamma_1\alpha_1f\mu K_2K_7\sigma_h
 +\alpha_1(1-f)\mu^2K_3K_4K_5K_7K_6
 +\gamma_1\alpha_1f\mu K_2K_5K_7K_6\sigma_h\\
&\quad -\gamma_1\alpha_1f\mu K_2K_5K_7K_6\sigma_hp
 +\alpha_1f\mu K_2K_4K_5K_7K_6\sigma_h
 -\alpha_1f\mu K_2K_4K_5K_7K_6\sigma_hp\\
&\quad +\alpha_1f\mu^2K_2K_4K_5K_7K_6
 +\alpha_1\phi_4K_2K_7 \psi_2\psi_1\gamma_1f\sigma_h
 -\alpha_1\phi_5 K_2K_5f\tau_2\gamma_1\sigma_hpK_6 \\
&\quad -\psi_2\psi_1\gamma_1\alpha_1f\mu K_2K_7\sigma_hp
 +\psi_2 \psi_1\gamma_1\alpha_1f\mu^2K_2K_7
 +\gamma_1\alpha_1f \mu^2K_2K_5K_7K_6\\
&\quad +\gamma_1\alpha_1f\mu K_2K_5K_7K_6\sigma_lp
 +\psi_1\gamma_1\alpha_1f\mu K_2K_7K_6\sigma_h
 -\psi_1\gamma_1\alpha_1f\mu K_2K_7K_6\sigma_hp\\
&\quad +\psi_1\gamma_1\alpha_1f\mu^2K_2K_7K_6
 +\psi_1\gamma_1\alpha_1f\mu K_2K_7K_6\sigma_lp
 +\alpha_1\phi_2K_2\gamma_1f\mu K_5K_7K_6\\
&\quad +\alpha_1\phi_2K_2\gamma_1f\sigma_hK_5K_7K_6
 +\alpha_1\phi_2K_2\gamma_1f\sigma_lpK_5K_7K_6
 +\alpha_1\phi_5K_2K_5f\tau_2\gamma_1\sigma_hK_6\\
&\quad  +\alpha_1\phi_4K_2K_7\psi_2\psi_1\gamma_1f\mu
 +\alpha_1\phi_4K_2K_7\psi_2\psi_1\gamma_1f\sigma_lp
 +\psi_2\psi_1\gamma_1\alpha_1f\mu K_2K_7\sigma_lp\\
&\quad +\alpha_1f\mu K_2K_5K_6\tau_2\gamma_1\sigma_h
 -\alpha_1f\mu K_2K_5K_6\tau_2\gamma_1\sigma_hp
 +\alpha_1f\mu^2K_2K_5K_6\tau_2\gamma_1\\
&\quad +\alpha_1f\mu K_2K_5K_6\tau_2\gamma_1\sigma_lp
 +\alpha_1f\mu K_2K_5K_6\tau_1\sigma_hK_4
 -\alpha_1f\mu K_2K_5K_6\tau_1\sigma_hpK_4\\
&\quad +\alpha_1f\mu^2K_2K_5K_6\tau_1K_4
 +\alpha_1f\mu K_2K_5K_6\tau_1\sigma_lpK_4
 +\alpha_1\phi_5K_2K_5f\tau_1\sigma_lpK_4K_6\\
&\quad +\alpha_1\phi_5K_2K_5f\tau_2\gamma_1\mu K_6
 +\alpha_1\phi_5K_2K_5f\tau_2\gamma_1\sigma_lpK_6\\
&\quad +\alpha_1\phi_1(1-f)\mu K_3K_4K_5K_7K_6
 -\alpha_1\phi_1(1-f)\sigma_hpK_3K_4K_5K_7K_6\\
&\quad
 -\alpha_1\phi_2K_2\gamma_1f\sigma_hpK_5K_7K_6
 +\alpha_1\phi_1(1-f)\sigma_hK_3K_4K_5K_7K_6\\
&\quad +\alpha_1\phi_3K_2\psi_1\gamma_1f\sigma_hK_7K_6
 -\alpha_1\phi_3K_2\psi_1\gamma_1f\sigma_hpK_7K_6\\
&\quad +\alpha_1\phi_3K_2\psi_1\gamma_1f\mu K_7K_6
  +\alpha_1\phi_5K_2K_5f\tau_1\sigma_hK_4K_6\\
&\quad -\alpha_1\phi_5K_2K_5f\tau_1\sigma_hpK_4K_6
  +\alpha_1\phi_5K_2K_5f\tau_1\mu K_4K_6\\
&\quad +\alpha_1\phi_3K_2\psi_1\gamma_1f\sigma_lpK_7K_6
 +\alpha_1\phi_1(1-f)\sigma_lpK_3K_4K_5K_7K_6\\
&\quad -\alpha_1\phi_4K_2K_7\psi_2\psi_1\gamma_1f\sigma_hp,
\end{align*}
\[
c_0= K_1K_2K_3K_4K_5K_6K_7(\mu+\sigma_l)(\mu+\sigma_h)(1-\mathcal{R}_c ).
\]

The positive endemic equilibria of the model
\eqref{eqn1} are then obtained by solving for $\lambda^{**}$ from the
quadratic \eqref{eq3} and substituting the results (positive values
of $\lambda^{**}$) into \eqref{eq2}. The coefficient $a_0$, of
\eqref{eq3}, is always positive, and $c_0$ is positive (negative)
if $\mathcal{R}_c $ is less than (greater than) unity, respectively. Thus, the
following result is established.

\begin{theorem}\label{thmEEP}
The model \eqref{eqn1} has:
\begin{itemize}
\item[(i)] a unique endemic equilibrium if
 $c_0<0\Leftrightarrow\mathcal{R}_c >1$;
\item[(ii)] a unique endemic equilibrium if $b_0<0$, and
 $c_0=0$ or $b_0^2-4a_0c_0=0$;
\item[(iii)] two endemic equilibrium if $c_0>0$, $b_0<0$ and
 $b_0^2-4a_0c_0>0$;
\item[(iv)] no endemic equilibrium otherwise.
\end{itemize}
\end{theorem}

It is clear from Theorem \ref{thmEEP} (Case (i)) that the model
has a unique endemic equilibrium whenever $\mathcal{R}_c >1$. Case (iii)
of Theorem \ref{thmEEP} suggests the possibility of backward
bifurcation in the model \eqref{eqn1} (where a stable DFE
co-exists with a stable EEP when $\mathcal{R}_c <1$). This phenomenon
is not considered in detail in the current study
(the reader may refer to
\cite{brauer,car,cast3,elbasha,feng,gomez,kribs,chandra,shar,
sharomi1,sharomi2,sharomi3,wang}, and some of the references
therein for discussions on backward bifurcation).


\subsubsection*{Global Stability of EEP: Special Case
($\theta_H=1$ and $\theta_P=0$)}
Consider, again, the model \eqref{eqn1} subject to the special
case where all the susceptible individuals (from both risk groups)
are equally likely to acquire infection (so that, $\theta_H=1$)
and the use of antiviral  prophylaxis gives perfect protection
against infection (i.e., $\theta_P=0$). Further, let
\begin{equation} \label{eqnstab}
  \begin{split}
\operatorname{sign}(S_L-S_L^{**})
&=\operatorname{sign}(S_H-S_H^{*})
=\operatorname{sign}(P-P^{**})\\
&=\operatorname{sign}(L-L^{**})
=\operatorname{sign}(A-A^{**})
=\operatorname{sign}(I_1-I_1^{**})\\
&=\operatorname{sign}(I_2-I_2^{**})
=\operatorname{sign}(H_1-H_1^{**})
=\operatorname{sign}(H_2-H_2^{**})\\
&=\operatorname{sign}(T-T^{**})
=\operatorname{sign}(R-R^{**})\}.
\end{split}
\end{equation}
It is convenient to define the region:
\begin{align*}
\mathcal{D}_0=\{& (S_L,S_H,P,L,A,I_1,I_2,H_1,H_2,T,R)
\in\mathcal{D}:\\
&L=A=I_1=I_2=H_1=H_2=T=R=0  \}.
\end{align*}
We claim the following result.

\begin{theorem}\label{th2}
The associated unique endemic equilibrium of the model \eqref{eqn1},
with $\theta_H=1$ and $\theta_P=0$, is GAS in
$\mathcal{D}\setminus\mathcal{D}_0$  whenever
$\tilde{\mathcal{R}_c} > 1$ and equation \eqref{eqnstab} holds.
\end{theorem}

\begin{proof}
Consider the model \eqref{eqn1} with $\theta_H=1$ and $\theta_P=0$.
Further, let $\tilde{\mathcal{R}_c} > 1$ (so that the
model \eqref{eqn1} has a unique EEP, as guaranteed by
Theorem \ref{thmEEP}). Furthermore, let the relations in
equation \eqref{eqnstab} hold.
Consider the Lyapunov function (Lyapunov functions of this
type have been used in the literature, such as in \cite{yang})
\begin{align*}
\mathcal{F}
&=|S_L-S_L^{**}|+|S_H-S_H^{**}|+|P-P^{**}|+|L-L^{**}|
 +|A-A^{**}|+|I_1-I_1^{**}|\\
&\quad +|I_2-I_2^{**}|+|H_1-H_1^{**}|+|H_2-H_2^{**}|+|T-T^{**}|+|R-R^{**}|.
\end{align*}
It follows that the right derivative, $\mathcal{D^+F}$, of
$\mathcal{F}$, along the solutions of \eqref{eqn1} with
$\theta_H=1$ and $\theta_P=0$, is given by:
\begin{align*}
&D^+\mathcal{F}\\
&=\operatorname{sign}(S_L-S_L^{**})\big\{ (1-p)\Pi-\lambda{S}_L
 -\sigma_LS_L-\mu{S_L}-[(1-p)\Pi -\lambda^{**}{S}_L^{**}\\
&\quad -\sigma_LS_L^{**}-\mu{S_L}^{**}]\big\}
 +\operatorname{sign}(S_H-S_H^{**})\big\{ p\Pi-\lambda{S}_H
 -\sigma_HS_H-\mu{S_H}\\
&\quad -[p\Pi-\lambda^{**}{S}_H^{**}
 -\sigma_HS_H^{**}-\mu{S_H}^{**}]\big\}\\
&\quad +\operatorname{sign}(P-P^{**})[\sigma_LS_L
 +\sigma_HS_H-\mu{P}-(\sigma_LS_L^{**}+\sigma_HS_H^{**}-\mu{P}^{**})]\\
&\quad +\operatorname{sign}(L-L^{**})\big\{[\lambda(S_L+S_H)
 -(\alpha_1+\mu){L}]-\big[\lambda^{**}(S_L^{**}+S_H^{**})\\
&\quad -(\alpha_1+\mu){L}^{**}\big]\big\} \\
&\quad +\operatorname{sign}(A-A^{**})\{\alpha_1(1-f)L
 -(\phi_1+\mu)A-[\alpha_1(1-f)L^{**}-(\phi_1+\mu)A^{**}]\}\\
&\quad +\operatorname{sign}(I_1-I_1^{**})\{\alpha_1f{L}
 -(\tau_1+\gamma_1+\mu)I_{1}-[\alpha_1f{L}^{**}
 -(\tau_1+\gamma_1+\mu)I_{1}^{**}]\}\\
&\quad +\operatorname{sign}(I_2-I_2^{**})\{\gamma_1I_1
 -(\tau_2+\psi_1+\phi_2+\mu+\delta)I_{2}\\
&\quad -[\gamma_1I_1^{**}
 -(\tau_2+\psi_1+\phi_2+\mu+\delta)I_{2}^{**}]\}\\
&\quad +\operatorname{sign}(H_1-H_1^{**})\{\psi_1I_2
 -(\psi_2+\phi_3+\mu+\theta_1\delta)H_1\\
&\quad -[\psi_1I_2^{**}
 -(\psi_2+\phi_3+\mu+\theta_1\delta)H_1^{**}]\}\\
&\quad +\operatorname{sign}(H_2-H_2^{**})\{\psi_2{H}_1
 -(\phi_4+\mu+\theta_2\delta)H_2\\
&\quad -[\psi_2{H}_1^{**}
 -(\phi_4+\mu+\theta_2\delta)H_2^{**}]\}\\
&\quad +\operatorname{sign}(T-T^{**})\{\tau_1{I}_1+\tau_2{I}_2
 -(\phi_5+\mu)T-[\tau_1{I}_1^{**}+\tau_2{I}_2^{**}
 -(\phi_5+\mu)T^{**}]\}\\
&\quad +\operatorname{sign}(R-R^{**})[\phi_1A+\phi_2I_2
 +\phi_3H_1+\phi_4H_2+\phi_5T-\mu{R}\\
&\hskip1in-(\phi_1A^{**}+\phi_2I_2^{**}+\phi_3H_1^{**}
 +\phi_4H_2^{**}+\phi_5T^{**}-\mu{R}^{**})] \\
&=\operatorname{sign}(S_L-S_L^{**})\{(1-p)\Pi
 -\lambda^{**}({S}_L-S_L^{**})-\sigma_L(S_L-S_L^{**})
 -\mu(S_L-S_L^{**})\}\\
&\quad +\operatorname{sign}(S_H-S_H^{**})\{p\Pi
 -\lambda^{**}({S}_H-S_H^{**})-\sigma_H(S_H-S_H^{**})
 -\mu(S_H-S_H^{**})\}\\
&\quad +\operatorname{sign}(P-P^{**})\{\sigma_L(S_L-S_L^{**})
 +\sigma_H(S_H-S_H^{**})-\mu(P-P^{**})\}\\
&\quad +\operatorname{sign}(L-L^{**})\{\lambda^{**}[(S_L-S_L^{**})
 +(S_H-S_H^{**})]-(\alpha_1+\mu)(L-L^{**})\}\\
&\quad +\operatorname{sign}(A-A^{**})\{\alpha_1(1-f)(L-L^{**})
 -(\phi_1+\mu)(A-A^{**})\}\\
&\quad +\operatorname{sign}(I_1-I_1^{**})\{\alpha_1f({L-L^{**}})
 -(\tau_1+\gamma_1+\mu)(I_{1}-I_1^{**})\}\\
&\quad +\operatorname{sign}(I_2-I_2^{**})\{\gamma_1(I_1-I_1^{**})
 -(\tau_2+\psi_1+\phi_2+\mu+\delta)(I_{2}-I_2^{**})\}\\
&\quad +\operatorname{sign}(H_1-H_1^{**})\{\psi_1(I_2-I_2^{**})
 -(\psi_2+\phi_3+\mu+\theta_1\delta)(H_1-H_1^{**})\}\\
&\quad +\operatorname{sign}(H_2-H_2^{**})\{\psi_2({H}_1-H_1^{**})
 -(\phi_4+\mu+\theta_2\delta)(H_2-H_2^{**})\}\\
&\quad +\operatorname{sign}(T-T^{**})\{\tau_1({I}_1-I_1^{**})
 +\tau_2({I}_2-{I}_2^{**})-(\phi_5+\mu)(T-T^{**})\}\\
&\quad +\operatorname{sign}(R-R^{**})\{\phi_1(A-A^{**})
 +\phi_2(I_2-I_2^{**})+\phi_3(H_1-H_1^{**})\\
&\quad +\phi_4(H_2-H_2^{**})+\phi_5(T-T^{**})
 -\mu{(R-R^{**})}\},\\
&=-\mu|S_L-S_L^{**}|-\mu|S_H-S_H^{**}|-\mu|P-P^{**}|
 -\mu|L-L^{**}|-\mu|I_1-I_1^{**}|\\
&\quad -(\mu+\delta)|I_2-I_2^{**}|
 -(\mu+\theta_1\delta)|H_1-H_1^{**}|
  -(\mu+\theta_2\delta)|H_2-H_2^{**}|\\
&\quad -\mu|R-R^{**}| \\
&\le-\mu|S_L-S_L^{**}|-\mu|S_H-S_H^{**}|-\mu|P-P^{**}|
 -\mu|L-L^{**}|-\mu|I_1-I_1^{**}|\\
&\quad -\mu|I_2-I_2^{**}|
 -\mu|H_1-H_1^{**}|-\mu|H_2-H_2^{**}|-\mu|R-R^{**}|\\
&=-\mu \mathcal{F}.
\end{align*}
Thus, $D^+{\mathcal{F}}\leq0$. Hence, $\mathcal{F}$ is a Lyapunov
function on $\mathcal{D}\setminus\mathcal{D}_0$. It follows, by
the LaSalle's Invariance Principle \cite{hale}, that every
solution to the equations of the model \eqref{eqn1}, with
$\theta_H=1$ and $\theta_P=0$, with initial conditions in
$\mathcal{D}\setminus\mathcal{D}_0$ approaches the associated
unique endemic equilibrium of the model \eqref{eqn1}, with
$\theta_P=0$ and $\theta_H=1$, as $t\to \infty$ if
$\tilde{\mathcal{R}}_c >1$ and equation \eqref{eqnstab} holds.
\end{proof}

\begin{figure}[tb]
\begin{center}
\includegraphics[width=0.6\textwidth]{fig3} %ee.pes
\end{center}
\caption{Simulations of \eqref{eqn1} showing the total number of
infected individuals as a function of time, using various initial conditions. 
Parameter values used are as in Table \ref{table1} with  
$p=0.5$, $\theta_H=1.2$, $\theta_P=0.9$, $\alpha_1=0.3$, $f=0.5$,
$\phi_4=\phi_5=1/7$, $c=1$, $\beta=0.4$, $\eta_1=0.1$ 
(so that $\mathcal{R}_c =1.15>1$)}
\label{ee}
\end{figure}


The sensitivity analysis of variables of the model
\eqref{det_model} using the sensitivity functions, $\partial
X/\partial q$, (where $X$ represents a variable, and $q$
represents a parameter), with respect to selective parameters
(given the size of the parameter space) was carried out. The
sensitivity analysis reveals that the sensitivity function of low
risk susceptible individuals, $S_L$, decreases (taking negative
values) with $\beta, \sigma_L,\alpha_1$, and $f$ initially,
attaining a minimum value in around 10 days. It then increases,
and finally becomes insensitive to the respective parameter. On
the other hand, $S_L$ increases with $\tau_1$ first, with a peak
within first 10 days. It then decreases, and finally becomes
insensitive to $\tau_1$. Similarly, the number of infectious
symptomatic individuals at early stage of infection ($I_1$) is
insensitive to $\pi,\eta_i,\beta$, and $\sigma_L$. However $I_1$
increases initially with $\alpha_1$, and $f$, attaining a peak in
less than 10 days, then decreases before becoming insensitive
asymptotically. On the other hand, it first decreases (taking
negative values), then increases with $\tau_1$ and $\gamma_1$, and
then becomes insensitive to the corresponding parameter. It should
be noted that the sensitivity functions of state variables take
the maximum magnitudes during transient time interval. However,
they become insensitive asymptotically. It should be mentioned
that this univariate approach does not account for a possible
influence of correlation between parameter estimates.

Figure \ref{ee} depicts the numerical simulation results obtained
for the case when $\tilde{\mathcal{R}}_c >1$, from which it is clear that all
initial solutions converged to the unique endemic equilibrium
(in line with Theorem \ref{th2}).

\begin{figure}[tb]
\begin{center}
\includegraphics[width=0.48\textwidth]{fig4a} % wave1fit.eps
\includegraphics[width=0.48\textwidth]{fig4b} % wave2fit.eps
\end{center}
\caption{Simulations of  model \eqref{eqn1} showing the time series 
 for the number of infected cases per day and its fit with the actual 
 confirmed new daily cases in Manitoba. 
 (A): First Wave. The first wave was assumed to begin in 
 April and end in early August. Here, we used $\theta_P=0.6$, 
 $\alpha_1=5$ and $\beta=.609$ in the beginning and it is assumed 
 that the combined use of antivirals (prophylaxis and therapeutic) 
 takes effect around the middle of June, (with $\theta_P=0.7$, 
 $\alpha_1=0.4$, $\beta=0.2$).
 Other parameter values are as in Table \ref{table1}. 
 (B): Second Wave. The second wave began in early October and ended 
 at the end of December. Here, $\theta_P=0.6$, $\alpha_1=3$, $\beta=0.6061$ 
 in the beginning of the wave}
\label{wave1wave2fit}
\end{figure}


\section{Numerical Simulations}

The model \eqref{eqn1} is further simulated using the parameter
values in Table \ref{table1} (unless otherwise stated) to quantitatively
assess the role of antivirals in curtailing the spread of the H1N1
pandemic. First of all, the model's output is compared with the
pandemic H1N1 data obtained from the province of Manitoba. The
results obtained, depicted in Figure \ref{wave1wave2fit}, show
that the model fits the observed data (for both the first and
second waves of the pandemic) reasonably well. It should be
mentioned that the model simulations for the first wave (Figure
\ref{wave1wave2fit}A) were based on the assumptions that 30\% of
Manitobans are in the high-risk (of infection) category, and that
the antivirals are available at the beginning of the pandemic. For
the second wave plot (Figure \ref{wave1wave2fit}B), it is
additionally assumed that 20\% of the Manitoban population have
pre-existing immunity against the H1N1 infection (due to their
assumed H1N1 infection during the first wave). We also point out
that insufficient H1N1 data at this point hinders the definition
of realistic ranges for the parameters, and a thorough sensitivity
analysis is not feasible.

The model is further simulated to assess the targeted use
of the antivirals as a preventive agent only (i.e.,
the drugs are only given to susceptible individuals,
and not to infected individuals with disease symptoms),
under the assumptions that all susceptible individuals are
equally likely to acquire the H1N1 infection ($\theta_H=1$)
and that all those who received prophylaxis are completely
protected against infection ($\theta_P=0$). It should be recalled
that, in such a case, Theorem \ref{thm_dfe} guarantees that
the disease will be eliminated if $\tilde{\mathcal{R}}_c \le 1$.

\begin{figure}[tb]
\begin{center}
\includegraphics[width=0.32\textwidth]{fig5a} % ProphylOnlyA.eps
\includegraphics[width=0.32\textwidth]{fig5b}
\includegraphics[width=0.32\textwidth]{fig5c}
\end{center}
 \caption{Simulations of  model \eqref{eqn1} in the absence of
antiviral treatment ($\tau_1=\tau_2=0$) showing the cumulative
number of new cases of infection for different effectiveness
levels of the prevention-only strategy. (A) Low effectiveness
level: $\sigma_L=\sigma_H=0.005$, $\beta=0.9$, $\theta_P=0$, 
$\theta_H=1$, $f=0.5$, $\alpha_1=0.9$ (so that $\mathcal{R}_P=0.0741$); 
(B) Moderate effectiveness level:
$\sigma_L=\sigma_H=0.02$, $\beta=0.9$, $\theta_P=0$, $\theta_H=1$, $f=0.5$,
$\alpha_1=0.9$ (so that $\mathcal{R}_P=0.0186$), 
(C) High effectiveness level:
$\sigma_L=\sigma_H=0.08$, $\beta=0.9$, $\theta_P=0$, $\theta_H=1$, $f=0.5$,
$\alpha_1=0.9$ (so tha, $\mathcal{R}_P=0.0047$). Other parameter
values  are as in Table \ref{table1}} 
\label{ProphylOnly}
\end{figure}

We consider three different levels of effectiveness for this
prevention-only targeted strategy, namely:
\begin{itemize}
\item[(i)] Low effectiveness level of the prevention-only strategy
 ($\sigma_L=\sigma_H=0.005$, $\tau_1=\tau_2=0$, $\theta_H=1$,
 $\theta_P=0$;
 so that $\mathcal{R}_P=\tilde{\mathcal{R}}_c |_{\tau_1=\tau_2=0}=0.0741)$;

\item[(ii)] Moderate effectiveness level of the prevention-only
 strategy ($\sigma_L=\sigma_H=0.02$,
  $\tau_1=\tau_2=0$, $\theta_H=1$, $\theta_P=0$; so that
  $\mathcal{R}_P=\tilde{\mathcal{R}}_c |_{\tau_1=\tau_2=0}=0.0186)$;

\item[(iii)] High effectiveness level of the prevention-only
strategy ($\sigma_L=\sigma_H=0.08$, $\tau_1=\tau_2=0$,
$\theta_H=1$, $\theta_P=0$; so that
$\mathcal{R}_P=\tilde{\mathcal{R}}_c |_{\tau_1=\tau_2=0}=0.0047)$.

\end{itemize}
In other words, the moderate effectiveness level of the
prevention-only strategy is assumed to be four times more
effective than the low effectiveness prevention-only strategy.
Similarly, the high effectiveness prevention-only strategy is
assumed to be four times more effective than the moderate
effectiveness prevention-only strategy. These effectiveness
levels are chosen arbitrarily.
The total cumulative number of new cases of infection is
computed over a span of one year. The results obtained,
depicted in Figure \ref{ProphylOnly}, show a decrease in
the cumulative number of new cases of infection with
increasing effectiveness level of the prevention-only strategy.
While the low effectiveness strategy resulted in close to a
million cumulative cases over one year (Figure \ref{ProphylOnly}A),
the moderate effectiveness level resulted in a decrease to
about 425,000 cases (Figure \ref{ProphylOnly}B).
Furthermore, the high effectiveness level strategy,
which is assumed to be sixteen times more effective than the
low effectiveness strategy, resulted in only about 120 new
cases over the same time period (Figure \ref{ProphylOnly}C).
Thus, these simulations suggest that the singular use
of antivirals as prophylaxis would have limited population-level
impact (in reducing disease burden) except if its effectiveness
level is very high.

Additional simulations are carried out to assess the effect  of
the combined use of the antivirals (both as prophylaxis and
therapeutic agents), under the assumptions that all susceptible
individuals are equally likely to acquire the H1N1 infection and
that all those who received prophylaxis are completely protected
against infection. Here, too, three effectiveness levels (of the
universal strategy) are considered as below:

\begin{itemize}
\item[(i)] Low effectiveness level of the universal strategy
 ($\sigma_L=\sigma_H=0.001$, $\tau_1=\tau_2=0.005$,
 $\theta_H=1$, $\theta_P=0$; so that $\tilde{\mathcal{R}}_c =0.0522)$;

\item[(ii)] Moderate effectiveness level of the universal strategy
($\sigma_L=\sigma_H=0.002$, $\tau_1=\tau_2=0.01$, $\theta_H=1$, $
\theta_P=0$; so that $\tilde{\mathcal{R}}_c =0.0248)$;

\item[(iii)] High effectiveness level of the universal strategy
($\sigma_L=\sigma_H=0.003$, $\tau_1=\tau_2=0.015$,
$\theta_H=1$, $\theta_P=0$; so that $\tilde{\mathcal{R}}_c =0.0157)$.
\end{itemize}


\begin{figure}[tb]
\begin{center}
\includegraphics[width=0.32\textwidth]{fig6a} % CombinA.eps
\includegraphics[width=0.32\textwidth]{fig6b}
\includegraphics[width=0.32\textwidth]{fig6c}
\end{center}
\caption{Simulations of  model \eqref{eqn1} showing the
cumulative number of new cases of infection for different
effectiveness levels of the prevention-treatment combined
strategy. 
(A) Low effectiveness level:
$\sigma_L=\sigma_H=0.001$, $\tau_1=\tau_2=0.005$, $\beta=0.14$, 
$\theta_P=0$, $\theta_H=1$, $f=0.5$, $\alpha_1=0.9$ 
(so that $\tilde{\mathcal{R}}_c =0.0522$); 
(B) Moderate effectiveness level:
$\sigma_L=\sigma_H=0.002$, $\tau_1=\tau_2=0.01$, $\beta=0.14$,
$\theta_P=0$, $\theta_H=1$, $f=0.5$, $\alpha_1=0.9$ (so that
 $\tilde{\mathcal{R}}_c =0.0248$), 
(C) High effectiveness level:
$\sigma_L=\sigma_H=0.003$, $\tau_1=\tau_2=0.015$, $\beta=0.14$,
$\theta_P=0$, $\theta_H=1$, $f=0.5$, $\alpha_1=0.9$ (so that
 $\tilde{\mathcal{R}}_c =0.0157$). Other
parameter values used are as given in Table \ref{table1}} 
\label{Combin}
\end{figure}

\begin{figure}[tb]
\begin{center}
\includegraphics[width=0.32\textwidth]{fig7a} % TimeToEliminateA.eps
\includegraphics[width=0.32\textwidth]{fig7b}
\includegraphics[width=0.32\textwidth]{fig7c}
\end{center}
\caption{Simulations of model \eqref{eqn1} showing the time
needed to eliminate the disease for different effectiveness levels
of the prevention-treatment combined strategy. 
(A) High effectiveness level:
$\sigma_L=\sigma_H=0.003$, $\tau_1=\tau_2=0.015$, $\beta=0.14$,
$\theta_P=0$, $\theta_H=1$, $f=0.2$ (so that $\tilde{\mathcal{R}}_c =0.0087$). 
(B) Moderate effectiveness level:
$\sigma_L=\sigma_H=0.002$, $\tau_1=\tau_2=0.01$, $\beta=0.14$,
$\theta_P=0$, $\theta_H=1$, $f=0.2$ (so that $\tilde{\mathcal{R}}_c =0.0136$), 
(C) Low effectiveness level:
$\sigma_L=\sigma_H=0.001$, $\tau_1=\tau_2=0.005$, $\beta=0.14$, 
$\theta_P=0$, $\theta_H=1$, $f=0.2$ (so that $\tilde{\mathcal{R}}_c =0.0281$); 
Other parameter values used
are as given in Table \ref{table1}} \label{TimeToEliminate}
\end{figure}


The moderate effectiveness level of the universal strategy is
assumed to be twice more effective than the low effectiveness
level. Similarly, the high effectiveness level is assumed to be
three times more effective than the low effectiveness level (it
should be noted that, in all these cases of the universal
strategy, the associated reproduction number $\tilde{\mathcal{R}}_c <1$ (so
that, by Theorem \ref{thm_dfe}, the disease will be eliminated
from the community). The simulations show that while the low
effectiveness level results in about 2800 cases (Figures
\ref{Combin}A), the moderate and high effectiveness levels of the
universal strategy resulted in about 75 and 15 cases, respectively
(Figures \ref{Combin}B, C). Thus, even the moderate effectiveness
level of the universal strategy will be extremely effective in
curtailing the spread of the disease. Figure \ref{TimeToEliminate}
shows the time to disease elimination using the three
effectiveness levels of the universal strategy. As depicted in
Figure \ref{TimeToEliminate}A, the disease can be eliminated in
about 150 days using the high effectiveness level of the universal
strategy. The time to eliminate the disease increases with
decreasing effectiveness level of the universal strategy (Figures
\ref{TimeToEliminate}B, C).

In summary, this study suggests that the use of antivirals as both
prophylaxis and therapeutic agents is more effective than their
targeted use as prophylaxis only. In other words, the prospect of
disease elimination from the community is greatly enhanced if the
universal strategy is used.

\subsection*{Conclusions} A deterministic model is designed and
rigorously analyzed to assess the impact of antiviral drugs in
curtailing the spread of disease of the 2009 swine influenza
pandemic. The analysis of the model, which consists of eleven
mutually-exclusive epidemiological compartments, shows the
following:
\begin{itemize}
\item[(i)] The disease-free equilibrium of the model is shown
to be globally-asymp\-totically stable under the following
conditions:
\begin{itemize}
    \item[(a)] the associated reproduction number $\tilde{\mathcal{R}}_c \le 1$;
     \item[(b)] all susceptible individuals are equally likely to acquire infection (so that, $\theta_H=1$);
     \item[(c)] susceptible individuals who received antivirals are fully protected (so that, $\theta_P=0$).
    \end{itemize}
\item[(ii)] The model has a unique endemic equilibrium when the associated reproduction threshold ($\mathcal{R}_c $) exceeds unity. The unique endemic equilibrium is shown to be globally-asymptotically stable for the special case where all susceptible individuals are equally likely to acquire infection and the use of antiviral prophylaxis gives perfect protection against infection.
\end{itemize}
Numerical simulations of the model, suggest the following:
\begin{itemize}
\item[(a)] The singular use of antivirals, as a preventive agent, has limited population-level impact in reducing disease burden, except if its effectiveness level is very high;
\item[(b)] The combined use of the antivirals, as preventive and therapeutic agents, offers great reduction in disease burden, and will result in disease elimination.
\end{itemize}

\subsection*{Acknowledgments} The authors are very grateful to Dr.
Abba Gumel for his valuable comments, and acknowledge the
assistance of the Analysis, Interpretation, and Research Group of
Health Information Management, Manitoba Health for providing us
with real data on the number of daily confirmed H1N1 cases in
Manitoba.

\begin{thebibliography}{99}

\bibitem{sty1} Bo\"elle, P. Y., Bernillon, P. and Desenclos,
J. C. (2000). A preliminary estimation of the reproduction ratio
for new influenza A(H1N1) from the outbreak in Mexico.
\emph{Euro Surveill.} \textbf{14}(19): pii=19205.

\bibitem{brauer} Brauer, F. (2004).
Backward bifurcations in simple vaccination models.
\emph{J. Math. Anal. Appl.} \textbf{298}(2): 418-431.

\bibitem{blower} Brian, J. C., Bradley, G. W. and Blower, S. (2009).
Modeling influenza epidemics and pandemics: insights into the
future of swine flu (H1N1). \emph{BMC Medicine}
doi:10.1186/1741-7015-7-30.

\bibitem{sty2} Carlos-Chavez, F., Peter, C. and Jose,  P. (2009).
The first influenza pandemic in the new millennium: lessons
learned hitherto for current control efforts and overall
pandemic preparedness.
\emph{Journal of Immune Based Therapies and Vaccines.}
doi:10.1186/1476-8518-7-2.

\bibitem{car} Carr, J. (1981).
Applications Centre Manifold Theory. Springer-Verlag, New York.

\bibitem{cast3} Castillo-Chavez, C. and Baojun, S. (2004).
Dynamical models of tuberculosis and their applications.
\emph{Mathematical Biosci. and Engrg.} \textbf{1}(2): 361-404.

\bibitem{cbc} CBC News - Health - Canada Enters 2nd Wave of H1N1 (2009).
(acc. Nov. 4, 2009)
\url{http://www.cbc.ca/health/story/2009/10/23/h1n1-second-wave-canada.html}.


\bibitem{resist} Centers for Disease Control and Prevention (2009).
Three reports of oseltamivir resistant novel influenza A (H1N1)
 viruses. \url{http://www.cdc.gov/h1n1flu/HAN/070909.htm}.
 (acc. Jan. 23, 2010).

\bibitem{cdc1} Centers for Disease Control and Prevention (2009).
(acc. Oct. 27, 2009).\\
  \url{http://www.cdc.gov/h1n1flu/background.htm}.
 

\bibitem{cdc4} Centers for Disease Control and Prevention (2009).
(acc. Oct. 27, 2009).\\
 \url{http://www.cdc.gov/media/pressrel/2009/r090729b.htm}.
 
 
\bibitem{cdc5} Centers for Disease Control and Prevention
 (CDC) (2009). Outbreak of swine-origin influenza A (H1N1)
 virus infection-Mexico, March-April 2009.
 \emph{MMWR Morb Mortal Wkly Rep}. 58 (Dispatch):1-3.

\bibitem{sty3} Christophe, F., et al (2009).
Pandemic potential of a strain of influenza A (H1N1): Early
findings. \emph{Science.} \textbf{324} (5934): 1557-1561.

\bibitem{vanden2002}  van den Driessche, P  and  Watmough, J. (2002).
Reproduction numbers and sub-threshold endemic equilibria for
compartmental models of disease transmission. \emph{Mathematical
Biosciences.} \textbf{180}: 29-48.

\bibitem{elbasha} Elbasha, E. H. and  Gumel, A. B. (2006).
Theoretical assessment of public health impact of imperfect
prophylactic HIV-1 vaccines with therapeutic benefits.
\emph{Bull. Math. Biol.} \textbf{68}: 577-614.

\bibitem{el} El Universal (2009). (acc. Oct. 27, 2009).\\
\url{http://www.eluniversal.com.mx/hemeroteca/edicion_impresa_20090406.html}.


\bibitem{feng} Feng, Z.,  Castillo-Chavez, C. and  Capurro, F. (2000).
A model for tuberculosis with exogenous reinfection.
\emph{Theor. Pop. Biol.} \textbf{57}: 235-247.

\bibitem{gen} Gen Bank Sequences From 2009 H1N1 Influenza Outbreak (2009).
(acc. Oct. 27, 2009). \url{http://www.ncbi.nlm.nih.gov/genomes/FLU/SwineFlu.html}.
    

\bibitem{marija} Gojovic, M. Z., Sanders, B., MEcDev, R. N., Fisman,
D., Krahn, M.D., and Bauch, C.T. (2009).
Modelling mitigation strategies for pandemic (H1N1).
\emph{CMAJ. DOI:10.1503/cmaj.091641}.

\bibitem{gomez} Gomez-Acevedo, H. and Li, M. (2005).
 Backward bifurcation in a model for HTLV-I infection of CD4+ T cells.
  \emph{Bull. Math. Biol.} \textbf{67}(1): 101-114.

\bibitem{hale} Hale, J. K. (1969).
Ordinary Differential Equations. John Wiley and Sons, New York.

\bibitem{hethcote2000} Hethcote, H. W. (2000).
The mathematics of infectious diseases.
\emph{SIAM Review}.\textbf{ 42}: 599-653.

\bibitem{sty4} Hiroshi, N., Don, K., Mick, R. and Johan,
A. P. H. (2009). Early Epidemiological Assessment of the
Virulence of Emerging Infectious Diseases: A Case Study of
an Influenza Pandemic. \emph{PLoS ONE}.  \textbf{4}(8) :e 6852.

\bibitem{jamieson} Jamieson, D. J., Honein M. A., Rasmussen,
S. A. et al. (2009).  H1N1 2009 influenza virus infection
during pregnancy in the USA.
{\it Lancet}. \textbf{374} (9688): 451-458.

\bibitem{kribs} Kribs-Zaleta, C. and Valesco-Hernandez,
J. (2000). A simple vaccination model with multiple endemic states.
\emph{Math Biosci.} \textbf{164}: 183-201.

\bibitem{kumar} Kumar, A, Zarychanski, R, Pinto, R, et al. (2009).
Critically ill patients with 2009 influenza A(H1N1) infection
in Canada. \emph{JAMA.} \textbf{302}(17): 1872-1879.

\bibitem{lasalle1976} LaSalle, J.P. (1976).
The Stability of Dynamical Systems. Regional Conference Series in
Applied Mathematics, SIAM, Philadelphia.

\bibitem{manitoba2} Manitoba Health:
Confirmed Cases of H1N1 Flu in Manitoba. (acc. Dec. 31, 2009).
\url{http://www.gov.mb.ca/health/publichealth/sri/stats1.html}.


\bibitem{Nuno2007} Nuno, M., Chowell, G.  and  Gumel, A. B. (2007).
Assessing transmission control measures, antivirals and vaccine
in curtailing pandemic influenza: scenarios for the US, UK,
and the Netherlands.
\emph{Proceedings of the Royal Society Interface.}
\textbf{4}(14): 505-521.

\bibitem{chandra} Podder, C. N. and  Gumel, A. B. (2009).
Qualitative dynamics of a vaccination model for HSV-2.
\emph{IMA Journal of Applied Mathematics.} \textbf{302}: 75-107.

\bibitem{sty5} Pourbohloul et al. (2009).
Initial human transmission dynamics of the pandemic (H1N1) 2009
virus in North America.
\emph{Influenza and Other Respiratory Viruses.} \textbf{3}(5): 215-222.

\bibitem{cancases} Public Health Agency of Canada (Week 10).  (2010).
(acc. Mar. 23, 2010). \url{http://www.phac-aspc.gc.ca/fluwatch/09-10}.
 

\bibitem{shar} Sharomi, O., Podder, C. N.,  Gumel, A. B., Mahmud,
S. M. and, E. Rubinstein. Modelling the Transmission Dynamics and
Control of the Novel 2009 Swine Influenza (H1N1) Pandemic.
\emph{Bull. Math. Biol.} To appear.

\bibitem{sharomi1} Sharomi, O and Gumel, A. B. (2009).
Re-infection-induced backward bifurcation in the transmission
dynamics of Chlamydia trachomatis.
\emph{J. Math. Anal. Appl.} \textbf{356}: 96-118.

\bibitem{sharomi2} Sharomi, O., Podder, C.N., Gumel, A. B.,
Elbasha, E. H. and Watmough, J. (2007). Role of incidence function
in vaccine-induced backward bifurcation in some HIV models.
\emph{Mathematical Biosciences.} \textbf{210}: 436-463.

\bibitem{sharomi3} Sharomi, O., Podder, C.N., Gumel, A.B., and
Song, B. (2008). Mathematical analysis of the transmission
dynamics of HIV/TB co-infection in the presence of treatment.
\emph{Math. Biosci. Engrg.} \textbf{5}(1): 145-174.

\bibitem{thieme} Thieme, H. R. (2003). Mathematics in Popluation
Biology. Princeton University Press.

\bibitem{uscdc1} United States Centers for Disease Control
and Prevention (2009). Pregnant women and novel influenza
 A (H1N1): Considerations for clinicians.
(acc. Nov. 5, 2009). \url{http://www.cdc.gov/h1n1flu/clinician_pregnant.htm}.


\bibitem{uscdc2} United States Centers for Disease Control (2009).
Information on people at high-risk of developing flu-related
complications.  \url{http://www.cdc.gov/h1n1flu/highrisk.htm.}
(acc. Nov. 5, 2009).

\bibitem{wang} Wang, W. and Zhao,  X.-Q. (2008).
Threshold dynamics for compartmental epidemic models
in periodic environments.
\emph{J. Dyn. Diff. Equat.} \textbf{20}: 699-717.

\bibitem{wrha} Winnipeg Regional Health Authority Report (2009).
Outbreak of novel H1N1 influenza A virus in the Winnipeg Health Region.
 \url{http://www.wrha.mb.ca/}.  (acc. Nov. 4, 2009).

\bibitem{who1} World Health Organization (2009).
Pandemic (H1N1) (2009)-update 71. (acc. Oct. 27, 2009).
\url{http://www.who.int/csr/don/2009_10_23/en/index.html}.


\bibitem{who2} World Health Organization (2009).
Influenza A (H1N1)-update 49. Global Alert and Response (GAR).
\url{http://www.who.int/csr/don/2009_06_15/en/index.html}.
(acc. Oct. 27, 2009).

\bibitem{who3} World Health Organization (2009).
Statement by Director-General. June 11, 2009.

\bibitem{whosalah1} World Health Organization (2009).
Human infection with new influenza A (H1N1) virus:
clinical observations from Mexico and other affected countries.
Weekly epidemiological record, May 2009; 84:185.
\url{http://www.who.int/wer/2009/wer8421.pdf}.
(acc. Nov. 5, 2009).

\bibitem{who9} World Health Organization (2009).
Pandemic (H1N1) 2009 - update 81 (acc. Mar. 5, 2010).
\url{http://www.who.int/csr/don/2010_03_05/en/index.html}.


\bibitem{yang} Yang, Y. and Xiao, Y. (2010).
Threshold dynamics for an HIV model in periodic environment.
\emph{JMAA}. \textbf{361}(1): 59-68

\end{thebibliography}

\end{document}