% Will be submitted to JMAA
% first draft @ 14 October 2008
% second draft @ 20 October 2008
% Third draft @ 4 November 2008
% Forth draft @ 14 November 2008
% Fifth draft @ 24 November 2008

\documentclass[11pt]{article}

\usepackage{amsfonts, amsmath, amssymb, amsthm}
\usepackage{newlfont, graphicx}
%\usepackage{colordvi}
%\usepackage[dvips,usenames]{color}

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
   % set the property of page
\setlength{\textwidth}{6in} \setlength{\textheight}{8.5in}
\setlength{\oddsidemargin}{.29in}
\renewcommand{\baselinestretch}{1.5}
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

%\makeatletter \@addtoreset{equation}{section}
%\@addtoreset{theo}{section} \makeatother

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
   % set special command
\newcommand{\qqed}{\quad \rule[-.5mm]{1.9mm}{3.2mm}}

   % set section and return equation label
%\newcommand{\newsection}[1]{\section{#1} \setcounter{equation}{0}}

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
   % set format of theorem, lemma etc.
\theoremstyle{plain}
\newtheorem{thm}{Theorem}[section]
\newtheorem{prop}[thm]{Proposition}
\newtheorem{lem}[thm]{Lemma}
\newtheorem{cor}[thm]{Corollary}
\newtheorem{conj}[thm]{Conjecture}
\newtheorem{remark}[thm]{Remark}

  % set format of definition, example etc.
\theoremstyle{definition}
\newtheorem{defn}{Definition}[section]
\newtheorem{ex}{Example}[section]


\newtheorem{involution}{Involution}[section]
\renewcommand{\theequation}{\arabic{section}.\arabic{equation}}
\renewcommand{\thesection}{\arabic{section}}
\renewcommand{\thethm}{\arabic{section}.\arabic{thm}}
\renewcommand{\thelem}{\arabic{section}.\arabic{lem}}
\renewcommand{\thefootnote}{\fnsymbol{footnote}}



%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
    % My definition
\newcommand{\f}{\frac}
\newcommand{\p}{\partial}
\newcommand{\disp}{\displaystyle}
\newcommand{\Om}{\Omega}
\newcommand{\om}{\omega}
\newcommand{\F}{\mathcal{F}}
\newcommand{\PS}{(\Om, \F, P)}
\newcommand{\R}{\mathbb{R}}
\newcommand{\N}{\mathbb{N}}
\newcommand{\de}{\triangle}
\newcommand{\Xn}{\{X_n\}}
\newcommand{\ep}{\epsilon}
\newcommand{\lt}{\mathcal{L}_T^2}
\newcommand{\mseq}{\overset{\textrm{m.s.}}{=}}
\newcommand{\msto}{\overset{\textrm{m.s.}}{\longrightarrow}}
\newcommand{\Peq}{\overset{\textrm{P}}{=}}
\newcommand{\Pto}{\overset{\textrm{P}}{\longrightarrow}}
\newcommand{\aseq}{\overset{\textrm{a.s.}}{=}}
\newcommand{\asto}{\overset{\textrm{a.s.}}{\longrightarrow}}
\newcommand{\Deq}{\overset{\textrm{D}}{=}}
\newcommand{\Dto}{\overset{\textrm{D}}{\longrightarrow}}
\newcommand{\del}{\delta}
\newcommand{\gam}{\gamma}
\newcommand{\eps}{\epsilon}
\newcommand{\Lam}{\Lambda}
\newcommand{\sut}{\subset}
\newcommand{\mto}{\mapsto}
\newcommand{\OGam}{\overset{\circ}{\Gamma}}
\newcommand{\Vdot}{\overset{\bullet}{V}}



\newcommand{\ba}{\begin{array}}
\newcommand{\ea}{\end{array}}


\newcommand{\tr}{\textrm{tr}}

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%


\title{Global Stability of Multi-Group 
Epidemic Models with Distributed Delays}

\author{Michael Y. Li, \quad\quad Zhisheng Shuai\\[8pt]
\small Department of Mathematical and Statistical Sciences,
University of Alberta\\
\small Edmonton, Alberta, T6G 2G1,
Canada\\[3pt]
\small \texttt{mli@math.ualberta.ca \quad zshuai@math.ualberta.ca}\\[5pt]
and\\[5pt]
Chuncheng Wang\\[5pt]
\small Department of Mathematics, Harbin Institute of Technology\\
\small Harbin, Heilongjiang, 150001, China\\[3pt]
\small \texttt{wangchuncheng@hit.edu.cn}}

\date{}

\begin{document}
\maketitle

\begin{abstract}
We investigate a class of multi-group SIR epidemic models with
arbitrarily distributed disease latency. The models are described by
systems of functional differential equations with infinitely
distributed delays.  We establish that the global dynamics are
completely determined by the basic reproduction number
$\mathcal{R}_0$. More specifically, we prove that, if $\mathcal{R}_0
\leq 1$, then the disease-free equilibrium is globally
asymptotically stable; if $\mathcal{R}_0
> 1$, then there exists a unique endemic equilibrium and it is
globally asymptotically stable. Our proof of global stability of the
endemic equilibrium utilizes a graph-theoretical approach to the
method of Lyapunov functionals.
\end{abstract}
\vskip 0.2cm

\noindent {\bf Keywords:} Multi-group SIR model, arbitrarily
distributed latency, infinite delay, global stability, Lyapunov
functional. \vskip 0.2cm

\noindent {\bf 2000 Mathematics Subject Classification}: 34K20,
92D30
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

\section{Introduction}\label{section-intro}\setcounter{equation}{0}

In the literature of mathematical epidemiology, multi-group
epidemic models have been proposed to describe the spread of many
infectious diseases in heterogeneous populations, such as measles,
mumps, gonorrhea, and HIV/AIDS. A heterogeneous host population
can be divided into several homogeneous groups according to modes
of transmission, contact patterns, or geographic distributions, so
that within-group and inter-group interactions could be modeled
separately. One of the earliest multi-group models is proposed by
Lajmanovich and Yorke \cite{ly} for the transmission of gonorrhea.
For a class of $n$-group SIS models, they have completely
established the global dynamics and proved the global stability of
a unique endemic equilibrium using a quadratic global Lyapunov
function. Various forms of multi-group models have subsequently
been studied. One of main mathematical challenges in the analysis
of multi-group models is the global stability of the endemic
equilibrium \cite{bc,he,ht,hcc,ls,th1,th2}. A complete resolution
of this problem has been elusive until recently. In
\cite{gls1,gls2}, for a class of multi-group SEIR models described
by ordinary differential equations, a graph-theoretic approach to
the method of global Lyapunov functions was proposed and used to
establish the global stability of a unique endemic equilibrium.

In the present paper, a multi-group epidemic model with more general
disease latency than that in \cite{gls2} is proposed to describe the
disease spread in a heterogeneous population. The total population
is divided into several homogeneous groups. Let $S_k, E_k, I_k$ and
$R_k$ denote the susceptible, infected but non-infectious,
infectious, and recovered populations in the $k$-th group,
respectively. Let $i_k(t,r)$ denote the population of infected
individuals with respect to the age of infection $r$ at the current
time $t$ and $I_k(t)=\int_{r=0}^\infty i_k(t,r)dr$. Let
$h_k(r)\geq 0$ be a continuous kernel function that represents
the infectivity at the age of infection $r$. The disease incidence
in the $k$-th group, assuming a bilinear incidence form, can be calculated as
\begin{equation}\label{incidence}\sum_{j=1}^n
\beta_{kj}S_k(t)\int_{r=0}^\infty h_j(r)i_j(t,r)dr,\end{equation}
where the sum takes into account of cross-infections from all groups. In
the special case $h_k(r)\equiv 1$, the incidence in
\eqref{incidence} becomes $\sum_{j=1}^n \beta_{kj}S_k(t)I_j(t)$.
Therefore, the model in \cite{gls2} can be generalized to the
following system of differential equations
\begin{equation}\label{model1}
\begin{aligned}
S_k^{\prime}&=\displaystyle\Lambda_k-
\sum_{j=1}^n\beta_{kj}S_k\int_{r=0}^{\infty}h_j(r)i_j(t, r)dr-d_k^SS_k,\\
E_k^{\prime}&=\displaystyle\sum_{j=1}^n\beta_{kj}S_k
\int_{r=0}^{\infty}h_j(r)i_j(t,r)dr - (d_k^E+\epsilon_k)E_k,\\
I_k^{\prime} &= \epsilon_k E_k -(d_k^I+\gamma_k)I_k, \vspace{0.4cm}\\
R_k^{\prime}&= \gamma_k I_k - d_k^R R_k,\quad \quad k=1,2,\cdots,n.
\end{aligned} 
\end{equation}
Here $\Lambda_k$ represents influx of individuals into the $k$-th
group, $\beta_{kj}$ represents the transmission coefficient between
compartments $S_k$ and $I_j$, $d_k^S, d_k^E, d_k^I$ and $d_k^R$
represent death rates of $S, E, I$ and $R$ compartments in the
$k$-th group, respectively, $\epsilon_k$ represents the rate of
becoming infectious after latent period in the $k$-th group, and
$\gamma_k$ represents the recovery rate of infectious individuals in
the $k$-th group. All parameter values are assumed to be nonnegative
and $\Lambda_k, d_k^S, d_k^I>0$ for all $k$. For detailed
discussions of the model, we refer the reader to \cite{gls2,rw} and
references therein.
%In the special case when $h_k(r)$ is an
%exponential function, model \eqref{model1} reduces to the SEIR model
%in \cite{gls2} described by a system of ordinary differential
%equations.
Note that
$$
\begin{array}{rcl} \displaystyle
\Big(\frac{\partial}{\partial t}+\frac{\partial}{\partial r}
\Big)i_k(t,r) &=&
-(d_k^I + \gamma_k) i_k(t,r),\\
i_k(t,0)&=&\epsilon_k E_k(t),
\end{array}
$$ 
whose solution is 
\begin{equation}\label{i}
i_k(t,r) = i_k(t-r,0)e^{-(d_k^I+\gamma_k)r} = \epsilon_k
E_k(t-r)e^{-(d_k^I+\gamma_k)r}.
\end{equation}
Substituting \eqref{i} into \eqref{model1} we obtain
\begin{equation}\label{model2}
\begin{aligned}
S_k^{\prime}&=\displaystyle\Lambda_k-
\sum_{j=1}^n\beta_{kj}S_k\int_{r=0}^{\infty}h_j(r)\epsilon_j E_j(t-r)e^{-(d_j^I+\gamma_j)r} dr-d_k^SS_k,\\
E_k^{\prime}&=\displaystyle\sum_{j=1}^n\beta_{kj}S_k\int_{r=0}^{\infty}h_j(r)\epsilon_j E_j(t-r)e^{-(d_j^I+\gamma_j)r} dr - (d_k^E+\epsilon_k)E_k,\\
I_k^{\prime} &= \epsilon_k E_k -(d_k^I+\gamma_k)I_k, \vspace{0.5cm}\\
R_k^{\prime}&= \gamma_k I_k - d_k^R R_k, \quad \quad k=1,2,\cdots,n.
\end{aligned} 
\end{equation}
Since the variables $I_k$ and $R_k$ do not appear in the first two
equations of \eqref{model2}, we can consider the
following reduced system with infinitely distributed time delays and
general kernel functions
\begin{equation}\label{model}
\begin{aligned}
S_k^{\prime}&=&\displaystyle\Lambda_k-
\sum_{j=1}^n\beta_{kj}S_k\int_{r=0}^{\infty}f_j(r)E_j(t-r) dr-d_k^SS_k,\\
E_k^{\prime}&=&\displaystyle\sum_{j=1}^n\beta_{kj}S_k\int_{r=0}^{\infty}f_j(r)
E_j(t-r) dr - (d_k^E+\epsilon_k)E_k, 
\end{aligned} \quad \quad k=1,2,\cdots,n.
\end{equation}
Here the kernel function $f_k(r)\geq 0$ is continuous and
$\int_{r=0}^\infty f_k(r) dr = a_k >0$. While system is derived from
an general age-structured model (\ref{model1}), it can also be interpreted   
as a multi-group model for an infectious disease whose latent period $r$ in hosts 
has a general probability density function $\frac 1{a_k}f_k(r)$, for the $k$-th group.
We will establish the global dynamics of system (\ref{model}).

The basic reproduction number
$\mathcal{R}_0$ is defined as the expected number of secondary
cases produced in an entirely susceptible population by a typical
infected individual during its entire infectious period
\cite{dhm}. Intuitively, if $\mathcal{R}_0 < 1$, the
disease dies out from the host population, and if $\mathcal{R}_0
>1$, the disease will persist. Let $S_k^0=\frac{\Lambda_k}{d_k^S}, \,
a_j = \int_{r=0}^\infty f_j(r) dr$. The next generation
matrix for system \eqref{model} is
\begin{equation}\label{M0}
M_0=\Big( \frac{\beta_{kj}S_k^0 a_j}{ d_k^E+\epsilon_k}
\Big)_{n\times n}.
\end{equation}
Motivated by \cite{dhm, vw, wz}, we define
the basic reproduction number as the spectral radius of $M_0$,
\begin{equation}\label{R0}
\mathcal{R}_0= \rho (M_0).
\end{equation}  
In the
special case of a single-group model, the definition of
$\mathcal{R}_0$ in \eqref{R0} agrees with that in \cite{rw}. In the
special case when $f_k(r)$ is an exponential function,
$\mathcal{R}_0$ reduces to that for the resulting ODE models as
given in \cite{gls2,vw}. For system \eqref{model}, we will establish
that the dynamical behaviors are completely determined by values of
$\mathcal{R}_0$. More specifically, if $\mathcal{R}_0\leq 1$, the
disease-free equilibrium is globally asymptotically stable and the
disease dies out; if $\mathcal{R}_0>1$, a unique endemic equilibrium
exists and is globally asymptotically stable, and the disease
persists at the endemic equilibrium. Our proof demonstrates that the
graph-theoretical approach developed in \cite{gls1,gls2} for systems
of ordinary differential equations is applicable to delay
differential systems like \eqref{model}.

The paper is organized as follows. In the next section, we prove
some preliminary results for system \eqref{model}. Our main results
are stated in Section~\ref{section-Main}. In
Section~\ref{section-DFE}, the global stability of the disease-free
equilibrium is proved. The global stability of the endemic
equilibrium  is proved in Section~\ref{section-EE}. For the
convenience of the reader, we include in Appendix results from graph
theory that are needed for our proof.









%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%


\section{Preliminaries}
\label{section-model}\setcounter{equation}{0}

We make the following assumption on the kernel function $f_k(r)$
in \eqref{incidence}
\begin{equation}\label{kernel}
\int_{r=0}^\infty f_k(r) e^{\lambda_k r} dr < \infty,
\end{equation} where $\lambda_k$ is a positive number, $
k=1,2,\ldots, n$. Define the following Banach space of fading memory
type (see e.g. \cite{ah} and references therein)
\begin{equation}\label{Ck}
\begin{array}{l}
C_{k} = \displaystyle\Big\{ \phi \in C((-\infty, 0], \mathbb{R}) :
\;  \phi(s)e^{\lambda_k s}
\;\mbox{is}\;\mbox{uniformly}\;\mbox{continuous}\;\mbox{on}\;
(-\infty,0],\\
\hspace{9cm}\displaystyle \mbox{and} \; \sup_{s\leq 0}
|\phi(s)|e^{\lambda_k s} < \infty  \Big\},
\end{array}\end{equation} with norm $||\phi||_k = \sup_{s\leq 0}
|\phi(s)| e^{\lambda_k s}$. For $\phi \in C_k$, let $\phi_{t}\in
C_{k}$ be such that $\phi_{t}(s)=\phi(t+s)$, $s \in (-\infty, 0]$.
Let $S_{k,0}\in \mathbb{R}_+$ and $\phi_k \in C_{k}$ such that
$\phi_k(s)\geq 0, s\in(-\infty,0]$. We consider solutions of system
\eqref{model}, $(S_1(t), E_{1\,t}, S_2(t),
E_{2\,t},\cdots,S_n(t),E_{n\,t})$, with initial conditions
\begin{equation}\label{initial}
S_k(0)=S_{k,0},  \quad E_{k\,0}=\phi_k, \quad k=1,2,\ldots, n.
\end{equation}
Standard theory of functional differential equations \cite{hmn}
implies $E_{k \, t} \in C_k$ for $t > 0$. We consider system
\eqref{model} in the phase space
\begin{equation}\label{X}
X=\prod_{k=1}^n \Big( \mathbb{R} \times C_{k} \Big).
\end{equation}
It can be verified that solutions of \eqref{model} in $X$ with
initial conditions \eqref{initial} remain nonnegative. In
particular, $S_k(t)>0$ for $t>0$. From the first equation of
\eqref{model}, we obtain $S_k^\prime(t)\leq\Lambda_k-d_k^SS_k(t)$.
Hence, $\limsup_{t\to\infty}S_k(t)\leq\frac{\Lambda_k}{d_k^S}$. For
each $k$, adding the two equations in \eqref{model} gives
$(S_k(t)+E_{k\,t}(0))^\prime\leq\Lambda_k-d_k^*(S_k(t)+E_{k\,t}(0))$,
where $d_k^*=\min\{d_k^S, d_k^E+\epsilon_k\}$. Hence,
$\limsup_{t\to\infty}(S_k(t)+E_{k\,t}(0))\leq\frac{\Lambda_k}{d_k^*}$.
Therefore, the following set is positively invariant for system
\eqref{model}.
\begin{equation}\label{theta}
\begin{aligned}
\Theta = \Big\{ (S_1, E_1(\cdot), \cdots, S_n,
E_n(\cdot)) \in X \;\Big|\; 0\leq &S_k \leq
\frac{\Lambda_k}{d_k^S},\; 0 \leq
S_k+E_k(0)\leq \frac{\Lambda_k}{d_k^*}, \\
& E_k(s) \geq 0, s\in (-\infty,0],
k=1,\ldots,n \Big\}.
\end{aligned}
\end{equation} 
All positive
semi-orbits in $\Theta$ are precompact in $X$ \cite{ah}, and thus
have non-empty $\omega$-limit sets. We have the following result.
\begin{lem}\label{lemma} All positive semi-orbits in $\Theta$ have
non-empty $\omega$-limit sets.
\end{lem}
\noindent Let
\begin{equation}\label{theta0}
\begin{aligned}
\overset{\circ}{\Theta} = \{ (S_1, E_1(\cdot),
\cdots, S_n, E_n(\cdot)) \in X \;\Big|\; 0< &S_k <
\frac{\Lambda_k}{d_k^S},\; 0 <
S_k+E_k(0) < \frac{\Lambda_k}{d_k^*}, \\
&E_k(s) > 0, s\in(-\infty,0], k=1,\ldots,n
\Big\}.
\end{aligned}
\end{equation} It can be shown that
$\overset{\circ}{\Theta}$ is the interior of $\Theta$.


The equilibria of \eqref{model} are the same as those of the
associated ODE system
\begin{equation}\label{SIR}
\begin{array}{rcl}
S_k^{\prime}&=&\displaystyle\Lambda_k-
\sum_{j=1}^n \beta_{kj}a_jS_k E_j-d_k^SS_k,\\
E_k^{\prime} &=& \displaystyle\sum_{j=1}^n\beta_{kj} a_j S_k E_j
-(d_k^E+\epsilon_k)E_k, \quad \quad
k=1,2,\ldots,n.\\
\end{array}
\end{equation} A similar system to \eqref{SIR} has been analyzed in
\cite{gls1}. We recall some results for system \eqref{SIR} which
have been proved in \cite{gls1}. In the positively invariant region
\begin{equation}\label{gamma}
\Gamma =\Big\{ (S_1, E_1, \cdots, S_n, E_n)\in
\mathbb{R}_+^{2n}\quad | \quad S_k\le \frac {\Lam_k}{d_k^S},\
S_k+E_k \leq \f{\Lambda_k}{d_k^*},\, 1\le k\le n\Big\},
\end{equation} system \eqref{SIR} has two possible equilibria: the
disease-free equilibrium $P_0=(S_1^0,0,\cdots,S_n^0,0)$, where
$S_k^0=\frac{\Lambda_k}{d_k^S}$, and the endemic equilibrium
$P^*=(S_1^*,E_1^*,\cdots,S_n^*,E_n^*)$, where $S_k^*,E_k^*>0$ and
satisfy
\begin{equation}\label{eq1}
\Lambda_k=\sum\limits_{j=1}^n\beta_{kj}a_jS_k^*E_j^*+d_k^SS_k^*,
\end{equation}
\begin{equation}\label{eq2}
\sum\limits_{j=1}^n\beta_{kj}a_j S_k^*E_j^*=(d_k^E+\epsilon_k)E_k^*.
\end{equation}
We assume that the transmission matrix $B=(\beta_{kj})$ is
irreducible. This is equivalent to assuming that for any two
distinct groups $k$ and $j$, individuals in $E_j$ can infect those
in $S_k$ directly or indirectly. The following result is proved in
\cite{gls1}.
\begin{prop}[Guo, Li,
Shuai \cite{gls1}]\label{GLS} Assume that $B=(\beta_{kj})$ is
irreducible.
\begin{enumerate}
    \item[$(1)$] If $\mathcal{R}_0\leq 1$,
    then $P_0$ is the only equilibrium for system \eqref{SIR} and it is
    globally asymptotically stable in $\Gamma$.
    \item[$(2)$] If $\mathcal{R}_0 > 1$, then $P_0$ is unstable and there exists a
    unique endemic equilibrium $P^*$ for
    system \eqref{SIR}. Furthermore, $P^*$ is globally asymptotically stable
    in the interior of  $\Gamma$.
\end{enumerate}
\end{prop}
Since the delay system \eqref{model} and the ODE system
\eqref{SIR} share the same equilibria, the following result
follows from Proposition~\ref{GLS}.
\begin{prop}\label{equilibria}
Assume that $B=(\beta_{kj})$ is irreducible.
\begin{enumerate}
    \item[$(1)$] If $\mathcal{R}_0\leq 1$, then $P_0$ is the only equilibrium
    for system \eqref{model} in $\Theta$.
    \item[$(2)$] If $\mathcal{R}_0>1$, then there exist two equilibria for
    system \eqref{model} in $\Theta$: the disease-free equilibrium $P_0$
    and a unique endemic equilibrium $P^*$ defined by equations \eqref{eq1}
    and \eqref{eq2}.
\end{enumerate}
\end{prop}



%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

\section{Main Result}
\label{section-Main}\setcounter{equation}{0}

The global dynamics of system \eqref{model} are completely
established in the following result.
\begin{thm}\label{main}
Assume that $B=(\beta_{kj})$ is irreducible.
\begin{enumerate}
\item[$(1)$] If $\mathcal{R}_0\leq1$, then the disease-free equilibrium
$P_0$ of system \eqref{model} is globally asymptotically stable in
$\Theta$. If $\mathcal{R}_0 > 1$, then $P_0$ is unstable.
\item[$(2)$] If $\mathcal{R}_0>1$, then the endemic equilibrium $P^*$ of system
\eqref{model} is globally asymptotically stable in
$\overset{\circ}{\Theta}$. \end{enumerate}\end{thm}

Biologically, Theorem~\ref{main} implies that, if the basic
reproduction number $\mathcal{R}_0\leq 1$, then the disease always
dies out from all groups; if $\mathcal{R}_0
> 1$, then the disease always persists in all groups at the unique endemic
equilibrium level, irrespective of the initial conditions. The
proof of the first part of Theorem~\ref{main} will be given in the
next section, and the second part in Section~\ref{section-EE}.

\vskip 10pt 
Theorem~\ref{main} includes several previous results.
Choose the kernel function as
$$
f_k(r)=\epsilon_k e^{-(d_k^I+\gamma_k) r}
$$
and let $\tilde{I}_k=\int_{r=0}^\infty f_k(r) E_k(t-r)dr$. Then 
$$ 
S_k^{\prime}= \displaystyle\Lambda_k-
\sum_{j=1}^n\beta_{kj}S_k\tilde{I}_j-d_k^SS_k,$$
$$E_k^\prime = \sum_{j=1}^n\beta_{kj}S_k
\tilde{I}_j-(d_k^E+\epsilon_k)E_k.
$$ 
Using integration by parts we obtain
$$
\tilde{I}_k^\prime =\int_{r=0}^\infty f_k(r) \frac{\partial E_k(t-r)}{\partial t}dr
= -\int_{r=0}^\infty f_k(r) \frac{\partial E_k(t-r)}{\partial r}dr
=\epsilon_k E_k - (d_k^I+\gamma_k)\tilde{I}_k.
$$ 
System \eqref{model} is thus reduced to a multi-group SEIR model governed by
the system of ordinary differential equations considered in
\cite{gls2}. Note that 
$$
\int_{r=0}^\infty f_k(r) dr =
 \frac{\epsilon_k}{d_k^I+\gamma_k},
$$ 
the basic reproduction number in (\ref{R0}) becomes
$$
\mathcal{R}_0= \rho \Big( \frac{\beta_{kj}\epsilon_k\Lambda_k}{(
d_k^E+\epsilon_k)(d_k^I+\gamma_k)d_k^S}\Big)_{1\le k, j\le n},
$$
which agrees with that given in \cite{gls1,gls2}. Thus the global stability
results in \cite{gls1,gls2} are special cases of Theorem~\ref{main}.

In the case $n=1$, system \eqref{model} reduces to a single-group
SEIR or SIR model with infinitely distributed delays studied in \cite{bcr,bhmt,bt,mthb,m1,m2,rw}.
Theorem~\ref{main} generalizes the global stability results in
\cite{m1,m2} to multi-group models.



%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%


\section{Proof of Theorem~\ref{main}-(1)}\label{section-DFE}\setcounter{equation}{0}

Since $B$ is irreducible, we know that $M_0$, as defined in
\eqref{M0}, is also irreducible, and has a positive left
eigenvector $(\omega_1,\omega_2,\cdots,\omega_n)$ corresponding to
the spectral radius $\rho(M_0)>0$. Let
$$
a_k=\frac{\omega_k}{d_k^I+\gamma_k} \quad \mbox{and} \quad
\alpha_k(r)=\int_{\sigma=r}^{\infty}f_k(\sigma)d\sigma.
$$
Consider a
Lyapunov functional
\begin{equation}\label{L}
L=\sum\limits_{k=1}^n a_k\Big(S_k-S_k^0-S_k^0\ln\frac{S_k}{S_k^0} +
E_k+
\sum\limits_{j=1}^n\beta_{kj}S_k^0\int_{r=0}^\infty\alpha_j(r)E_j(t-r)dr\Big).
\end{equation}
Note that $\Lambda_k=d_k^S S_k^0$,
$\alpha_k(0)=\int_{\sigma=0}^\infty f_k(\sigma)d\sigma=a_k$, and
$\frac{S_k}{S_k^0}+\frac{S_k^0}{S_k}\geq 2$ with equality holding if
and only if $S_k=S_k^0$. Differentiating $L$ along the solution of
system \eqref{model} and using integration by parts, we obtain
\begin{equation}\label{L'}
\begin{array}{rcl}
L^\prime &=&\displaystyle\sum_{k=1}^n a_k\Big(\Lambda_k -d_k^S S_k
-\Lambda_k \frac{S_k^0}{S_k}+\sum_{j=1}^n\beta_{kj}S_k^0
\int_{r=0}^{\infty}f_j(r)E_j(t-r)\,dr +d_k^S S_k^0 \\
&& \quad \displaystyle -(d_k^E+\epsilon_k)E_k
+\sum_{j=1}^n\beta_{kj}S_k^0\int_{r=0}^\infty\alpha_j(r)\frac{\partial
E_j(t-r)}{\partial t}\;dr\Big)\\
&=&\displaystyle\sum_{k=1}^n a_k\Big[d_k^SS_k^0 \Big(2-
\frac{S_k}{S_k^0}-\frac{S_k^0}{S_k}\Big) +
\sum_{j=1}^n\beta_{kj}S_k^0
\int_{r=0}^{\infty}f_j(r)E_j(t-r)\,dr\\
&&\quad\displaystyle - (d_k^E+\epsilon_k)E_k
+\sum_{j=1}^n\beta_{kj}S_k^0\int_{r=0}^\infty\alpha_j(r)\Big( -
\frac{\partial
E_j(t-r)}{\partial r} \Big) \;dr\Big]\\
%%%%
&=&\displaystyle\sum_{k=1}^n a_k\Big[d_k^SS_k^0 \Big(2-
\frac{S_k}{S_k^0}-\frac{S_k^0}{S_k}\Big) +
\sum_{j=1}^n\beta_{kj}S_k^0
\int_{r=0}^{\infty}f_j(r)E_j(t-r)\,dr\\
&&\quad\displaystyle - (d_k^E+\epsilon_k)E_k
+\sum_{j=1}^n\beta_{kj}S_k^0 \Big(a_j \,E_j-\int_{r=0}^\infty f_j(r) E_j(t-r)\;dr \Big)\Big]\\
%%%%
&\leq &\displaystyle\sum\limits_{k=1}^n
\frac{\omega_k}{d_k^I+\gamma_k}\Big(\sum_{j=1}^n\beta_{kj}a_jS_k^0E_j-
(d_k^E+\epsilon_k)E_k\Big)\;= \;(\omega_1,\omega_2,\cdots,\omega_n)(M_0 E-E)\\
&=&(\rho(M_0)-1)(\omega_1,\omega_2,\cdots,\omega_n)E \,\leq\, 0,
\quad\quad\quad \mbox{if}\;\mathcal{R}_0\leq 1.
\end{array}
\end{equation}
Here $E(t)=(E_1(t),E_2(t),\cdots, E_n(t))^T$.  Denote
$$Y=\{(S_1, E_1(\cdot), \cdots, S_n, E_n(\cdot)) \in \Theta
\quad|\quad L^\prime=0\},$$ and $Z$ be the largest compact
invariant set in $Y$. We will show $Z=\{P_0\}$. From \eqref{L'},
$L^\prime=0$ implies $S_k(t)\equiv S_k^0=\frac{\Lambda_k}{d_k^S}$
for each $k$. Hence, from the first equation of \eqref{model}, we
obtain
$$\sum_{j=1}^n\beta_{kj} \int_{r=0}^{\infty}f_j(r)E_j(t-r)dr=0,$$
and thus
$$\beta_{kj} \int_{r=0}^{\infty}f_j(r)E_j(t-r)dr=0,$$
for $1 \leq k, j \leq n$. Then, by irreducibility of $B$, for each
$j$, there exists $k\not = j$ such that $\beta_{kj}\not=0$, thus
$$\int_{r=0}^{\infty}f_j(r)E_j(t-r)dr=0.$$
This implies that in $Z$, $E_{j\,t}(s)= 0, s\in(-\infty,0]$,
$j=1,2,\ldots,n$. Therefore, $Z=\{P_0\}$.

Using Lemma~\ref{lemma} and the LaSalle-Lyapunov Theorem (see
\cite[Theorem 3.4.7]{la} or \cite[Theorem 5.3.1]{hv}), we conclude
that $P_0$ is globally attractive in $\Theta$ if $\mathcal{R}_0\leq
1$. Furthermore, it can be verified that $P_0$ is locally stable
using the same proof as one for Corollary~5.3.1 in \cite{hv}. In
fact, we can show that there exists a nonnegative monotone
increasing continuous function $a(r)$ such that $a(|(S_1(t), E_1(t),
\cdots, S_n(t), E_n(t))|) \leq L(S_1(t), E_{1\,t}, \cdots, S_n(t),
E_{n\,t}) \leq L(S_1(0), E_{1\,0}, \cdots, S_n(0), E_{n\,0})$.
Therefore, $P_0$ is globally asymptotically stable if
$\mathcal{R}_0\leq 1$. On the other hand, if $\mathcal{R}_0
>1$, then $-L$ serves as a Lyapunov functional for system
\eqref{model}. The same proof as in Theorem~5.3.3 of \cite{hv} can
be used to show that $P_0$ is unstable.


%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

\section{Proof of Theorem~\ref{main}-(2)}
\label{section-EE}\setcounter{equation}{0}


The global stability of the endemic equilibrium of the single-group
model with delays has been proved in \cite{m1,m2}. In the following,
we consider the case $n\geq 2$. Let
$P^*=(S_1^*,E_1^*,\cdots,S_n^*,E_n^*)$ denote the unique endemic
equilibrium of system \eqref{model}. Set
\begin{equation}\label{barbeta}
\bar{\beta}_{kj}=\beta_{kj}S_k^* E_j^*, \quad 1 \leq k,j \leq n,\;
n\geq 2,\end{equation} and
\begin{equation}\label{barB}
\overline{B} = \left[\ba{cccc}\sum_{l\not=1}\bar{\beta}_{1l}
& -\bar{\beta}_{21} & \cdots & -\bar{\beta}_{n1}\\[5pt]
  -\bar{\beta}_{12} & \sum_{l\not=2}\bar{\beta}_{2l} & \cdots
  & -\bar{\beta}_{n2}\\[5pt]
  \vdots & \vdots &\ddots &\vdots\\[5pt]
  -\bar{\beta}_{1n} & -\bar{\beta}_{2n} & \cdots &
  \sum_{l\not=n}\bar{\beta}_{nl}
  \ea\right].
\end{equation}
Note that $\overline{B}$ is the Laplacian matrix of the matrix
$(\bar{\beta}_{kj})$ (see Appendix). Since $(\beta_{kj})$ is
irreducible, matrices $(\bar{\beta}_{kj})$ and $\overline{B}$ are
also irreducible. Let $C_{kj}$ denote the cofactor of the $(k,j)$
entry of $\overline{B}$. We know that system $\overline{B}v=0$ has
a positive solution $v=(v_1, v_2, \cdots, v_n)$, where
$v_k=C_{kk}>0$ for $k=1,2,\ldots,n$, by Theorem~A of Appendix.
Consider a Lyapunov functional
\begin{equation}\label{V}
V=V_1+V_2,
\end{equation}
where
\begin{equation}\label{V1}
V_1 = \displaystyle\sum\limits_{k=1}^nv_k\Big(S_k-S_k^*-S_k^*\ln
\frac{S_k}{S_k^*}+E_k-E_k^*-E_k^*\ln \frac{E_k}{E_k^*}\Big),
\end{equation}
\begin{equation}\label{V2}
V_2 =
\displaystyle\sum_{k,j=1}^nv_k\beta_{kj}S_k^{*}\int_{r=0}^{\infty}
\alpha_j(r) \Big( E_j(t-r) -E_j^* - E_j^* \ln \frac{E_j(t-r)}{E_j^*}
\Big) \;dr,
\end{equation}
and
$$\alpha_j(r)=\int_{\sigma=r}^\infty f_j(\sigma) d\sigma.$$
Differentiating $V_1$ along the solution of system \eqref{model}
and using equilibrium equations \eqref{eq1} and \eqref{eq2}, we
obtain
\begin{equation}\label{V1'}
\begin{array}{rcl}
V_1^\prime &=&\displaystyle\sum_{k=1}^n v_k \Big(\Lambda_k-d_k^SS_k-
\frac{\Lambda_kS_k^*}{S_k} +\sum\limits_{j=1}^n\beta_{kj}S_k^*\int_{r=0}^{\infty}f_j(r)E_j(t-r)dr +d_k^SS_k^* \\
& &\displaystyle \quad
-(d_k^E+\epsilon_k)E_k-\frac{E_k^*}{E_k}\sum_{j=1}^n \beta_{kj}S_k
\int_{r=0}^\infty f_j(r)E_j(t-r)dr +
(d_k^E+\epsilon_k)E_k^* \Big)\\
&=&\displaystyle\sum_{k=1}^n v_k \Big[
d_k^SS_k^*\Big(2-\frac{S_k^*}{S_k}-\frac{S_k}{S_k^*}\Big)+
\sum_{j=1}^n \beta_{kj}S_k^*E_j^* \Big(a_j\Big(2-\frac{S_k^*}{S_k}-\frac{E_k}{E_k^*}\Big)\\
&& \displaystyle\quad +\frac{1}{E_j^*}
\int_{r=0}^{\infty}f_j(r)E_j(t-r)dr -\frac{S_k E_k^* }{S_k^* E_k
E_j^*} \int_{r=0}^\infty f_j(r)E_j(t-r)dr\Big)\Big].
\end{array}
\end{equation}
Differentiating $V_2$ along the solution of system \eqref{model} and
using integration by parts, we obtain
\begin{equation}\label{V2'}
\begin{array}{rcl}
V_2^\prime &=&\displaystyle \sum_{k,j=1}^n v_k \beta_{kj} S_k^*
\int_{r=0}^{\infty} \alpha_j(r) \frac{\partial}{\partial t}\Big(
E_j(t-r) -E_j^* - E_j^* \ln \frac{E_j(t-r)}{E_j^*} \Big) \;dr\\
&=& \displaystyle \sum_{k,j=1}^n v_k  \beta_{kj} S_k^*
\int_{r=0}^{\infty} \alpha_j(r) \Big[-\frac{\partial}{\partial
r}\Big(E_j(t-r) -E_j^* - E_j^* \ln \frac{E_j(t-r)}{E_j^*} \Big)\Big] \;dr\\
&=&\displaystyle \sum_{k,j=1}^n v_k \beta_{kj} S_k^* E_j^{*}\Big(
\frac{a_j E_j}{E_j^{*}}- \frac{1}{E_j^{*}}
\int_{r=0}^{\infty}f_j(r)E_j(t-r)dr
-\int_{r=0}^{\infty}f_j(r)\ln\frac{ E_j(t)}{E_j(t-r)}\;dr\Big).
\end{array}
\end{equation}
Combining \eqref{V1'} and \eqref{V2'} and using expression
\eqref{barbeta}, we have
\begin{equation}\label{V'}
\begin{array}{rcl}
V^\prime &=&\displaystyle\sum_{k=1}^n v_k
d_k^SS_k^*\Big(2-\frac{S_k^*}{S_k}-\frac{S_k}{S_k^*}\Big) +
\sum_{k,j=1}^n v_k \bar{\beta}_{kj}\Big[ a_j\Big(2-\frac{S_k^*}{S_k}
-\frac{E_k}{E_k^*} +
\frac{E_j}{E_j^{*}}\Big) \\
& & \displaystyle \quad -\frac{S_k E_k^* }{S_k^* E_k E_j^*}
\int_{r=0}^\infty f_j(r)E_j(t-r)dr
-\int_{r=0}^{\infty}f_j(r)\ln\frac{
E_j(t)}{E_j(t-r)}dr\Big]\\
& \leq & \displaystyle \sum_{k,j=1}^n v_k \bar{\beta}_{kj}
\int_{r=0}^{\infty}f_j(r) \Big[
\frac{E_j}{E_j^{*}}-\frac{E_k}{E_k^*}- \ln \frac{S_k^*}{S_k} \cdot
\frac{S_k E_k^* E_j(t-r)}{S_k^* E_k E_j^*}
\cdot \frac{ E_j(t)}{E_j(t-r)}\\
& & \displaystyle \quad
+\Big(1-\frac{S_k^*}{S_k}+\ln\frac{S_k^*}{S_k}\Big) + \Big( 1 -
\frac{S_k E_k^* E_j(t-r)}{S_k^* E_k E_j^*} + \ln\frac{S_k E_k^*
E_j(t-r)}{S_k^* E_k E_j^*}\Big)\Big]
 \;dr\\
&\leq&\displaystyle \sum_{k,j=1}^n v_k \bar{\beta}_{kj} \Big(
\frac{E_j}{E_j^{*}}-\frac{E_k}{E_k^*}\Big)
-\sum_{k,j=1}^n v_k \bar{\beta}_{kj} \ln\frac{E_k^*E_j}{E_kE_j^*}\\
&=:&H_1 - H_2.
\end{array}
\end{equation}
In the above derivation, we have used two facts:
$\frac{S_k^*}{S_k}+\frac{S_k}{S_k^*}\geq 2$ with equality holding
if and only if $S_k=S_k^*$, and $1-x+\ln x \leq 0$ for all $x>0$
with equality holding if and only if $x=1$.

We first show $H_1\equiv 0$ for all $E_1, E_2, \cdots, E_n > 0$. It
follows from $\overline{B}v=0$ that
$$\sum_{j=1}^n\bar{\beta}_{jk}v_j=\sum_{i=1}^n\bar{\beta}_{ki}v_k,$$
or, using $\bar{\beta}_{jk}=\beta_{jk}S_j^*E_k^*$,
$$\sum_{j=1}^n \beta_{jk}S_j^*E_k^*v_j=\sum_{i=1}^n
\beta_{ki}S_k^*E_i^*v_k, \quad k=1,2,\ldots,n.$$ This implies that
$$\sum_{k,j=1}^n v_k \beta_{kj}S_k^*E_j
=\sum_{k=1}^n E_k \sum_{j=1}^n \beta_{jk}S_j^*v_j =\sum_{k=1}^n
\frac{E_k}{E_k^*} \sum_{i=1}^n
\beta_{ki}S_k^*E_i^*v_k=\sum_{k,j=1}^n v_k
\beta_{kj}S_k^*E_j^*\frac{E_k}{E_k^*},$$ and thus $H_1\equiv 0$ for
all $E_1, E_2, \cdots, E_n >0$.

Next we show $H_2\equiv 0$ for all $E_1, E_2, \cdots, E_n >0$. Let
$G$ denote the directed graph associated with matrix
$(\bar{\beta}_{kj})$. $G$ has vertices $\{1,2,\cdots, n\}$ with a
directed arc $(k,j)$ from $k$ to $j$ iff $\bar{\beta}_{kj}\not=0$.
$E(G)$ denotes the set of all directed arcs of $G$. Using
Kirchhoff's Matrix-Tree Theorem (see Theorem A in Appendix), we know
that $v_k=C_{kk}$ can be interpreted as a sum of weights of all
directed spanning subtrees $T$ of $G$ that are rooted at vertex $k$.
Consequently, each term in $v_k\bar{\beta}_{kj}$ is the weight
$w(Q)$ of a unicyclic subgraph $Q$ of $G$, obtained from such a tree
$T$ by adding a directed arc $(k,j)$ from the root $k$ to vertex
$j$. Note that the arc $(k,j)$ is part of the unique cycle $CQ$ of
$Q$, and that the same unicyclic graph $Q$ can be formed when each
arc of $CQ$ is added to a corresponding rooted tree $T$. Therefore,
the double sum in $H_2$ can be reorganized as a sum over all
unicyclic subgraphs $Q$ containing vertices $\{1,2,\cdots,n\}$, that
is,
\begin{equation}\label{H2}
H_2=\sum_{Q}H_Q,
\end{equation} where
\begin{equation}\label{HQ}\begin{array}{rcl}
H_Q &=& \displaystyle w(Q) \;\cdot\;
 \sum_{(k,j)\in E(CQ)} \ln\frac{E_k^*E_j}{E_kE_j^*}\\
&=&\displaystyle w(Q) \;\cdot \;\ln\; \Big(\prod_{(k,j)\in E(CQ)}
\frac{E_k^*E_j}{E_kE_j^*}\Big).
\end{array}
\end{equation}
Since $E(CQ)$ is the set of arcs of a cycle $CQ$, we have
$$\prod_{(k,j)\in E(CQ)}
\frac{E_k^*E_j}{E_kE_j^*} = 1, \quad \mbox{and\; thus} \quad \ln
\Big(\prod_{(k,j)\in E(CQ)} \frac{E_k^*E_j}{E_kE_j^*}\Big)=0.$$ This
implies $H_Q=0$ for each $Q$, and hence $H_2\equiv 0$ for all $E_1,
E_2, \cdots, E_n >0$. Therefore, we obtain $V^\prime \leq 0$ for all
$(S_1, E_1(\cdot), \cdots, S_n, E_n(\cdot))\in
\overset{\circ}{\Theta}$. Furthermore, if $\beta_{kj}\not =0$, then,
by \eqref{V'}, $V^\prime=0$ implies
$$\Big(1-\frac{S_k^*}{S_k}+\ln\frac{S_k^*}{S_k}\Big)
+\Big(1-\frac{S_k E_k^* E_j(t-r)}{S_k^* E_k E_j^*}+\ln\frac{S_k
E_k^* E_j(t-r)}{S_k^* E_k E_j^*}\Big) =0,$$ for all $r\in [0,
\infty)$ and $t>0$. Using the fact that $1-x+\ln x \leq 0$ for all
$x>0$ with equality holding iff $x=1$, we obtain
\begin{equation}\label{ratio}
S_k=S_k^*, \quad \frac{E_k}{E_k^*}=\frac{E_j(t-r)}{E_j^*}, \quad
r\in [0, \infty), \quad t>0.\end{equation} Now let $p$ and $q$
denote any two distinct groups, namely, two distinct vertices of the
directed graph $G$ associated with the irreducible matrix
$(\bar{\beta}_{kj})$. Then, by the strong connectivity of $G$, there
exists oriented path between $p$ and $q$. Applying \eqref{ratio} to
each arc $(k,j)$ of such a path, we can see that
$$S_p = S_p^*, \quad \frac{E_p(t-r)}{E_p^*} =
\frac{E_q(t-r)}{E_q^*} = c, \quad r\in [0,\infty), \quad t>0,\quad
\mbox{\;for\, all\;} p,\, q.$$ Therefore, $V^\prime =0$ if and only
if
\begin{equation}\label{ratio2}
S_k=S_k^*,\quad E_k(t-r)\equiv c E_k^*, \quad k=1, 2, \cdots, n,
\quad r\in[0,\infty),\quad t>0, \end{equation} where $c>0$ is an
arbitrary number. Substituting \eqref{ratio2} into the first
equation of system \eqref{model}, we obtain
\begin{equation}\label{V'0}
0 = \Lambda_k-c \sum_{j=1}^n\beta_{kj}a_jS_k^*E_j^*-d_k^SS_k^*.
\end{equation}
The right-hand-side of \eqref{V'0} is strictly decreasing in $a$. By
(\ref{eq1}), we know that \eqref{V'0} holds iff $c=1$, namely at
$P^*.$ Therefore, the only compact invariant subset of the set where
$V^\prime=0$ is the singleton $\{P^*\}$. By a similar argument as in
Section~\ref{section-DFE}, $P^*$ is globally asymptotically stable
in $\overset{\circ}{\Theta}$ if $R_0>1$. This establishes
Theorem~\ref{main}-(2).
%\end{proof}




%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

\vspace{10pt} \noindent{\bf Acknowledgments.}

This research was supported in part by grants from the Natural
Science and Engineering Research Council of Canada (NSERC) and
Canada Foundation for Innovation (CFI). Z. Shuai acknowledges the
support of an Izaak Walton Killam Memorial Scholarship. C. Wang
acknowledges the support of the State Scholarship Fund from China
Scholarship Council. The authors would like to thank Dr. C.
Connell McCluskey for helpful discussions and thank an anonymous
referee for valuable suggestions.


\vspace{20pt}



%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
\begin{appendix}
\setcounter{equation}{0} \setcounter{thm}{0}
\renewcommand{\thesection}{\Alph{section}}
\renewcommand{\theequation}{\Alph{section}.\arabic{equation}}
\setcounter{section}{1}
\renewcommand{\thethm}{\Alph{thm}}
\noindent\textbf{Appendix: Kirchhoff's Matrix-Tree
Theorem}\vspace{10pt}

%For $n\geq 2$, an $n\times n$ matrix $A=(a_{ij})$ is
%\textit{reducible} if, for some permutation matrix $P$,
%$$PAP^T=\left[
%          \begin{array}{cc}
%            A_1 & 0 \\
%            A_2 & A_3 \\
%          \end{array}
%        \right],
%$$
%where $A_1$ and $A_3$ are square matrices. The following
%properties of irreducible matrices are standard \cite{bp}.
%\begin{enumerate}
%    \item[P1.] If $A$ is nonnegative (namely, $a_{ij}\geq 0$) and irreducible,
%    then $\rho(A)$ is a simple eigenvalue,
%    and $A$ has a positive eigenvector $x$
%    corresponding to $\rho(A)$.
%    \item[P2.] If $A$ is nonnegative and irreducible,
%    and $D$ is diagonal and positive, then $AD$ is irreducible.
%\end{enumerate}
%
%Irreducibility of $A$ can be verified using the associated
%directed graph.

Given a nonnegative matrix $A=(a_{ij})$, the \textit{directed
graph} $G(A)$ associated with $A=(a_{ij})$ has vertices $\{1, 2,
\cdots, n\}$ with a directed arc $(i,j)$ from $i$ to $j$ iff
$a_{ij}\not=0$. It is \textit{strongly connected} if any two
distinct vertices are joined by an oriented path. Matrix $A$ is
irreducible if and only if $G(A)$ is strongly connected \cite{bp}.

A \textit{tree} is a connected graph with no cycles. A subtree $T$
of a graph $G$ is said to be \textit{spanning} if $T$ contains all
the vertices of $G$. A \textit{directed tree} is a tree in which
each edge has been replaced by an arc directed one way or the
other. A directed tree is said to be \textit{rooted} at a vertex,
called the root, if every arc is oriented in the direction towards
to the root. An \textit{oriented cycle} in a directed graph is a
simple closed oriented path. A \textit{unicyclic graph} is a
directed graph consisting of a collection of disjoint rooted
directed trees whose root are on an oriented cycle. We refer the
reader to \cite{mo} for more details of these concepts.

For a given nonnegative matrix $A=(a_{ij})$, let
\begin{equation}\label{K} L =
\left[\ba{cccc}\sum_{l\not=1} a_{1l}
& -a_{21} & \cdots & -a_{n1}\\[5pt]
  -a_{12} & \sum_{l\not=2}a_{2l} & \cdots
  & -a_{n2}\\[5pt]
  \vdots & \vdots &\ddots &\vdots\\[5pt]
  -a_{1n} & -a_{2n} & \cdots &
  \sum_{l\not=n}a_{nl}
  \ea\right]
\end{equation} be the Laplacian matrix of the directed graph $G(A)$
and $C_{ij}$ denote the cofactor of the $(i,j)$ entry of $L$. For
the linear system
\begin{equation}\label{linear}
L v=0,
\end{equation} the following result holds, see \cite{gls2}.

\begin{thm}[Kirchhoff's Matrix-Tree Theorem]\label{MTT}
Assume that $n\geq 2$ and that $A$ is irreducible. Then following
results hold.
\begin{enumerate}
\item[$(1)$] The solution space of system \eqref{linear} has dimension 1,
with a basis $(v_1, v_2, \cdots, v_n) = (C_{11}, C_{22}, \cdots,
C_{nn}).$
\item[$(2)$] For $1\leq k \leq n$,
\begin{equation}\label{Ckk}
C_{kk}=\sum_{T \in \mathbb{T}_k} w(T) = \sum_{T \in
\mathbb{T}_k}\prod_{(r,m)\in E(T)}a_{rm}>0,
\end{equation}
where $\mathbb{T}_k$ is the set of all directed spanning subtrees
of $G(A)$ that are rooted at vertex $k$, $w(T)$ is the weight of a
directed tree $T$, and $E(T)$ denotes the set of directed arcs in
$T$.
\end{enumerate}
\end{thm}
\end{appendix}



%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
\begin{thebibliography}{9999}

\bibitem{ah}
F. V. Atkinson and J. R. Haddock, { On determining phase spaces
for functional differential equations}, Funkcial. Ekvac. {31}
(1988) 331--347.

\bibitem{bc}
E. Beretta and V. Capasso, { Global stability results for a
multigroup SIR epidemic model}, in: T. G. Hallam, L. J. Gross, and
S. A. Levin (Eds.), Mathematical Ecology,  World Scientific,
Teaneck, NJ, 1988, pp. 317--342.

\bibitem{bcr}
E. Beretta, V. Capasso, and F. Rinaldi, { Global stability results
for a generalized Lotka-Volterra system with distributed delays:
applications to predator-prey and to epidemic systems}, J. Math.
Biol. {26} (1988) 661--688.

\bibitem{bhmt}
E. Beretta, T. Hara, W. Ma, and Y. Takeuchi, { Global asymptotic
stability of an $SIR$ epidemic model with distributed time delay},
Nonlinear Anal. {47} (2001) 4107--4115.

\bibitem{bt}
E. Beretta and Y. Takeuchi, { Global stability of an SIR epidemic
model with time delays}, J. Math. Biol. {33} (1995) 250--260.


\bibitem{bp}
A. Berman and R. J. Plemmons, { Nonnegative Matrices in the
Mathematical Sciences}, Academic Press, New York, 1979.

%\bibitem{c}
%
%K. L. Cooke, {Stability analysis for a vector disease model},
%Rocky Mount. J. Math. {9} (1979) 31--42.


\bibitem{dhm}
O. Diekmann, J. A. P. Heesterbeek, and J. A. J. Metz, { On the
definition and the computation of the basic reproduction ratio
$R_0$ in models for infectious diseases in heterogeneous
populations}, J. Math. Biol. {28} (1990) 365--382.


\bibitem{gls1}
H. Guo, M. Y. Li, and Z. Shuai, { Global stability of the endemic
equilibrium of multigroup SIR epidemic models}, Canadian Appl.
Math. Quart. {14} (2006) 259--284.

\bibitem{gls2}
H. Guo, M. Y. Li, and Z. Shuai, { A graph-theoretic approach to
the method of global Lyapunov functions}, Proc. Amer. Math. Soc.
{136} (2008) 2793--2802.


\bibitem{hv}
J. K. Hale and S. M. V. Lunel, { Introduction to Functional
Differential Equations}, Applied Mathematical Sciences, Vol. 99,
Springer, New York, 1993.


\bibitem{he}
H. W. Hethcote, { An immunization model for a heterogeneous
population}, Theor. Popu. Biol. {14} (1978) 338--349.

\bibitem{ht}
H. W. Hethcote and H. R. Thieme, { Stability of the endemic
equilibrium in epidemic models with subpopulations}, Math. Biosci.
{75} (1985) 205--227.

\bibitem{hmn}
Y. Hino, S. Murakami, and T. Naito, { Functional Differential
Equations with Infinite Delay}, Lecture Notes in Mathematics, Vol.
1473, Springer, Berlin, 1991.


\bibitem{hcc}
W. Huang, K. L. Cooke, and C. Castillo-Chavez, { Stability and
bifurcation for a multiple-group model for the dynamics of
HIV/AIDS transmission}, SIAM J. Appl. Math. {52} (1992) 835--854.


\bibitem{ly}
A. Lajmanovich and J. A. York, { A deterministic model for
gonorrhea in a nonhomogeneous population}, Math. Biosci. {28}
(1976) 221--236.

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


\bibitem{ls}
X. Lin and J. W.-H. So, { Global stability of the endemic
equilibrium and uniform persistence in epidemic models with
subpopulations}, J. Austral. Math. Soc. Ser. B {34} (1993)
282--295.

\bibitem{mthb}
W. Ma, Y. Takeuchi, T. Hara, and E. Beretta, { Permanence of an
SIR epidemic model with distributed time delays}, Tohoku Math. J.
{54} (2002) 581--591.


\bibitem{m1}
C. C. McCluskey, { Complete global stability for an SIR epidemic
model with delay - distributed or discrete}, to appear in
Nonlinear Anal. RWA.

\bibitem{m2}
C. C. McCluskey, { Global stability for an SEIR epidemiological
model with varying infectivity and infinity delay}, Math. Biosci.
Eng. {6} (2009) 603--610.

\bibitem{mo}
J. W. Moon, { Counting Labelled Trees}, Canadian Mathematical
Congress, Montreal, 1970.


\bibitem{rw}
G. R\"{o}st and J. Wu, { SEIR epidemiological model with varying
infectivity and finite delay}, Math. Biosci. Eng. {5} (2008)
389--402.


\bibitem{th1}
H. R. Thieme, { Local stability in epidemic models for
heterogeneous populations}, in: V. Capasso, E. Grosso,  and S. L.
Paveri-Fontana (Eds.), Mathematics in Biology and Medicine,
Lecture Notes in Biomathematics, Vol. 57, Springer, Berlin, 1985,
pp. 185--189.

\bibitem{th2}
H. R. Thieme, { Mathematics in Population Biology}, Princeton
University Press, Princeton, 2003.

\bibitem{vw}
P. van den Driessche and J. Watmough, { Reproduction numbers and
sub-threshold endemic equilibria for compartmental models of
disease transmission}, Math. Biosci. {180} (2002) 29--48.

\bibitem{wz}
W. Wang and X.-Q. Zhao, {An epidemic model with population
dispersal and infection period}, SIAM J. Appl. Math. {66} (2006)
1454--1472.


\end{thebibliography}

\end{document}
