Global Dynamic Behavior of a Multigroup Cholera Model with Indirect Transmission

Abstract

For a multigroup cholera model with indirect transmission, the infection for a susceptible person is almost invariably transmitted by drinking contaminated water in which pathogens, V. cholerae, are present. The basic reproduction number ℛ₀ is identified and global dynamics are completely determined by ℛ₀. It shows that ℛ₀ is a globally threshold parameter in the sense that if it is less than one, the disease-free equilibrium is globally asymptotically stable; whereas if it is larger than one, there is a unique endemic equilibrium which is global asymptotically stable. For the proof of global stability with the disease-free equilibrium, we use the comparison principle; and for the endemic equilibrium we use the classical method of Lyapunov function and the graph-theoretic approach.

1. Introduction

Cholera, a waterborne gastroenteric infection, remains a significant threat to public health in the developing world. Outbreaks of cholera occur cyclically, usually twice per year in endemic areas, and the intensity of these outbreaks varies over longer periods [1]. Hence, in the last few decades, enormous attention is being paid to the cholera disease and several mathematical dynamic models have been developed to study the transmission of cholera [1–7]. In these papers, they consider the population is uniformly mixed, but many factors can lead to heterogeneity in a host population. So in this paper we divide different population into different groups, which can be divided geographically into communities, cities, and countries, to incorporate differential infectivity of multiple strains of the disease agent.

In the case of cholera, the transmission usually occurs through ingestion of contaminated water or feces rather than through casual human-human contact [1]. Therefore, direct contact of healthy person with an infected person is not a risk for contracting infection, whereas a healthy person may contract infection by drinking contaminated water in which pathogens, V. cholerae, are present [2]. The members of this bacterial genus (V. cholerae) naturally colonize in lakes, rivers, and estuaries. Therefore we consider that cholera transmits to other individuals via bacteria in the aquatic environment and formulates a multi-group epidemic model for cholera. Let𝑁be the total population which is divided into four epidemiological compartments, susceptible compartment𝑆, infectious compartment𝐼, recovered compartment 𝑅, and vaccinated compartment𝑉. Let𝑊be the density of V. cholerae in the aquatic environment. As a consequence of the increase in the density of virulent V. cholerae in the aquatic environment, humans become infected and begin to shed increasing numbers of bacteria into the aquatic environment, further elevating bacterial density and exacerbating the outbreak [1]. The growth rate of density of bacteria in the aquatic environment is assumed to be proportional to the number of infectious individuals. We assume that the immunity induced by vaccination is perfect; therefore individuals in vaccinated individuals𝑉cannot be infected. The model is called a multi-group cholera SIRVW epidemic model.

In recent years, multi-group epidemic models have been used to describe the transmission dynamics of many infectious disease in heterogeneous individuals, such as HIV/AIDS [8], dengue [9], West-Nile virus [10], sexually transmitted diseases [11]; and so on. It is well known that global dynamics of multi-group models with higher dimensions, especially the global stability of the endemic equilibrium, is a very challenging problem. Lajmanovich and York [12] proved global stability of the unique endemic equilibrium by using a quadratic global Lyapunov function on a class of 𝑛-group SIS models for gonorrhea; Hethcote [13] proved global stability of the endemic equilibrium for multi-group SIR model without vital dynamics; Thieme [14] proved global stability of the endemic equilibrium of multi-group SEIRS model under certain restrictions. However, they only proved global stability of the endemic equilibrium for multi-group model under some special conditions. In 2006, Guo et al. [15] have first succeeded to establish the complete global dynamics for a multi-group SIR model, by making use of the theory of non-negative matrices, Lyapunov functions and a subtle grouping technique in estimating the derivatives of Lyapunov functions guided by graph theory. By using the results or ideas of [15], the papers [16, 17] proved the global stability of the endemic equilibrium for multi-group model with nonlinear incidence rates and the papers [18, 19] proved the global stability of the endemic equilibrium for multi-group model with distributed delays.

Distinguishing from these multi-group models with direct transmission from person to person, a multi-group cholera model with indirect transmission from the bacteria of the aquatic environment to person is proposed in this paper. We prove that the disease-free equilibrium is globally asymptotically stable if ℛ₀<1, while an endemic equilibrium exists uniquely and is globally asymptotically stable if ℛ₀>1.

The organization of this paper is as follows. In Section 2, we construct a multi-group cholera epidemiological and give some dynamic analysis on the disease-free equilibrium and the endemic equilibrium. An example is given in Section 3 and some conclusions are included in Section 4.

2. Mathematical Modeling and Analysis

For a multi-group epidemic model with cholera, the population of human is divided into𝑛discrete groups, where𝑛∈ℕ. Let 𝑆_𝑖(𝑡), 𝐼_𝑖(𝑡), 𝑅_𝑖(𝑡), and 𝑉_𝑖(𝑡) be the numbers of susceptible, infectious, recovered, and vaccinated individuals in group 𝑖=1,2,…,𝑛at time𝑡, respectively. Let𝑊_𝑖(𝑡)be the density of bacteria in the aquatic environment in group𝑖=1,2,…,𝑛at time𝑡. Based on the assumptions in Section 1, the disease transmission rate of cholera between compartments𝑆_𝑖and𝑊_𝑗is denoted by𝛽_𝑖𝑗, which means the susceptible individuals in the𝑖th group can contact the bacteria of the aquatic environment in the𝑗th (𝑗=1,2,…,𝑛) group. So the new infection that occurred in the𝑖th group is given by∑^𝑛_𝑗=1𝛽_𝑖𝑗𝑆_𝑖𝑊_𝑗. The recruitment rate of individuals into𝑆_𝑖(𝑡)compartment with the𝑖th group is given by a constant𝐴_𝑖. Within the𝑖th group, it is assumed that natural death of human is𝑑_𝑖. A simple immunization policy is considered where the vaccination rate in𝑆_𝑖(𝑡)compartment is given by a constant𝛾_𝑖 and the losing immunity rate from vaccination individuals is𝜆_𝑖. We assume that individuals in𝐼_𝑖(𝑡)compartment recover with a rate constant𝑟_𝑖. In𝑊_𝑖(𝑡)compartment, the brucella shedding rate from𝐼_𝑖(𝑡)compartment is𝑘_𝑖, and the decaying rate of brucella is𝛿_𝑖. So a general multi-group SIRVW epidemic model is described by the following system of differential equations: 𝑑𝑆_𝑖𝑑𝑡=𝐴_𝑖−𝑑_𝑖+𝛾_𝑖𝑆_𝑖+𝜆_𝑖𝑉_𝑖−^𝑛_𝑗=1𝛽_𝑖𝑗𝑆_𝑖𝑊_𝑗,𝑑𝐼_𝑖=𝑑𝑡^𝑛_𝑗=1𝛽_𝑖𝑗𝑆_𝑖𝑊_𝑗,−𝑑_𝑖+𝑟_𝑖𝐼_𝑖,𝑑𝑅_𝑖𝑑𝑡=𝑟_𝑖𝐼_𝑖−𝑑_𝑖𝑅_𝑖,𝑑𝑉_𝑖𝑑𝑡=𝛾_𝑖𝑆_𝑖−𝜆_𝑖+𝑑_𝑖𝑉_𝑖,𝑑𝑊_𝑖𝑑𝑡=𝑘_𝑖𝐼_𝑖−𝛿_𝑖𝑊_𝑖,𝑖=1,2,…,𝑛. (1)

The parameters 𝐴_𝑖, 𝑑_𝑖, 𝜆_𝑖, 𝛾_𝑖, 𝑘_𝑖, and𝛿_𝑖are positive for all𝑖=1,2,…,𝑛, which is made for the biological justification. And we assume that𝛽_𝑖𝑗is nonnegative for all𝑖,𝑗=1,2,…,𝑛and𝑛-square matrix𝐵=(𝛽_𝑖𝑗)_{1≤𝑖,𝑗≤𝑛}is irreducible, which implies that every pair of groups is joined by an infectious path so that the presence of an infectious individual in the first group can cause infection in the second group.

