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

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 number 𝑅₀. 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.

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 into 𝑚 distinct groups (𝑚≥1). For 1≤𝑖≤𝑚, the 𝑖th group is further partitioned into four disjoint classes: the susceptible individuals, exposed individuals, infectious individuals, and recovered individuals, whose numbers of individuals at time 𝑡 are denoted by 𝑆_𝑖(𝑡),𝐸_𝑖(𝑡),𝐼_𝑖(𝑡), and 𝑅_𝑖(𝑡), 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 𝑃(𝑡) denote the probability that an exposed individual remains in the time 𝑡 after entering the exposed class. For 1≤𝑖,𝑗≤𝑚, 𝛽_𝑖𝑗≥0 denotes the coefficient of transmission between compartments 𝑆_𝑖 and 𝐼_𝑗. It is assumed that 𝑚-square matrix (𝛽_𝑖𝑗)_{1≤𝑖,𝑗≤𝑚} is irreducible [7]. So the proportion of exposed individuals can be expressed by the integral 𝐸_𝑖(𝑡)=^𝑚_𝑗=1𝛽_𝑖𝑗^𝑡₀𝑓_𝑖𝑗𝑆_𝑖(𝑢),𝐼_𝑗(𝑢)𝑒^{−𝛿𝑗(𝑡−𝑢)}𝑃_𝑗(𝑡−𝑢)𝑑𝑢, (1) where the sum takes into account cross-infections from all groups. Integrals in (1) are in the Riemann-Stieltjes sense. 𝑃_𝑗(𝑡) satisfies the following reasonable properties: (𝐴)𝑃_𝑗∶[0,∞)→[0,1] is nonincreasing, piecewise continuous with possibly finitely many jumps and satisfies 𝑃_𝑗(0⁺)=1, and lim_𝑡→∞𝑃_𝑗(𝑡)=0 with ∫^∞₀𝑃_𝑗(𝑡)𝑑𝑡 is positive and finite.

Differentiating (1) gives 𝐸^′_𝑖(𝑡)=^𝑚_𝑗=1𝛽_𝑖𝑗𝑓_𝑖𝑗𝑆_𝑖(𝑡),𝐼_𝑗(𝑡+_𝑚𝑗=1)𝛽_𝑖𝑗^𝑡₀𝑓_𝑖𝑗𝑆_𝑖(𝑢),𝐼_𝑗(𝑢)𝑒^{−𝛿𝑗(𝑡−𝑢)}×𝑃^′_𝑗(𝑡−𝑢)𝑑𝑢−𝛿_𝑖𝐸_𝑖(𝑡). (2) 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 𝑃^′_𝑗(𝑡−𝑢)≤0 due to the aformentioned property) from the exposed class; hence 𝐼^′_𝑖(𝑡)=−^𝑚_𝑗=1𝛽_𝑖𝑗^𝑡₀𝑓_𝑖𝑗𝑆_𝑖(𝑢),𝐼_𝑗(𝑢)𝑒^{−𝛿𝑗(𝑡−𝑢)}×𝑃^′_𝑗(𝑡−𝑢)𝑑𝑢−𝛿_𝑖+𝜀_𝑖+𝛾_𝑖𝐼_𝑖(𝑡). (3) Let ℎ_𝑗(𝑡)=−𝑃^′_𝑗(𝑡) be the probability density function for the time (a random variable) it takes for an infected individual in the 𝑖th group to become infectious. Then (4) becomes 𝐼^′_𝑖(𝑡)=^𝑚_𝑗=1𝛽_𝑖𝑗^𝑡₀𝑓_𝑖𝑗𝑆_𝑖(𝑢),𝐼_𝑗(𝑢)𝑒^{−𝛿𝑗(𝑡−𝑢)}×ℎ_𝑗(𝑡−𝑢)𝑑𝑢−𝛿_𝑖+𝜀_𝑖+𝛾_𝑖𝐼_𝑖(𝑡). (4)

Within the 𝑖th group, 𝜑_𝑖(𝑆_𝑖) denotes the growth rate of 𝑆_𝑖, 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: 𝑆^′_𝑖(𝑡)=𝜑_𝑖𝑆_𝑖(𝑡)−^𝑚_𝑗=1𝛽_𝑖𝑗𝑓_𝑖𝑗𝑆_𝑖(𝑡),𝐼_𝑗(𝑡,𝐸_′𝑖)(𝑡)=^𝑚_𝑗=1𝛽_𝑖𝑗𝑓_𝑖𝑗𝑆_𝑖(𝑡),𝐼_𝑗(𝑡)−^𝑚_𝑗=1𝛽_𝑖𝑗^𝑡₀𝑓_𝑖𝑗𝑆_𝑖(𝑢),𝐼_𝑗(𝑢)𝑒^{−𝛿𝑗(𝑡−𝑢)}×ℎ_𝑗(𝑡−𝑢)𝑑𝑢−𝛿_𝑖𝐸_𝑖(𝑡),𝐼^′_𝑖(𝑡)=^𝑚_𝑗=1𝛽_𝑖𝑗^𝑡₀𝑓_𝑖𝑗𝑆_𝑖(𝑢),𝐼_𝑗(𝑢)𝑒^{−𝛿𝑗(𝑡−𝑢)}×ℎ_𝑗(𝑡−𝑢)𝑑𝑢−𝛿_𝑖+𝜀_𝑖+𝛾_𝑖𝐼_𝑖(𝑡),𝑅^′_𝑖(𝑡)=𝛾_𝑖𝐼_𝑖(𝑡)−𝛿_𝑖𝑅_𝑖(𝑡). (5)

Since the variables 𝐸_𝑖 and 𝑅_𝑖 do not appear in the first and third equations of model (5), 𝐸_𝑖(𝑡) and 𝑅_𝑖(𝑡), 𝑖=1,…,𝑚, can be decoupled from the 𝑆_𝑖(𝑡) and 𝐼_𝑖(𝑡) equations; we only need to consider the subsystem of (5) consisting of only the 𝑆_𝑖 and 𝐼_𝑖 equations: 𝑆^′_𝑖(𝑡)=𝜑_𝑖𝑆_𝑖(𝑡)−^𝑚_𝑗=1𝛽_𝑖𝑗𝑓_𝑖𝑗𝑆_𝑖(𝑡),𝐼_𝑗(𝑡,𝐼_′𝑖)(𝑡)=^𝑚_𝑗=1𝛽_𝑖𝑗^𝑡₀𝑓_𝑖𝑗𝑆_𝑖(𝑢),𝐼_𝑗(𝑢)𝑒^{−𝛿𝑗(𝑡−𝑢)}×ℎ_𝑗(𝑡−𝑢)𝑑𝑢−𝛿_𝑖+𝜀_𝑖+𝛾_𝑖𝐼_𝑖(𝑡), (6) where 𝛿_𝑖 denotes the natural death rates of 𝐼_𝑖 compartments in the 𝑖th group, 𝜀_𝑖 is the death rate caused by disease in the 𝑖th group, and 𝛾_𝑖 is the rate of recovery of infectious individuals in the 𝑖th group. In what follows we investigate the global stability of system (5).

