# Stability Analysis of a Multigroup Epidemic Model with General Exposed Distribution and Nonlinear Incidence Rates

Authors: Ling Zhang; Jingmei Pang; Jinliang Wang
Journal: Abstract and Applied Analysis
 (2013)

---

## Abstract

We investigate a class of multigroup epidemic models with general exposed distribution and nonlinear incidence rates. For a simpler case that assumes an identical natural death rate for all groups, and with a gamma distribution for exposed distribution is considered. Some sufficient conditions are obtained to ensure that the global dynamics are completely determined by the basic production numberR0. The proofs of the main results exploit the method of constructing Lyapunov functionals and a graph-theoretical technique in estimating the derivatives of Lyapunov functionals.


---

## Body

## 1. Introduction

Multigroup epidemic models have been used in the literature to describe the transmission dynamics of many different infectious diseases such as mumps, measles, gonorrhea, HIV/AIDS and vector borne diseases such as Malaria [1]. In the models, 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 intergroup interactions can be modeled separately. It is well known that global dynamics of multigroup models with higher dimensions, especially the global stability of the endemic equilibrium, are a very challenging problem. Guo et al. [2] proposed a graph-theoretic approach to the method of global Lyapunov functions and used it to resolve the open problem on the uniqueness and global stability of the endemic equilibrium of a multigroup SIR model with varying subpopulation sizes. Subsequently, a series of studies on the global stability of multigroup epidemic models were produced in the literature (see e.g., [2–5]).In the present paper, a more general multigroup epidemic model is proposed and studied to describe the disease spread in a heterogeneous host population with general exposed distribution and nonlinear incidence rate. The host population is divided intom distinct groups (m≥1). For 1≤i≤m, the ith group is further partitioned into four disjoint classes: the susceptible individuals, exposed individuals, infectious individuals, and recovered individuals, whose numbers of individuals at time t are denoted by Si(t),Ei(t),Ii(t), and Ri(t), respectively. Susceptible individuals infected with the disease but not yet infective are in the exposed (latent) class.It is pointed in [6] that a fixed latent period can be considered as an approximation of the mean latent period, and this would be appropriate for those diseases whose latent periods vary only relatively slightly. For example, poliomyelitis has a latent period of 1–3 days (comparing to its much longer infectious period of 14–20 days). However disease such as tuberculosis, including bovine tuberculosis (a disease spread from animal to animal mainly by direct contact), may take months to develop to the infectious stage and also can relapse. Since the time it takes from the moment of new infection to the moment of becoming infectious may differ from disease to disease, even for the same disease, it differs from individual to individual, and it is indeed a random variable. It is thus of interest from both mathematical and biological viewpoints to investigate whether sustained oscillations are the result of general exposed distribution.Following the method of [6], we also assume that the disease does not cause deaths during the latent period, taking the natural death rate into consideration. Let P(t) denote the probability that an exposed individual remains in the time t after entering the exposed class. For 1≤i,j≤m, βij≥0 denotes the coefficient of transmission between compartments Si and Ij. It is assumed that m-square matrix (βij)1≤i,j≤m is irreducible [7]. So the proportion of exposed individuals can be expressed by the integral
(1)Ei(t)=∑j=1m‍βij∫0t‍fij(Si(u),Ij(u))e-δj(t-u)Pj(t-u)du,
where the sum takes into account cross-infections from all groups. Integrals in (1) are in the Riemann-Stieltjes sense. Pj(t) satisfies the following reasonable properties: (A)
Pj:[0,∞)→[0,1] is nonincreasing, piecewise continuous with possibly finitely many jumps and satisfies Pj(0+)=1, and limt→∞Pj(t)=0 with ∫0∞‍Pj(t)dt is positive and finite.Differentiating (1) gives
(2)Ei′(t)=∑j=1m‍βijfij(Si(t),Ij(t))+∑j=1m‍βij∫0t‍fij(Si(u),Ij(u))e-δj(t-u)×Pj′(t-u)du-δiEi(t).
The first term on the right hand side in (2) is the rate at which new infected individuals come into the exposed class, and the last term explains the natural deaths. The second term accounts for the rate at which the individuals move to the infectious class (noting that Pj′(t-u)≤0 due to the aformentioned property) from the exposed class; hence
(3)Ii′(t)=-∑j=1m‍βij∫0t‍fij(Si(u),Ij(u))e-δj(t-u)×Pj′(t-u)du-(δi+εi+γi)Ii(t).
Let hj(t)=-Pj′(t) be the probability density function for the time (a random variable) it takes for an infected individual in the ith group to become infectious. Then (4) becomes
(4)Ii′(t)=∑j=1m‍βij∫0t‍fij(Si(u),Ij(u))e-δj(t-u)×hj(t-u)du-(δi+εi+γi)Ii(t).Within theith group, φi(Si) denotes the growth rate of Si, which includes both the production and the natural death of susceptible individuals. Therefore, under the assumptions, the model to be studied takes the following differential and integral equations form:
(5)Si′(t)=φi(Si(t))-∑j=1m‍βijfij(Si(t),Ij(t)),Ei′(t)=∑j=1m‍βijfij(Si(t),Ij(t))-∑j=1m‍βij∫0t‍fij(Si(u),Ij(u))e-δj(t-u)×hj(t-u)du-δiEi(t),Ii′(t)=∑j=1m‍βij∫0t‍fij(Si(u),Ij(u))e-δj(t-u)×hj(t-u)du-(δi+εi+γi)Ii(t),Ri′(t)=γiIi(t)-δiRi(t).Since the variablesEi and Ri do not appear in the first and third equations of model (5), Ei(t) and Ri(t), i=1,…,m, can be decoupled from the Si(t) and Ii(t) equations; we only need to consider the subsystem of (5) consisting of only the Si and Ii equations:
(6)Si′(t)=φi(Si(t))-∑j=1m‍βijfij(Si(t),Ij(t)),Ii′(t)=∑j=1m‍βij∫0t‍fij(Si(u),Ij(u))e-δj(t-u)×hj(t-u)du-(δi+εi+γi)Ii(t),
where δi denotes the natural death rates of Ii compartments in the ith group, εi is the death rate caused by disease in the ith group, and γi is the rate of recovery of infectious individuals in the ith group. In what follows we investigate the global stability of system (5).Whenm=1, P(t)=eϵt, and with bilinear incidence rate, system (5) will reduce to the standard SEIR ordinary differential equation (ODE) model studied in [8, 9], and if we further assume that P(t) is a step function, system (5) becomes the SEIR model with a discrete delay studied in [10]. Recently, a model of this type, but including the possibility of disease relapse, has been proposed in [11, 12] to investigate the transmission of herpes, and its global dynamics have been completely investigated in [5, 13].To express the main idea and the approaches more clearly, we consider a simpler case in which all groups share the same natural death rate:δj=δ for j=1,2,…,m. Further, we assume that the functions hj(u) are disease specific only, implying that hj(u)=h(u) for j=1,2,…,m. We choose the gamma distribution:
(7)h(u)=hn,b(u)=un-1(n-1)!bne-u/b,
where b>0 is a real number and n>1 is an integer, which is widely used and can approximate several frequently used distributions. For example, when b→0+, hn,b(s) will approach the Dirac delta function, and when n=1, hn,b(s) is an exponentially decaying function.The main object of this paper is to carry out the well-known “linear chain trick” to system (6), transfer system into an equivalent ordinary differential equations system, and establish its global dynamics. We derive the basic reproductive number R0 and show that R0 completely determines the global dynamics of system (6). More specifically, if R0≤1, the disease-free equilibrium is globally asymptotically stable and the disease dies out; if R0>1, a unique endemic equilibrium exists and is globally asymptotically stable, and the disease persists at the endemic equilibrium. The global stability of P* rules out any possibility for Hopf bifurcations and the existence of sustained oscillations. We should point out here that this work is motivated by Yuan and Zou [11, 12, 14].  In the proof we demonstrate that the graph-theoretic approach developed in [2, 3] can be successfully applied to construct suitable Lyapunov functionals and thus prove the global stability of the endemic equilibrium for model (6) with general exposed distribution and nonlinear incidence rate. Our work is also based on a recent work by Sun and Shi [15], which resolved the dynamics of multigroup SEIR epidemic models with nonlinear incidence of infection and nonlinear removal functions between compartments.In Section2, we first give the model, preliminaries and the basic reproduction number R0. The global stability of the corresponding equilibria for R0≤1 and R0>1 is shown, respectively, in Section 3—the key results of this paper. And in Section 4, some numerical simulations are shown to illustrate the effectiveness of the proposed result.