Observe that the variable𝑅_𝑖does not appear in the first and last two equations of system (1); this allows us to consider the following reduced system: 𝑑𝑆_𝑖𝑑𝑡=𝐴_𝑖−𝑑_𝑖+𝛾_𝑖𝑆_𝑖+𝜆_𝑖𝑉_𝑖−^𝑛_𝑗=1𝛽_𝑖𝑗𝑆_𝑖𝑊_𝑗,𝑑𝐼_𝑖=𝑑𝑡^𝑛_𝑗=1𝛽_𝑖𝑗𝑆_𝑖𝑊_𝑗−𝑑_𝑖+𝑟_𝑖𝐼_𝑖,𝑑𝑉_𝑖𝑑𝑡=𝛾_𝑖𝑆_𝑖−𝜆_𝑖+𝑑_𝑖𝑉_𝑖,𝑑𝑊_𝑖𝑑𝑡=𝑘_𝑖𝐼_𝑖−𝛿_𝑖𝑊_𝑖,𝑖=1,2,…,𝑛. (2)

For each group𝑖, adding the four equations in system (2) gives 𝑑𝑆_𝑖+𝐼_𝑖+𝑉_𝑖𝑑𝑡=𝐴_𝑖−𝑑_𝑖𝑆_𝑖+𝐼_𝑖+𝑉_𝑖−𝑟_𝑖𝐼_𝑖≤𝐴_𝑖−𝑑_𝑖𝑆_𝑖+𝐼_𝑖+𝑉_𝑖, (3)

then it follows that lim_𝑡→∞sup𝑆_𝑖+𝐼_𝑖+𝑉_𝑖≤𝐴_𝑖𝑑_𝑖,lim_𝑡→∞sup𝑊_𝑖≤𝑘_𝑖𝐴_𝑖𝑑_𝑖𝛿_𝑖. (4)

Therefore, omega limit sets of system (2) are contained in the following bounded region in the nonnegative cone ofℝ^4𝑛: 𝑋=𝑆_𝑖,𝐼_𝑖,𝑉_𝑖,𝑊_𝑖∣𝑆_𝑖,𝐼_𝑖,𝑉_𝑖,𝑊_𝑖≥0,0≤𝑆_𝑖+𝐼_𝑖+𝑉_𝑖≤𝐴_𝑖𝑑_𝑖,𝑊_𝑖≤𝑘_𝑖𝐴_𝑖𝑑_𝑖𝛿_𝑖,𝑖=1,2,…,𝑛. (5)

It can be verified that region𝑋is positively invariant with respect to system (2). System (2) always has a disease-free equilibrium 𝑃₀=𝑆⁰₁,0,𝑉⁰₁,0,…,𝑆⁰_𝑖,0,𝑉⁰_𝑖,0,…,𝑆⁰_𝑛,0,𝑉⁰_𝑛,0 (6)

on the boundary of𝑋, where 𝑆⁰_𝑖=𝐴_𝑖𝜆_𝑖+𝑑_𝑖𝑑_𝑖𝜆_𝑖+𝑑_𝑖+𝛾_𝑖,𝑉⁰_𝑖=𝐴_𝑖𝛾_𝑖𝑑_𝑖𝜆_𝑖+𝑑_𝑖+𝛾_𝑖. (7)

2.1. The Basic Reproduction Number

According to the next generation matrix formulated in papers [20–22], we define the basic reproduction numberℛ₀of system (2). In order to formulateℛ₀, we order the infected variables first by disease state and then by group, that is, 𝐼₁,𝑊₁,𝐼₂,𝑊₂,…,𝐼_𝑛,𝑊_𝑛. (8)

Consider the following auxiliary system: 𝑑𝐼_𝑖=𝑑𝑡^𝑛_𝑗=1𝛽_𝑖𝑗𝑆_𝑖𝑊_𝑗−𝑑_𝑖+𝑟_𝑖𝐼_𝑖,𝑑𝑊_𝑖𝑑𝑡=𝑘_𝑖𝐼_𝑖−𝛿_𝑖𝑊_𝑖,𝑖=1,2,…,𝑛. (9)

Follow the recipe from van den Driessche and Watmough [21] to obtain 𝐹=0𝛽₁₁𝑆⁰₁0𝛽₁₂𝑆⁰₁⋯0𝛽_1𝑛𝑆⁰₁0000⋯000𝛽₂₁𝑆⁰₂0𝛽₂₂𝑆⁰₂⋯0𝛽_2𝑛𝑆⁰₂0000⋯00⋮⋮⋮⋮⋮⋮0𝛽_𝑛1𝑆⁰_𝑛0𝛽_𝑛2𝑆⁰_𝑛⋯0𝛽_𝑛𝑛𝑆⁰_𝑛0000⋯00_2𝑛×2𝑛,𝑑𝑉=¹+𝑟₁000⋯00−𝑘₁𝛿₁00⋯0000𝑑₂+𝑟₂0⋯0000−𝑘₂𝛿₂⋯00⋮⋮⋮⋮⋮⋮0000⋯𝑑_𝑛+𝑟_𝑛00000⋯−𝑘₁𝛿_{𝑛2𝑛×2𝑛}. (10)

We can get the inverse of𝑉, which equals𝑉⁻¹=1𝑑₁+𝑟₁𝑘000⋯00₁𝛿₁(𝑑₁+𝑟₁)1𝛿₁100⋯0000𝑑₂+𝑟₂𝑘0⋯0000²𝛿₂(𝑑₂+𝑟₂)1𝛿₂1⋯00⋮⋮⋮⋮⋮⋮0000⋯𝑑_𝑛+𝑟_𝑛0𝑘0000⋯^𝑛𝛿_𝑛𝑑_𝑛+𝑟_𝑛1𝛿_{𝑛2𝑛×2𝑛}. (11)

Thus, the next generation matrix is𝐹𝑉⁻¹, 𝐹𝑉⁻¹=𝐴₁₁⋯𝐴_1𝑛𝐵₁₁⋯𝐵_1𝑛𝐴⋮⋮⋮⋮_𝑛1⋯𝐴_𝑛𝑛𝐵_𝑛1⋯𝐵_𝑛𝑛0⋯00⋯0⋮⋮⋮⋮0⋯00⋯0_2𝑛×2𝑛,𝐴𝐴=¹¹𝐴₁₂⋯𝐴_1𝑛𝐴₂₁𝐴₂₂⋯𝐴_2𝑛𝐴⋮⋮⋮_𝑛1𝐴_𝑛2⋯𝐴_{𝑛𝑛𝑛×𝑛}. (12)

So we can calculate the basic reproduction number of system (2), ℛ₀=𝜌𝐹𝑉⁻¹=𝜌(𝐴), (13)

where 𝐴_𝑖𝑗=𝛽_𝑖𝑗𝑘_𝑗𝑆⁰_𝑖𝛿_𝑗𝑑_𝑗+𝑟_𝑗,𝑆⁰_𝑖=𝐴_𝑖𝜆_𝑖+𝑑_𝑖𝑑_𝑖𝜆_𝑖+𝑑_𝑖+𝛾_𝑖,𝑖=1,2,…,𝑛, (14)

and𝜌denotes the spectral radius. As we will show,ℛ₀is the key threshold parameters whose values completely characterize the global dynamics of system (2).

2.2. Global Stability of the Disease-Free Equilibrium of System (2)

For the disease-free equilibrium𝑃₀of system (2), we have the following property.

Theorem 1. Ifℛ₀<1, the disease-free equilibrium𝑃₀of system (2) is globally asymptotically stable in the region𝑋.