When 𝑚=1, 𝑃(𝑡)=𝑒^𝜖𝑡, 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 𝑃(𝑡) 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: 𝛿_𝑗=𝛿 for 𝑗=1,2,…,𝑚. Further, we assume that the functions ℎ_𝑗(𝑢) are disease specific only, implying that ℎ_𝑗(𝑢)=ℎ(𝑢) for 𝑗=1,2,…,𝑚. We choose the gamma distribution: ℎ(𝑢)=ℎ_𝑛,𝑏(𝑢)=𝑢^𝑛−1(𝑛−1)!𝑏^𝑛𝑒^−𝑢/𝑏, (7) where 𝑏>0 is a real number and 𝑛>1 is an integer, which is widely used and can approximate several frequently used distributions. For example, when 𝑏→0⁺, ℎ_𝑛,𝑏(𝑠) will approach the Dirac delta function, and when 𝑛=1, ℎ_𝑛,𝑏(𝑠) 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 𝑅₀ and show that 𝑅₀ completely determines the global dynamics of system (6). More specifically, if 𝑅₀≤1, the disease-free equilibrium is globally asymptotically stable and the disease dies out; if 𝑅₀>1, a unique endemic equilibrium exists and is globally asymptotically stable, and the disease persists at the endemic equilibrium. The global stability of 𝑃^∗ 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 Section 2, we first give the model, preliminaries and the basic reproduction number 𝑅₀. The global stability of the corresponding equilibria for 𝑅₀≤1 and 𝑅₀>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 the 𝑖th group 𝜑_𝑖(𝑆_𝑖) and the transmission functions 𝑓_𝑖𝑗(𝑆_𝑖,𝐼_𝑗).(𝐴₁)𝜑_𝑖 are 𝐶¹ non-increasing functions on [0,∞) with 𝜑_𝑖(0)>0, and there is a unique positive solution 𝜉=𝑆⁰_𝑖 for the equation 𝜑_𝑖(𝜉)=0. 𝜑_𝑖(𝑆)>0 for 0≤𝑆<𝑆⁰_𝑖, and 𝜑_𝑖(𝑆)<0 for 𝑆>𝑆⁰_𝑖; that is 𝜑_𝑖𝑆_𝑖−𝜑_𝑖𝑆⁰_𝑖𝑆_𝑖−𝑆⁰_𝑖<0,for𝑆_𝑖≠𝑆⁰_𝑖,𝑖=1,2,…,𝑚. (8)(𝐴₂)𝑓_𝑖𝑗(𝑆_𝑖,𝐼_𝑗)≤𝑐_𝑖𝑗(𝑆_𝑖)𝐼_𝑗 for all 𝐼_𝑗>0. (𝐴₃)𝑐_𝑖𝑗(𝑆_𝑖)≤𝑐_𝑖𝑗(𝑆⁰_𝑖), 0<𝑆_𝑖<𝑆⁰_𝑖, 𝑖,𝑗=1,…,𝑚.

Following the technique and method in [14], define 𝑏̂≡𝑏,1+𝛿𝑏 (9) which can absorb the exponential term 𝑒^−𝛿𝑢 into the delay kernel. The second equation in (6) can be rewritten as 𝐼^′_𝑖(𝑡)=^𝑚_𝑗=1𝛽_𝑖𝑗(1+𝛿𝑏)^𝑛𝑡₀𝑓_𝑖𝑗𝑆_𝑖,𝐼_𝑗ℎ_𝑛,̂𝑏()−𝑡−𝑢𝑑𝑢𝛿+𝜀_𝑖+𝛾_𝑖𝐼_𝑖. (10) For 𝑙=1,…,𝑛, let𝑦_𝑖,𝑙(𝑡)=^𝑚_𝑗=1𝛽_𝑖𝑗𝑏̂(1+𝛿𝑏)^𝑛𝑡₀𝑓_𝑖𝑗𝑆_𝑖,𝐼_𝑗ℎ_𝑙,̂𝑏(𝑡−𝑢)𝑑𝑢,𝑖=1,2,…,𝑚. (11) Thus, for 𝑙∈{2,…,𝑛}, we obtain 𝑦̇_𝑖,𝑙=ℎ_𝑙,̂𝑏(0)^𝑚_𝑗=1𝛽_𝑖𝑗𝑏̂(1+𝛿𝑏)^𝑛𝑓_𝑖𝑗𝑆_𝑖,𝐼_𝑗+^𝑚_𝑗=1𝛽_𝑖𝑗𝑏̂(1+𝛿𝑏)^𝑛𝑡_−∞(𝑙−1)(𝑡−𝑢)^𝑙−2(𝑙−1)!𝑏̂^𝑙𝑒^{−(𝑡−𝑢)/}𝑏̂𝑓_𝑖𝑗𝑆_𝑖,𝐼_𝑗−𝑑𝑢_𝑚𝑗=1𝛽_𝑖𝑗𝑏̂(1+𝛿𝑏)^𝑛𝑡_−∞(𝑡−𝑢)^𝑙−1(𝑙−1)!𝑏̂^𝑙+1𝑒^{−(𝑡−𝑢)/}𝑏̂𝑓_𝑖𝑗𝑆_𝑖,𝐼_𝑗=𝑦𝑑𝑢_{𝑖,𝑙−1}−𝑦_𝑖,𝑙𝑏̂. (12) For 𝑙=1, we have 𝑦_𝑖,1=^𝑚_𝑗=1𝛽_𝑖𝑗𝑏̂(1+𝛿𝑏)^𝑛𝑡_−∞𝑒^{−(𝑡−𝑢)/}𝑏̂𝑏̂𝑓_𝑖𝑗𝑆_𝑖,𝐼_𝑗𝑑𝑢,𝑖=1,…,𝑚. (13) It follows that 𝑦̇_𝑖,1=^𝑚_𝑗=1𝛽_𝑖𝑗(1+𝛿𝑏)^𝑛𝑓_𝑖𝑗𝑆_𝑖,𝐼_𝑗−^𝑚_𝑗=1𝛽_𝑖𝑗(1+𝛿𝑏)^𝑛𝑡_−∞𝑒^{−(𝑡−𝑢)/}𝑏̂𝑏̂𝑓_𝑖𝑗𝑆_𝑖,𝐼_𝑗=𝑑𝑢_𝑚𝑗=1𝛽_𝑖𝑗(1+𝛿𝑏)^𝑛𝑓_𝑖𝑗𝑆_𝑖,𝐼_𝑗−1𝑏̂𝑦_𝑖,1,𝑖=1,…,𝑚. (14)