## 2. Preliminaries

We make the following basic assumptions for the intrinsic growth rate of susceptible individuals in theith group φi(Si) and the transmission functions fij(Si,Ij).(A1)
φi are C1 non-increasing functions on [0,∞) with φi(0)>0, and there is a unique positive solution ξ=Si0 for the equation φi(ξ)=0. φi(S)>0 for 0≤S<Si0, and φi(S)<0 for S>Si0; that is
(8)[φi(Si)-φi(Si0)](Si-Si0)<0,forSi≠Si0,i=1,2,…,m.(A2)
fij(Si,Ij)≤cij(Si)Ij for all Ij>0.(A3)
cij(Si)≤cij(Si0), 0<Si<Si0, i,j=1,…,m.Following the technique and method in [14], define
(9)b^≡b1+δb,
which can absorb the exponential term e-δu into the delay kernel. The second equation in (6) can be rewritten as
(10)Ii′(t)=∑j=1m‍βij(1+δb)n∫0t‍fij(Si,Ij)hn,b^(t-u)du-(δ+εi+γi)Ii.
For l=1,…,n, let(11)yi,l(t)=∑j=1m‍βijb^(1+δb)n∫0t‍fij(Si,Ij)hl,b^(t-u)du,i=1,2,…,m.
Thus, for l∈{2,…,n}, we obtain
(12)y˙i,l=hl,b^(0)∑j=1m‍βijb^(1+δb)nfij(Si,Ij)+∑j=1m‍βijb^(1+δb)n∫-∞t‍(l-1)(t-u)l-2(l-1)!b^le-(t-u)/b^fij(Si,Ij)du-∑j=1m‍βijb^(1+δb)n∫-∞t‍(t-u)l-1(l-1)!b^l+1e-(t-u)/b^fij(Si,Ij)du=[yi,l-1-yi,l]b^.
For l=1, we have
(13)yi,1=∑j=1mβijb^(1+δb)n∫-∞te-(t-u)/b^b^fij(Si,Ij)du,i=1,…,m.
It follows that
(14)y˙i,1=∑j=1m‍βij(1+δb)nfij(Si,Ij)-∑j=1m‍βij(1+δb)n∫-∞t‍e-(t-u)/b^b^fij(Si,Ij)du=∑j=1m‍βij(1+δb)nfij(Si,Ij)-1b^yi,1,i=1,…,m.Thus the integro-differential system (6) is equivalent to the ordinary differential equations
(15)Si′(t)=φi(Si(t))-∑j=1m‍βijfij(Si(t),Ij(t)),yi,1′(t)=1(1+δb)n∑j=1m‍βijfij(Si(t),Ij(t))-1b^yi,1(t),yi,2′(t)=1b^(yi,1(t)-yi,2(t)),i=1,2,…,m,⋮yi,n′(t)=1b^(yi,n-1(t)-yi,n(t)),Ii′(t)=1b^yi,n(t)-(δ+εi+γi)Ii(t).For initial condition(16)(S1(0),y1,1(0),…,y1,n(0),I1(0),S2(0),y2,1(0),…,y1,n(0),I2(0),…,Sm(0),ym,1(0),…,ym,n(0),Im(0))∈Rm(n+2),
the existence, uniqueness, and continuity of the solution (Si,yi,1,yi,2,…,yi,n,Ii) of system (15) follow from the standard theory of Volterra integro-differential equation [16]. It can also be verified that every solution of (15) with nonnegative initial condition remains nonnegative.It follows from (A1) and the first equation of (15) that limsupt→∞Si(t)≤Si0 for all i=1,2,…,m. Let Nφi be the maximum of the function φi on R+ and let q be a positive real number such that q>b^Nφi. Denote by Υi the ith tube for system (15); that is,
(17)Υi=(Si,yi,1,yi,2,…,yi,n,Ii).
It follows from a similar argument to that in [14] that  we can show that the set Dϵ defined by
(18)Γϵ={(Si,yi,1,yi,2,…,yi,n,Ii)∈R+m(n+2)∣Si≤Si0+ϵ,Si+(1+δb)nyi,1≤q+Si0,yi,l≤q+Si0+lϵ(1+δb)n,Ii≤q+Si0+(n+1)ϵb^(1+δb)n(δ+ϵi+γi),i=1,2,…,m,l=2,3,…,n{R+m(n+2)}}
is a forward invariant compact absorbing set for system (15) for ϵ>0 and that the set Γ0 (i.e., when ϵ=0) is a forward invariant compact set.Under the assumption(A1), we know that system (15) always has the disease-free equilibrium
(19)P0=(S10,0,…,0,I10,S20,0,…,0,I20,…,Sm0,0,…,0,Im0)∈Rm(n+2).
An equilibrium P* of (6) has the form P*=(S1*,I1*,S2*,I2*,…,Sm*,Im*)∈R2m with Si*>0, Ii*>0, i=1,…,m. Translating to the equivalent system (15), P* is corresponding to
(20)P¯*=(S1*,y1,1*,…,y1,n*,I1*,S2*,y2,1*,…,y2,n*,I2*,…,Sm*,ym,1*,…,ym,n*,Im*)∈Rm(n+2).P¯* in the interior of Γ0 is called an endemic equilibrium, and it satisfies the following equilibrium equations:
(21)0=φi(Si*)-∑j=1m‍βijfij(Si*,Ij*),0=1(1+δb)n∑j=1m‍βijfij(Si*,Ij*)-1b^yi,1*,0=1b^(yi,1*-yi,2*),⋮0=1b^(yi,n-1*-yi,n*),0=1b^yi,n*-(δ+εi+γi)Ii*.
The basic reproduction number R0 is defined as the expected number of secondary cases produced by single infectious individual during its entire period of infectiousness in a completely susceptible population. For system (15), we can calculate it as the spectral radius of a matrix called the next generation matrix. Let
(22)ℱ=(c11(S10)β11(1+δb)n⋯c1m(S10)β1m(1+δb)n⋮⋱⋮cm1(Sm0)βm1(1+δb)n⋯cmm(Sm0)βmm(1+δb)n),𝒱=diag(δ+εi+γi)=(δ+ε1+γ10⋯00δ+ε2+γ2⋯0⋮⋮⋱⋮00⋯δ+εm+γm).
Then the next generation matrix is
(23)ℱ𝒱-1=(c11(S10)β11(1+δb)n(δ+ε1+γ1)⋯c1m(S10)β1m(1+δb)n(δ+εm+γm)⋮⋱⋮cm1(S10)βm1(1+δb)n(δ+ε1+γ1)⋯cmm(S10)βmm(1+δb)n(δ+εm+γm)),
and hence, the basic reproduction number R0 is
(24)R0=ρ(ℱ𝒱-1)=max{|λ|;λ∈σ(ℱ𝒱-1)},
where ρ(·) and σ(·) denote the spectral radius and the set of eigenvalues of a matrix, respectively. Since it can be verified that system (15) satisfies conditions (A1)–(A5) of Theorem 2 of [17], we have the following proposition.Lemma 1.
For system (15), the disease-free equilibrium P0 is locally asymptotically stable if R0<1, while it is unstable if R0>1.Following the method of [2], one defines a matrix
(25)M0=𝒱-1ℱ=(c11(S10)β11(1+δb)n(δ+ε1+γ1)⋯c1m(S10)β1m(1+δb)n(δ+ε1+γ1)⋮⋱⋮cm1(S10)βm1(1+δb)n(δ+εm+γm)⋯cmm(S10)βmm(1+δb)n(δ+εm+γm)),
whose spectral radius has a similar threshold property to that of R0, since both of the nonnegative matrices ℱ𝒱-1 and M0 are irreducible, and hence from the Perron-Frobenius theorem [7] that their spectral radii are given by each of their simple eigenvalues. Thus, we obtain R0=ρ(ℱ𝒱-1)=ρ(𝒱-1ℱ)=ρ(M0). Then the following lemma immediately follows.Lemma 2.
ρ
(
M
0
)
≤
1 if and only if R0≤1.