Proof. Let𝑀=𝐹−𝑉, and define𝑠(𝑀)=max{Re𝜆∶𝜆is an eigenvalue of𝑀}, so𝑠(𝑀)is a simple eigenvalue of𝑀with a positive eigenvector [23]. By Theorem  2 in [21], there hold two equivalences: ℛ₀>1⟺𝑠(𝑀)>0,ℛ₀<1⟺𝑠(𝑀)<0. (15)
To prove the locally stability of disease-free equilibrium, we check the hypotheses (A1)–(A5) in [21]. Hypotheses (A1)–(A4) are easily verified, while (A5) is satisfied if all eigenvalues of the4𝑛×4𝑛matrix 𝐽|_𝑃0=𝐽𝑀0₃𝐽_{44𝑛×4𝑛} (16)
have negative real parts, where𝐽₃=−𝐹,𝐽₄=−𝑑₁+𝛾₁𝜆₁𝛾00⋯00₁−𝑑₁+𝜆₁𝑑00⋯0000−₂+𝛾₂𝜆₂⋯0000𝛾₂−𝑑₂+𝜆₂𝑑⋯00⋮⋮⋮⋮⋮⋮0000⋯−_𝑛+𝛾_𝑛𝜆_𝑛0000⋯𝛾_𝑛−(𝑑_𝑛+𝜆_𝑛)_2𝑛×2𝑛. (17)
Calculate the eigenvalues of𝐽₄: 𝑠𝐽₄=max−𝑑₁,…,−𝑑_𝑛,−𝑑₁+𝜆₁+𝛾₁−𝑑,…,_𝑛+𝜆_𝑛+𝛾_𝑛<0. (18)
Ifℛ₀<1, then𝑠(𝑀)<0and𝑠(𝐽|_𝑃0)<0, and the disease-free equilibrium𝑃₀of (2) is locally asymptotically stable.
Now we will prove that the disease-free equilibrium𝑃₀of system (2) is globally attractive whenℛ₀<1. From the third equation of system (2), we have 𝑑𝑉_𝑖𝑑𝑡=𝛾_𝑖𝑆_𝑖−𝜆_𝑖+𝑑_𝑖𝑉_𝑖=𝛾_𝑖𝑁_𝑖−𝐼_𝑖+𝑉_𝑖−𝜆_𝑖+𝑑_𝑖𝑉_𝑖≤𝛾_𝑖𝐴_𝑖𝑑_𝑖−𝜆_𝑖+𝛾_𝑖+𝑑_𝑖𝑉_𝑖. (19)
So we can have that, for a small enough positive number𝜖₁, there exists 𝑡_𝑖>0, 𝑖=1,2,…,𝑛 such that for all𝑡>𝑡_𝑖, 𝑉_𝑖≤𝐴_𝑖𝛾_𝑖𝑑_𝑖𝜆_𝑖+𝛾_𝑖+𝑑_𝑖+𝜖₁=𝑉⁰_𝑖+𝜖₁. (20)
Also from the equations of system (2), we have 𝑑𝑆_𝑖𝑑𝑡=𝐴_𝑖+𝜆_𝑖𝑉_𝑖−𝛾_𝑖+𝑑_𝑖𝑆_𝑖−^𝑛_𝑗=1𝛽_𝑖𝑗𝑆_𝑖𝑊_𝑗≤𝐴_𝑖+𝜆_𝑖𝑉⁰_𝑖+𝜖₁−𝑑_𝑖+𝛾_𝑖𝑆_𝑜. (21)
Then lim_𝑡→∞sup𝑆_𝑖=𝐴_𝑖+𝜆_𝑖𝑉⁰_𝑖+𝜖₁𝑑_𝑖+𝛾_𝑖=𝑆⁰_𝑖+𝜖₂,𝜖₂=𝜆_𝑖𝜖₁𝑑_𝑖+𝛾_𝑖. (22)
From system (9) and𝑆_𝑖≤𝑆⁰_𝑖+𝜖₂with all𝑡>𝑡_𝑖. Thus, when𝑡>𝑡_𝑖, we derive 𝑑𝐼_𝑖=𝑆𝑑𝑡^0𝑖+𝜖₂^𝑛_𝑗=1𝛽_𝑖𝑗𝑊_𝑗−𝑑_𝑖+𝑟_𝑖𝐼_𝑖,𝑑𝑊_𝑖𝑑𝑡=𝑘_𝑖𝐼_𝑖−𝛿_𝑖𝑊_𝑖,𝑖=1,2,…,𝑛. (23)
Consider the following auxiliary system 𝑑𝐼^′_𝑖=𝑆𝑑𝑡^0𝑖+𝜖₂^𝑛_𝑗=1𝛽_𝑖𝑗𝑊^′_𝑗−𝑑_𝑖+𝑟_𝑖𝐼^′_𝑖,𝑑𝑊^′_𝑖𝑑𝑡=𝑘_𝑖𝐼^′_𝑖−𝛿_𝑖𝑊^′_𝑖,𝑖=1,2,…,𝑛. (24)
Let𝑀₀be the matrix defined by 𝑀₀=0𝛽₁₁0𝛽₁₂⋯0𝛽_1𝑛0000⋯000𝛽₂₁0𝛽₂₂⋯0𝛽_2𝑛0000⋯00⋮⋮⋮⋮⋮⋮0𝛽_𝑛10𝛽_𝑛2⋯0𝛽_𝑛𝑛0000⋯00_2𝑛×2𝑛, (25)
and set𝑀₁=𝑀+𝜖₂𝑀₀. It follows from Theorem  2 in [21] thatℛ₀<1if and only if𝑠(𝑀)<0. Thus, there exists an𝜖₂>0small enough such that𝑠(𝑀₁)<0. Using the Perron-Frobenius theorem, all eigenvalues of the matrix𝑀₁have negative real parts when𝑠(𝑀₁)<0. Therefore it has 𝐼^′₁(𝑡),𝑊^′₁(𝑡),𝐼^′₂(𝑡),𝑊^′₂(𝑡),…,𝐼^′_𝑛(𝑡),𝑊^′_𝑛(𝑡)⟶(0,0,0,0,…,0,0),𝑡⟶∞, (26)
which implies that the zero solution of system (24) is globally asymptotically stable. Using the comparison principle of Smith and Waltman [23], we know that 𝐼₁(𝑡),𝑊₁(𝑡),𝐼₂(𝑡),𝑊₂(𝑡),…,𝐼_𝑛(𝑡),𝑊_𝑛(𝑡)⟶(0,0,0,0,…,0,0),𝑡⟶∞. (27)
By the theory of asymptotic autonomous system of Thieme [24], it is also known that 𝑆₁(𝑡),𝑉₁(𝑡),…,𝑆_𝑛(𝑡),𝑉_𝑛(𝑡)⟶𝑆₁(0),𝑉₁(0),…,𝑆_𝑛(0),𝑉_𝑛(0),𝑡⟶∞. (28)
So𝑃₀is globally attractive whenℛ₀<1. It follows that the disease-free equilibrium𝑃₀of (2) is globally asymptotically stable whenℛ₀<1. This completes the proof.

2.3. The Uniform Persistence and Unique Positive Solution of System (2)

In this section, we give the proof of the uniform persistence and the unique positive solution of system (2). Define 𝑋₀=𝑆_𝑖,𝐼_𝑖,𝑉_𝑖,𝑊_𝑖∈𝑋∣𝐼_𝑖,𝑊_𝑖>0,𝑖=1,2,…,𝑛,𝜕𝑋₀=𝑋∣𝑋₀. (29)

Theorem 2. Whenℛ₀>1, there exists a positive constant𝜀₁such that when‖𝐼_𝑖(0)‖<𝜀₁,‖𝑊_𝑖(0)‖<𝜀₁for(𝑆_𝑖(0),𝐼_𝑖(0),𝑉_𝑖(0),𝑊_𝑖(0))∈𝑋₀limsup_𝑡→∞max𝐼_𝑖(𝑡),𝑊_𝑖(𝑡)>𝜀₁,𝑖=1,2,…,𝑛. (30)