Thus the integro-differential system (6) is equivalent to the ordinary differential equations 𝑆^′_𝑖(𝑡)=𝜑_𝑖𝑆_𝑖(𝑡)−^𝑚_𝑗=1𝛽_𝑖𝑗𝑓_𝑖𝑗𝑆_𝑖(𝑡),𝐼_𝑗(𝑡,𝑦_′𝑖,1)(𝑡)=1(1+𝛿𝑏)^𝑛𝑚_𝑗=1𝛽_𝑖𝑗𝑓_𝑖𝑗𝑆_𝑖(𝑡),𝐼_𝑗(𝑡)−1𝑏̂𝑦_𝑖,1(𝑡),𝑦^′_𝑖,2(𝑡)=1𝑏̂𝑦_𝑖,1(𝑡)−𝑦_𝑖,2(𝑡)⋮𝑦,𝑖=1,2,…,𝑚,_′𝑖,𝑛(𝑡)=1𝑏̂𝑦_{𝑖,𝑛−1}(𝑡)−𝑦_𝑖,𝑛(𝑡),𝐼^′_𝑖(𝑡)=1𝑏̂𝑦_𝑖,𝑛(𝑡)−𝛿+𝜀_𝑖+𝛾_𝑖𝐼_i(𝑡). (15)

For initial condition 𝑆₁(0),𝑦_1,1(0),…,𝑦_1,𝑛(0),𝐼₁(0),𝑆₂(0),𝑦_2,1(0),…,𝑦_1,𝑛(0),𝐼₂(0)𝑆,…,_𝑚(0),𝑦_𝑚,1(0),…,𝑦_𝑚,𝑛(0),𝐼_𝑚(0)∈R^𝑚(𝑛+2), (16) the existence, uniqueness, and continuity of the solution (𝑆_𝑖,𝑦_𝑖,1,𝑦_𝑖,2,…,𝑦_𝑖,𝑛,𝐼_𝑖) 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 (𝐴₁) and the first equation of (15) that limsup_𝑡→∞𝑆_𝑖(𝑡)≤𝑆⁰_𝑖 for all 𝑖=1,2,…,𝑚. Let 𝑁_𝜑𝑖 be the maximum of the function 𝜑_𝑖 on R₊ and let 𝑞 be a positive real number such that 𝑞>𝑏̂𝑁_𝜑𝑖. Denote by Υ_𝑖 the 𝑖th tube for system (15); that is, Υ_𝑖=𝑆_𝑖,𝑦_𝑖,1,𝑦_𝑖,2,…,𝑦_𝑖,𝑛,𝐼_𝑖. (17) It follows from a similar argument to that in [14] that we can show that the set 𝐷_𝜖 defined by Γ_𝜖=𝑆_𝑖,𝑦_𝑖,1,𝑦_𝑖,2,…,𝑦_𝑖,𝑛,𝐼_𝑖∈R₊^𝑚(𝑛+2)∣𝑆_𝑖≤𝑆⁰_𝑖+𝜖,𝑆_𝑖+(1+𝛿𝑏)^𝑛𝑦_𝑖,1≤𝑞+𝑆⁰_𝑖,𝑦_𝑖,𝑙≤𝑞+𝑆⁰_𝑖+𝑙𝜖(1+𝛿𝑏)^𝑛,𝐼_𝑖≤𝑞+𝑆⁰_𝑖+(𝑛+1)𝜖𝑏̂(1+𝛿𝑏)^𝑛𝛿+𝜖_𝑖+𝛾_𝑖,𝑖=1,2,…,𝑚,𝑙=2,3,…,𝑛 (18) is a forward invariant compact absorbing set for system (15) for 𝜖>0 and that the set Γ₀ (i.e., when 𝜖=0) is a forward invariant compact set.