## 3. Main Results

The following main theorems are summarized in terms of system (15).Theorem 3.
Assume that the functionsφi and fij satisfy assumptions (A1)–(A3), and the matrix B=(βij)m×m is irreducible and R0 is defined by (24). (i)
IfR0≤1, then P0 is the unique equilibrium of system (15), and P0 is globally asymptotically stable in Γ0.(ii)
IfR0>1, then P0 is unstable, and system (15) is uniformly persistent in Γ0.Proof.
Let us defineM(S)=(βijcij(Si)/(1+δb)n(δ+εi+γi))m×m, where S=(S1,S2,…,Sm)T. Note that M(S0)=M0. Since B=(βij)m×m is irreducible, the matrix M0 is also irreducible.
First we claim that there does not exist any endemic equilibriumP¯* in Ω. Suppose that S≠S0. Then we have 0<M(S)<M0. Since nonnegative matrix M(S)+M0 is irreducible, it follows from the Perron-Frobenius theorem (see Corollary  2.1.5 of [7]) that ρ(M(S))<ρ(M0)≤1. This implies that M(S)I=I has only the trivial solution I=0, where I=(I1,…,Im)T. Hence the claim is true. Next we claim that the disease-free equilibrium P0 is globally asymptotically stable in Γ0. From the Perron-Frobenius (see Theorem  2.1.4 of [7]), the nonnegative irreducible matrix M0 has a strictly positive left eigenvector (ω1,ω2,…,ωm) associated with the eigenvalue ρ(M0) such that
(26)(ω1,ω2,…,ωm)ρ(M0)=(ω1,ω2,…,ωm)M0.
Using this positive eigenvector, we construct the following Lyapunov function:
(27)VDFE=∑i=1m‍ωiδ+εi+γi(∑j=1n‍yi,j+Ii).
Computing the derivative of VDFE along the solutions of (15) in Γ0, we get
(28)VDFE′=∑i=1m‍[∑j=1m‍ωiβij(1+δb)n(δ+εi+γi)fij(Si,Ij)-ωiIi]≤∑i=1m‍[∑j=1m‍ωiβijcij(Si)(1+δb)n(δ+εi+γi)Ij-ωiIi]≤∑i=1m‍[∑j=1m‍ωiβijcij(Si0)(1+δb)n(δ+εi+γi)Ij-ωiIi]=(ω1,…,ωm)[M0I-I]=[ρ(M0)-1](ω1,…,ωm)I.
Thus, under the assumption R0=ρ(M0)<1, VDFE′≤0, and VDFE′=0 if and only if I=0 and S=S0=(S10,S20,…,Sm0). Suppose that ρ(M0)=1. Then it follows from the previous that VDFE′=0 implies
(29)(ω1,…,ωm)M0I=(ω1,…,ωm)I.
Hence, if S≠S0, then (ω1,…,ωm)M(S)<(ω1,…,ωm)M0=ρ(M0)(ω1,…,ωm)=(ω1,…,ωm) and thus I=0 is the only solution of (29). Summarizing the statements, we see that VDFE′=0 if and only if I=0 or S=S0, which implies that the compact invariant subset of the set where VDFE′=0 is only the singleton P0. Thus, by LaSalle’s invariance principle [18], it follows that the disease-free equilibrium E0 is globally asymptotically stable in Γ0.
IfR0=ρ(M0)>1, then
(30)(ω1,ω2,…,ωm)M0-(ω1,ω2,…,ωm)=[ρ(M0)-1](ω1,ω2,…,ωm)>0,
and then, by continuity, we can obtain
(31)VDFE′=(ω1,…,ωm)[M0I-I]>0,
in a neighborhood of P0 in Γ0; then P0 is unstable.
AssumeR0=ρ(M0)>1. By the uniform persistence result from [19] and a similar argument as in the proof of [2], the instability of P0 implies the uniform persistence of (15). This together with the dissipativity of (15) resulted from the forward invariant and compact property of Γ0 stated previously, implies which that  (15) has an equilibrium in Γ0, denoted by P¯* (see, e.g., Theorem D.3 in [20]).In what follows we prove that the endemic equilibriumP¯* of system (15) is globally asymptotically stable when R0>1.Throughout the paper, we denote(32)H(z)=z-1-lnz.
Then H(z)≥0 for z>0 and has global minimum at z=1.For convenience of notations, setβ-ij=βijfij(Si*,Ij*),1≤i,j≤m, and
(33)B-=[∑l≠1‍β-1l-β-21⋯-β-ml-β-12∑l≠2‍β-2l⋯-β-m2⋮⋮⋱⋮-β-1m-β-2m⋯∑l≠m‍β-ml].
Then, B- is also irreducible. One knows that the solution space of the linear system
(34)B-v=0
has dimension 1 and
(35)(v1,v2,…vm)=(C11,…,Cmm)
gives a base of this space, where Ckk>0,k=1,2,…,m, is the cofactor of the kth diagonal entry of B-. To get the global stability of P¯*, the following assumptions in [15] are proposed: (A4)
:[φi(Si)-φi(Si*)](Si-Si*)<0 for Si≠Si*, Si∈[0,Si0],(A5)
: ForSi≠Si*, [φi(Si)-φi(Si*)]·[fii(Si,Ii*)-fii(Si*,Ii*)]<0.(A6)
: ForSi,Ij>0,
(36)(fii(Si*,Ii*)fij(Si,Ij)-fii(Si,Ii*)fij(Si*,Ij*))·(fii(Si*,Ii*)fij(Si,Ij)Ij-fii(Si,Ii*)fij(Si*,Ij*)Ij*)≤0.A difficult mathematical question for system (15) is that of whether the endemic equilibrium P¯* is unique when R0>1 and whether P¯* is globally asymptotically stable when it is unique. Our main global stability result is given.Theorem 4.
Consider system (15). Assume that (A4)–(A6) hold and the matrix B=(βij)m×m is irreducible. If R0>1, then there is a unique endemic equilibrium P¯* for system (15), and P¯* is globally asymptotically stable in Γ0.Proof.
We show thatP¯* is globally asymptotically stable in Γ0, which implies that there exists a unique endemic equilibrium.
Consider a Lyapunov function as(37)VEE=Si-fii(Si*,Ii*)∫Si*Si‍dξfii(ξ,Ii*)+(1+δb)n[∑j=1n‍(yi,j-yi,j*-yi,j*lnyi,jyi,j*)+Ii-Ii*-Ii*lnIiIi*{∑j=1n‍}].
This function has a linear partVEE expressed by
(38)LEE=Si+(1+δb)n[∑j=1n‍(yi,j-yi,j*)+Ii-Ii*].
First, calculating the derivatives of LEE, we obtain
(39)LEE′=φi(Si)-(1+δb)n(δ+εi+γi)Ii.
Calculating the time derivative ofVEE along the solutions of system (15) and using equilibrium equation (21), we have
(40)VEE′=LEE′-fii(Si*,Ii*)fii(Si,Ii*)S˙i+(1+δb)n[∑j=1n‍yi,j*yi,jy˙i,j+Ii*IiI˙i]=φi(Si)-(1+δb)n(δ+εi+γi)Ii-{(φi(Si)-∑j=1m‍βijfij(Si,Ij))fii(Si*,Ii*)fii(Si,Ii*)+(1+δb)n[∑j=1m‍βijfij(Si,Ij)yi,1*(1+δb)nyi,1-yi,1*b^+1b^∑k=2n‍yi,k*(yi,k-1yi,k-1)+yi,nIi*b^Ii-(δ+εi+γi)Ii*{(∑j=1m‍)}]}=φi(Si)(1-fii(Si*,Ii*)fii(Si,Ii*))-∑j=1m‍βijfij(Si*,Ij*)yi,1*fij(Si,Ij)yi,1fij(Si*,Ij*)-(1+δb)nb^∑k=2n‍yi,k*yi,k-1yi,k+(1+δb)nb^nyi*-(1+δb)nb^yi,n*yi,nIi*yi,n*Ii+(1+δb)n(δ+εi+γi)Ii*+∑j=1m‍βijfij(Si*,Ij*)fij(Si,Ij)fii(Si*,Ii*)fij(Si*,Ij*)fii(Si,Ii*)-(1+δb)n(δ+εi+γi)Ii=(φi(Si)-φi(Si*))(1-fii(Si*,Ii*)fii(Si,Ii*))+∑j=1m‍β-ij{{fij(Si,Ij)fii(Si*,Ii*)fij(Si*,Ij*)fii(Si,Ii*)}n+2-fii(Si*,Ii*)fii(Si,Ii*)-∑k=2n‍yi,k*yi,k-1yi,kyi,k-1*-yi,nIi*yi,n*Ii-IiIi*-fij(Si,Ij)yi,1*fij(Si*,Ij*)yi,1+fij(Si,Ij)fii(Si*,Ii*)fij(Si*,Ij*)fii(Si,Ii*){∑k=2n‍}}.
It follows from the assumptions (A4)-(A5) that VEE′ can be estimated by
(41)VEE′≤∑i,j=1m‍β-ij{{[1-Ijfii(Si,Ii*)fij(Si*,Ij*)Ij*fii(Si*,Ii*)fij(Si,Ij)]}Gi(Ii)-Gj(Ij)+H(fii(Si*,Ii*)fii(Si,Ii*))+H(fij(Si,Ij)yi,1*fij(Si*,Ij*)yi,1)+∑k=2n‍H(yi,k*yi,k-1yi,kyi,k-1*)+H(yi,nIi*yi,n*Ii)+H(Ijfii(Si,Ii*)fij(Si*,Ij*)Ij*fii(Si*,Ii*)fij(Si,Ij))+[fii(Si*,Ii*)fij(Si,Ij)fii(Si,Ii*)fij(Si*,Ij*)-1]·[1-Ijfii(Si,Ii*)fij(Si*,Ij*)Ij*fii(Si*,Ii*)fij(Si,Ij)]}.
From the assumption (A6) and (32), we know that
(42)VEE′≤∑i,j=1m‍β-ij{Gi(Ii)-Gj(Ij)},
where Gi(Ii)=-Ii/Ii*+ln(Ii/Ii*).
Obviously, the equalities in (41) and (42) hold if and only if
(43)fii(Si*,Ii*)fii(Si,Ii*)=1,(1-fii(Si*,Ii*)fii(Si,Ii*))[φi(Si)-φi(Si*)]=0,[fij(Si,Ij)fii(Si*,Ii*)fij(Si*,Ij*)fii(Si,Ii*)-1]×[1-fii(Si,Ii*)fij(Si*,Ij*)Ijfii(Si,Ii*)fij(Si,Ij)Ij*]=0.
That is, Si=Si*,Ii=Ii*,i=1,2,…,m. We can show that VEE and β-ij satisfy the assumptions of Theorem 3.1 and Corollary  3.3 in [21]. Therefore, the function
(44)L=∑i=1n‍viVEE
is a Lyapunov function for system (15); namely, L′|(15)≤0 for P¯*∈Γ0. One can only show that the largest invariant subset, where L′|(15)=0, is the singleton {P¯*} by the same argument as in [2–5, 13, 21]. By LaSalle’s invariance principle, P¯* is globally asymptotically stable in Γ0. This completes the proof of Theorem 4.Remark 5.
We show a complete proof for global asymptotic stability of unique endemic equilibrium of system (15). In the case of fij(Si,Ij)=SiIj, system (15) will reduce to the system studied in [14, 22]. Here Theorem 4 extends related results in [14, 22] to a result to a more general case allowing a nonlinear incidence rate. Our result also cover the related results of single group model in [13] for the case of f(S,I)=f(S)I.