Proof. Consider the following system: 𝑑𝑆_𝑖𝑑𝑡=𝐴_𝑖+𝜆_𝑖𝑉_𝑖−𝛾_𝑖+𝑑_𝑖𝑆_𝑖,𝑑𝑉_𝑖𝑑𝑡=𝛾_𝑖𝑆_𝑖−𝜆_𝑖+𝑑_𝑖𝑉_𝑖,𝑖=1,2,…,𝑛. (31)
Using Corollary  3.2 in Zhao and Jing [25], it then follows that system (31) has a unique positive equilibrium(𝑆⁰₁,𝑉⁰₁,…,𝑆⁰_𝑛,𝑉⁰_𝑛) which is globally asymptotically stable.
As toℛ₀>1⇔𝑠(𝑀)>0, choose𝜀>0small enough such that𝑠(𝑀₂)>0, where𝑀₂=𝑀−𝜀𝑀₀. Let us consider a perturbed system 𝑑𝑆_𝑖𝑑𝑡=𝐴_𝑖−𝑑_𝑖+𝛾_𝑖𝑆_𝑖+𝜆_𝑖𝑉_𝑖−𝜀₁𝑆_𝑖^𝑛_𝑗=1𝛽_𝑖𝑗,𝑑𝑉_𝑖𝑑𝑡=𝛾_𝑖𝑆_𝑖−𝜆_𝑖+𝑑_𝑖𝑉_𝑖,𝑖=1,2,…,𝑛. (32)
From our previous analysis of system (32), we can restrict𝜀₁>0small enough such that (32) admits a unique positive equilibrium(𝑆⁰_𝑖(𝜀₁),𝑉⁰_𝑖(𝜀₁),𝑖=1,2,…,𝑛)which is globally asymptotically stable.𝑆⁰_𝑖(𝜀₁)is continuous in𝜀₁, so we can further restrict𝜀₁small enough such that𝑆⁰_𝑖(𝜀₁)>𝑆⁰_𝑖−𝜀, 𝑖=1,2,…,𝑛.
For the sake of contradiction of Theorem 2, there is a𝑇>0such that𝐼_𝑖(𝑡)<𝜀₁, 𝑊_𝑖(𝑡)<𝜀₁, 𝑖=1,2,…,𝑛, for all𝑡≥𝑇. Then for𝑡≥𝑇, we have 𝑑𝑆_𝑖𝑑𝑡≥𝐴_𝑖−𝑑_𝑖+𝛾_𝑖𝑆_𝑖+𝜆_𝑖𝑅_𝑖−𝜀₁𝑆_𝑖^𝑛_𝑗=1𝛽_𝑖𝑗,𝑑𝑅_𝑖𝑑𝑡=𝛾_𝑖𝑆_𝑖−𝜆_𝑖+𝑑_𝑖𝑅_𝑖,𝑖=1,2,…,𝑛. (33)
Since the equilibrium (𝑆⁰_𝑖(𝜀₁), 𝑉⁰_𝑖(𝜀₁), 𝑖=1,2,…,𝑛) of (32) is globally asymptotically stable and𝑆⁰_𝑖(𝜀₁)>𝑆⁰_𝑖−𝜀, 𝑖=1,2,…,𝑛. There exists a𝑇₁>𝑇>0such that𝑆_𝑖(𝑡)>𝑆⁰_𝑖−𝜀, 𝑖=1,2,…,𝑛 for 𝑡>𝑇₁. As a consequence, for 𝑡>𝑇₁, there holds 𝑑𝐼_𝑖≥𝑆𝑑𝑡^0𝑖−𝜀^𝑛_𝑗=1𝛽_𝑖𝑗𝑊_𝑗−𝑑_𝑖+𝑟_𝑖𝐼_𝑖,𝑑𝑊_𝑖𝑑𝑡=𝑘_𝑖𝐼_𝑖−𝛿_𝑖𝑊_𝑖,𝑖=1,2,…,𝑛. (34)
Consider the following system 𝑑𝐼^′_𝑖=𝑆𝑑𝑡^0𝑖−𝜀^𝑛_𝑗=1𝛽_𝑖𝑗𝑊^′_𝑗−𝑑_𝑖+𝑟_𝑖𝐼^′_𝑖,𝑑𝑊^′_𝑖𝑑𝑡=𝑘_𝑖𝐼^′_𝑖−𝛿_𝑖𝑊^′_𝑖,𝑖=1,2,…,𝑛. (35)
Since the matrix𝑀₂has positive eigenvalue𝑠(𝑀₂)with a positive eigenvector. It is easy to see that 𝐼^′₁(𝑡),𝑊^′₁(𝑡),𝐼^′₂(𝑡),𝑊^′₂(𝑡),…,𝐼^′_𝑛(𝑡),𝑊^′_𝑛(𝑡)⟶(∞,∞,∞,∞,…,∞,∞),𝑡⟶∞. (36)
Using the comparison principle of Smith and Waltman [23], we also know that 𝐼₁(𝑡),𝑊₁(𝑡),𝐼₂(𝑡),𝑊₂(𝑡),…,𝐼_𝑛(𝑡),𝑊_𝑛(𝑡)⟶(∞,∞,∞,∞,…,∞,∞),𝑡⟶∞, (37)
which leads to a contradiction, therefore we claim that limsup_𝑡→∞max𝐼_𝑖(𝑡),𝑊_𝑖(𝑡)>𝜀₁,𝑖=1,2,…,𝑛. (38)
This completes the proof.

We also have the following result of system (2).

Theorem 3. Ifℛ₀>1, then system (2) admits at least one positive equilibrium, and there is a positive constant𝜀such that every solution(𝑆_𝑖(𝑡),𝐼_𝑖(𝑡),𝑉_𝑖(𝑡),𝑊_𝑖(𝑡))of the system (2) with(𝑆_𝑖(0),𝐼_𝑖(0),𝑉_𝑖(0),𝑊_𝑖(0))∈𝑋₀satisfies minliminf_𝑡→∞𝐼_𝑖(𝑡),liminf_𝑡→∞𝑊_𝑖(𝑡)≥𝜀,𝑖=1,2,…,𝑛, (39)
which implies that system (2) is uniformly persistent.