Under the assumption (𝐴₁), we know that system (15) always has the disease-free equilibrium 𝑃₀=𝑆⁰₁,0,…,0,𝐼⁰₁,𝑆⁰₂,0,…,0,𝐼⁰₂,…,𝑆⁰_𝑚,0,…,0,𝐼⁰_𝑚∈R^𝑚(𝑛+2). (19) An equilibrium 𝑃^∗ of (6) has the form 𝑃^∗=(𝑆^∗₁,𝐼^∗₁,𝑆^∗₂,𝐼^∗₂,…,𝑆^∗_𝑚,𝐼^∗_𝑚)∈R^2𝑚 with 𝑆^∗_𝑖>0, 𝐼^∗_𝑖>0, 𝑖=1,…,𝑚. Translating to the equivalent system (15), 𝑃^∗ is corresponding to 𝑃^∗=𝑆^∗₁,𝑦^∗_1,1,…,𝑦^∗_1,𝑛,𝐼^∗₁,𝑆^∗₂,𝑦^∗_2,1,…,𝑦^∗_2,𝑛,𝐼^∗₂𝑆,…,_∗𝑚,𝑦^∗_𝑚,1,…,𝑦^∗_𝑚,𝑛,𝐼^∗_𝑚∈R^𝑚(𝑛+2). (20)𝑃^∗ in the interior of Γ₀ is called an endemic equilibrium, and it satisfies the following equilibrium equations: 0=𝜑_𝑖𝑆^∗_𝑖−^𝑚_𝑗=1𝛽_𝑖𝑗𝑓_𝑖𝑗𝑆^∗_𝑖,𝐼^∗_𝑗,10=(1+𝛿𝑏)^𝑛𝑚_𝑗=1𝛽_𝑖𝑗𝑓_𝑖𝑗𝑆^∗_𝑖,𝐼^∗_𝑗−1𝑏̂𝑦^∗_𝑖,1,10=𝑏̂𝑦^∗_𝑖,1−𝑦^∗_𝑖,2,⋮10=𝑏̂𝑦^∗_{𝑖,𝑛−1}−𝑦^∗_𝑖,𝑛,10=𝑏̂𝑦^∗_𝑖,𝑛−𝛿+𝜀_𝑖+𝛾_𝑖𝐼^∗_𝑖. (21) The basic reproduction number 𝑅₀ 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 𝑐ℱ=¹¹𝑆⁰₁𝛽₁₁(1+𝛿𝑏)^𝑛⋯𝑐_1𝑚𝑆⁰₁𝛽_1𝑚(1+𝛿𝑏)^𝑛𝑐⋮⋱⋮_𝑚1𝑆⁰_𝑚𝛽_𝑚1(1+𝛿𝑏)^𝑛⋯𝑐_𝑚𝑚𝑆⁰_𝑚𝛽_𝑚𝑚(1+𝛿𝑏)^𝑛,𝒱=diag𝛿+𝜀_𝑖+𝛾_𝑖=𝛿+𝜀₁+𝛾₁0⋯00𝛿+𝜀₂+𝛾₂⋯0⋮⋮⋱⋮00⋯𝛿+𝜀_𝑚+𝛾_𝑚. (22) Then the next generation matrix is ℱ𝒱⁻¹=𝑐₁₁𝑆⁰₁𝛽₁₁(1+𝛿𝑏)^𝑛𝛿+𝜀₁+𝛾₁⋯𝑐_1𝑚𝑆⁰₁𝛽_1𝑚(1+𝛿𝑏)^𝑛𝛿+𝜀_𝑚+𝛾_𝑚𝑐⋮⋱⋮_𝑚1𝑆⁰₁𝛽_𝑚1(1+𝛿𝑏)^𝑛𝛿+𝜀₁+𝛾₁⋯𝑐_𝑚𝑚𝑆⁰₁𝛽_𝑚𝑚(1+𝛿𝑏)^𝑛𝛿+𝜀_𝑚+𝛾_𝑚, (23) and hence, the basic reproduction number 𝑅₀ is 𝑅₀=𝜌ℱ𝒱⁻¹=max|𝜆|;𝜆∈𝜎ℱ𝒱⁻¹, (24) 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 (𝐴₁)--(𝐴₅) of Theorem 2 of [17], we have the following proposition.

Lemma 1. For system (15), the disease-free equilibrium 𝑃₀ is locally asymptotically stable if 𝑅₀<1, while it is unstable if 𝑅₀>1.

Following the method of [2], one defines a matrix 𝑀⁰=𝒱⁻¹ℱ=𝑐₁₁𝑆⁰₁𝛽₁₁(1+𝛿𝑏)^𝑛𝛿+𝜀₁+𝛾₁⋯𝑐_1𝑚𝑆⁰₁𝛽_1𝑚(1+𝛿𝑏)^𝑛𝛿+𝜀₁+𝛾₁𝑐⋮⋱⋮_𝑚1𝑆⁰₁𝛽_𝑚1(1+𝛿𝑏)^𝑛𝛿+𝜀_𝑚+𝛾_𝑚⋯𝑐_𝑚𝑚𝑆⁰₁𝛽_𝑚𝑚(1+𝛿𝑏)^𝑛𝛿+𝜀_𝑚+𝛾_𝑚, (25) whose spectral radius has a similar threshold property to that of 𝑅₀, since both of the nonnegative matrices ℱ𝒱⁻¹ and 𝑀⁰ 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 𝑅₀=𝜌(ℱ𝒱⁻¹)=𝜌(𝒱⁻¹ℱ)=𝜌(𝑀⁰). Then the following lemma immediately follows.

Lemma 2. 𝜌(𝑀⁰)≤1 if and only if 𝑅₀≤1.

3. Main Results

The following main theorems are summarized in terms of system (15).

Theorem 3. Assume that the functions 𝜑_𝑖 and 𝑓_𝑖𝑗 satisfy assumptions (𝐴₁)--(𝐴₃), and the matrix 𝐵=(𝛽_𝑖𝑗)_𝑚×𝑚 is irreducible and 𝑅₀ is defined by (24). (i)If 𝑅₀≤1, then 𝑃₀ is the unique equilibrium of system (15), and 𝑃₀ is globally asymptotically stable in Γ₀. (ii)If 𝑅₀>1, then 𝑃₀ is unstable, and system (15) is uniformly persistent in Γ₀.