## 4. Numerical Example

Consider the system (15) when m=2, n=2, φi(Si(t))=3-Si, and fij(Si,Ij)=SiIj, i,j=1,2. One then has a two-group model as follows:
(45)S1′(t)=3-S1-[β11S1(t)I1(t)+β12S1(t)I2(t)],y1,1′(t)=1(1+δb)n[β11S1(t)I1(t)+β12S1(t)I2(t)]-1b^y1,1(t),y1,2′(t)=1b^(y1,1(t)-y1,2(t)),I1′(t)=1b^y1,2(t)-(δ+ε1+γ1)I1(t),S2′(t)=3-S2-[β21S2(t)I1(t)+β22S2(t)I2(t)],y2,1′(t)=1(1+δb)n[β21S2(t)I1(t)+β22S2(t)I2(t)]-1b^y2,1(t),y2,2′(t)=1b^(y2,1(t)-y2,2(t)),I2′(t)=1b^y2,2(t)-(δ+ε2+γ2)I2(t).
If we choose parameters as β11=5/24,β12=1,β21=1/36,β22=1/2,δ=0.8,ε1=2, ε2=2, γ1=1/4, and γ2=1/4, we can compute R0=0.051<1, and hence P0=(3,0,0,0,3,0,0,0) is the unique equilibrium of system (45) and it is globally stable from Theorem 4 (see Figures 1 and 2).Trajectories ofS1(t),I1(t),S2(t), and I2(t) for R0=0.051<1, and P0=(3,0,0,0,3,0,0,0) is globally stable. S1(t),S2(t),I1(t), and I2(t) versus t are illustrated by (a), (b), (c), and (d). Initial values are S1(0)=9,S2(0)=1,y1,1(0)=2,y1,2(0)=2,y2,1(0)=0,y2,2(0)=0,I1(0)=6, and I2(0)=2.
(a)
(b)
(c)
(d)Numerical simulation of (45) with R0=0.051<1; hence P0=(3,0,0,0,3,0,0,0) is globally stable. Graphs (a) and (b) illustrate that S1(t),y1,1(t),y1,2(t) and I1(t) will eventually towards to steady state. Graphs (c) and (d) illustrate that S2(t),y2,1(t),y2,2(t), and I2(t) will eventually towards to steady state. Initial values are S1(0)=9,S2(0)=1,y1,1(0)=2,y1,2(0)=2,y2,1(0)=0,y2,2(0)=0,I1(0)=6, and I2(0)=2.
(a)
(b)
(c)
(d)On the other hand, ifβij are chosen as β11=0.7,β12=1,β21=0.8,β22=1,δ=0.5,ε1=0.02, ε2=0.03, γ1=0.05, and γ2=0,05, we can compute R0=1.67355>1, and hence P¯*= (0.347644, 0.0760948, 0.0760948, 4.51674, 0.330353, 0.0765909, 0.0765909, 4.4678) is the unique equilibrium of system (45) and it is globally stable from Theorem 4 (see Figures 3 and 4).Trajectories ofS1(t),I1(t),S2(t), and I2(t) for R0=1.67355>1, and P¯*=  (0.347644, 0.0760948, 0.0760948, 4.51674, 0.330353, 0.0765909, 0.0765909, 4.4678) is globally stable. S1(t),S2(t),I1(t), and I2(t) versus t are illustrated by (a), (b), (c), and (d). Initial values are S1(0)=6,S2(0)=2,y1,1(0)=3,y1,2(0)=3,y2,1(0)=0.1,y2,2(0)=0.1,I1(0)=1.5, and I2(0)=0.5.
(a)
(b)
(c)
(d)Numerical simulation of (45) with R0=1.67355>1; hence P*= (0.347644,  0.0760948,  0.0760948,  4.51674,  0.330353,  0.0765909,  0.0765909,  4.4678) is globally stable. Graphs (a) and (b) illustrate that S1(t),y1,1(t),y1,2(t), and I1(t) will eventually towards to steady state. Graphs (c) and (d) illustrate that S2(t),y2,1(t),y2,2(t), and I2(t) will eventually towards to steady state. Initial values are S1(0)=6,S2(0)=2,y1,1(0)=3,y1,2(0)=3,y2,1(0)=0.1,y2,2(0)=0.1,I1(0)=1.5, and I2(0)=0.5.
(a)
(b)
(c)
(d)



---

*Source: 354287-2013-09-15.xml*