Proof. Now we prove that system (2) is uniformly persistent with respect to(𝑋₀,𝜕𝑋₀). By the form of (2), it is easy to see that both𝑋and𝑋₀are positively invariant and𝜕𝑋₀is relatively closed in𝑋. Furthermore system (2) is point dissipative. Let 𝑀_𝜕=𝑆_𝑖(0),𝐼_𝑖(0),𝑉_𝑖(0),𝑊_𝑖(0)∣𝑆_𝑖(𝑡),𝐼_𝑖(𝑡),𝑉_𝑖(𝑡),𝑊_𝑖(𝑡)∈𝜕𝑋₀,∀𝑡≥0,𝑖=1,2,…,𝑛. (40)
It is easy to show that 𝑀_𝜕=𝑆_𝑖(𝑡),0,𝑉_𝑖(𝑡),0∣𝑆_𝑖(𝑡),𝑉_𝑖(𝑡)≥0,𝑖=1,2,…,𝑛. (41)
Noting that{(𝑆_𝑖(𝑡),0,𝑉_𝑖(𝑡),0)∣𝑆_𝑖(𝑡),𝑉_𝑖(𝑡)≥0,𝑖=1,2,…,𝑛}⊆𝑀_𝜕. We only need to prove𝑀_𝜕⊆{(𝑆_𝑖(𝑡),0,𝑉_𝑖(𝑡),0)∣𝑆_𝑖(𝑡),𝑉_𝑖(𝑡)≥0,𝑖=1,2,…,𝑛}. Assume(𝑆_𝑖(0),𝐼_𝑖(0),𝑉_𝑖(0),𝑊_𝑖(0),𝑖=1,2,…,𝑛)∈𝑀_𝜕. It suffices to show that𝐼_𝑖(𝑡)=0, 𝑊_𝑖(𝑡)=0, for all 𝑡≥0, 𝑖=1,2,…,𝑛. Suppose not. Then there exist an 𝑖₀, 1≤𝑖₀≤𝑛, and𝑡₀≥0such that𝐼_𝑖0(𝑡₀)>0, 𝑊_𝑖0(𝑡₀)>0 and we partition{1,2,…,𝑛}into two sets𝑄₁and𝑄₂such that 𝐼_𝑖𝑡₀,𝑊_𝑖𝑡₀^𝑇=0,∀𝑖∈𝑄₁,𝐼_𝑖𝑡₀,𝑊_𝑖𝑡₀^𝑇>0,∀𝑖∈𝑄₂. (42)
𝑄₁is nonempty due to the definition of𝑀_𝜕.𝑄₂is non-empty since𝐼_𝑖0(𝑡₀)>0,𝑊_𝑖0(𝑡₀)>0. For any𝑖∈𝑄₂ and we have that 𝑑𝑊_𝑖𝑡₀𝑑𝑡₀=𝑘_𝑖𝐼_𝑖𝑡₀−𝛿_𝑖𝑊_𝑖𝑡₀>𝑘_𝑖𝐼_𝑖𝑡₀,𝑖∈𝑄₂. (43)
It follows that there is an𝜂>0such that𝐼_𝑖(𝑡)>0, for𝑡₀<𝑡<𝑡₀+𝜂, 𝑖∈𝑄₂.This means that(𝑆_𝑖(𝑡),𝐼_𝑖(𝑡),𝑉_𝑖(𝑡),𝑊_𝑖(𝑡),𝑖=1,2,…,𝑛)does not belong to 𝜕𝑋₀for𝑡₀<𝑡<𝑡₀+𝜂, which contradicts the assumption that (𝑆_𝑖(0),𝐼_𝑖(0),𝑉_𝑖(0),𝑊_𝑖(0),𝑖=1,2,…,𝑛)∈𝑀_𝜕. This proves the system (41).
𝑃₀is globally asymptotically stable for system (2). It is clear that there is only an equilibriaum𝑃₀in𝑀_𝜕 and by aforementioned claim, it then follows that𝑃₀is isolated invariant set in𝑋,𝑊^𝑠(𝑃₀)∩𝑋₀=∅. Clearly, every orbit in𝑀_𝜕converges to𝑃₀,𝑃₀is acyclic in𝑀_𝜕. Using Theorem  4.6 in Thieme [26], we conclude that the system (2) is uniformly persistent with respect to(𝑋₀,𝜕𝑋₀). By the result of [27, 28], system (2) has an equilibrium(𝑆^∗₁,𝐼^∗₁,𝑉^∗₁,𝑊^∗₁,…,𝑆^∗_𝑛,𝐼^∗_𝑛,𝑉^∗_𝑛,𝑊^∗_𝑛)∈𝑋₀. We further claim that𝑆^∗_𝑖,𝑉^∗_𝑖>0, 𝑖=1,2,…,𝑛. Suppose that𝑆^∗_𝑖=𝑉^∗_𝑖=0, 𝑖=1,2,…,𝑛; from of (2), we can get𝐼^∗_𝑖=𝑊^∗_𝑖=0, 𝑖=1,2,…,𝑛. It is a contradiction. Then (𝑆^∗₁,𝐼^∗₁,𝑉^∗₁,𝑊^∗₁,…,𝑆^∗_𝑛,𝐼^∗_𝑛,𝑉^∗_𝑛,𝑊^∗_𝑛)∈𝑋₀ is a componentwise positive equilibrium of system (2). This completes the proof.

The following theorem shows that there exists a unique positive solution for system (2) whenℛ₀>1.

Theorem 4. Ifℛ₀>1, then there only exists a unique positive equilibrium 𝑃^∗ for system (2).

Proof. Consider the following system: 𝐴_𝑖−𝑑_𝑖+𝛾_𝑖𝑆_𝑖+𝜆_𝑖𝑉_𝑖−^𝑛_𝑗=1𝛽_𝑖𝑗𝑆_𝑖𝑊_𝑗=0,_𝑛𝑗=1𝛽_𝑖𝑗𝑆_𝑖𝑊_𝑗−𝑑_𝑖+𝑟_𝑖𝐼_𝑖𝛾=0,_𝑖𝑆_𝑖−𝜆_𝑖+𝑑_𝑖𝑉_𝑖𝑘=0,_𝑖𝐼_𝑖−𝛿_𝑖𝑊_𝑖=0,𝑖=1,2,…,𝑛. (44)
We have that 𝑆_𝑖=𝑑_𝑖+𝜆_𝑖𝑑_𝑖𝑑_𝑖+𝜆_𝑖+𝛾_𝑖𝐴_𝑖−𝑑_𝑖+𝑟_𝑖𝐼_𝑖,𝑊_𝑖=𝑘_𝑖𝐼_𝑖𝛿_𝑖,𝑉_𝑖=𝛾_𝑖𝑆_𝑖𝑑_𝑖+𝜆_𝑖,𝑖=1,2,…,𝑛. (45)
Hence, the equilibrium of system (2) is equal to the following system: 𝐵_𝑖𝐴_𝑖−𝑛_𝑖𝐼_𝑖^𝑛_𝑗=1𝛽_𝑖𝑗𝐼_𝑗−𝑛_𝑖𝐼_𝑖=0,𝑖=1,2,…,𝑛, (46)
where 𝐵_𝑖=𝑘_𝑖𝑑_𝑖+𝜆_𝑖𝑑_𝑖𝛿_𝑖𝑑_𝑖+𝜆_𝑖+𝛾_𝑖,𝑛_𝑖=𝑑_𝑖+𝑟_𝑖,𝑖=1,2,…,𝑛. (47)
Therefore, we only need to prove that (46) has a unique positive equilibrium whenℛ₀>1. Use the method in [12] to demonstrate the unique positive equilibrium of (46). First we prove that 𝐼^∗_𝑖=ℎ, 𝑖=1,2,…,𝑛, is the only positive solution of (46). Assume that 𝐼^∗_𝑖=ℎand𝐼^∗_𝑖=𝑘are two positive solutions of (46), both nonzero. Ifℎ≠𝑘, thenℎ_𝑖≠𝑘_𝑖for some 𝑖 (𝑖=1,2,…,𝑛). Assume without loss of generality thatℎ₁>𝑘₁ and moreover thatℎ₁/𝑘₁≥ℎ_𝑖/𝑘_𝑖for all 𝑖 (𝑖=1,2,…,𝑛). Sinceℎand𝑘are positive solutions of (46), we substitute them into (46). We obtain 0=𝐵₁𝐴₁−𝑛₁ℎ₁^𝑛_𝑗=1𝛽_1𝑗ℎ_𝑗−𝑛₁ℎ₁=𝐵₁𝐴₁−𝑛₁𝑘₁^𝑛_𝑗=1𝛽_1𝑗𝑘_𝑗−𝑛₁𝑘₁, (48)
so 0=𝐵₁𝐴₁−𝑛₁ℎ₁^𝑛_𝑗=1𝛽_1𝑗ℎ_𝑗𝑘₁ℎ₁−𝑛₁𝑘₁=𝐵₁𝐴₁−𝑛₁𝑘₁^𝑛_𝑗=1𝛽_1𝑗𝑘_𝑗−𝑛₁𝑘₁,𝐵₁𝐴₁−𝑛₁ℎ₁^𝑛_𝑗=1𝛽_1𝑗ℎ_𝑗𝑘₁ℎ₁=𝐵₁𝐴₁−𝑛₁𝑘₁^𝑛_𝑗=1𝛽_1𝑗𝑘_𝑗. (49)
But(ℎ_𝑖/ℎ₁)𝑘₁≤𝑘_𝑖and𝐵₁(𝐴₁−𝑛₁ℎ₁)<𝐵₁(𝐴₁−𝑛₁𝑘₁); thus from the above equalities we get 𝐵₁𝐴₁−𝑛₁ℎ₁^𝑛_𝑗=1𝛽_1𝑗ℎ_𝑗𝑘₁ℎ₁≤𝐵₁𝐴₁−𝑛₁ℎ₁^𝑛_𝑗=1𝛽_1𝑗𝑘_𝑗<𝐵₁𝐴₁−𝑛₁𝑘₁^𝑛_𝑗=1𝛽_1𝑗𝑘_𝑗. (50)
This is a contradiction, so there is only one positive solution𝐼^∗_𝑖=ℎ, 𝑖=1,2,…,𝑛, of (46). So when ℛ₀>1, there only exists a unique positive equilibrium for system (2).