Proof. Let us define 𝑀(𝑆)=(𝛽_𝑖𝑗𝑐_𝑖𝑗(𝑆_𝑖)/(1+𝛿𝑏)^𝑛(𝛿+𝜀_𝑖+𝛾_𝑖))_𝑚×𝑚, where 𝑆=(𝑆₁,𝑆₂,…,𝑆_𝑚)^𝑇. Note that 𝑀(𝑆₀)=𝑀⁰. Since 𝐵=(𝛽_𝑖𝑗)_𝑚×𝑚 is irreducible, the matrix 𝑀⁰ is also irreducible.
First we claim that there does not exist any endemic equilibrium 𝑃^∗ in Ω. Suppose that 𝑆≠𝑆₀. Then we have 0<𝑀(𝑆)<𝑀⁰. Since nonnegative matrix 𝑀(𝑆)+𝑀⁰ is irreducible, it follows from the Perron-Frobenius theorem (see Corollary  2.1.5 of [7]) that 𝜌(𝑀(𝑆))<𝜌(𝑀⁰)≤1. This implies that 𝑀(𝑆)𝐼=𝐼 has only the trivial solution 𝐼=0, where 𝐼=(𝐼₁,…,𝐼_𝑚)^𝑇. Hence the claim is true. Next we claim that the disease-free equilibrium 𝑃₀ is globally asymptotically stable in Γ₀. From the Perron-Frobenius (see Theorem  2.1.4 of [7]), the nonnegative irreducible matrix 𝑀⁰ has a strictly positive left eigenvector (𝜔₁,𝜔₂,…,𝜔_𝑚) associated with the eigenvalue 𝜌(𝑀⁰) such that 𝜔₁,𝜔₂,…,𝜔_𝑚𝜌𝑀⁰=𝜔₁,𝜔₂,…,𝜔_𝑚𝑀⁰. (26) Using this positive eigenvector, we construct the following Lyapunov function: 𝑉_DFE=^𝑚_𝑖=1𝜔_𝑖𝛿+𝜀_𝑖+𝛾_𝑖^𝑛_𝑗=1𝑦_𝑖,𝑗+𝐼_𝑖. (27) Computing the derivative of 𝑉_DFE along the solutions of (15) in Γ₀, we get 𝑉^′_DFE=^𝑚_𝑖=1^𝑚_𝑗=1𝜔_𝑖𝛽_𝑖𝑗(1+𝛿𝑏)^𝑛𝛿+𝜀_𝑖+𝛾_𝑖𝑓_𝑖𝑗𝑆_𝑖,𝐼_𝑗−𝜔_𝑖𝐼_𝑖≤^𝑚_𝑖=1^𝑚_𝑗=1𝜔_𝑖𝛽_𝑖𝑗𝑐_𝑖𝑗𝑆_𝑖(1+𝛿𝑏)^𝑛𝛿+𝜀_𝑖+𝛾_𝑖𝐼_𝑗−𝜔_𝑖𝐼_𝑖≤^𝑚_𝑖=1^𝑚_𝑗=1𝜔_𝑖𝛽_𝑖𝑗𝑐_𝑖𝑗𝑆⁰_𝑖(1+𝛿𝑏)^𝑛𝛿+𝜀_𝑖+𝛾_𝑖𝐼_𝑗−𝜔_𝑖𝐼_𝑖=𝜔₁,…,𝜔_𝑚𝑀⁰=𝜌𝑀𝐼−𝐼₀−1𝜔₁,…,𝜔_𝑚𝐼. (28) Thus, under the assumption 𝑅₀=𝜌(𝑀⁰)<1, 𝑉^′_DFE≤0, and 𝑉^′_DFE=0 if and only if 𝐼=0 and 𝑆=𝑆⁰=(𝑆⁰₁,𝑆⁰₂,…,𝑆⁰_𝑚). Suppose that 𝜌(𝑀⁰)=1. Then it follows from the previous that 𝑉^′_DFE=0 implies 𝜔₁,…,𝜔_𝑚𝑀⁰𝐼=𝜔₁,…,𝜔_𝑚𝐼. (29) Hence, if 𝑆≠𝑆₀, then (𝜔₁,…,𝜔_𝑚)𝑀(𝑆)<(𝜔₁,…,𝜔_𝑚)𝑀⁰=𝜌(𝑀⁰)(𝜔₁,…,𝜔_𝑚)=(𝜔₁,…,𝜔_𝑚) and thus 𝐼=0 is the only solution of (29). Summarizing the statements, we see that 𝑉^′_DFE=0 if and only if 𝐼=0 or 𝑆=𝑆₀, which implies that the compact invariant subset of the set where 𝑉^′_DFE=0 is only the singleton 𝑃₀. Thus, by LaSalle’s invariance principle [18], it follows that the disease-free equilibrium 𝐸⁰ is globally asymptotically stable in Γ₀.
If 𝑅₀=𝜌(𝑀⁰)>1, then 𝜔₁,𝜔₂,…,𝜔_𝑚𝑀⁰−𝜔₁,𝜔₂,…,𝜔_𝑚=𝜌𝑀⁰−1𝜔₁,𝜔₂,…,𝜔_𝑚>0, (30) and then, by continuity, we can obtain 𝑉^′_DFE=𝜔₁,…,𝜔_𝑚𝑀⁰𝐼−𝐼>0, (31) in a neighborhood of 𝑃₀ in Γ₀; then 𝑃₀ is unstable.
Assume 𝑅₀=𝜌(𝑀⁰)>1. By the uniform persistence result from [19] and a similar argument as in the proof of [2], the instability of 𝑃₀ implies the uniform persistence of (15). This together with the dissipativity of (15) resulted from the forward invariant and compact property of Γ₀ stated previously, implies which that (15) has an equilibrium in Γ₀, denoted by 𝑃^∗ (see, e.g., Theorem D.3 in [20]).

In what follows we prove that the endemic equilibrium 𝑃^∗ of system (15) is globally asymptotically stable when 𝑅₀>1.

Throughout the paper, we denote 𝐻(𝑧)=𝑧−1−ln𝑧. (32) Then 𝐻(𝑧)≥0 for 𝑧>0 and has global minimum at 𝑧=1.

For convenience of notations, set 𝛽_𝑖𝑗=𝛽_𝑖𝑗𝑓_𝑖𝑗(𝑆^∗_𝑖,𝐼^∗_𝑗),1≤𝑖,𝑗≤𝑚, and 𝐵=^𝑙≠1𝛽_1𝑙−𝛽₂₁⋯−𝛽_𝑚𝑙−𝛽_12𝑙≠2𝛽_2𝑙⋯−𝛽_𝑚2−⋮⋮⋱⋮𝛽_1𝑚−𝛽_2𝑚⋯_𝑙≠𝑚𝛽_𝑚𝑙. (33) Then, 𝐵 is also irreducible. One knows that the solution space of the linear system 𝐵𝑣=0 (34) has dimension 1 and 𝑣₁,𝑣₂,…𝑣_𝑚=𝐶₁₁,…,𝐶_𝑚𝑚 (35) gives a base of this space, where 𝐶_𝑘𝑘>0,𝑘=1,2,…,𝑚, is the cofactor of the 𝑘th diagonal entry of 𝐵. To get the global stability of 𝑃^∗, the following assumptions in [15] are proposed: (𝐴₄): [𝜑_𝑖(𝑆_𝑖)−𝜑_𝑖(𝑆^∗_𝑖)](𝑆_𝑖−𝑆^∗_𝑖)<0 for 𝑆_𝑖≠𝑆^∗_𝑖, 𝑆_𝑖∈[0,𝑆⁰_𝑖], (𝐴₅): For 𝑆_𝑖≠𝑆^∗_𝑖, [𝜑_𝑖(𝑆_𝑖)−𝜑_𝑖(𝑆^∗_𝑖)]⋅[𝑓_𝑖𝑖(𝑆_𝑖,𝐼^∗_𝑖)−𝑓_𝑖𝑖(𝑆^∗_𝑖,𝐼^∗_𝑖)]<0. (𝐴₆): For 𝑆_𝑖,𝐼_𝑗>0, 𝑓_𝑖𝑖𝑆^∗_𝑖,𝐼^∗_𝑖𝑓_𝑖𝑗𝑆_𝑖,𝐼_𝑗−𝑓_𝑖𝑖𝑆_𝑖,𝐼^∗_𝑖𝑓_𝑖𝑗𝑆^∗_𝑖,𝐼^∗_𝑗⋅𝑓_𝑖𝑖𝑆^∗_𝑖,𝐼^∗_𝑖𝑓_𝑖𝑗𝑆_𝑖,𝐼_𝑗𝐼_𝑗−𝑓_𝑖𝑖𝑆_𝑖,𝐼^∗_𝑖𝑓_𝑖𝑗𝑆^∗_𝑖,𝐼^∗_𝑗𝐼^∗_𝑗≤0. (36)A difficult mathematical question for system (15) is that of whether the endemic equilibrium 𝑃^∗ is unique when 𝑅₀>1 and whether 𝑃^∗ is globally asymptotically stable when it is unique. Our main global stability result is given.