2.4. Global Stability of the Unique Endemic Solution of System (2)

In this section, we prove that the unique endemic equilibrium of system (2) is globally asymptotically stable in𝑋₀. In order to prove global stability of the endemic equilibrium, the Lyapunov function will be used. In the following, we also use a Lyapunov function to prove global stability of the endemic equilibrium.

Theorem 5. Ifℛ₀>1, the unique positive equilibrium𝑃^∗of system (2) is globally asymptotically stable in𝑋₀.

Proof. Following [15] we define 𝜉_𝑖𝑗=𝛽_𝑖𝑗𝑆^∗_𝑖𝑊^∗_𝑗,1≤𝑖,𝑗≤𝑛,𝑛≥2,𝐁=^𝑛𝑗≠1𝜉_1𝑗−𝜉₂₁⋯−𝜉_𝑛1−𝜉^𝑛_12𝑗≠2𝜉_2𝑗⋯−𝜉_𝑛2⋮⋮⋱⋮−𝜉_1𝑛−𝜉_2𝑛⋯^𝑛_𝑗≠𝑛𝜉_{𝑛𝑗𝑛×𝑛}, (51)(52)
which is a Laplacian matrix whose column sums are zero and which is irreducible. Therefore, it follows from Lemma  2.1 of [15] that the solution space of linear system 𝐁𝜁=0, (53)
has dimension 1 with a basis 𝜁=𝜁₁,𝜁₂,…,𝜁_𝑛^𝑇=𝑐₁,𝑐₂,…,𝑐_𝑛^𝑇, (54)
where𝑐_𝑖denotes the cofactor of the𝑖th diagonal entry of𝐁. Note that from (53) we have that ^𝑛_𝑗=1𝜁_𝑖𝜉_𝑖𝑗=^𝑛_𝑗=1𝜁_𝑗𝜉_𝑗𝑖,𝑖=1,2,…,𝑛. (55)
For such𝜁=(𝜁₁,𝜁₂,…,𝜁_𝑛)we define a Lyapunov function 𝐿()=𝐒,𝐈,𝐕,𝐖_𝑛𝑖=1𝜁_𝑖𝑆_𝑖−𝑆^∗_𝑖−𝑆^∗_𝑖𝑆ln^∗𝑖𝑆_𝑖+𝐼_𝑖−𝐼^∗_𝑖−𝐼^∗_𝑖𝐼ln^∗𝑖𝐼_𝑖+𝑉_𝑖−𝑉^∗_𝑖−𝑉^∗_𝑖𝑉ln^∗𝑖𝑉_𝑖+𝑑_𝑖+𝑟_𝑖𝑘_𝑖𝑊_𝑖−𝑊^∗_𝑖−𝑊^∗_𝑖𝑊ln^∗𝑖𝑊_𝑖, (56)
where𝐒=(𝑆₁,𝑆₂,…,𝑆_𝑛), 𝐈=(𝐼₁,𝐼₂,…,𝐼_𝑛), 𝐕=(𝑉₁,𝑉₂,…,𝑉_𝑛), and𝐖=(𝑊₁,𝑊₂,…,𝑊_𝑛). It is easy to see that𝐿(𝐒,𝐈,𝐕,𝐖)≥0for all(𝐒,𝐈,𝐕,𝐖)≥0and the equality𝐿(𝐒,𝐈,𝐕,𝐖)=0holds if and only if(𝐒,𝐈,𝐕,𝐖)=(𝐒^∗,𝐈^∗,𝐕^∗,𝐖^∗). The derivative along the trajectories of system (2) is 𝐿^′()=𝐒,𝐈,𝐕,𝐖_𝑛𝑖=1𝜁_𝑖𝐴_𝑖−𝑑_𝑖+𝛾_𝑖𝑆_𝑖+𝜆_𝑖𝑉_𝑖−^𝑛_𝑗=1𝛽_𝑖𝑗𝑆_𝑖𝑊_𝑗−𝑆^∗_𝑖𝑆_𝑖𝐴_𝑖−𝑑_𝑖+𝛾_𝑖𝑆_𝑖+𝜆_𝑖𝑉_𝑖−^𝑛_𝑗=1𝛽_𝑖𝑗𝑆_𝑖𝑊_𝑗+^𝑛_𝑗=1𝛽_𝑖𝑗𝑆_𝑖𝑊_𝑗−𝑑_𝑖+𝑟_𝑖𝐼_𝑖−𝐼^∗_𝑖𝐼_𝑖^𝑛_𝑗=1𝛽_𝑖𝑗𝑆_𝑖𝑊_𝑗−𝑑_𝑖+𝑟_𝑖𝐼_𝑖+𝛾_𝑖𝑆_𝑖−𝜆_𝑖+𝑑_𝑖𝑉_𝑖−𝑉^∗_𝑖𝑉_𝑖𝛾_𝑖𝑆_𝑖−𝜆_𝑖+𝑑_𝑖𝑉_𝑖+𝑑_𝑖+𝑟_𝑖𝑘_𝑖𝑘_𝑖𝐼_𝑖−𝛿_𝑖𝑊_𝑖−𝑊^∗_𝑖𝑊_𝑖𝑘_𝑖𝐼_𝑖−𝛿_𝑖𝑊_𝑖=𝐿₁+𝐿₂+𝐿₃. (57)
From system (44), we have 𝐴_𝑖=𝑑_𝑖+𝛾_𝑖𝑆^∗_𝑖−𝜆_𝑖𝑉^∗_𝑖+^𝑛_𝑗=1𝛽_𝑖𝑗𝑆^∗_𝑖𝑊^∗_𝑗,_𝑛𝑗=1𝛽_𝑖𝑗𝑆^∗_𝑖𝑊^∗_𝑗=𝑑_𝑖+𝑟_𝑖𝐼^∗_𝑖=𝛿_𝑖𝑑_𝑖+𝑟_𝑖𝑊^∗_𝑖𝑘_𝑖. (58)(59)
So 𝐿₁=^𝑛_𝑖=1𝜁_𝑖^𝑛_𝑗=1𝛽_𝑖𝑗𝑆^∗_𝑖𝑊_𝑗−𝛿_𝑖𝑑_𝑖+𝑟_𝑖𝑊_𝑖𝑘_𝑖,𝐿₂=^𝑛_𝑖=1𝜁_𝑖𝑑_𝑖+𝛾_𝑖𝑆^∗_𝑖−𝜆_𝑖𝑉^∗_𝑖−𝑑_𝑖+𝛾_𝑖𝑆_𝑖+𝜆_𝑖𝑉_𝑖+𝑆^∗_𝑖𝑆_𝑖𝑑_𝑖+𝛾_𝑖𝑆^∗_𝑖−𝜆_𝑖𝑉^∗−𝑑_𝑖+𝛾_𝑖𝑆_𝑖+𝜆_𝑖𝑉_𝑖+𝛾_𝑖𝑆_𝑖−𝜆_𝑖+𝑑_𝑖𝑉_𝑖+𝑉^∗_𝑖𝑉_𝑖𝛾_𝑖𝑆_𝑖−𝜆_𝑖+𝑑_𝑖𝑉_𝑖=^𝑛_𝑖=1𝜁_𝑖𝑑_𝑖𝑆^∗_𝑖𝑆2−^𝑖𝑆^∗_𝑖−𝑆^∗_𝑖𝑆_i+𝜆_𝑖𝑉^∗_𝑖𝑆2−^𝑖𝑉^∗_𝑖𝑆^∗_𝑖𝑉_𝑖−𝑆^∗_𝑖𝑉_𝑖𝑆_𝑖𝑉^∗_𝑖+𝑑_𝑖𝑉^∗_𝑖𝑉3−^𝑖𝑉^∗_𝑖−𝑆^∗_𝑖𝑆_𝑖−𝑆_𝑖𝑉^∗_𝑖𝑆^∗_𝑖𝑉_𝑖𝐿≤0,₃=^𝑛_𝑖=1𝜁_𝑖3^𝑛_𝑗=1𝛽_𝑖𝑗𝑆^∗_𝑖𝑊^∗_𝑗−^𝑛_𝑗=1𝛽_𝑖𝑗𝑆^∗_𝑖𝑊^∗_𝑗𝑆^∗_𝑖𝑆_𝑖−^𝑛_𝑗=1𝛽_𝑖𝑗𝑆_𝑖𝑊_𝑗𝐼^∗_𝑖𝐼_𝑖−𝑑_𝑖+𝑟_𝑖𝐼_𝑖𝑊^∗_𝑖𝑊_𝑖. (60)
Now we claim that ^𝑛_𝑖=1𝜁_𝑖^𝑛_𝑗=1𝛽_𝑖𝑗𝑆^∗_𝑖𝑊_𝑗=^𝑛_𝑖=1𝜁_𝑖𝛿_𝑖𝑑_𝑖+𝑟_𝑖𝑊_𝑖𝑘_𝑖. (61)
Appealing to (51), (55) and (59), ^𝑛𝑛_𝑖=1𝑗=1𝜁_𝑖𝛽_𝑖𝑗𝑆^∗_𝑖𝑊_𝑗=^𝑛𝑛_𝑖=1𝑗=1𝜁_𝑗𝛽_𝑗𝑖𝑆^∗_𝑗𝑊_𝑖=^𝑛𝑛_𝑖=1𝑗=1𝑊_𝑖𝑊^∗_𝑖𝜁_𝑗𝛽_𝑗𝑖𝑆^∗_𝑗𝑊^∗_𝑖=^𝑛_𝑖=1𝑊_𝑖𝑊^∗_𝑖^𝑛_𝑗=1𝜁_𝑗𝜉_𝑗𝑖=^𝑛_𝑖=1𝑊_𝑖𝑊^∗_𝑖^𝑛_𝑗=1𝜁_𝑖𝜉_𝑖𝑗=^𝑛_𝑖=1𝜁_𝑖𝛿_𝑖𝑑_𝑖+𝑟_𝑖𝑊_𝑖𝑘_𝑖. (62)
From (61) we have 𝐿^′()≤𝐒,𝐈,𝐕,𝐖_𝑛𝑖=1𝜁_𝑖3^𝑛_𝑗=1𝛽_𝑖𝑗𝑆^∗_𝑖𝑊^∗_𝑗−^𝑛_𝑗=1𝛽_𝑖𝑗𝑆^∗_𝑖𝑊^∗_𝑗𝑆^∗_𝑖𝑆_𝑖−^𝑛_𝑗=1𝛽_𝑖𝑗𝑆_𝑖𝑊_𝑗𝐼^∗_𝑖𝐼_𝑖−𝑑_𝑖+𝑟_𝑖𝐼_𝑖𝑊^∗_𝑖𝑊_𝑖=^𝑛_𝑖,𝑗=1𝜁_𝑖𝜉_𝑖𝑗𝑆3−^∗𝑖𝑆_𝑖−𝑆_𝑖𝑊_𝑗𝐼^∗_𝑖𝑆^∗_𝑖𝑊^∗_𝑗𝐼_𝑖−𝑊^∗_𝑖𝐼_𝑖𝑊_𝑖𝐼^∗_𝑖=∶𝐻_𝑛𝑆₁,𝐼₁,𝑊₁,…,𝑆_𝑛,𝐼_𝑛,𝑊_𝑛. (63)
Next we show that𝐻_𝑛≤0for all(𝑆₁,𝐼₁,𝑊₁,…,𝑆_𝑛,𝐼_𝑛,𝑊_𝑛)∈𝑋₀by applying the graph-theoretic approach developed in [29–31]. As in [29],𝐿=𝐺(𝐵)denotes the directed graph associated with matrix𝐁,𝑄presents a subgraph of𝐿,𝐶𝑄denotes the unique elementary cycle of𝑄,𝐸(𝐶𝑄)presents the set of directed arcs in 𝐶𝑄, and𝑙=𝑙(𝑄)denotes the number of arcs in𝐶𝑄. Then𝐻_𝑛can be rewritten as 𝐻_𝑛=_𝑄𝐻_𝑛,𝑄, (64)
where 𝐻_𝑛,𝑄=_{()(𝑄)𝑟,𝑚∈𝐸}𝜉_𝑟𝑚×3𝑙−_{()𝑖,𝑗∈𝐸𝐶𝑄}𝑆^∗_𝑖𝑆_𝑖+𝑆_𝑖𝑊_𝑗𝐼^∗_𝑖𝑆^∗_𝑖𝑊^∗_𝑗𝐼_𝑖+𝑊^∗_𝑖𝐼_𝑖𝑊_𝑖𝐼^∗_𝑖. (65)
For instance, 𝐻₁=𝐻₁𝑆₁,𝐼₁,𝑊₁=_𝑖=𝑗=1𝜁₁𝜉₁₁𝑆3−^∗1𝑆₁−𝑆₁𝑊₁𝐼^∗₁𝑆^∗₁𝑊^∗₁𝐼₁−𝑊^∗₁𝐼₁𝑊₁𝐼^∗₁𝐻≤0,₂=𝐻₂𝑆₁,𝐼₁,𝑊₁,𝑆₂,𝐼₂,𝑊₂=²_𝑖,𝑗=1𝜁_𝑖𝜉_𝑖𝑗𝑆3−^∗𝑖𝑆_𝑖−𝑆_𝑖𝑊_𝑗𝐼^∗_𝑖𝑆^∗_𝑖𝑊^∗_𝑗𝐼_𝑖−𝑊^∗_𝑖𝐼_𝑖𝑊_𝑗𝐼^∗_𝑖=𝜉₁₁𝜉₂₁𝑆3−^∗1𝑆₁−𝑆₁𝑊₁𝐼^∗₁𝑆^∗₁𝑊^∗₁𝐼₁−𝑊^∗₁𝐼₁𝑊₁𝐼^∗₁+𝜉₂₂𝜉₁₂𝑆3−^∗2𝑆₂−𝑆₂𝑊₂𝐼^∗₂𝑆^∗₂𝑊^∗₂𝐼₂−𝑊^∗₂𝐼₂𝑊₂𝐼^∗₂+𝜉₁₂𝜉₂₁𝑆6−^∗1𝑆₁−𝑆₁𝑊₁𝐼^∗₁𝑆^∗₁𝑊^∗₁𝐼₁−𝑊^∗₁𝐼₁𝑊₁𝐼^∗₁−𝑆^∗₂𝑆₂−𝑆₂𝑊₂𝐼^∗₂𝑆^∗₂𝑊^∗₂𝐼₂−𝑊^∗₂𝐼₂𝑊₂𝐼^∗₂≤0. (66)
Note that for each unicycle graph𝑄, it is easy to see that _{(𝑖,𝑗)∈𝐸(𝐶𝑄)}𝑆^∗_𝑖𝑆_𝑖⋅𝑆_𝑖𝑊_𝑗𝐼^∗_𝑖𝑆^∗_𝑖𝑊^∗_𝑗𝐼_𝑖⋅𝑊^∗_𝑖𝐼_𝑖𝑊_𝑖𝐼^∗_𝑖=_{(𝑖,𝑗)∈𝐸(𝐶𝑄)}𝑊^∗_𝑖𝑊_𝑗𝑊_𝑖𝑊^∗_𝑗=1. (67)
Therefore, _{()𝑖,𝑗∈𝐸𝐶𝑄}𝑆^∗_𝑖𝑆_𝑖+𝑆_𝑖𝑊_𝑗𝐼^∗_𝑖𝑆^∗_𝑖𝑊^∗_𝑗𝐼_𝑖+𝑊^∗_𝑖𝐼_𝑖𝑊_𝑖𝐼^∗_𝑖≥3𝑙, (68)
and hence𝐻_𝑛,𝑄≤0for each𝑄, and𝐻_𝑛,𝑄=0if and only if 𝑆^∗_𝑖𝑆_𝑖=𝑆_𝑖𝑊_𝑗𝐼^∗_𝑖𝑆^∗_𝑖𝑊^∗_𝑗𝐼_𝑖=𝑊^∗_𝑖𝐼_𝑖𝑊_𝑖𝐼^∗_𝑖,𝑖,𝑗∈𝐸(𝐶𝑄). (69)
Thus 𝐿^′(𝐒,𝐈,𝐕,𝐖)≤𝐻_𝑛≤0. (70)
The equality 𝐿^′(𝐒,𝐈,𝐕,𝐖)=0 holds if and only if 𝑆_𝑖=𝑆^∗_𝑖, 𝐼_𝑖=𝐼^∗_𝑖, 𝑉_𝑖=𝑉^∗_𝑖, and 𝑊_𝑖=𝑊^∗_𝑖 for all𝑖=1,2,…,𝑛. Therefore, following from LaSalle’s Invariance Principle [32], the unique endemic equilibrium𝑃^∗of system (2) is globally asymptotically stable. This completes the proof.

3. A Numerical Example

Consider the system (1), when𝑖=2, one has the two-group model as follows: 𝑑𝑆₁𝑑𝑡=𝐴₁−𝑑₁+𝛾₁𝑆₁+𝜆₁𝑉₁−𝛽₁₁𝑆₁𝑊₁+𝛽₁₂𝑆₁𝑊₂,𝑑𝐼₁𝑑𝑡=𝛽₁₁𝑆₁𝑊₁+𝛽₁₂𝑆₁𝑊₂−𝑑₁+𝑟₁𝐼₁,𝑑𝑅₁𝑑𝑡=𝑟₁𝐼_𝑖−𝑑₁𝑅₁,𝑑𝑉₁𝑑𝑡=𝛾₁𝑆₁−𝜆₁+𝑑₁𝑉₁,𝑑𝑊₁𝑑𝑡=𝑘₁𝐼₁−𝛿₁𝑊₁,𝑑𝑆₂𝑑𝑡=𝐴₂−𝑑₂+𝛾₂𝑆₂+𝜆₂𝑉₂−𝛽₂₁𝑆₂𝑊₁+𝛽₂₂𝑆₂𝑊₂,𝑑𝐼₂𝑑𝑡=𝛽₂₁𝑆₂𝑊₁+𝛽₂₂𝑆₂𝑊₂−𝑑₂+𝑟₂𝐼₂,𝑑𝑅₂𝑑𝑡=𝑟₂𝐼₂−𝑑₂𝑅₂,𝑑𝑉₂𝑑𝑡=𝛾₂𝑆₂−𝜆₂+𝑑₂𝑉₂,𝑑𝑊₂𝑑𝑡=𝑘₂𝐼₂−𝛿₂𝑊₂. (71)