Theorem 4. Consider system (15). Assume that (𝐴₄)--(𝐴₆) hold and the matrix 𝐵=(𝛽_𝑖𝑗)_𝑚×𝑚 is irreducible. If 𝑅₀>1, then there is a unique endemic equilibrium 𝑃^∗ for system (15), and 𝑃^∗ is globally asymptotically stable in Γ₀.

Proof. We show that 𝑃^∗ is globally asymptotically stable in Γ₀, which implies that there exists a unique endemic equilibrium.
Consider a Lyapunov function as 𝑉_EE=𝑆_𝑖−𝑓_𝑖𝑖𝑆^∗_𝑖,𝐼^∗_𝑖^𝑆𝑖_𝑆∗𝑖𝑑𝜉𝑓_𝑖𝑖𝜉,𝐼^∗_𝑖+(1+𝛿𝑏)^𝑛𝑛_𝑗=1𝑦_𝑖,𝑗−𝑦^∗_𝑖,𝑗−𝑦^∗_𝑖,𝑗𝑦ln^𝑖,𝑗𝑦^∗_𝑖,𝑗+𝐼_𝑖−𝐼^∗_𝑖−𝐼^∗_𝑖𝐼ln^𝑖𝐼^∗_𝑖. (37)
This function has a linear part 𝑉_EE expressed by 𝐿_EE=𝑆_𝑖+(1+𝛿𝑏)^𝑛𝑛_𝑗=1𝑦_𝑖,𝑗−𝑦^∗_𝑖,𝑗+𝐼_𝑖−𝐼^∗_𝑖. (38) First, calculating the derivatives of 𝐿_EE, we obtain 𝐿^′_EE=𝜑_𝑖𝑆_𝑖−(1+𝛿𝑏)^𝑛𝛿+𝜀_𝑖+𝛾_𝑖𝐼_𝑖. (39)
Calculating the time derivative of 𝑉_EE along the solutions of system (15) and using equilibrium equation (21), we have 𝑉^′_EE=𝐿^′_EE−𝑓_𝑖𝑖𝑆^∗_𝑖,𝐼^∗_𝑖𝑓_𝑖𝑖𝑆_𝑖,𝐼^∗_𝑖𝑆̇_𝑖+(1+𝛿𝑏)^𝑛𝑛_𝑗=1𝑦^∗_𝑖,𝑗𝑦_𝑖,𝑗𝑦̇_𝑖,𝑗+𝐼^∗_𝑖𝐼_𝑖𝐼̇_𝑖=𝜑_𝑖𝑆_𝑖−(1+𝛿𝑏)^𝑛𝛿+𝜀_𝑖+𝛾_𝑖𝐼_𝑖−𝜑_𝑖𝑆_𝑖−^𝑚_𝑗=1𝛽_𝑖𝑗𝑓_𝑖𝑗𝑆_𝑖,𝐼_𝑗𝑓_𝑖𝑖𝑆^∗_𝑖,𝐼^∗_𝑖𝑓_𝑖𝑖𝑆_𝑖,𝐼^∗_𝑖+(1+𝛿𝑏)^𝑛𝑚_𝑗=1𝛽_𝑖𝑗𝑓_𝑖𝑗𝑆_𝑖,𝐼_𝑗𝑦^∗_𝑖,1(1+𝛿𝑏)^𝑛𝑦_𝑖,1−𝑦^∗_𝑖,1𝑏̂+1𝑏̂^𝑛_𝑘=2𝑦^∗_𝑖,𝑘𝑦_{𝑖,𝑘−1}𝑦_𝑖,𝑘+𝑦−1_𝑖,𝑛𝐼^∗_i𝑏̂𝐼_𝑖−𝛿+𝜀_𝑖+𝛾_𝑖𝐼^∗_𝑖=𝜑_𝑖𝑆_𝑖𝑓1−^𝑖𝑖𝑆^∗_𝑖,𝐼^∗_𝑖𝑓_𝑖𝑖𝑆_𝑖,𝐼^∗_𝑖−^𝑚_𝑗=1𝛽_𝑖𝑗𝑓_𝑖𝑗𝑆^∗_𝑖,𝐼^∗_𝑗𝑦^∗_𝑖,1𝑓_𝑖𝑗𝑆_𝑖,𝐼_𝑗𝑦_𝑖,1𝑓_𝑖𝑗𝑆^∗_𝑖,𝐼^∗_𝑗−(1+𝛿𝑏)^𝑛𝑏̂^𝑛_𝑘=2𝑦^∗_𝑖,𝑘𝑦_{𝑖,𝑘−1}𝑦_𝑖,𝑘+(1+𝛿𝑏)^𝑛𝑏̂𝑛𝑦^∗_𝑖−(1+𝛿𝑏)^𝑛𝑏̂𝑦^∗_𝑖,𝑛𝑦_𝑖,𝑛𝐼^∗_𝑖𝑦^∗_𝑖,𝑛𝐼_𝑖+(1+𝛿𝑏)^𝑛𝛿+𝜀_𝑖+𝛾_𝑖𝐼^∗_𝑖+^𝑚_𝑗=1𝛽_𝑖𝑗𝑓_𝑖𝑗𝑆^∗_𝑖,𝐼^∗_𝑗𝑓_𝑖𝑗𝑆_𝑖,𝐼_𝑗𝑓_𝑖𝑖𝑆^∗_𝑖,𝐼^∗_𝑖𝑓_𝑖𝑗𝑆^∗_𝑖,𝐼^∗_𝑗𝑓_𝑖𝑖𝑆_𝑖,𝐼^∗_𝑖−(1+𝛿𝑏)^𝑛𝛿+𝜀_𝑖+𝛾_𝑖𝐼_𝑖=𝜑_𝑖𝑆_𝑖−𝜑_𝑖𝑆^∗_𝑖𝑓1−^𝑖𝑖𝑆^∗_𝑖,𝐼^∗_𝑖𝑓_𝑖𝑖𝑆_𝑖,𝐼^∗_𝑖+^𝑚_𝑗=1𝛽_𝑖𝑗𝑓𝑛+2−^𝑖𝑖𝑆^∗_𝑖,𝐼^∗_𝑖𝑓_𝑖𝑖𝑆_𝑖,𝐼^∗_𝑖−^𝑛_𝑘=2𝑦^∗_𝑖,𝑘𝑦_{𝑖,𝑘−1}𝑦_𝑖,𝑘𝑦^∗_{𝑖,𝑘−1}−𝑦_𝑖,𝑛𝐼^∗_𝑖𝑦^∗_𝑖,𝑛𝐼_𝑖−𝐼_𝑖𝐼^∗_𝑖−𝑓_𝑖𝑗𝑆_𝑖,𝐼_𝑗𝑦^∗_𝑖,1𝑓_𝑖𝑗𝑆^∗_𝑖,𝐼^∗_𝑗𝑦_𝑖,1+𝑓_𝑖𝑗𝑆_𝑖,𝐼_𝑗𝑓_𝑖𝑖𝑆^∗_𝑖,𝐼^∗_𝑖𝑓_𝑖𝑗𝑆^∗_𝑖,𝐼^∗_𝑗𝑓_𝑖𝑖𝑆_𝑖,𝐼^∗_𝑖. (40) It follows from the assumptions (𝐴₄)-(𝐴₅) that 𝑉^′_EE can be estimated by 𝑉^′_EE≤^𝑚_𝑖,𝑗=1𝛽_𝑖𝑗𝐺_𝑖𝐼_𝑖−𝐺_𝑗𝐼_𝑗𝑓+𝐻^𝑖𝑖𝑆^∗_𝑖,𝐼^∗_𝑖𝑓_𝑖𝑖𝑆_𝑖,𝐼^∗_𝑖𝑓+𝐻^𝑖𝑗𝑆_𝑖,𝐼_𝑗𝑦^∗_𝑖,1𝑓_𝑖𝑗𝑆^∗_𝑖,𝐼^∗_𝑗𝑦_𝑖,1+^𝑛_𝑘=2𝐻𝑦^∗_𝑖,𝑘𝑦_{𝑖,𝑘−1}𝑦_𝑖,𝑘𝑦^∗_{𝑖,𝑘−1}𝑦+𝐻^𝑖,𝑛𝐼^∗_𝑖𝑦^∗_𝑖,𝑛𝐼_𝑖𝐼+𝐻^𝑗𝑓_𝑖𝑖𝑆_𝑖,𝐼^∗_𝑖𝑓_𝑖𝑗𝑆^∗_𝑖,𝐼^∗_𝑗𝐼^∗_𝑗𝑓_𝑖𝑖𝑆^∗_𝑖,𝐼^∗_𝑖𝑓_𝑖𝑗𝑆_𝑖,𝐼_𝑗+𝑓_𝑖𝑖𝑆^∗_𝑖,𝐼^∗_𝑖𝑓_𝑖𝑗𝑆_𝑖,𝐼_𝑗𝑓_𝑖𝑖𝑆_𝑖,𝐼^∗_𝑖𝑓_𝑖𝑗𝑆^∗_𝑖,𝐼^∗_𝑗⋅𝐼−11−^𝑗𝑓_𝑖𝑖𝑆_𝑖,𝐼^∗_𝑖𝑓_𝑖𝑗𝑆^∗_𝑖,𝐼^∗_𝑗𝐼^∗_𝑗𝑓_𝑖𝑖𝑆^∗_𝑖,𝐼^∗_𝑖𝑓_𝑖𝑗𝑆_𝑖,𝐼_𝑗. (41) From the assumption (𝐴₆) and (32), we know that 𝑉^′_EE≤^𝑚_𝑖,𝑗=1𝛽_𝑖𝑗𝐺_𝑖𝐼_𝑖−𝐺_𝑗𝐼_𝑗, (42) where 𝐺_𝑖(𝐼_𝑖)=−𝐼_𝑖/𝐼^∗_𝑖+ln(𝐼_𝑖/𝐼^∗_𝑖).
Obviously, the equalities in (41) and (42) hold if and only if 𝑓_𝑖𝑖𝑆^∗_𝑖,𝐼^∗_𝑖𝑓_𝑖𝑖𝑆_𝑖,𝐼^∗_𝑖𝑓=1,1−^𝑖𝑖𝑆^∗_𝑖,𝐼^∗_𝑖𝑓_𝑖𝑖𝑆_𝑖,𝐼^∗_𝑖𝜑_𝑖𝑆_𝑖−𝜑_𝑖𝑆^∗_𝑖𝑓=0,_𝑖𝑗𝑆_𝑖,𝐼_𝑗𝑓_𝑖𝑖𝑆^∗_𝑖,𝐼^∗_𝑖𝑓_𝑖𝑗𝑆^∗_𝑖,𝐼^∗_𝑗𝑓_𝑖𝑖𝑆_𝑖,𝐼^∗_𝑖×𝑓−11−^𝑖𝑖𝑆_𝑖,𝐼^∗_𝑖𝑓_𝑖𝑗𝑆^∗_𝑖,𝐼^∗_𝑗𝐼_𝑗𝑓_𝑖𝑖𝑆_𝑖,𝐼^∗_𝑖𝑓_𝑖𝑗𝑆_𝑖,𝐼_𝑗𝐼^∗_𝑗=0. (43) That is, 𝑆_𝑖=𝑆^∗_𝑖,𝐼_𝑖=𝐼^∗_𝑖,𝑖=1,2,…,𝑚. We can show that 𝑉_EE and 𝛽_𝑖𝑗 satisfy the assumptions of Theorem 3.1 and Corollary  3.3 in [21]. Therefore, the function 𝐿=^𝑛_𝑖=1𝑣_𝑖𝑉_EE (44) is a Lyapunov function for system (15); namely, 𝐿^′|₍₁₅₎≤0 for 𝑃^∗∈Γ₀. One can only show that the largest invariant subset, where 𝐿^′|₍₁₅₎=0, is the singleton {𝑃^∗} by the same argument as in [2–5, 13, 21]. By LaSalle’s invariance principle, 𝑃^∗ is globally asymptotically stable in Γ₀. 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 𝑓_𝑖𝑗(𝑆_𝑖,𝐼_𝑗)=𝑆_𝑖𝐼_𝑗, 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 𝑓(𝑆,𝐼)=𝑓(𝑆)𝐼.