We can give the basic reproduction number of system (71), which is ℛ^′₀=𝐴₁₁+𝐴₂₂+𝐴₁₁−𝐴₂₂²+4𝐴₁₂𝐴₂₁2, (72)

where 𝐴_𝑖𝑗=𝛽_𝑖𝑗𝑘_𝑗𝑆⁰_𝑖𝛿_𝑗𝑑_𝑗+𝑟_𝑗,𝑆⁰_𝑖=𝐴_𝑖𝜆_𝑖+𝑑_𝑖𝑑_𝑖𝜆_𝑖+𝑑_𝑖+𝛾_𝑖,𝑖=1,2. (73)

Taking𝐴₁=150, 𝐴₂=220, 𝑑₁=0.1, 𝑑₂=0.1, 𝜆₁=0.4, 𝜆₂=0.6, 𝜆₁=0.5, 𝜆₂=0.5, 𝑟₁=1, 𝑟₂=1, 𝑘₁=10, 𝑘₂=10, 𝛿₁=8, 𝛿₂=8 and using Matlab ODE solver, we run numerical simulations for two cases.

If 𝛽₁₁=0.00048, 𝛽₁₂=0.0004, 𝛽₂₁=0.0004, and 𝛽₂₂=0.00045, we have ℛ^′₀≈0.9804<1. Hence the disease-free equilibrium of system (71) is globally asymptotically stable (see Figure 1(a)). If𝛽₁₁=0.0025, 𝛽₁₂=0.001, 𝛽₂₁=0.001, and 𝛽₂₂=0.0020, we have ℛ^′₀≈3.6594>1. Hence the endemic equilibrium of system (71) is globally asymptotically stable (see Figure 1(b)).

(a)

(b)

Figure 1 

(a) The disease dies out in both groups. (b) The disease persists in both groups. Initial conditions are: 𝑆1(0)=280, 𝐼1(0)=40, 𝑅1(0)=10, 𝑉1(0)=130, 𝑊1(0)=250, 𝑆2(0)=260, 𝐼2(0)=20, 𝑅2(0)=10, 𝑉2(0)=130, 𝑊2(0)=300.

4. Conclusion

Cholera epidemic has become a major health problem for many developing countries. From good understanding of the transmission dynamics of cholera in many emergent epidemic regions, the heterogeneous host population can be divided into several homogeneous groups according to modes of transmission, contact patterns, or geographic distributions. Hence, in this paper, we proposed a multi-group cholera SIRVW epidemiological model. In order to distinguish many multi-group models with direct transmission from person to person, we only considered this multi-group cholera model with indirect transmission from the bacteria of the aquatic environment to person. Firstly, the basic reproduction numberℛ₀of this model is given. Then, it is found that the model has two non-negative equilibria, the disease-free equilibrium and the endemic equilibrium. The disease-free equilibrium exists without any condition whereas the endemic equilibrium exists providedℛ₀>1. Finally, through the analysis of the model it has been found that the global asymptotic behavior of multi-group SIRVW model is completely determined by the size ofℛ₀. That is, the disease-free equilibrium is globally asymptotically stable ifℛ₀<1while an endemic equilibrium exists uniquely and is globally asymptotically stable ifℛ₀>1. By running numerical simulations for the cases of two-groups model, we can see that the disease-free equilibrium of system (71) is globally stable whenℛ^′₀<1and the unique endemic equilibrium of system (71) is globally stable whenℛ^′₀>1.

Acknowledgments

This research is supported by the National Natural Science Foundation of China under Grants 11301490, 11301491, 11331009, 11171314, and 11147015, Natural Science Foundation of Shan’Xi Province Grant no. 2012021002-1, the specialized research fund for the doctoral program of higher education preferential development, no. 20121420130001, China Postdoctoral Science Foundation under Grant no. 2012M520814, Shanghai Postdoctoral Science Foundation under Grants no. 13R21410100 and IDRC104519-010.