4. Numerical Example

Consider the system (15) when 𝑚=2, 𝑛=2, 𝜑_𝑖(𝑆_𝑖(𝑡))=3−𝑆_𝑖, and 𝑓_𝑖𝑗(𝑆_𝑖,𝐼_𝑗)=𝑆_𝑖𝐼_𝑗, 𝑖,𝑗=1,2. One then has a two-group model as follows: 𝑆^′₁(𝑡)=3−𝑆₁−𝛽₁₁𝑆₁(𝑡)𝐼₁(𝑡)+𝛽₁₂𝑆₁(𝑡)𝐼₂(𝑡,𝑦_′1,1)(𝑡)=1(1+𝛿𝑏)^𝑛𝛽₁₁𝑆₁(𝑡)𝐼₁(𝑡)+𝛽₁₂𝑆₁(𝑡)𝐼₂(𝑡−1)𝑏̂𝑦_1,1(𝑡),𝑦^′_1,2(𝑡)=1𝑏̂𝑦_1,1(𝑡)−𝑦_1,2(𝑡),𝐼^′₁(𝑡)=1𝑏̂𝑦_1,2(𝑡)−𝛿+𝜀₁+𝛾₁𝐼₁(𝑡),𝑆^′₂(𝑡)=3−𝑆₂−𝛽₂₁𝑆₂(𝑡)𝐼₁(𝑡)+𝛽₂₂𝑆₂(𝑡)𝐼₂(𝑡,𝑦_′2,1)(𝑡)=1(1+𝛿𝑏)^𝑛𝛽₂₁𝑆₂(𝑡)𝐼₁(𝑡)+𝛽₂₂𝑆₂(𝑡)𝐼₂(𝑡−1)𝑏̂𝑦_2,1(𝑡),𝑦^′_2,2(𝑡)=1𝑏̂𝑦_2,1(𝑡)−𝑦_2,2(𝑡),𝐼^′₂(𝑡)=1𝑏̂𝑦_2,2(𝑡)−𝛿+𝜀₂+𝛾₂𝐼₂(𝑡). (45) If we choose parameters as 𝛽₁₁=5/24,𝛽₁₂=1,𝛽₂₁=1/36,𝛽₂₂=1/2,𝛿=0.8,𝜀₁=2, 𝜀₂=2, 𝛾₁=1/4, and 𝛾₂=1/4, we can compute 𝑅₀=0.051<1, and hence 𝑃₀=(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).

(a)

(b)

(c)

(d)

Figure 1 

Trajectories of 𝑆1(𝑡),𝐼1(𝑡),𝑆2(𝑡), and 𝐼2(𝑡) for 𝑅0=0.051<1, and 𝑃0=(3,0,0,0,3,0,0,0) is globally stable. 𝑆1(𝑡),𝑆2(𝑡),𝐼1(𝑡), and 𝐼2(𝑡) versus 𝑡 are illustrated by (a), (b), (c), and (d). Initial values are 𝑆1(0)=9,𝑆2(0)=1,𝑦1,1(0)=2,𝑦1,2(0)=2,𝑦2,1(0)=0,𝑦2,2(0)=0,𝐼1(0)=6, and 𝐼2(0)=2.

(a)

(b)

(c)

(d)

Figure 2 

Numerical simulation of (45) with 𝑅0=0.051<1; hence 𝑃0=(3,0,0,0,3,0,0,0) is globally stable. Graphs (a) and (b) illustrate that 𝑆1(t),𝑦1,1(𝑡),𝑦1,2(𝑡) and 𝐼1(𝑡) will eventually towards to steady state. Graphs (c) and (d) illustrate that 𝑆2(𝑡),𝑦2,1(𝑡),𝑦2,2(𝑡), and 𝐼2(𝑡) will eventually towards to steady state. Initial values are 𝑆1(0)=9,𝑆2(0)=1,𝑦1,1(0)=2,𝑦1,2(0)=2,𝑦2,1(0)=0,𝑦2,2(0)=0,𝐼1(0)=6, and 𝐼2(0)=2.

On the other hand, if 𝛽_𝑖𝑗 are chosen as 𝛽₁₁=0.7,𝛽₁₂=1,𝛽₂₁=0.8,𝛽₂₂=1,𝛿=0.5,𝜀₁=0.02, 𝜀₂=0.03, 𝛾₁=0.05, and 𝛾₂=0,05, we can compute 𝑅₀=1.67355>1, and hence 𝑃^∗= (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).

(a)

(b)

(c)

(d)

Figure 3 

Trajectories of 𝑆1(𝑡),𝐼1(𝑡),𝑆2(𝑡), and 𝐼2(𝑡) for 𝑅0=1.67355>1, and 𝑃∗= (0.347644, 0.0760948, 0.0760948, 4.51674, 0.330353, 0.0765909, 0.0765909, 4.4678) is globally stable. 𝑆1(𝑡),𝑆2(𝑡),𝐼1(𝑡), and 𝐼2(𝑡) versus 𝑡 are illustrated by (a), (b), (c), and (d). Initial values are 𝑆1(0)=6,𝑆2(0)=2,𝑦1,1(0)=3,𝑦1,2(0)=3,𝑦2,1(0)=0.1,𝑦2,2(0)=0.1,𝐼1(0)=1.5, and 𝐼2(0)=0.5.

(a)

(b)

(c)

(d)

Figure 4 

Numerical simulation of (45) with 𝑅0=1.67355>1; hence 𝑃∗= (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 𝑆1(𝑡),𝑦1,1(𝑡),𝑦1,2(𝑡), and 𝐼1(𝑡) will eventually towards to steady state. Graphs (c) and (d) illustrate that 𝑆2(𝑡),𝑦2,1(𝑡),𝑦2,2(𝑡), and 𝐼2(𝑡) will eventually towards to steady state. Initial values are 𝑆1(0)=6,𝑆2(0)=2,𝑦1,1(0)=3,𝑦1,2(0)=3,𝑦2,1(0)=0.1,𝑦2,2(0)=0.1,𝐼1(0)=1.5, and 𝐼2(0)=0.5.

Acknowledgments

The authors wish to thank the reviewers for their valuable comments and suggestions that led to truly significant improvement of the paper. J. Wang is supported by the Science and Technology Research Project of the Department of Education of Heilongjiang Province (no. 12531495).