Stability analysis of an HIV/AIDS epidemic model with sexual transmission in a patchy environment

ABSTRACT

A multi-patch HIV/AIDS model with heterosexual transmission is formulated to investigate the impact of population migration on the spread of HIV/AIDS. We derive the basic reproduction number $R_0$ and prove that if $R_0<1$, the disease-free equilibrium is globally asymptotically stable. If $R_0>1$ and certain conditions are satisfied, the endemic equilibrium is globally asymptotically stable. We apply the model to two patches and conduct numerical simulations. If HIV/AIDS becomes extinct in each patch when two patches are isolated, the disease remains extinct in two patches when the population migration occurs; if HIV/AIDS spreads in each patch when two patches are isolated, the disease remains persistent in two patches when the population migration occurs; if the disease disappears in one patch and spreads in the other patch when they are isolated, the disease can spread or disappear in two patches if migration rates of individuals are suitably chosen.

KEYWORDS:

• Multi-patch
• HIV/AIDS
• migration
• basic reproduction number
• equilibrium

MATHEMATICS SUBJECT CLASSIFICATIONS:

• 34D20
• 34M04

1. Introduction

HIV/AIDS remains one of the world's most significant public health challenges [1]. AIDS is mainly transmitted through three ways: sexual contact, blood and mother and infant, and sexual transmission is the most important way to spread AIDS in China. Since the discovery of HIV in the early 1980s, the disease has spread in successive waves to most regions around the global [37]. More than 95% of the newly diagnosed HIV infected people in China are sexually transmitted, and about 70% of them are heterosexual transmission [4]. With the growth of travel among cities and countries, new infectious diseases are developing regionally and globally much faster than before [19]. In the process of population migration, HIV/AIDS can easily spread from one region to another.

Many disease transmission models in multiple patches environment have been proposed and investigated. See , for example [6, 7]. Li et al. [22] established the SIR model of population migration in n patches, analysed the stability of disease-free equilibrium, and proved endemic equilibrium is globally asymptotically stable by constructing Lyapunov function and combining graph theory. Guo et al. [16] also studied global stability of the endemic equilibrium by constructing Lyapunov function and combining graph theory. Arino et al. [5] established an epidemic model of n cities to study the impact of human movement between cities on disease transmission. Wang et al. [33] established an SIS model of n patches, which assumes that the susceptible individuals and infected individuals of each patch have the same diffusion rate. Under the same assumptions, Jin et al. [18] further proved the uniqueness and global stability of endemic equilibrium by using monotone dynamic system theory. Wang et al. [34] established a multi-patch SEIQR model to discuss the impact of entry–exit inspection on disease control. Gao et al. [15] established a multi-patch SIS model to study the impact of media reports and human movement on disease transmission, and conducted numerical simulation on the two patches. The results show that when the two patches are isolated, the disease disappears in the first patch and prevails in the second patch. When the two patches are connected by travel, the disease may extinct in both patches. Cosner et al. [12] established a patch model to study the impact of human movement on media transmission and took two patches as examples for numerical simulation. Mccormack et al. [25] studied the multi-patch deterministic SIS model and random SIR model for wild animal diseases and took three patches as examples to conduct numerical simulation to study the influence of three different movement rates on disease transmission. Gao et al. [14] established the SIS epidemic model for n patches with a standard incidence rate. Finally, two patches are taken as examples for numerical simulation. Zhang et al. [35] established a two-patch model for the spread of West Nile virus. Wang et al. [32] proposed an epidemic model to describe the dynamics of disease spread between two patches due to population dispersal. As far as we know, few people have established a model of HIV/AIDS transmission in a patchy environment.

In the research on HIV/AIDS, Bacaër et al. [8] proposed an HIV/AIDS epidemic model, considering both injecting drug use and heterosexual transmission and give out some numerical results. Zhao et al. [37] established a model of AIDS transmission including three age groups, analysed the global dynamics of the system and further demonstrated the difference of AIDS transmission in different age groups by numerical simulation. Cai et al. [10] established an HIV/AIDS epidemic model with treatment. However, there have been relatively few studies linking HIV/AIDS to patches. Kheiri et al. [20] formulated a fractional order model for the HIV/AIDS epidemic in a patchy environment to investigate the effect of human movement on the spread of HIV/AIDS epidemic among patches. The model was established in combination with Euler model and adopted bilinear incidence. The basic reproduction number $R_0$ was derived and numerical simulation is done in two patches. The simulation showed that the disease died out in patch 2 and became endemic in patch 1 when there was no migration between two patches, that is, they are isolated. The disease became endemic in both patches when there was migration between two patches. When two sexes are considered, there is still less research on the impact of population migration on the spread of HIV/AIDS.

Cosner et al. [12] and Sattenspiel et al. [27] introduced two models of human movement, including Euler movement model and Lagrange movement model. The former only considers the individuals' migration from one place to another without tracking the individuals ' behaviour. It considers that the individuals will not return after migration. After reaching a place, the individuals have no memory of the previous position. The latter, in contrast, tracks individuals ' behaviour, which considers the individuals temporary travel to another location and return to the previous location. Citron et al. [11] further introduced the two models and pointed out that the flow model is the Euler movement model and the simple trip model is the Lagrangian movement model. Examples are given to illustrate how to combine the two models with the research problems. Combined with the Euler movement model, our goal is establish a patch model of HIV/AIDS to study the impact of population migration on the spread of HIV/AIDS among regions.

The rest of this paper is organized as follows. In the next section, combined with the Euler movement model, we develop a new HIV/AIDS model with heterosexual transmission and constant recruitment rate in a patchy environment. In Section 3, we derive the basic reproduction number $R_0$ which serves as a threshold for the disease extinction and persistence. The global stability of the disease-free equilibrium is proved. When certain conditions are satisfied, the global stability of the endemic equilibrium is also proved. In Section 4, we apply the model to two patches. In Section 5, we conduct numerical simulations using the data from Beijing City and Yunnan Province in China. First, we present time series diagrams under different conditions. Next, we use contour plots to show the relationship between $R_0$ and migration rates. Finally, we make a comparison between no migration and migration in two patches. In Section 6, we conduct sensitivity analysis of $R_0$. In Section 7, we give a brief conclusion.

2. Model formulation

In this section, we propose a susceptible-infected-AIDS (SIA) model with standard incidence for the heterosexual transmission of HIV/AIDS with population migration between patches. A patch can be interpreted as a city, a province or a country, etc. Let n be the total number of patches. In patch i, the population is divided into six subclasses: the susceptible female, the susceptible male, the infective female (no clinical symptoms of AIDS), the infective male, the AIDS female and the AIDS male whose numbers are denoted by $S_{if}(t),S_{im}(t),I_{if}(t),I_{im}(t),A_{if}(t),A_{im}(t)$ respectively. Susceptible individuals are transformed into infected individuals through contact with heterosexual infected individuals. The total population in patch i at time t is given by $N_i(t)=S_{if}(t)+S_{im}(t)+I_{if}(t)+I_{im}(t)+A_{if}(t)+A_{im}(t)$. The flow diagram of HIV/AIDS transmission in patch i is illustrated in Figure 1.

We assume that susceptible individuals become infected individuals through contact with heterosexual infected individuals and they have the same transmission rate. Moreover, we assume that susceptible individuals, infected individuals and AIDS individuals of each patch have the same migration rate. Based on these assumptions, a multi-patch SIA epidemic model with standard incidence is expressed in terms of the following system of ordinary differential equations: (1) $\begin{array}{rl}\frac{\mathrm{d}{S}_{if}}{\mathrm{d}t}& ={\epsilon }_i{\mathrm{\Lambda }}_i-{\beta }_i\frac{S_{if}{I}_{im}}{N_{im}}-{\mu }_{if}{S}_{if}+\sum _{j=1,j\ne i}^n(m_{ij}{S}_{jf}-m_{ji}{S}_{if}),\\ \frac{\mathrm{d}{I}_{if}}{\mathrm{d}t}& ={\beta }_i\frac{S_{if}{I}_{im}}{N_{im}}-({\mu }_{if}+{\gamma }_{if})I_{if}+\sum _{j=1,j\ne i}^n(m_{ij}{I}_{jf}-m_{ji}{I}_{if}),\\ \frac{\mathrm{d}{A}_{if}}{\mathrm{d}t}& ={\gamma }_{if}{I}_{if}-({\mu }_{if}+d_{if})A_{if}+\sum _{j=1,j\ne i}^n(m_{ij}{A}_{jf}-m_{ji}{A}_{if}),\\ \frac{\mathrm{d}{S}_{im}}{\mathrm{d}t}& =(1-{\epsilon }_i){\mathrm{\Lambda }}_i-{\beta }_i\frac{S_{im}{I}_{if}}{N_{if}}-{\mu }_{im}{S}_{im}+\sum _{j=1,j\ne i}^n(m_{ij}{S}_{jm}-m_{ji}{S}_{im}),\\ \frac{\mathrm{d}{I}_{im}}{\mathrm{d}t}& ={\beta }_i\frac{S_{im}{I}_{if}}{N_{if}}-({\mu }_{im}+{\gamma }_{im})I_{im}+\sum _{j=1,j\ne i}^n(m_{ij}{I}_{jm}-m_{ji}{I}_{im}),\\ \frac{\mathrm{d}{A}_{im}}{\mathrm{d}t}& ={\gamma }_{im}{I}_{im}-({\mu }_{im}+d_{im})A_{im}+\sum _{j=1,j\ne i}^n(m_{ij}{A}_{jm}-m_{ji}{A}_{im}),\end{array}$(1) where $i=1,2,\dots ,n$, $N_{if}=S_{if}+I_{if}+A_{if}$, $N_{im}=S_{im}+I_{im}+A_{im}$. Then $N_i(t)=N_{if}(t)+N_{im}(t)$. $N(t)$ is the total population over all patches at time t, then $N(t)=\sum _{i=1}^n{N}_i(t)$. The parameters in this model are all non-negative constants. We summarize the state variables and the model parameters in Table 1. The migration matrix $M=(m_{ij})$ is not required to be symmetric, namely, the migration rate from the patch i to j may not be the same as that from the patch j to i, $0\le m_{ij},m_{ji}\le 1$, where, $i,j=1,2,\dots ,n$. Here the movement of humans between patches is governed by the Euler approach [11], that is, humans change their residences when they move from one patch to another. Furthermore, it is assumed that all parameters in the model are strictly positive with the exception of the migration rates.

Figure 1. The flow diagram of HIV/AIDS transmission in patch i.

The initial conditions of system (1(1) $\begin{array}{rl}\frac{\mathrm{d}{S}_{if}}{\mathrm{d}t}& ={\epsilon }_i{\mathrm{\Lambda }}_i-{\beta }_i\frac{S_{if}{I}_{im}}{N_{im}}-{\mu }_{if}{S}_{if}+\sum _{j=1,j\ne i}^n(m_{ij}{S}_{jf}-m_{ji}{S}_{if}),\\ \frac{\mathrm{d}{I}_{if}}{\mathrm{d}t}& ={\beta }_i\frac{S_{if}{I}_{im}}{N_{im}}-({\mu }_{if}+{\gamma }_{if})I_{if}+\sum _{j=1,j\ne i}^n(m_{ij}{I}_{jf}-m_{ji}{I}_{if}),\\ \frac{\mathrm{d}{A}_{if}}{\mathrm{d}t}& ={\gamma }_{if}{I}_{if}-({\mu }_{if}+d_{if})A_{if}+\sum _{j=1,j\ne i}^n(m_{ij}{A}_{jf}-m_{ji}{A}_{if}),\\ \frac{\mathrm{d}{S}_{im}}{\mathrm{d}t}& =(1-{\epsilon }_i){\mathrm{\Lambda }}_i-{\beta }_i\frac{S_{im}{I}_{if}}{N_{if}}-{\mu }_{im}{S}_{im}+\sum _{j=1,j\ne i}^n(m_{ij}{S}_{jm}-m_{ji}{S}_{im}),\\ \frac{\mathrm{d}{I}_{im}}{\mathrm{d}t}& ={\beta }_i\frac{S_{im}{I}_{if}}{N_{if}}-({\mu }_{im}+{\gamma }_{im})I_{im}+\sum _{j=1,j\ne i}^n(m_{ij}{I}_{jm}-m_{ji}{I}_{im}),\\ \frac{\mathrm{d}{A}_{im}}{\mathrm{d}t}& ={\gamma }_{im}{I}_{im}-({\mu }_{im}+d_{im})A_{im}+\sum _{j=1,j\ne i}^n(m_{ij}{A}_{jm}-m_{ji}{A}_{im}),\end{array}$(1) ) is given by $S_{if}\ge 0,S_{im}\ge 0,I_{if}\ge 0,I_{im}\ge 0,A_{if}\ge 0,A_{im}\ge 0,i=1,2,\dots ,n.$

Table 1. State variables and parameters interpretation of system (1(1) $\begin{array}{rl}\frac{\mathrm{d}{S}_{if}}{\mathrm{d}t}& ={\epsilon }_i{\mathrm{\Lambda }}_i-{\beta }_i\frac{S_{if}{I}_{im}}{N_{im}}-{\mu }_{if}{S}_{if}+\sum _{j=1,j\ne i}^n(m_{ij}{S}_{jf}-m_{ji}{S}_{if}),\\ \frac{\mathrm{d}{I}_{if}}{\mathrm{d}t}& ={\beta }_i\frac{S_{if}{I}_{im}}{N_{im}}-({\mu }_{if}+{\gamma }_{if})I_{if}+\sum _{j=1,j\ne i}^n(m_{ij}{I}_{jf}-m_{ji}{I}_{if}),\\ \frac{\mathrm{d}{A}_{if}}{\mathrm{d}t}& ={\gamma }_{if}{I}_{if}-({\mu }_{if}+d_{if})A_{if}+\sum _{j=1,j\ne i}^n(m_{ij}{A}_{jf}-m_{ji}{A}_{if}),\\ \frac{\mathrm{d}{S}_{im}}{\mathrm{d}t}& =(1-{\epsilon }_i){\mathrm{\Lambda }}_i-{\beta }_i\frac{S_{im}{I}_{if}}{N_{if}}-{\mu }_{im}{S}_{im}+\sum _{j=1,j\ne i}^n(m_{ij}{S}_{jm}-m_{ji}{S}_{im}),\\ \frac{\mathrm{d}{I}_{im}}{\mathrm{d}t}& ={\beta }_i\frac{S_{im}{I}_{if}}{N_{if}}-({\mu }_{im}+{\gamma }_{im})I_{im}+\sum _{j=1,j\ne i}^n(m_{ij}{I}_{jm}-m_{ji}{I}_{im}),\\ \frac{\mathrm{d}{A}_{im}}{\mathrm{d}t}& ={\gamma }_{im}{I}_{im}-({\mu }_{im}+d_{im})A_{im}+\sum _{j=1,j\ne i}^n(m_{ij}{A}_{jm}-m_{ji}{A}_{im}),\end{array}$(1) ).

State variables and parameters	Interpretation
$S_{if}(t)$	The number of susceptible female in patch i at time t
$S_{im}(t)$	The number of susceptible male in patch i at time t
$I_{if}(t)$	The number of infective female in patch i at time t
$I_{im}(t)$	The number of infective male in patch i at time t
$A_{if}(t)$	The number of AIDS female in patch i at time t
$A_{im}(t)$	The number of AIDS male in patch i at time t
${\mathrm{\Lambda }}_i$	The recruitment rate of the population in patch i per unit time
${\epsilon }_i$	The fraction of recruited susceptible $S_{if}$ in patch i
${\beta }_i$	The transmission rate between one susceptible individual and one heterosexual infected individual in patch i
${\mu }_{if}$	The natural death rate of female in patch i
${\mu }_{im}$	The natural death rate of male in patch i
$d_{if}$	The death rate due to AIDS of female in patch i
$d_{im}$	The death rate due to AIDS of male in patch i
${\gamma }_{if}$	The transfer rate from stage $I_{if}$ to $A_{if}$ in patch i
${\gamma }_{im}$	The transfer rate from stage $I_{im}$ to $A_{im}$ in patch i
$m_{ji}$	The migration rate from patch i to j
$m_{ij}$	The migration rate from patch j to i

Theorem 2.1

Solutions of system (1(1) $\begin{array}{rl}\frac{\mathrm{d}{S}_{if}}{\mathrm{d}t}& ={\epsilon }_i{\mathrm{\Lambda }}_i-{\beta }_i\frac{S_{if}{I}_{im}}{N_{im}}-{\mu }_{if}{S}_{if}+\sum _{j=1,j\ne i}^n(m_{ij}{S}_{jf}-m_{ji}{S}_{if}),\\ \frac{\mathrm{d}{I}_{if}}{\mathrm{d}t}& ={\beta }_i\frac{S_{if}{I}_{im}}{N_{im}}-({\mu }_{if}+{\gamma }_{if})I_{if}+\sum _{j=1,j\ne i}^n(m_{ij}{I}_{jf}-m_{ji}{I}_{if}),\\ \frac{\mathrm{d}{A}_{if}}{\mathrm{d}t}& ={\gamma }_{if}{I}_{if}-({\mu }_{if}+d_{if})A_{if}+\sum _{j=1,j\ne i}^n(m_{ij}{A}_{jf}-m_{ji}{A}_{if}),\\ \frac{\mathrm{d}{S}_{im}}{\mathrm{d}t}& =(1-{\epsilon }_i){\mathrm{\Lambda }}_i-{\beta }_i\frac{S_{im}{I}_{if}}{N_{if}}-{\mu }_{im}{S}_{im}+\sum _{j=1,j\ne i}^n(m_{ij}{S}_{jm}-m_{ji}{S}_{im}),\\ \frac{\mathrm{d}{I}_{im}}{\mathrm{d}t}& ={\beta }_i\frac{S_{im}{I}_{if}}{N_{if}}-({\mu }_{im}+{\gamma }_{im})I_{im}+\sum _{j=1,j\ne i}^n(m_{ij}{I}_{jm}-m_{ji}{I}_{im}),\\ \frac{\mathrm{d}{A}_{im}}{\mathrm{d}t}& ={\gamma }_{im}{I}_{im}-({\mu }_{im}+d_{im})A_{im}+\sum _{j=1,j\ne i}^n(m_{ij}{A}_{jm}-m_{ji}{A}_{im}),\end{array}$(1) ) with nonnegative initial conditions uniquely exist and remain nonnegative for all time $t\ge$0. Moreover, the system (1(1) $\begin{array}{rl}\frac{\mathrm{d}{S}_{if}}{\mathrm{d}t}& ={\epsilon }_i{\mathrm{\Lambda }}_i-{\beta }_i\frac{S_{if}{I}_{im}}{N_{im}}-{\mu }_{if}{S}_{if}+\sum _{j=1,j\ne i}^n(m_{ij}{S}_{jf}-m_{ji}{S}_{if}),\\ \frac{\mathrm{d}{I}_{if}}{\mathrm{d}t}& ={\beta }_i\frac{S_{if}{I}_{im}}{N_{im}}-({\mu }_{if}+{\gamma }_{if})I_{if}+\sum _{j=1,j\ne i}^n(m_{ij}{I}_{jf}-m_{ji}{I}_{if}),\\ \frac{\mathrm{d}{A}_{if}}{\mathrm{d}t}& ={\gamma }_{if}{I}_{if}-({\mu }_{if}+d_{if})A_{if}+\sum _{j=1,j\ne i}^n(m_{ij}{A}_{jf}-m_{ji}{A}_{if}),\\ \frac{\mathrm{d}{S}_{im}}{\mathrm{d}t}& =(1-{\epsilon }_i){\mathrm{\Lambda }}_i-{\beta }_i\frac{S_{im}{I}_{if}}{N_{if}}-{\mu }_{im}{S}_{im}+\sum _{j=1,j\ne i}^n(m_{ij}{S}_{jm}-m_{ji}{S}_{im}),\\ \frac{\mathrm{d}{I}_{im}}{\mathrm{d}t}& ={\beta }_i\frac{S_{im}{I}_{if}}{N_{if}}-({\mu }_{im}+{\gamma }_{im})I_{im}+\sum _{j=1,j\ne i}^n(m_{ij}{I}_{jm}-m_{ji}{I}_{im}),\\ \frac{\mathrm{d}{A}_{im}}{\mathrm{d}t}& ={\gamma }_{im}{I}_{im}-({\mu }_{im}+d_{im})A_{im}+\sum _{j=1,j\ne i}^n(m_{ij}{A}_{jm}-m_{ji}{A}_{im}),\end{array}$(1) ) is point dissipative.

Proof.

It is simple to verify that the vector field defined by the right-hand side of (1(1) $\begin{array}{rl}\frac{\mathrm{d}{S}_{if}}{\mathrm{d}t}& ={\epsilon }_i{\mathrm{\Lambda }}_i-{\beta }_i\frac{S_{if}{I}_{im}}{N_{im}}-{\mu }_{if}{S}_{if}+\sum _{j=1,j\ne i}^n(m_{ij}{S}_{jf}-m_{ji}{S}_{if}),\\ \frac{\mathrm{d}{I}_{if}}{\mathrm{d}t}& ={\beta }_i\frac{S_{if}{I}_{im}}{N_{im}}-({\mu }_{if}+{\gamma }_{if})I_{if}+\sum _{j=1,j\ne i}^n(m_{ij}{I}_{jf}-m_{ji}{I}_{if}),\\ \frac{\mathrm{d}{A}_{if}}{\mathrm{d}t}& ={\gamma }_{if}{I}_{if}-({\mu }_{if}+d_{if})A_{if}+\sum _{j=1,j\ne i}^n(m_{ij}{A}_{jf}-m_{ji}{A}_{if}),\\ \frac{\mathrm{d}{S}_{im}}{\mathrm{d}t}& =(1-{\epsilon }_i){\mathrm{\Lambda }}_i-{\beta }_i\frac{S_{im}{I}_{if}}{N_{if}}-{\mu }_{im}{S}_{im}+\sum _{j=1,j\ne i}^n(m_{ij}{S}_{jm}-m_{ji}{S}_{im}),\\ \frac{\mathrm{d}{I}_{im}}{\mathrm{d}t}& ={\beta }_i\frac{S_{im}{I}_{if}}{N_{if}}-({\mu }_{im}+{\gamma }_{im})I_{im}+\sum _{j=1,j\ne i}^n(m_{ij}{I}_{jm}-m_{ji}{I}_{im}),\\ \frac{\mathrm{d}{A}_{im}}{\mathrm{d}t}& ={\gamma }_{im}{I}_{im}-({\mu }_{im}+d_{im})A_{im}+\sum _{j=1,j\ne i}^n(m_{ij}{A}_{jm}-m_{ji}{A}_{im}),\end{array}$(1) ) is locally Lipschitz on $R_{+}^{6n}$, so there exists a unique solution of system (1(1) $\begin{array}{rl}\frac{\mathrm{d}{S}_{if}}{\mathrm{d}t}& ={\epsilon }_i{\mathrm{\Lambda }}_i-{\beta }_i\frac{S_{if}{I}_{im}}{N_{im}}-{\mu }_{if}{S}_{if}+\sum _{j=1,j\ne i}^n(m_{ij}{S}_{jf}-m_{ji}{S}_{if}),\\ \frac{\mathrm{d}{I}_{if}}{\mathrm{d}t}& ={\beta }_i\frac{S_{if}{I}_{im}}{N_{im}}-({\mu }_{if}+{\gamma }_{if})I_{if}+\sum _{j=1,j\ne i}^n(m_{ij}{I}_{jf}-m_{ji}{I}_{if}),\\ \frac{\mathrm{d}{A}_{if}}{\mathrm{d}t}& ={\gamma }_{if}{I}_{if}-({\mu }_{if}+d_{if})A_{if}+\sum _{j=1,j\ne i}^n(m_{ij}{A}_{jf}-m_{ji}{A}_{if}),\\ \frac{\mathrm{d}{S}_{im}}{\mathrm{d}t}& =(1-{\epsilon }_i){\mathrm{\Lambda }}_i-{\beta }_i\frac{S_{im}{I}_{if}}{N_{if}}-{\mu }_{im}{S}_{im}+\sum _{j=1,j\ne i}^n(m_{ij}{S}_{jm}-m_{ji}{S}_{im}),\\ \frac{\mathrm{d}{I}_{im}}{\mathrm{d}t}& ={\beta }_i\frac{S_{im}{I}_{if}}{N_{if}}-({\mu }_{im}+{\gamma }_{im})I_{im}+\sum _{j=1,j\ne i}^n(m_{ij}{I}_{jm}-m_{ji}{I}_{im}),\\ \frac{\mathrm{d}{A}_{im}}{\mathrm{d}t}& ={\gamma }_{im}{I}_{im}-({\mu }_{im}+d_{im})A_{im}+\sum _{j=1,j\ne i}^n(m_{ij}{A}_{jm}-m_{ji}{A}_{im}),\end{array}$(1) ) $t\ge 0$. Next, we prove the nonnegativity of the solution. By the first equation (1(1) $\begin{array}{rl}\frac{\mathrm{d}{S}_{if}}{\mathrm{d}t}& ={\epsilon }_i{\mathrm{\Lambda }}_i-{\beta }_i\frac{S_{if}{I}_{im}}{N_{im}}-{\mu }_{if}{S}_{if}+\sum _{j=1,j\ne i}^n(m_{ij}{S}_{jf}-m_{ji}{S}_{if}),\\ \frac{\mathrm{d}{I}_{if}}{\mathrm{d}t}& ={\beta }_i\frac{S_{if}{I}_{im}}{N_{im}}-({\mu }_{if}+{\gamma }_{if})I_{if}+\sum _{j=1,j\ne i}^n(m_{ij}{I}_{jf}-m_{ji}{I}_{if}),\\ \frac{\mathrm{d}{A}_{if}}{\mathrm{d}t}& ={\gamma }_{if}{I}_{if}-({\mu }_{if}+d_{if})A_{if}+\sum _{j=1,j\ne i}^n(m_{ij}{A}_{jf}-m_{ji}{A}_{if}),\\ \frac{\mathrm{d}{S}_{im}}{\mathrm{d}t}& =(1-{\epsilon }_i){\mathrm{\Lambda }}_i-{\beta }_i\frac{S_{im}{I}_{if}}{N_{if}}-{\mu }_{im}{S}_{im}+\sum _{j=1,j\ne i}^n(m_{ij}{S}_{jm}-m_{ji}{S}_{im}),\\ \frac{\mathrm{d}{I}_{im}}{\mathrm{d}t}& ={\beta }_i\frac{S_{im}{I}_{if}}{N_{if}}-({\mu }_{im}+{\gamma }_{im})I_{im}+\sum _{j=1,j\ne i}^n(m_{ij}{I}_{jm}-m_{ji}{I}_{im}),\\ \frac{\mathrm{d}{A}_{im}}{\mathrm{d}t}& ={\gamma }_{im}{I}_{im}-({\mu }_{im}+d_{im})A_{im}+\sum _{j=1,j\ne i}^n(m_{ij}{A}_{jm}-m_{ji}{A}_{im}),\end{array}$(1) ), we have $\frac{\mathrm{d}{S}_{if}}{\mathrm{d}t}>{\epsilon }_i{\mathrm{\Lambda }}_i-{\beta }_i\frac{S_{if}{I}_{im}}{N_{im}}-{\mu }_{if}{S}_{if}-m_i{S}_{if},$where $m_i=\sum _{j=1,j\ne i}^n{m}_{ji}.$Then $\frac{\mathrm{d}{S}_{if}}{\mathrm{d}t}+{\beta }_i\frac{S_{if}{I}_{im}}{N_{im}}+{\mu }_{if}{S}_{if}+m_i{S}_{if}>{\epsilon }_i{\mathrm{\Lambda }}_i.$Let $\lambda (t)={\beta }_i\frac{I_{im}(t)}{N_{im}(t)}$, then $\begin{array}{rl}& \frac{\mathrm{d}{S}_{if}(t)}{\mathrm{d}t}{\mathrm{e}}^{{\int }_0^t\lambda (u)\mathrm{d}u+({\mu }_{if}+m_i)t}+S_{if}(t)(\lambda (t)+{\mu }_{if}+m_i){\mathrm{e}}^{{\int }_0^t\lambda (u)\mathrm{d}u+({\mu }_{if}+m_i)t}\\ & >{\epsilon }_i{\mathrm{\Lambda }}_i{\mathrm{e}}^{{\int }_0^t\lambda (u)\mathrm{d}u+({\mu }_{if}+m_i)t},\end{array}$namely $\frac{\mathrm{d}}{\mathrm{d}t}[S_{if}(t){\mathrm{e}}^{{\int }_0^t\lambda (u)\mathrm{d}u+({\mu }_{if}+m_i)t}]>{\epsilon }_i{\mathrm{\Lambda }}_i{\mathrm{e}}^{{\int }_0^t\lambda (u)\mathrm{d}u+({\mu }_{if}+m_i)t}.$Then $S_{if}(t){\mathrm{e}}^{{\int }_0^t\lambda (u)\mathrm{d}u+({\mu }_{if}+m_i)t}-S_{if}(0)>{\int }_0^t{\epsilon }_i{\mathrm{\Lambda }}_i{\mathrm{e}}^{{\int }_0^u\lambda (v)dv+({\mu }_{if}+m_i)u}\mathrm{d}u.$Hence, we have $\begin{array}{rl}{S}_{if}(t)& =S_{if}(0){\mathrm{e}}^{-{\int }_0^t\lambda (u)\mathrm{d}u+({\mu }_{if}+m_i)t}+{\mathrm{e}}^{-{\int }_0^t\lambda (u)\mathrm{d}u+({\mu }_{if}+m_i)t}\\ & \times {\int }_0^t{\epsilon }_i{\mathrm{\Lambda }}_i{\mathrm{e}}^{{\int }_0^u\lambda (v)dv+({\mu }_{if}+m_i)u}\mathrm{d}u\ge 0.\end{array}$Similarly, $I_{if}(t)\ge 0,A_{if}(t)\ge 0,S_{im}(t)\ge 0,I_{im}(t)\ge 0,A_{im}(t)\ge 0$.

Let $\mu =\min \{{\mu }_{1f},{\mu }_{2f},\cdots ,{\mu }_{nf},{\mu }_{1m},{\mu }_{2m},\cdots ,{\mu }_{nm}\}$ and $\overline{\mathrm{\Lambda }}=\sum _{i=1}^n{\mathrm{\Lambda }}_i$. The total population over all patches, denoted by $N=\sum _{i=1}^n{N}_i$, satisfies $\frac{\mathrm{d}N}{\mathrm{d}t}=\sum _{i=1}^n({\mathrm{\Lambda }}_i-{\mu }_{if}{N}_{if}-{\mu }_{im}{N}_{im}-d_{if}{A}_{if}-d_{im}{A}_{im})\le \overline{\mathrm{\Lambda }}-\mu N,$which implies that $\underset{t\to \mathrm{\infty }}{lim}supN\le \frac{\overline{\mathrm{\Lambda }}}{\mu }$. Thus the system (1(1) $\begin{array}{rl}\frac{\mathrm{d}{S}_{if}}{\mathrm{d}t}& ={\epsilon }_i{\mathrm{\Lambda }}_i-{\beta }_i\frac{S_{if}{I}_{im}}{N_{im}}-{\mu }_{if}{S}_{if}+\sum _{j=1,j\ne i}^n(m_{ij}{S}_{jf}-m_{ji}{S}_{if}),\\ \frac{\mathrm{d}{I}_{if}}{\mathrm{d}t}& ={\beta }_i\frac{S_{if}{I}_{im}}{N_{im}}-({\mu }_{if}+{\gamma }_{if})I_{if}+\sum _{j=1,j\ne i}^n(m_{ij}{I}_{jf}-m_{ji}{I}_{if}),\\ \frac{\mathrm{d}{A}_{if}}{\mathrm{d}t}& ={\gamma }_{if}{I}_{if}-({\mu }_{if}+d_{if})A_{if}+\sum _{j=1,j\ne i}^n(m_{ij}{A}_{jf}-m_{ji}{A}_{if}),\\ \frac{\mathrm{d}{S}_{im}}{\mathrm{d}t}& =(1-{\epsilon }_i){\mathrm{\Lambda }}_i-{\beta }_i\frac{S_{im}{I}_{if}}{N_{if}}-{\mu }_{im}{S}_{im}+\sum _{j=1,j\ne i}^n(m_{ij}{S}_{jm}-m_{ji}{S}_{im}),\\ \frac{\mathrm{d}{I}_{im}}{\mathrm{d}t}& ={\beta }_i\frac{S_{im}{I}_{if}}{N_{if}}-({\mu }_{im}+{\gamma }_{im})I_{im}+\sum _{j=1,j\ne i}^n(m_{ij}{I}_{jm}-m_{ji}{I}_{im}),\\ \frac{\mathrm{d}{A}_{im}}{\mathrm{d}t}& ={\gamma }_{im}{I}_{im}-({\mu }_{im}+d_{im})A_{im}+\sum _{j=1,j\ne i}^n(m_{ij}{A}_{jm}-m_{ji}{A}_{im}),\end{array}$(1) ) is point dissipative.

Let $\begin{array}{rl}A& =\left[\begin{array}{cccc}{\mu }_{1f}+{\displaystyle \sum _{j\ne 1}{m}_{j1}}& -m_{12}& \cdots & -m_{1n}\\ -m_{21}& {\mu }_{2f}+{\displaystyle \sum _{j\ne 2}{m}_{j2}}& \cdots & -m_{2n}\\ ⋮& ⋮& \ddots & ⋮\\ -m_{n1}& -m_{n2}& \cdots & {\mu }_{nf}+{\displaystyle \sum _{j\ne n}{m}_{jn}}\end{array}\right],\\ B& =\left[\begin{array}{cccc}{\mu }_{1m}+{\displaystyle \sum _{j\ne 1}{m}_{j1}}& -m_{12}& \cdots & -m_{1n}\\ -m_{21}& {\mu }_{2m}+{\displaystyle \sum _{j\ne 2}{m}_{j2}}& \cdots & -m_{2n}\\ ⋮& ⋮& \ddots & ⋮\\ -m_{n1}& -m_{n2}& \cdots & {\mu }_{nm}+{\displaystyle \sum _{j\ne n}{m}_{jn}}\end{array}\right],\end{array}$since all off-diagonal entries of A and B are non-positive, it follows from the first and the forth equations of (1(1) $\begin{array}{rl}\frac{\mathrm{d}{S}_{if}}{\mathrm{d}t}& ={\epsilon }_i{\mathrm{\Lambda }}_i-{\beta }_i\frac{S_{if}{I}_{im}}{N_{im}}-{\mu }_{if}{S}_{if}+\sum _{j=1,j\ne i}^n(m_{ij}{S}_{jf}-m_{ji}{S}_{if}),\\ \frac{\mathrm{d}{I}_{if}}{\mathrm{d}t}& ={\beta }_i\frac{S_{if}{I}_{im}}{N_{im}}-({\mu }_{if}+{\gamma }_{if})I_{if}+\sum _{j=1,j\ne i}^n(m_{ij}{I}_{jf}-m_{ji}{I}_{if}),\\ \frac{\mathrm{d}{A}_{if}}{\mathrm{d}t}& ={\gamma }_{if}{I}_{if}-({\mu }_{if}+d_{if})A_{if}+\sum _{j=1,j\ne i}^n(m_{ij}{A}_{jf}-m_{ji}{A}_{if}),\\ \frac{\mathrm{d}{S}_{im}}{\mathrm{d}t}& =(1-{\epsilon }_i){\mathrm{\Lambda }}_i-{\beta }_i\frac{S_{im}{I}_{if}}{N_{if}}-{\mu }_{im}{S}_{im}+\sum _{j=1,j\ne i}^n(m_{ij}{S}_{jm}-m_{ji}{S}_{im}),\\ \frac{\mathrm{d}{I}_{im}}{\mathrm{d}t}& ={\beta }_i\frac{S_{im}{I}_{if}}{N_{if}}-({\mu }_{im}+{\gamma }_{im})I_{im}+\sum _{j=1,j\ne i}^n(m_{ij}{I}_{jm}-m_{ji}{I}_{im}),\\ \frac{\mathrm{d}{A}_{im}}{\mathrm{d}t}& ={\gamma }_{im}{I}_{im}-({\mu }_{im}+d_{im})A_{im}+\sum _{j=1,j\ne i}^n(m_{ij}{A}_{jm}-m_{ji}{A}_{im}),\end{array}$(1) ) that $\begin{array}{rl}\frac{\mathrm{d}{S}_{if}}{\mathrm{d}t}& \le {\epsilon }_i{\mathrm{\Lambda }}_i-{\mu }_{if}{S}_{if}+\sum _{j=1,j\ne i}^n(m_{ij}{S}_{jf}-m_{ji}{S}_{if})=(A{S}_f^0-A{S}_f{)}_i\le 0,\\ \frac{\mathrm{d}{S}_{im}}{\mathrm{d}t}& \le (1-{\epsilon }_i){\mathrm{\Lambda }}_i-{\mu }_{im}{S}_{im}+\sum _{j=1,j\ne i}^n(m_{ij}{S}_{jm}-m_{ji}{S}_{im})=(A{S}_m^0-A{S}_m{)}_i\le 0,\end{array}$when $S_{if}=S_{if}^0,S_{jf}\le S_{jf}^0,S_{im}=S_{im}^0,S_{jm}\le S_{jm}^0$ for $j\ne i$. Thus the feasible region of (1(1) $\begin{array}{rl}\frac{\mathrm{d}{S}_{if}}{\mathrm{d}t}& ={\epsilon }_i{\mathrm{\Lambda }}_i-{\beta }_i\frac{S_{if}{I}_{im}}{N_{im}}-{\mu }_{if}{S}_{if}+\sum _{j=1,j\ne i}^n(m_{ij}{S}_{jf}-m_{ji}{S}_{if}),\\ \frac{\mathrm{d}{I}_{if}}{\mathrm{d}t}& ={\beta }_i\frac{S_{if}{I}_{im}}{N_{im}}-({\mu }_{if}+{\gamma }_{if})I_{if}+\sum _{j=1,j\ne i}^n(m_{ij}{I}_{jf}-m_{ji}{I}_{if}),\\ \frac{\mathrm{d}{A}_{if}}{\mathrm{d}t}& ={\gamma }_{if}{I}_{if}-({\mu }_{if}+d_{if})A_{if}+\sum _{j=1,j\ne i}^n(m_{ij}{A}_{jf}-m_{ji}{A}_{if}),\\ \frac{\mathrm{d}{S}_{im}}{\mathrm{d}t}& =(1-{\epsilon }_i){\mathrm{\Lambda }}_i-{\beta }_i\frac{S_{im}{I}_{if}}{N_{if}}-{\mu }_{im}{S}_{im}+\sum _{j=1,j\ne i}^n(m_{ij}{S}_{jm}-m_{ji}{S}_{im}),\\ \frac{\mathrm{d}{I}_{im}}{\mathrm{d}t}& ={\beta }_i\frac{S_{im}{I}_{if}}{N_{if}}-({\mu }_{im}+{\gamma }_{im})I_{im}+\sum _{j=1,j\ne i}^n(m_{ij}{I}_{jm}-m_{ji}{I}_{im}),\\ \frac{\mathrm{d}{A}_{im}}{\mathrm{d}t}& ={\gamma }_{im}{I}_{im}-({\mu }_{im}+d_{im})A_{im}+\sum _{j=1,j\ne i}^n(m_{ij}{A}_{jm}-m_{ji}{A}_{im}),\end{array}$(1) ) can be chosen as $\begin{array}{rl}& X=\{S_{1f},I_{1f},A_{1f},\dots ,S_{nf},I_{nf},A_{nf},S_{1m},I_{1m},A_{1m},\dots ,S_{nm},I_{nm},A_{nm}\in {\mathbb{R}}_{+}^{6n}|N\\ & \le \frac{\overline{\mathrm{\Lambda }}}{\mu },S_{if}\le S_{if}^0,S_{im}\le S_{im}^0,1\le i\le n\}.\end{array}$It can be verified that X is positively invariant with respect to (1(1) $\begin{array}{rl}\frac{\mathrm{d}{S}_{if}}{\mathrm{d}t}& ={\epsilon }_i{\mathrm{\Lambda }}_i-{\beta }_i\frac{S_{if}{I}_{im}}{N_{im}}-{\mu }_{if}{S}_{if}+\sum _{j=1,j\ne i}^n(m_{ij}{S}_{jf}-m_{ji}{S}_{if}),\\ \frac{\mathrm{d}{I}_{if}}{\mathrm{d}t}& ={\beta }_i\frac{S_{if}{I}_{im}}{N_{im}}-({\mu }_{if}+{\gamma }_{if})I_{if}+\sum _{j=1,j\ne i}^n(m_{ij}{I}_{jf}-m_{ji}{I}_{if}),\\ \frac{\mathrm{d}{A}_{if}}{\mathrm{d}t}& ={\gamma }_{if}{I}_{if}-({\mu }_{if}+d_{if})A_{if}+\sum _{j=1,j\ne i}^n(m_{ij}{A}_{jf}-m_{ji}{A}_{if}),\\ \frac{\mathrm{d}{S}_{im}}{\mathrm{d}t}& =(1-{\epsilon }_i){\mathrm{\Lambda }}_i-{\beta }_i\frac{S_{im}{I}_{if}}{N_{if}}-{\mu }_{im}{S}_{im}+\sum _{j=1,j\ne i}^n(m_{ij}{S}_{jm}-m_{ji}{S}_{im}),\\ \frac{\mathrm{d}{I}_{im}}{\mathrm{d}t}& ={\beta }_i\frac{S_{im}{I}_{if}}{N_{if}}-({\mu }_{im}+{\gamma }_{im})I_{im}+\sum _{j=1,j\ne i}^n(m_{ij}{I}_{jm}-m_{ji}{I}_{im}),\\ \frac{\mathrm{d}{A}_{im}}{\mathrm{d}t}& ={\gamma }_{im}{I}_{im}-({\mu }_{im}+d_{im})A_{im}+\sum _{j=1,j\ne i}^n(m_{ij}{A}_{jm}-m_{ji}{A}_{im}),\end{array}$(1) ).

3. Mathematical analysis

In this section, we first study the existence and uniqueness of the disease-free equilibrium for (1(1) $\begin{array}{rl}\frac{\mathrm{d}{S}_{if}}{\mathrm{d}t}& ={\epsilon }_i{\mathrm{\Lambda }}_i-{\beta }_i\frac{S_{if}{I}_{im}}{N_{im}}-{\mu }_{if}{S}_{if}+\sum _{j=1,j\ne i}^n(m_{ij}{S}_{jf}-m_{ji}{S}_{if}),\\ \frac{\mathrm{d}{I}_{if}}{\mathrm{d}t}& ={\beta }_i\frac{S_{if}{I}_{im}}{N_{im}}-({\mu }_{if}+{\gamma }_{if})I_{if}+\sum _{j=1,j\ne i}^n(m_{ij}{I}_{jf}-m_{ji}{I}_{if}),\\ \frac{\mathrm{d}{A}_{if}}{\mathrm{d}t}& ={\gamma }_{if}{I}_{if}-({\mu }_{if}+d_{if})A_{if}+\sum _{j=1,j\ne i}^n(m_{ij}{A}_{jf}-m_{ji}{A}_{if}),\\ \frac{\mathrm{d}{S}_{im}}{\mathrm{d}t}& =(1-{\epsilon }_i){\mathrm{\Lambda }}_i-{\beta }_i\frac{S_{im}{I}_{if}}{N_{if}}-{\mu }_{im}{S}_{im}+\sum _{j=1,j\ne i}^n(m_{ij}{S}_{jm}-m_{ji}{S}_{im}),\\ \frac{\mathrm{d}{I}_{im}}{\mathrm{d}t}& ={\beta }_i\frac{S_{im}{I}_{if}}{N_{if}}-({\mu }_{im}+{\gamma }_{im})I_{im}+\sum _{j=1,j\ne i}^n(m_{ij}{I}_{jm}-m_{ji}{I}_{im}),\\ \frac{\mathrm{d}{A}_{im}}{\mathrm{d}t}& ={\gamma }_{im}{I}_{im}-({\mu }_{im}+d_{im})A_{im}+\sum _{j=1,j\ne i}^n(m_{ij}{A}_{jm}-m_{ji}{A}_{im}),\end{array}$(1) ) and calculate the basic reproduction number $R_0$. Then, we study the stability of the disease-free equilibrium. Finally, we study the existence, uniqueness and stability of the endemic equilibrium.

3.1. Disease-free equilibrium and the basic reproduction number

To find the disease-free equilibrium of (1(1) $\begin{array}{rl}\frac{\mathrm{d}{S}_{if}}{\mathrm{d}t}& ={\epsilon }_i{\mathrm{\Lambda }}_i-{\beta }_i\frac{S_{if}{I}_{im}}{N_{im}}-{\mu }_{if}{S}_{if}+\sum _{j=1,j\ne i}^n(m_{ij}{S}_{jf}-m_{ji}{S}_{if}),\\ \frac{\mathrm{d}{I}_{if}}{\mathrm{d}t}& ={\beta }_i\frac{S_{if}{I}_{im}}{N_{im}}-({\mu }_{if}+{\gamma }_{if})I_{if}+\sum _{j=1,j\ne i}^n(m_{ij}{I}_{jf}-m_{ji}{I}_{if}),\\ \frac{\mathrm{d}{A}_{if}}{\mathrm{d}t}& ={\gamma }_{if}{I}_{if}-({\mu }_{if}+d_{if})A_{if}+\sum _{j=1,j\ne i}^n(m_{ij}{A}_{jf}-m_{ji}{A}_{if}),\\ \frac{\mathrm{d}{S}_{im}}{\mathrm{d}t}& =(1-{\epsilon }_i){\mathrm{\Lambda }}_i-{\beta }_i\frac{S_{im}{I}_{if}}{N_{if}}-{\mu }_{im}{S}_{im}+\sum _{j=1,j\ne i}^n(m_{ij}{S}_{jm}-m_{ji}{S}_{im}),\\ \frac{\mathrm{d}{I}_{im}}{\mathrm{d}t}& ={\beta }_i\frac{S_{im}{I}_{if}}{N_{if}}-({\mu }_{im}+{\gamma }_{im})I_{im}+\sum _{j=1,j\ne i}^n(m_{ij}{I}_{jm}-m_{ji}{I}_{im}),\\ \frac{\mathrm{d}{A}_{im}}{\mathrm{d}t}& ={\gamma }_{im}{I}_{im}-({\mu }_{im}+d_{im})A_{im}+\sum _{j=1,j\ne i}^n(m_{ij}{A}_{jm}-m_{ji}{A}_{im}),\end{array}$(1) ), we consider the following linear system: $\begin{array}{rl}A{S}_f& =\epsilon \mathrm{\Lambda },\\ B{S}_m& =(1-\epsilon )\mathrm{\Lambda },\end{array}$where $S_f=(S_{1f},S_{2f},\dots ,S_{nf}{)}^{\mathrm{T}}$, $S_m=(S_{1m},S_{2m},\dots ,S_{nm}{)}^{\mathrm{T}}$, $\epsilon =\mathrm{diag}({\epsilon }_1,{\epsilon }_2,\dots ,{\epsilon }_n)$, $\mathrm{\Lambda }=({\mathrm{\Lambda }}_1,{\mathrm{\Lambda }}_2,\dots ,{\mathrm{\Lambda }}_n{)}^{\mathrm{T}}$. Since all off-diagonal entries of A are non-positive and the sum of the entries in each column of A is positive, so A is a nonsingular M-matrix and $A^{-1}\ge 0$ [9]. Hence, linear system $A{S}_f=\epsilon \mathrm{\Lambda }$ has a unique positive solution $S_f^0=(S_{1f}^0,S_{2f}^0,\dots ,S_{nf}^0{)}^{\mathrm{T}}=A^{-1}\epsilon \mathrm{\Lambda }>0$. Mathematically, if $I_{if}=0$ for all i at a steady state, then by summing the third equation of (1(1) $\begin{array}{rl}\frac{\mathrm{d}{S}_{if}}{\mathrm{d}t}& ={\epsilon }_i{\mathrm{\Lambda }}_i-{\beta }_i\frac{S_{if}{I}_{im}}{N_{im}}-{\mu }_{if}{S}_{if}+\sum _{j=1,j\ne i}^n(m_{ij}{S}_{jf}-m_{ji}{S}_{if}),\\ \frac{\mathrm{d}{I}_{if}}{\mathrm{d}t}& ={\beta }_i\frac{S_{if}{I}_{im}}{N_{im}}-({\mu }_{if}+{\gamma }_{if})I_{if}+\sum _{j=1,j\ne i}^n(m_{ij}{I}_{jf}-m_{ji}{I}_{if}),\\ \frac{\mathrm{d}{A}_{if}}{\mathrm{d}t}& ={\gamma }_{if}{I}_{if}-({\mu }_{if}+d_{if})A_{if}+\sum _{j=1,j\ne i}^n(m_{ij}{A}_{jf}-m_{ji}{A}_{if}),\\ \frac{\mathrm{d}{S}_{im}}{\mathrm{d}t}& =(1-{\epsilon }_i){\mathrm{\Lambda }}_i-{\beta }_i\frac{S_{im}{I}_{if}}{N_{if}}-{\mu }_{im}{S}_{im}+\sum _{j=1,j\ne i}^n(m_{ij}{S}_{jm}-m_{ji}{S}_{im}),\\ \frac{\mathrm{d}{I}_{im}}{\mathrm{d}t}& ={\beta }_i\frac{S_{im}{I}_{if}}{N_{if}}-({\mu }_{im}+{\gamma }_{im})I_{im}+\sum _{j=1,j\ne i}^n(m_{ij}{I}_{jm}-m_{ji}{I}_{im}),\\ \frac{\mathrm{d}{A}_{im}}{\mathrm{d}t}& ={\gamma }_{im}{I}_{im}-({\mu }_{im}+d_{im})A_{im}+\sum _{j=1,j\ne i}^n(m_{ij}{A}_{jm}-m_{ji}{A}_{im}),\end{array}$(1) ) up from 1 to n, we have $-\sum _{i=1}^n({\mu }_{if}+d_{if})A_{if}+\sum _{i=1}^n\sum _{j=1}^n(m_{ij}{A}_{jf}-m_{ji}{A}_{if})=0.$Hence, $-\sum _{i=1}^n({\mu }_{if}+d_{if})A_{if}=0$. This implies $A_{if}=0$ for $i=1,2,\dots ,n$. In the same way, linear system $B{S}_m=(1-\epsilon )\mathrm{\Lambda }$ has also a unique positive solution $S_m^0=(S_{1m}^0,S_{2m}^0,\dots ,S_{nm}^0{)}^{\mathrm{T}}=B^{-1}(1-\epsilon )\mathrm{\Lambda }>0$ and $A_{im}=0$ for $i=1,2,\dots ,n$. As a consequence, system (1(1) $\begin{array}{rl}\frac{\mathrm{d}{S}_{if}}{\mathrm{d}t}& ={\epsilon }_i{\mathrm{\Lambda }}_i-{\beta }_i\frac{S_{if}{I}_{im}}{N_{im}}-{\mu }_{if}{S}_{if}+\sum _{j=1,j\ne i}^n(m_{ij}{S}_{jf}-m_{ji}{S}_{if}),\\ \frac{\mathrm{d}{I}_{if}}{\mathrm{d}t}& ={\beta }_i\frac{S_{if}{I}_{im}}{N_{im}}-({\mu }_{if}+{\gamma }_{if})I_{if}+\sum _{j=1,j\ne i}^n(m_{ij}{I}_{jf}-m_{ji}{I}_{if}),\\ \frac{\mathrm{d}{A}_{if}}{\mathrm{d}t}& ={\gamma }_{if}{I}_{if}-({\mu }_{if}+d_{if})A_{if}+\sum _{j=1,j\ne i}^n(m_{ij}{A}_{jf}-m_{ji}{A}_{if}),\\ \frac{\mathrm{d}{S}_{im}}{\mathrm{d}t}& =(1-{\epsilon }_i){\mathrm{\Lambda }}_i-{\beta }_i\frac{S_{im}{I}_{if}}{N_{if}}-{\mu }_{im}{S}_{im}+\sum _{j=1,j\ne i}^n(m_{ij}{S}_{jm}-m_{ji}{S}_{im}),\\ \frac{\mathrm{d}{I}_{im}}{\mathrm{d}t}& ={\beta }_i\frac{S_{im}{I}_{if}}{N_{if}}-({\mu }_{im}+{\gamma }_{im})I_{im}+\sum _{j=1,j\ne i}^n(m_{ij}{I}_{jm}-m_{ji}{I}_{im}),\\ \frac{\mathrm{d}{A}_{im}}{\mathrm{d}t}& ={\gamma }_{im}{I}_{im}-({\mu }_{im}+d_{im})A_{im}+\sum _{j=1,j\ne i}^n(m_{ij}{A}_{jm}-m_{ji}{A}_{im}),\end{array}$(1) ) has a unique disease-free equilibrium $E_0=(S_{1f}^0,0,0,S_{1m}^0,0,0,\dots ,S_{nf}^0,0,0,S_{nm}^0,0,0)$.

To derive the basic reproduction number for (1(1) $\begin{array}{rl}\frac{\mathrm{d}{S}_{if}}{\mathrm{d}t}& ={\epsilon }_i{\mathrm{\Lambda }}_i-{\beta }_i\frac{S_{if}{I}_{im}}{N_{im}}-{\mu }_{if}{S}_{if}+\sum _{j=1,j\ne i}^n(m_{ij}{S}_{jf}-m_{ji}{S}_{if}),\\ \frac{\mathrm{d}{I}_{if}}{\mathrm{d}t}& ={\beta }_i\frac{S_{if}{I}_{im}}{N_{im}}-({\mu }_{if}+{\gamma }_{if})I_{if}+\sum _{j=1,j\ne i}^n(m_{ij}{I}_{jf}-m_{ji}{I}_{if}),\\ \frac{\mathrm{d}{A}_{if}}{\mathrm{d}t}& ={\gamma }_{if}{I}_{if}-({\mu }_{if}+d_{if})A_{if}+\sum _{j=1,j\ne i}^n(m_{ij}{A}_{jf}-m_{ji}{A}_{if}),\\ \frac{\mathrm{d}{S}_{im}}{\mathrm{d}t}& =(1-{\epsilon }_i){\mathrm{\Lambda }}_i-{\beta }_i\frac{S_{im}{I}_{if}}{N_{if}}-{\mu }_{im}{S}_{im}+\sum _{j=1,j\ne i}^n(m_{ij}{S}_{jm}-m_{ji}{S}_{im}),\\ \frac{\mathrm{d}{I}_{im}}{\mathrm{d}t}& ={\beta }_i\frac{S_{im}{I}_{if}}{N_{if}}-({\mu }_{im}+{\gamma }_{im})I_{im}+\sum _{j=1,j\ne i}^n(m_{ij}{I}_{jm}-m_{ji}{I}_{im}),\\ \frac{\mathrm{d}{A}_{im}}{\mathrm{d}t}& ={\gamma }_{im}{I}_{im}-({\mu }_{im}+d_{im})A_{im}+\sum _{j=1,j\ne i}^n(m_{ij}{A}_{jm}-m_{ji}{A}_{im}),\end{array}$(1) ), we use the next generation operator approach as described by Diekmann [13] and van den Driessche and Watmough [31]. Let $x=(I_{1f},\dots ,I_{nf},I_{1m},\dots ,I_{nm},A_{1f},\dots ,A_{nf},A_{1m},\dots ,A_{nm}),$then system (1(1) $\begin{array}{rl}\frac{\mathrm{d}{S}_{if}}{\mathrm{d}t}& ={\epsilon }_i{\mathrm{\Lambda }}_i-{\beta }_i\frac{S_{if}{I}_{im}}{N_{im}}-{\mu }_{if}{S}_{if}+\sum _{j=1,j\ne i}^n(m_{ij}{S}_{jf}-m_{ji}{S}_{if}),\\ \frac{\mathrm{d}{I}_{if}}{\mathrm{d}t}& ={\beta }_i\frac{S_{if}{I}_{im}}{N_{im}}-({\mu }_{if}+{\gamma }_{if})I_{if}+\sum _{j=1,j\ne i}^n(m_{ij}{I}_{jf}-m_{ji}{I}_{if}),\\ \frac{\mathrm{d}{A}_{if}}{\mathrm{d}t}& ={\gamma }_{if}{I}_{if}-({\mu }_{if}+d_{if})A_{if}+\sum _{j=1,j\ne i}^n(m_{ij}{A}_{jf}-m_{ji}{A}_{if}),\\ \frac{\mathrm{d}{S}_{im}}{\mathrm{d}t}& =(1-{\epsilon }_i){\mathrm{\Lambda }}_i-{\beta }_i\frac{S_{im}{I}_{if}}{N_{if}}-{\mu }_{im}{S}_{im}+\sum _{j=1,j\ne i}^n(m_{ij}{S}_{jm}-m_{ji}{S}_{im}),\\ \frac{\mathrm{d}{I}_{im}}{\mathrm{d}t}& ={\beta }_i\frac{S_{im}{I}_{if}}{N_{if}}-({\mu }_{im}+{\gamma }_{im})I_{im}+\sum _{j=1,j\ne i}^n(m_{ij}{I}_{jm}-m_{ji}{I}_{im}),\\ \frac{\mathrm{d}{A}_{im}}{\mathrm{d}t}& ={\gamma }_{im}{I}_{im}-({\mu }_{im}+d_{im})A_{im}+\sum _{j=1,j\ne i}^n(m_{ij}{A}_{jm}-m_{ji}{A}_{im}),\end{array}$(1) ) can be written as $\frac{\mathrm{d}x}{\mathrm{d}t}=\mathcal{F}(x)-\mathcal{V}(x),$where $\begin{array}{rl}\mathcal{F}(x)& ={({\beta }_1\frac{S_{1f}{I}_{1m}}{N_{1m}},\dots ,{\beta }_n\frac{S_{nf}{I}_{nm}}{N_{nm}},{\beta }_1\frac{S_{1m}{I}_{1f}}{N_{1f}},\dots ,{\beta }_n\frac{S_{nm}{I}_{nf}}{N_{nf}},0,\dots ,0,0,\dots ,0)}^{\mathrm{T}},\\ \mathcal{V}(x)& =\left[\begin{array}{c}({\mu }_{1f}+{\gamma }_{1f}+{\displaystyle \sum _{j\ne 1}{m}_{j1})I_{1f}-{\displaystyle \sum _{j\ne 1}{m}_{1j}{I}_{jf}}}\\ ⋮\\ ({\mu }_{nf}+{\gamma }_{nf}+{\displaystyle \sum _{j\ne n}{m}_{jn})I_{nf}-{\displaystyle \sum _{j\ne n}{m}_{nj}{I}_{jf}}}\\ ({\mu }_{1m}+{\gamma }_{1m}+{\displaystyle \sum _{j\ne 1}{m}_{j1})I_{1m}-{\displaystyle \sum _{j\ne 1}{m}_{1j}{I}_{jm}}}\\ ⋮\\ ({\mu }_{nm}+{\gamma }_{nm}+{\displaystyle \sum _{j\ne n}{m}_{jn})I_{nm}-{\displaystyle \sum _{j\ne n}{m}_{nj}{I}_{jm}}}\\ -{\gamma }_{1f}{I}_{1f}+({\mu }_{1f}+d_{1f}+{\displaystyle \sum _{j\ne 1}{m}_{j1})A_{1f}-{\displaystyle \sum _{j\ne 1}{m}_{1j}{A}_{jf}}}\\ ⋮\\ -{\gamma }_{nf}{I}_{1f}+({\mu }_{nf}+d_{nf}+{\displaystyle \sum _{j\ne n}{m}_{jn})A_{nf}-{\displaystyle \sum _{j\ne n}{m}_{nj}{A}_{jf}}}\\ -{\gamma }_{1m}{I}_{1m}+({\mu }_{1m}+d_{1m}+{\displaystyle \sum _{j\ne 1}{m}_{j1})A_{1m}-{\displaystyle \sum _{j\ne 1}{m}_{1j}{A}_{jm}}}\\ ⋮\\ -{\gamma }_{nm}{I}_{1m}+({\mu }_{nm}+d_{nm}+{\displaystyle \sum _{j\ne n}{m}_{jn})A_{nm}-{\displaystyle \sum _{j\ne n}{m}_{nj}{A}_{jm}}}\end{array}\right].\end{array}$Then the new infection and transition matrices are, respectively, given by $F=\left[\begin{array}{cccc}0& F_{12}& 0& 0\\ F_{21}& 0& 0& 0\\ 0& 0& 0& 0\\ 0& 0& 0& 0\end{array}\right],V=\left[\begin{array}{cccc}{V}_{11}& 0& 0& 0\\ 0& V_{22}& 0& 0\\ V_{31}& 0& V_{33}& 0\\ 0& V_{42}& 0& V_{44}\end{array}\right],$where $\begin{array}{rl}{F}_{12}& =\mathrm{diag}({\beta }_1\frac{S_{1f}^0}{S_{1m}^0},{\beta }_2\frac{S_{2f}^0}{S_{2m}^0},\dots ,{\beta }_n\frac{S_{nf}^0}{S_{nm}^0}),F_{21}=\mathrm{diag}({\beta }_1\frac{S_{1m}^0}{S_{1f}^0},{\beta }_2\frac{S_{2m}^0}{S_{2f}^0},\dots ,{\beta }_n\frac{S_{nm}^0}{S_{nf}^0}),\\ V_{11}& =\left[\begin{array}{cccc}{\mu }_{1f}+{\gamma }_{1f}+\sum _{j\ne 1}{m}_{j1}& -m_{12}& \cdots & -m_{1n}\\ -m_{21}& {\mu }_{2f}+{\gamma }_{2f}+\sum _{j\ne 2}{m}_{j2}& \cdots & -m_{2n}\\ ⋮& ⋮& \ddots & ⋮\\ -m_{n1}& -m_{n2}& \cdots & {\mu }_{nf}+{\gamma }_{nf}+\sum _{j\ne n}{m}_{jn}\end{array}\right],\\ V_{22}& =\left[\begin{array}{cccc}{\mu }_{1m}+{\gamma }_{1m}+\sum _{j\ne 1}{m}_{j1}& -m_{12}& \cdots & -m_{1n}\\ -m_{21}& {\mu }_{2m}+{\gamma }_{2m}+\sum _{j\ne 2}{m}_{j2}& \cdots & -m_{2n}\\ ⋮& ⋮& \ddots & ⋮\\ -m_{n1}& -m_{n2}& \cdots & {\mu }_{nm}+{\gamma }_{nm}+\sum _{j\ne n}{m}_{jn}\end{array}\right],\\ V_{33}& =\left[\begin{array}{cccc}{\mu }_{1f}+d_{1f}+\sum _{j\ne 1}{m}_{j1}& -m_{12}& \cdots & -m_{1n}\\ -m_{21}& {\mu }_{2f}+d_{2f}+\sum _{j\ne 2}{m}_{j2}& \cdots & -m_{2n}\\ ⋮& ⋮& \ddots & ⋮\\ -m_{n1}& -m_{n2}& \cdots & {\mu }_{nf}+d_{nf}+\sum _{j\ne n}{m}_{jn}\end{array}\right],\\ V_{44}& =\left[\begin{array}{cccc}{\mu }_{1m}+d_{1m}+\sum _{j\ne 1}{m}_{j1}& -m_{12}& \cdots & -m_{1n}\\ -m_{21}& {\mu }_{2m}+d_{2m}+\sum _{j\ne 2}{m}_{j2}& \cdots & -m_{2n}\\ ⋮& ⋮& \ddots & ⋮\\ -m_{n1}& -m_{n2}& \cdots & {\mu }_{nm}+d_{nm}+\sum _{j\ne n}{m}_{jn}\end{array}\right],\end{array}$

and $V_{31}=\mathrm{diag}(-{\gamma }_{1f},-{\gamma }_{2f},\dots ,-{\gamma }_{nf}),V_{42}=\mathrm{diag}(-{\gamma }_{1m},-{\gamma }_{2m},\dots ,-{\gamma }_{nm}).$Since $V_{ii}$, for i = 1, 2, 3, 4, is a strictly diagonally dominant matrix, by the Gershgorin circle theorem, the real parts of its eigenvalues are positive and therefore $V_{ii}^{-1}$ exists. So the inverse of V exists and equals $V^{-1}=\left[\begin{array}{cccc}{V}_{11}^{-1}& 0& 0& 0\\ 0& V_{22}^{-1}& 0& 0\\ -V_{33}^{-1}{V}_{31}{V}_{11}^{-1}& 0& V_{33}^{-1}& 0\\ 0& -V_{44}^{-1}{V}_{42}{V}_{22}^{-1}& 0& V_{44}^{-1}\end{array}\right].$Thus the next generation matrix of system (1(1) $\begin{array}{rl}\frac{\mathrm{d}{S}_{if}}{\mathrm{d}t}& ={\epsilon }_i{\mathrm{\Lambda }}_i-{\beta }_i\frac{S_{if}{I}_{im}}{N_{im}}-{\mu }_{if}{S}_{if}+\sum _{j=1,j\ne i}^n(m_{ij}{S}_{jf}-m_{ji}{S}_{if}),\\ \frac{\mathrm{d}{I}_{if}}{\mathrm{d}t}& ={\beta }_i\frac{S_{if}{I}_{im}}{N_{im}}-({\mu }_{if}+{\gamma }_{if})I_{if}+\sum _{j=1,j\ne i}^n(m_{ij}{I}_{jf}-m_{ji}{I}_{if}),\\ \frac{\mathrm{d}{A}_{if}}{\mathrm{d}t}& ={\gamma }_{if}{I}_{if}-({\mu }_{if}+d_{if})A_{if}+\sum _{j=1,j\ne i}^n(m_{ij}{A}_{jf}-m_{ji}{A}_{if}),\\ \frac{\mathrm{d}{S}_{im}}{\mathrm{d}t}& =(1-{\epsilon }_i){\mathrm{\Lambda }}_i-{\beta }_i\frac{S_{im}{I}_{if}}{N_{if}}-{\mu }_{im}{S}_{im}+\sum _{j=1,j\ne i}^n(m_{ij}{S}_{jm}-m_{ji}{S}_{im}),\\ \frac{\mathrm{d}{I}_{im}}{\mathrm{d}t}& ={\beta }_i\frac{S_{im}{I}_{if}}{N_{if}}-({\mu }_{im}+{\gamma }_{im})I_{im}+\sum _{j=1,j\ne i}^n(m_{ij}{I}_{jm}-m_{ji}{I}_{im}),\\ \frac{\mathrm{d}{A}_{im}}{\mathrm{d}t}& ={\gamma }_{im}{I}_{im}-({\mu }_{im}+d_{im})A_{im}+\sum _{j=1,j\ne i}^n(m_{ij}{A}_{jm}-m_{ji}{A}_{im}),\end{array}$(1) ) is $F{V}^{-1}=\left[\begin{array}{cccc}0& F_{12}{V}_{22}^{-1}& 0& 0\\ F_{21}{V}_{11}^{-1}& 0& 0& 0\\ 0& 0& 0& 0\\ 0& 0& 0& 0\end{array}\right].$By calculating $(F{V}^{-1}{)}^2$, we find the basic reproduction number $R_0=\sqrt{\rho (F_{12}{V}_{22}^{-1}{F}_{21}{V}_{11}^{-1})},$where ρ represents the spectral radius of the matrix.

The basic reproduction number for the ith patch in isolation (i.e. there is no migration between patch i and other patches) is given by $R_0^{(i)}=\frac{{\beta }_i}{\sqrt{({\mu }_{if}+{\gamma }_{if})({\mu }_{im}+{\gamma }_{im})}}.$

3.2. Global stability of the disease-free equilibrium

Following Theorem 2 of [31], we have the following result on the local stability of $E_0$. $F-V=\left[\begin{array}{cccc}-V_{11}& F_{12}& 0& 0\\ F_{21}& -V_{22}& 0& 0\\ -V_{31}& 0& -V_{33}& 0\\ 0& -V_{42}& 0& -V_{44}\end{array}\right].$Define $s(F-V)$=max{$Re\lambda :\lambda$ is an eigenvalue of F−V}, so $s(F-V)$ is a simple eigenvalue of F−V with a positive eigenvector [28]. By Theorem 2 of van den Driessche and Watmough [31], there hold two equivalences: $s(F-V)>0⟺R_0>1,s(F-V)<0⟺R_0<1.$

Theorem 3.1

The disease-free equilibrium $E_0$ of system (1(1) $\begin{array}{rl}\frac{\mathrm{d}{S}_{if}}{\mathrm{d}t}& ={\epsilon }_i{\mathrm{\Lambda }}_i-{\beta }_i\frac{S_{if}{I}_{im}}{N_{im}}-{\mu }_{if}{S}_{if}+\sum _{j=1,j\ne i}^n(m_{ij}{S}_{jf}-m_{ji}{S}_{if}),\\ \frac{\mathrm{d}{I}_{if}}{\mathrm{d}t}& ={\beta }_i\frac{S_{if}{I}_{im}}{N_{im}}-({\mu }_{if}+{\gamma }_{if})I_{if}+\sum _{j=1,j\ne i}^n(m_{ij}{I}_{jf}-m_{ji}{I}_{if}),\\ \frac{\mathrm{d}{A}_{if}}{\mathrm{d}t}& ={\gamma }_{if}{I}_{if}-({\mu }_{if}+d_{if})A_{if}+\sum _{j=1,j\ne i}^n(m_{ij}{A}_{jf}-m_{ji}{A}_{if}),\\ \frac{\mathrm{d}{S}_{im}}{\mathrm{d}t}& =(1-{\epsilon }_i){\mathrm{\Lambda }}_i-{\beta }_i\frac{S_{im}{I}_{if}}{N_{if}}-{\mu }_{im}{S}_{im}+\sum _{j=1,j\ne i}^n(m_{ij}{S}_{jm}-m_{ji}{S}_{im}),\\ \frac{\mathrm{d}{I}_{im}}{\mathrm{d}t}& ={\beta }_i\frac{S_{im}{I}_{if}}{N_{if}}-({\mu }_{im}+{\gamma }_{im})I_{im}+\sum _{j=1,j\ne i}^n(m_{ij}{I}_{jm}-m_{ji}{I}_{im}),\\ \frac{\mathrm{d}{A}_{im}}{\mathrm{d}t}& ={\gamma }_{im}{I}_{im}-({\mu }_{im}+d_{im})A_{im}+\sum _{j=1,j\ne i}^n(m_{ij}{A}_{jm}-m_{ji}{A}_{im}),\end{array}$(1) ) is locally asymptotically stable if $R_0<1$ and unstable if $R_0>1$.

Proof.

To prove the local stability of disease-free equilibrium $E_0$ of system (1(1) $\begin{array}{rl}\frac{\mathrm{d}{S}_{if}}{\mathrm{d}t}& ={\epsilon }_i{\mathrm{\Lambda }}_i-{\beta }_i\frac{S_{if}{I}_{im}}{N_{im}}-{\mu }_{if}{S}_{if}+\sum _{j=1,j\ne i}^n(m_{ij}{S}_{jf}-m_{ji}{S}_{if}),\\ \frac{\mathrm{d}{I}_{if}}{\mathrm{d}t}& ={\beta }_i\frac{S_{if}{I}_{im}}{N_{im}}-({\mu }_{if}+{\gamma }_{if})I_{if}+\sum _{j=1,j\ne i}^n(m_{ij}{I}_{jf}-m_{ji}{I}_{if}),\\ \frac{\mathrm{d}{A}_{if}}{\mathrm{d}t}& ={\gamma }_{if}{I}_{if}-({\mu }_{if}+d_{if})A_{if}+\sum _{j=1,j\ne i}^n(m_{ij}{A}_{jf}-m_{ji}{A}_{if}),\\ \frac{\mathrm{d}{S}_{im}}{\mathrm{d}t}& =(1-{\epsilon }_i){\mathrm{\Lambda }}_i-{\beta }_i\frac{S_{im}{I}_{if}}{N_{if}}-{\mu }_{im}{S}_{im}+\sum _{j=1,j\ne i}^n(m_{ij}{S}_{jm}-m_{ji}{S}_{im}),\\ \frac{\mathrm{d}{I}_{im}}{\mathrm{d}t}& ={\beta }_i\frac{S_{im}{I}_{if}}{N_{if}}-({\mu }_{im}+{\gamma }_{im})I_{im}+\sum _{j=1,j\ne i}^n(m_{ij}{I}_{jm}-m_{ji}{I}_{im}),\\ \frac{\mathrm{d}{A}_{im}}{\mathrm{d}t}& ={\gamma }_{im}{I}_{im}-({\mu }_{im}+d_{im})A_{im}+\sum _{j=1,j\ne i}^n(m_{ij}{A}_{jm}-m_{ji}{A}_{im}),\end{array}$(1) ), we check the hypotheses (A1)–(A5) in [31]. Hypotheses (A1)–(A4) are easily verified, while (A5) is satisfied if all eigenvalues of the $6n\times 6n$ matrix $J{\mid }_{E_0}=\left[\begin{array}{cc}F-V& 0\\ J_3& J_4\end{array}\right]$ have negative real parts, where $J_3=\left[\begin{array}{cc}0& -F_{12}\\ -F_{21}& 0\end{array}\right],J_4=\left[\begin{array}{cc}-A& 0\\ 0& -B\end{array}\right].$The eigenvalues of $J{\mid }_{E_0}$ are the eigenvalues of F−V and $J_4$. The matrices A and B are two nonsingular M-matrices. Thus $-A$ and $-B$ have all eigenvalues with negative real parts, so $J_4$ has all eigenvalues with negative real parts. Consequently, the local stability of the disease-free equilibrium depends only on eigenvalues of F−V. Thus if $R_0<1$, then $s(F-V)<0$ and $s(J{\mid }_{E_0})<0$, the disease-free equilibrium $E_0$ of system (1(1) $\begin{array}{rl}\frac{\mathrm{d}{S}_{if}}{\mathrm{d}t}& ={\epsilon }_i{\mathrm{\Lambda }}_i-{\beta }_i\frac{S_{if}{I}_{im}}{N_{im}}-{\mu }_{if}{S}_{if}+\sum _{j=1,j\ne i}^n(m_{ij}{S}_{jf}-m_{ji}{S}_{if}),\\ \frac{\mathrm{d}{I}_{if}}{\mathrm{d}t}& ={\beta }_i\frac{S_{if}{I}_{im}}{N_{im}}-({\mu }_{if}+{\gamma }_{if})I_{if}+\sum _{j=1,j\ne i}^n(m_{ij}{I}_{jf}-m_{ji}{I}_{if}),\\ \frac{\mathrm{d}{A}_{if}}{\mathrm{d}t}& ={\gamma }_{if}{I}_{if}-({\mu }_{if}+d_{if})A_{if}+\sum _{j=1,j\ne i}^n(m_{ij}{A}_{jf}-m_{ji}{A}_{if}),\\ \frac{\mathrm{d}{S}_{im}}{\mathrm{d}t}& =(1-{\epsilon }_i){\mathrm{\Lambda }}_i-{\beta }_i\frac{S_{im}{I}_{if}}{N_{if}}-{\mu }_{im}{S}_{im}+\sum _{j=1,j\ne i}^n(m_{ij}{S}_{jm}-m_{ji}{S}_{im}),\\ \frac{\mathrm{d}{I}_{im}}{\mathrm{d}t}& ={\beta }_i\frac{S_{im}{I}_{if}}{N_{if}}-({\mu }_{im}+{\gamma }_{im})I_{im}+\sum _{j=1,j\ne i}^n(m_{ij}{I}_{jm}-m_{ji}{I}_{im}),\\ \frac{\mathrm{d}{A}_{im}}{\mathrm{d}t}& ={\gamma }_{im}{I}_{im}-({\mu }_{im}+d_{im})A_{im}+\sum _{j=1,j\ne i}^n(m_{ij}{A}_{jm}-m_{ji}{A}_{im}),\end{array}$(1) ) is locally asymptotically stable.

Theorem 3.2

The disease-free equilibrium $E_0$ of system (1(1) $\begin{array}{rl}\frac{\mathrm{d}{S}_{if}}{\mathrm{d}t}& ={\epsilon }_i{\mathrm{\Lambda }}_i-{\beta }_i\frac{S_{if}{I}_{im}}{N_{im}}-{\mu }_{if}{S}_{if}+\sum _{j=1,j\ne i}^n(m_{ij}{S}_{jf}-m_{ji}{S}_{if}),\\ \frac{\mathrm{d}{I}_{if}}{\mathrm{d}t}& ={\beta }_i\frac{S_{if}{I}_{im}}{N_{im}}-({\mu }_{if}+{\gamma }_{if})I_{if}+\sum _{j=1,j\ne i}^n(m_{ij}{I}_{jf}-m_{ji}{I}_{if}),\\ \frac{\mathrm{d}{A}_{if}}{\mathrm{d}t}& ={\gamma }_{if}{I}_{if}-({\mu }_{if}+d_{if})A_{if}+\sum _{j=1,j\ne i}^n(m_{ij}{A}_{jf}-m_{ji}{A}_{if}),\\ \frac{\mathrm{d}{S}_{im}}{\mathrm{d}t}& =(1-{\epsilon }_i){\mathrm{\Lambda }}_i-{\beta }_i\frac{S_{im}{I}_{if}}{N_{if}}-{\mu }_{im}{S}_{im}+\sum _{j=1,j\ne i}^n(m_{ij}{S}_{jm}-m_{ji}{S}_{im}),\\ \frac{\mathrm{d}{I}_{im}}{\mathrm{d}t}& ={\beta }_i\frac{S_{im}{I}_{if}}{N_{if}}-({\mu }_{im}+{\gamma }_{im})I_{im}+\sum _{j=1,j\ne i}^n(m_{ij}{I}_{jm}-m_{ji}{I}_{im}),\\ \frac{\mathrm{d}{A}_{im}}{\mathrm{d}t}& ={\gamma }_{im}{I}_{im}-({\mu }_{im}+d_{im})A_{im}+\sum _{j=1,j\ne i}^n(m_{ij}{A}_{jm}-m_{ji}{A}_{im}),\end{array}$(1) ) is globally asymptotically stable if $R_0<1$.

Proof.

Since $S_{if}\le N_{if}$, $S_{im}\le N_{im}$, thus (2) $\begin{array}{rl}\frac{\mathrm{d}{I}_{if}}{\mathrm{d}t}& \le {\beta }_i\frac{S_{if}^0}{S_{im}^0}{I}_{im}-({\mu }_{if}+{\gamma }_{if})I_{if}+\sum _{j=1}^n(m_{ij}{I}_{jf}-m_{ji}{I}_{if}),\\ \frac{\mathrm{d}{I}_{im}}{\mathrm{d}t}& \le {\beta }_i\frac{S_{im}^0}{S_{if}^0}{I}_{if}-({\mu }_{im}+{\gamma }_{im})I_{im}+\sum _{j=1}^n(m_{ij}{I}_{jm}-m_{ji}{I}_{im}),\\ \frac{\mathrm{d}{A}_{if}}{\mathrm{d}t}& ={\gamma }_{if}{I}_{if}-({\mu }_{if}+d_{if})A_{if}+\sum _{j=1}^n(m_{ij}{A}_{jf}-m_{ji}{A}_{if}),\\ \frac{\mathrm{d}{A}_{im}}{\mathrm{d}t}& ={\gamma }_{im}{I}_{im}-({\mu }_{im}+d_{im})A_{im}+\sum _{j=1}^n(m_{ij}{A}_{jm}-m_{ji}{A}_{im}).\end{array}$(2) Consider the following auxiliary system: (3) $\begin{array}{rl}\frac{\mathrm{d}{\overline{I}}_{if}}{\mathrm{d}t}& ={\beta }_i\frac{S_{if}^0}{S_{im}^0}{\overline{I}}_{im}-({\mu }_{if}+{\gamma }_{if}){\overline{I}}_{if}+\sum _{j=1}^n(m_{ij}{\overline{I}}_{jf}-m_{ji}{\overline{I}}_{if}),\\ \frac{\mathrm{d}{\overline{I}}_{im}}{\mathrm{d}t}& ={\beta }_i\frac{S_{im}^0}{S_{if}^0}{\overline{I}}_{if}-({\mu }_{im}+{\gamma }_{im}){\overline{I}}_{im}+\sum _{j=1}^n(m_{ij}{\overline{I}}_{jm}-m_{ji}{\overline{I}}_{im}),\\ \frac{\mathrm{d}{\overline{A}}_{if}}{\mathrm{d}t}& ={\gamma }_{if}{\overline{I}}_{if}-({\mu }_{if}+d_{if}){\overline{A}}_{if}+\sum _{j=1}^n(m_{ij}{\overline{A}}_{jf}-m_{ji}{\overline{A}}_{if}),\\ \frac{\mathrm{d}{\overline{A}}_{im}}{\mathrm{d}t}& ={\gamma }_{im}{\overline{I}}_{im}-({\mu }_{im}+d_{im}){\overline{A}}_{im}+\sum _{j=1}^n(m_{ij}{\overline{A}}_{jm}-m_{ji}{\overline{A}}_{im}).\end{array}$(3) The right side has the coefficient matrix F−V. Since $R_0<1$ if and only if $s(F-V)<0$. By the Perron Frobenius theorem, all eigenvalues of the matrix F−V have negative real parts when $s(F-V)<0$. Therefore it has $\begin{array}{rl}& ({\overline{I}}_{1f},\dots ,{\overline{I}}_{nf},{\overline{I}}_{1m},\dots ,{\overline{I}}_{nm},{\overline{A}}_{1f},\dots ,{\overline{A}}_{nf},{\overline{A}}_{1m},\dots ,{\overline{A}}_{nm})\\ & \to (0,\dots ,0,0,\dots ,0,0,\dots ,0,0,\dots ,0)\end{array}$as $t\to \mathrm{\infty }.$ Using the comparison principle of Smith and Waltman, namely Theorem B.1[28], we know that $\begin{array}{rl}& (I_{1f},\dots ,I_{nf},I_{1m},\dots ,I_{nm},A_{1f},\dots ,A_{nf},A_{1m},\dots ,A_{nm})\\ & \to (0,\dots ,0,0,\dots ,0,0,\dots ,0,0,\dots ,0)\end{array}$as $t\to \mathrm{\infty }$. So, (1(1) $\begin{array}{rl}\frac{\mathrm{d}{S}_{if}}{\mathrm{d}t}& ={\epsilon }_i{\mathrm{\Lambda }}_i-{\beta }_i\frac{S_{if}{I}_{im}}{N_{im}}-{\mu }_{if}{S}_{if}+\sum _{j=1,j\ne i}^n(m_{ij}{S}_{jf}-m_{ji}{S}_{if}),\\ \frac{\mathrm{d}{I}_{if}}{\mathrm{d}t}& ={\beta }_i\frac{S_{if}{I}_{im}}{N_{im}}-({\mu }_{if}+{\gamma }_{if})I_{if}+\sum _{j=1,j\ne i}^n(m_{ij}{I}_{jf}-m_{ji}{I}_{if}),\\ \frac{\mathrm{d}{A}_{if}}{\mathrm{d}t}& ={\gamma }_{if}{I}_{if}-({\mu }_{if}+d_{if})A_{if}+\sum _{j=1,j\ne i}^n(m_{ij}{A}_{jf}-m_{ji}{A}_{if}),\\ \frac{\mathrm{d}{S}_{im}}{\mathrm{d}t}& =(1-{\epsilon }_i){\mathrm{\Lambda }}_i-{\beta }_i\frac{S_{im}{I}_{if}}{N_{if}}-{\mu }_{im}{S}_{im}+\sum _{j=1,j\ne i}^n(m_{ij}{S}_{jm}-m_{ji}{S}_{im}),\\ \frac{\mathrm{d}{I}_{im}}{\mathrm{d}t}& ={\beta }_i\frac{S_{im}{I}_{if}}{N_{if}}-({\mu }_{im}+{\gamma }_{im})I_{im}+\sum _{j=1,j\ne i}^n(m_{ij}{I}_{jm}-m_{ji}{I}_{im}),\\ \frac{\mathrm{d}{A}_{im}}{\mathrm{d}t}& ={\gamma }_{im}{I}_{im}-({\mu }_{im}+d_{im})A_{im}+\sum _{j=1,j\ne i}^n(m_{ij}{A}_{jm}-m_{ji}{A}_{im}),\end{array}$(1) ) is an asymptotically autonomous system with limit affine system $\begin{array}{rl}\frac{\mathrm{d}{S}_{if}}{\mathrm{d}t}& ={\epsilon }_i{\mathrm{\Lambda }}_i-({\mu }_{if}+\sum _{j=1,j\ne i}^n{m}_{ji})S_{if}+\sum _{j=1,j\ne i}^n{m}_{ij}{S}_{jf},\\ \frac{\mathrm{d}{S}_{im}}{\mathrm{d}t}& =(1-{\epsilon }_i){\mathrm{\Lambda }}_i-({\mu }_{im}+\sum _{j=1,j\ne i}^n{m}_{ji})S_{im}+\sum _{j=1,j\ne i}^n{m}_{ij}{S}_{jm}.\end{array}$By the theory of asymptotic autonomous system of Thieme [29], it is also known that $(S_{1f},\dots ,S_{nf},S_{1m},\cdots ,S_{nm})\to (S_{1f}^0,\dots ,S_{nf}^0,S_{1m}^0,\dots ,S_{nm}^0)$, $t\to \mathrm{\infty }$. So $E_0$ is globally attractive when $R_0<1$. It follows that the disease-free equilibrium $E_0$ of system (1(1) $\begin{array}{rl}\frac{\mathrm{d}{S}_{if}}{\mathrm{d}t}& ={\epsilon }_i{\mathrm{\Lambda }}_i-{\beta }_i\frac{S_{if}{I}_{im}}{N_{im}}-{\mu }_{if}{S}_{if}+\sum _{j=1,j\ne i}^n(m_{ij}{S}_{jf}-m_{ji}{S}_{if}),\\ \frac{\mathrm{d}{I}_{if}}{\mathrm{d}t}& ={\beta }_i\frac{S_{if}{I}_{im}}{N_{im}}-({\mu }_{if}+{\gamma }_{if})I_{if}+\sum _{j=1,j\ne i}^n(m_{ij}{I}_{jf}-m_{ji}{I}_{if}),\\ \frac{\mathrm{d}{A}_{if}}{\mathrm{d}t}& ={\gamma }_{if}{I}_{if}-({\mu }_{if}+d_{if})A_{if}+\sum _{j=1,j\ne i}^n(m_{ij}{A}_{jf}-m_{ji}{A}_{if}),\\ \frac{\mathrm{d}{S}_{im}}{\mathrm{d}t}& =(1-{\epsilon }_i){\mathrm{\Lambda }}_i-{\beta }_i\frac{S_{im}{I}_{if}}{N_{if}}-{\mu }_{im}{S}_{im}+\sum _{j=1,j\ne i}^n(m_{ij}{S}_{jm}-m_{ji}{S}_{im}),\\ \frac{\mathrm{d}{I}_{im}}{\mathrm{d}t}& ={\beta }_i\frac{S_{im}{I}_{if}}{N_{if}}-({\mu }_{im}+{\gamma }_{im})I_{im}+\sum _{j=1,j\ne i}^n(m_{ij}{I}_{jm}-m_{ji}{I}_{im}),\\ \frac{\mathrm{d}{A}_{im}}{\mathrm{d}t}& ={\gamma }_{im}{I}_{im}-({\mu }_{im}+d_{im})A_{im}+\sum _{j=1,j\ne i}^n(m_{ij}{A}_{jm}-m_{ji}{A}_{im}),\end{array}$(1) ) is globally asymptotically stable when $R_0<1$.

3.3. The uniform persistence of system (1)

Using the techniques of persistence theory [36], we can show the uniform persistence of the disease and the existence of at least one endemic equilibrium when $R_0>1$. Thus the basic reproduction number $R_0$ is a threshold parameter of the disease dynamics. The proof below is analogous to those of Theorem 2.3 in Wang and Zhao [33] and Theorem 3.7 in Gao and Ruan [15].

Theorem 3.3

For system (1(1) $\begin{array}{rl}\frac{\mathrm{d}{S}_{if}}{\mathrm{d}t}& ={\epsilon }_i{\mathrm{\Lambda }}_i-{\beta }_i\frac{S_{if}{I}_{im}}{N_{im}}-{\mu }_{if}{S}_{if}+\sum _{j=1,j\ne i}^n(m_{ij}{S}_{jf}-m_{ji}{S}_{if}),\\ \frac{\mathrm{d}{I}_{if}}{\mathrm{d}t}& ={\beta }_i\frac{S_{if}{I}_{im}}{N_{im}}-({\mu }_{if}+{\gamma }_{if})I_{if}+\sum _{j=1,j\ne i}^n(m_{ij}{I}_{jf}-m_{ji}{I}_{if}),\\ \frac{\mathrm{d}{A}_{if}}{\mathrm{d}t}& ={\gamma }_{if}{I}_{if}-({\mu }_{if}+d_{if})A_{if}+\sum _{j=1,j\ne i}^n(m_{ij}{A}_{jf}-m_{ji}{A}_{if}),\\ \frac{\mathrm{d}{S}_{im}}{\mathrm{d}t}& =(1-{\epsilon }_i){\mathrm{\Lambda }}_i-{\beta }_i\frac{S_{im}{I}_{if}}{N_{if}}-{\mu }_{im}{S}_{im}+\sum _{j=1,j\ne i}^n(m_{ij}{S}_{jm}-m_{ji}{S}_{im}),\\ \frac{\mathrm{d}{I}_{im}}{\mathrm{d}t}& ={\beta }_i\frac{S_{im}{I}_{if}}{N_{if}}-({\mu }_{im}+{\gamma }_{im})I_{im}+\sum _{j=1,j\ne i}^n(m_{ij}{I}_{jm}-m_{ji}{I}_{im}),\\ \frac{\mathrm{d}{A}_{im}}{\mathrm{d}t}& ={\gamma }_{im}{I}_{im}-({\mu }_{im}+d_{im})A_{im}+\sum _{j=1,j\ne i}^n(m_{ij}{A}_{jm}-m_{ji}{A}_{im}),\end{array}$(1) ), assume that $M=(m_{ij})$ is irreducible, then if $R_0>1$, the disease is uniformly persistent, i.e. there exists a constant k>0, such that every solution ${\mathrm{\Phi }}_t(x_0)=(S_{1f}(t),\dots ,S_{nf}(t),S_{1m}(t),\dots ,S_{nm}(t)$, $I_{1f}(t),\cdots ,I_{nf}(t),I_{1m}(t),\cdots$, $I_{nm}(t),A_{1f}(t),\cdots ,A_{nf}(t),A_{1m}(t),\cdots ,A_{nm}(t))$ of system (1(1) $\begin{array}{rl}\frac{\mathrm{d}{S}_{if}}{\mathrm{d}t}& ={\epsilon }_i{\mathrm{\Lambda }}_i-{\beta }_i\frac{S_{if}{I}_{im}}{N_{im}}-{\mu }_{if}{S}_{if}+\sum _{j=1,j\ne i}^n(m_{ij}{S}_{jf}-m_{ji}{S}_{if}),\\ \frac{\mathrm{d}{I}_{if}}{\mathrm{d}t}& ={\beta }_i\frac{S_{if}{I}_{im}}{N_{im}}-({\mu }_{if}+{\gamma }_{if})I_{if}+\sum _{j=1,j\ne i}^n(m_{ij}{I}_{jf}-m_{ji}{I}_{if}),\\ \frac{\mathrm{d}{A}_{if}}{\mathrm{d}t}& ={\gamma }_{if}{I}_{if}-({\mu }_{if}+d_{if})A_{if}+\sum _{j=1,j\ne i}^n(m_{ij}{A}_{jf}-m_{ji}{A}_{if}),\\ \frac{\mathrm{d}{S}_{im}}{\mathrm{d}t}& =(1-{\epsilon }_i){\mathrm{\Lambda }}_i-{\beta }_i\frac{S_{im}{I}_{if}}{N_{if}}-{\mu }_{im}{S}_{im}+\sum _{j=1,j\ne i}^n(m_{ij}{S}_{jm}-m_{ji}{S}_{im}),\\ \frac{\mathrm{d}{I}_{im}}{\mathrm{d}t}& ={\beta }_i\frac{S_{im}{I}_{if}}{N_{if}}-({\mu }_{im}+{\gamma }_{im})I_{im}+\sum _{j=1,j\ne i}^n(m_{ij}{I}_{jm}-m_{ji}{I}_{im}),\\ \frac{\mathrm{d}{A}_{im}}{\mathrm{d}t}& ={\gamma }_{im}{I}_{im}-({\mu }_{im}+d_{im})A_{im}+\sum _{j=1,j\ne i}^n(m_{ij}{A}_{jm}-m_{ji}{A}_{im}),\end{array}$(1) ) with $x_0=(S_{1f}(0),\cdots$, $S_{nf}(0)S_{1m}(0),\cdots ,S_{nm}(0)$, $I_{1f}(0),\cdots ,I_{nf}(0),I_{1m}(0)$, $\cdots ,I_{nm}(0),A_{1f}(0),\dots ,A_{nf}(0),A_{1m}(0),\dots ,A_{nm}(0))\in X_0$ satisfies $\begin{array}{rl}& \underset{t\to \mathrm{\infty }}{lim}\inf (I_{1f}(t),\dots ,I_{nf}(t),I_{1m}(t),\dots ,I_{nm}(t),A_{1f}(t),\dots ,A_{nf}(t),A_{1m}(t),\dots ,A_{nm}(t))\\ & >(k,\dots ,k{)}_{1\times 4n},\end{array}$and (1(1) $\begin{array}{rl}\frac{\mathrm{d}{S}_{if}}{\mathrm{d}t}& ={\epsilon }_i{\mathrm{\Lambda }}_i-{\beta }_i\frac{S_{if}{I}_{im}}{N_{im}}-{\mu }_{if}{S}_{if}+\sum _{j=1,j\ne i}^n(m_{ij}{S}_{jf}-m_{ji}{S}_{if}),\\ \frac{\mathrm{d}{I}_{if}}{\mathrm{d}t}& ={\beta }_i\frac{S_{if}{I}_{im}}{N_{im}}-({\mu }_{if}+{\gamma }_{if})I_{if}+\sum _{j=1,j\ne i}^n(m_{ij}{I}_{jf}-m_{ji}{I}_{if}),\\ \frac{\mathrm{d}{A}_{if}}{\mathrm{d}t}& ={\gamma }_{if}{I}_{if}-({\mu }_{if}+d_{if})A_{if}+\sum _{j=1,j\ne i}^n(m_{ij}{A}_{jf}-m_{ji}{A}_{if}),\\ \frac{\mathrm{d}{S}_{im}}{\mathrm{d}t}& =(1-{\epsilon }_i){\mathrm{\Lambda }}_i-{\beta }_i\frac{S_{im}{I}_{if}}{N_{if}}-{\mu }_{im}{S}_{im}+\sum _{j=1,j\ne i}^n(m_{ij}{S}_{jm}-m_{ji}{S}_{im}),\\ \frac{\mathrm{d}{I}_{im}}{\mathrm{d}t}& ={\beta }_i\frac{S_{im}{I}_{if}}{N_{if}}-({\mu }_{im}+{\gamma }_{im})I_{im}+\sum _{j=1,j\ne i}^n(m_{ij}{I}_{jm}-m_{ji}{I}_{im}),\\ \frac{\mathrm{d}{A}_{im}}{\mathrm{d}t}& ={\gamma }_{im}{I}_{im}-({\mu }_{im}+d_{im})A_{im}+\sum _{j=1,j\ne i}^n(m_{ij}{A}_{jm}-m_{ji}{A}_{im}),\end{array}$(1) ) admits at least one endemic equilibrium.

Proof.

Let $\begin{array}{rl}X& =\{(S_{1f},\dots ,S_{nf},S_{1m},\dots ,S_{nm},I_{1f},\dots ,I_{nf},I_{1m},\dots ,I_{nm},A_{1f},\dots ,A_{nf},A_{1m},\dots ,A_{nm}):\\ & S_{if}\ge 0,S_{im}\ge 0,I_{if}\ge 0,I_{im}\ge 0,A_{if}\ge 0,A_{im}\ge 0,i=1,2,\dots ,n\},\\ X_0& =\{(S_{1f},\cdots ,S_{nf},S_{1m},\cdots ,S_{nm},I_{1f},\cdots ,I_{nf},I_{1m},\cdots ,I_{nm},A_{1f},\cdots ,A_{nf},A_{1m},\cdots ,A_{nm})\\ & \in X:I_{if}>0,I_{im}>0,A_{if}>0,A_{im}>0,i=1,2,\dots ,n\},\\ \mathrm{\partial }{X}_0& =X\mathrm{\setminus }{X}_0=\{(S_{1f},\cdots ,S_{nf},S_{1m},\cdots ,S_{nm},I_{1f},\cdots ,I_{nf},I_{1m},\cdots ,I_{nm},A_{1f},\cdots ,A_{nf},A_{1m},\cdots ,A_{nm})\\ & \in X:I_{if}=0,I_{im}=0,A_{if}=0,A_{im}=0,\mathrm{forsome}i\in \{1,2,\dots ,n\}\}.\end{array}$

It suffices to prove that $\mathrm{\partial }{X}_0$ repels uniformly the solutions of system (1(1) $\begin{array}{rl}\frac{\mathrm{d}{S}_{if}}{\mathrm{d}t}& ={\epsilon }_i{\mathrm{\Lambda }}_i-{\beta }_i\frac{S_{if}{I}_{im}}{N_{im}}-{\mu }_{if}{S}_{if}+\sum _{j=1,j\ne i}^n(m_{ij}{S}_{jf}-m_{ji}{S}_{if}),\\ \frac{\mathrm{d}{I}_{if}}{\mathrm{d}t}& ={\beta }_i\frac{S_{if}{I}_{im}}{N_{im}}-({\mu }_{if}+{\gamma }_{if})I_{if}+\sum _{j=1,j\ne i}^n(m_{ij}{I}_{jf}-m_{ji}{I}_{if}),\\ \frac{\mathrm{d}{A}_{if}}{\mathrm{d}t}& ={\gamma }_{if}{I}_{if}-({\mu }_{if}+d_{if})A_{if}+\sum _{j=1,j\ne i}^n(m_{ij}{A}_{jf}-m_{ji}{A}_{if}),\\ \frac{\mathrm{d}{S}_{im}}{\mathrm{d}t}& =(1-{\epsilon }_i){\mathrm{\Lambda }}_i-{\beta }_i\frac{S_{im}{I}_{if}}{N_{if}}-{\mu }_{im}{S}_{im}+\sum _{j=1,j\ne i}^n(m_{ij}{S}_{jm}-m_{ji}{S}_{im}),\\ \frac{\mathrm{d}{I}_{im}}{\mathrm{d}t}& ={\beta }_i\frac{S_{im}{I}_{if}}{N_{if}}-({\mu }_{im}+{\gamma }_{im})I_{im}+\sum _{j=1,j\ne i}^n(m_{ij}{I}_{jm}-m_{ji}{I}_{im}),\\ \frac{\mathrm{d}{A}_{im}}{\mathrm{d}t}& ={\gamma }_{im}{I}_{im}-({\mu }_{im}+d_{im})A_{im}+\sum _{j=1,j\ne i}^n(m_{ij}{A}_{jm}-m_{ji}{A}_{im}),\end{array}$(1) ) in $X_0$. Clearly, $\mathrm{\partial }{X}_0$ is relatively closed in X. It is immediate that X and $X_0$ are positively invariant. Theorem 2.1 implies that system (1(1) $\begin{array}{rl}\frac{\mathrm{d}{S}_{if}}{\mathrm{d}t}& ={\epsilon }_i{\mathrm{\Lambda }}_i-{\beta }_i\frac{S_{if}{I}_{im}}{N_{im}}-{\mu }_{if}{S}_{if}+\sum _{j=1,j\ne i}^n(m_{ij}{S}_{jf}-m_{ji}{S}_{if}),\\ \frac{\mathrm{d}{I}_{if}}{\mathrm{d}t}& ={\beta }_i\frac{S_{if}{I}_{im}}{N_{im}}-({\mu }_{if}+{\gamma }_{if})I_{if}+\sum _{j=1,j\ne i}^n(m_{ij}{I}_{jf}-m_{ji}{I}_{if}),\\ \frac{\mathrm{d}{A}_{if}}{\mathrm{d}t}& ={\gamma }_{if}{I}_{if}-({\mu }_{if}+d_{if})A_{if}+\sum _{j=1,j\ne i}^n(m_{ij}{A}_{jf}-m_{ji}{A}_{if}),\\ \frac{\mathrm{d}{S}_{im}}{\mathrm{d}t}& =(1-{\epsilon }_i){\mathrm{\Lambda }}_i-{\beta }_i\frac{S_{im}{I}_{if}}{N_{if}}-{\mu }_{im}{S}_{im}+\sum _{j=1,j\ne i}^n(m_{ij}{S}_{jm}-m_{ji}{S}_{im}),\\ \frac{\mathrm{d}{I}_{im}}{\mathrm{d}t}& ={\beta }_i\frac{S_{im}{I}_{if}}{N_{if}}-({\mu }_{im}+{\gamma }_{im})I_{im}+\sum _{j=1,j\ne i}^n(m_{ij}{I}_{jm}-m_{ji}{I}_{im}),\\ \frac{\mathrm{d}{A}_{im}}{\mathrm{d}t}& ={\gamma }_{im}{I}_{im}-({\mu }_{im}+d_{im})A_{im}+\sum _{j=1,j\ne i}^n(m_{ij}{A}_{jm}-m_{ji}{A}_{im}),\end{array}$(1) ) is point dissipative.

Denote $M_{\alpha }=\{x_0\in \mathrm{\partial }{X}_0:{\mathrm{\Phi }}_t(x_0)\in \mathrm{\partial }{X}_{0}fort\ge 0,i=1,2,\dots ,n\},D=\{x_0\in X:I_{if}=0,I_{im}=0,A_{if}=0,A_{im}=0,i=1,2,\dots ,n\}$. Obviously, $D\subset M_{\alpha }$. On the other hand, we have $I_{1f}(0)+\cdots +I_{nf}(0)+I_{1m}(0)+\cdots +I_{nm}(0)+A_{1f}(0)+\cdots +A_{nf}(0)+A_{1m}(0)+\cdots +A_{nm}(0)>0$ for any $x_0\in \mathrm{\partial }{X}_0\setminus D$. By the irreducibility of the migration matrix $M=(m_{ij}{)}_{n\times n}$, we know that ${\mathrm{\Phi }}_t(x_0)\in X_0$ for all t>0. Therefore, $x_0\notin M_{\alpha }$ and $M_{\alpha }\subset D$, which implies that $M_{\alpha }=D$.

The disease-free equilibrium $E_0$ is the unique equilibrium in $M_{\alpha }$. Let $W^s(E_0)$ be the stable manifold of $E_0$. We now show that $W^s(E_0)\bigcap X_0=\varnothing$ when $R_0>1$.

Let $M_{\epsilon }=F-V-ϵF$. Since $s(F-V)>0$ if and only if $R_0>1$, there is an ${ϵ}_1>0$ such that $s(M_{ϵ})>0$ for $ϵ\in [0,{ϵ}_1]$. Choose η small enough such that $\frac{S_{if}(0)}{N_{im}(0)}\ge \frac{S_{if}^0}{S_{im}^0}-{ϵ}_1$ and $\frac{S_{im}(0)}{N_{if}(0)}\ge \frac{S_{im}^0}{S_{if}^0}-{ϵ}_1$ for $i=1,2,\dots ,n$, $\parallel x_0-E_0\parallel \le \eta$.

We now show that $\underset{t\to \mathrm{\infty }}{lim}\sup \parallel x_0-E_0\parallel >\eta$ for $x_0\in X_0$, where $\parallel \cdot \parallel$ is the usual Euclidean norm. Suppose not, then $\parallel x_0-E_0\parallel \le \eta$ for $t\ge 0$ and $\begin{array}{rl}\frac{\mathrm{d}{I}_{if}}{\mathrm{d}t}& \ge ({\beta }_i\frac{S_{if}^0}{S_{im}^0}-{ϵ}_1)I_{im}-({\mu }_{if}+{\gamma }_{if}+\sum _{j\ne i}{m}_{ji})I_{if}+\sum _{j\ne i}{m}_{ij}{I}_{jf},\\ \frac{\mathrm{d}{I}_{im}}{\mathrm{d}t}& \ge ({\beta }_i\frac{S_{im}^0}{S_{if}^0}-{ϵ}_1)I_{if}-({\mu }_{im}+{\gamma }_{im}+\sum _{j\ne i}{m}_{ji})I_{im}+\sum _{j\ne i}{m}_{ij}{I}_{jm},\\ \frac{\mathrm{d}{A}_{if}}{\mathrm{d}t}& ={\gamma }_{if}{I}_{if}-({\mu }_{if}+d_{if}+\sum _{j\ne i}{m}_{ji})A_{if}+\sum _{j\ne i}{m}_{ij}{A}_{jf},\\ \frac{\mathrm{d}{A}_{im}}{\mathrm{d}t}& ={\gamma }_{im}{I}_{im}-({\mu }_{im}+d_{im}+\sum _{j\ne i}{m}_{ji})A_{im}+\sum _{j\ne i}{m}_{ij}{A}_{jm}.\end{array}$Consider an auxiliary system $\frac{\mathrm{d}u}{\mathrm{d}t}=M_{{ϵ}_1}u.$Note that $M_{{ϵ}_1}$ has a positive eigenvalue $s(M_{{ϵ}_1})$ associated to a positive eigenvector. By the comparison theorem, we have $\underset{t\to \mathrm{\infty }}{lim}(I_{if}(t),I_{im}(t),A_{if}(t),A_{im}(t))=(\mathrm{\infty },\mathrm{\infty },\mathrm{\infty },\mathrm{\infty }),i=1,2,\dots ,n.$The above conclusions contradict the assumptions. Since $E_0$ is globally stable in $M_{\alpha }$, it follows that $E_0$ is an isolated invariant set and $W^s(E_0)\bigcap X_0=\varnothing$. By Theorem 4.6 in Thieme [30], system (1(1) $\begin{array}{rl}\frac{\mathrm{d}{S}_{if}}{\mathrm{d}t}& ={\epsilon }_i{\mathrm{\Lambda }}_i-{\beta }_i\frac{S_{if}{I}_{im}}{N_{im}}-{\mu }_{if}{S}_{if}+\sum _{j=1,j\ne i}^n(m_{ij}{S}_{jf}-m_{ji}{S}_{if}),\\ \frac{\mathrm{d}{I}_{if}}{\mathrm{d}t}& ={\beta }_i\frac{S_{if}{I}_{im}}{N_{im}}-({\mu }_{if}+{\gamma }_{if})I_{if}+\sum _{j=1,j\ne i}^n(m_{ij}{I}_{jf}-m_{ji}{I}_{if}),\\ \frac{\mathrm{d}{A}_{if}}{\mathrm{d}t}& ={\gamma }_{if}{I}_{if}-({\mu }_{if}+d_{if})A_{if}+\sum _{j=1,j\ne i}^n(m_{ij}{A}_{jf}-m_{ji}{A}_{if}),\\ \frac{\mathrm{d}{S}_{im}}{\mathrm{d}t}& =(1-{\epsilon }_i){\mathrm{\Lambda }}_i-{\beta }_i\frac{S_{im}{I}_{if}}{N_{if}}-{\mu }_{im}{S}_{im}+\sum _{j=1,j\ne i}^n(m_{ij}{S}_{jm}-m_{ji}{S}_{im}),\\ \frac{\mathrm{d}{I}_{im}}{\mathrm{d}t}& ={\beta }_i\frac{S_{im}{I}_{if}}{N_{if}}-({\mu }_{im}+{\gamma }_{im})I_{im}+\sum _{j=1,j\ne i}^n(m_{ij}{I}_{jm}-m_{ji}{I}_{im}),\\ \frac{\mathrm{d}{A}_{im}}{\mathrm{d}t}& ={\gamma }_{im}{I}_{im}-({\mu }_{im}+d_{im})A_{im}+\sum _{j=1,j\ne i}^n(m_{ij}{A}_{jm}-m_{ji}{A}_{im}),\end{array}$(1) ) is uniformly persistent with respect to $(X_0,\mathrm{\partial }{X}_0)$.

3.4. Uniqueness and global stability of the endemic equilibrium

In this section, if $R_0>1$, we derive a sufficient condition that the endemic equilibrium is unique and globally asymptotically stable. Our proof utilizes the graph-theoretical approach developed in [16, 17, 23].

Theorem 3.4

Assume that $M=(m_{ij})$ is irreducible, and there exist ${\eta }_{j1},{\eta }_{j2}>0$ such that $S_{jf}^{\ast }={\eta }_{j1}{I}_{jf}^{\ast },S_{jm}^{\ast }={\eta }_{j2}{I}_{jm}^{\ast }$ for all $1\le i,j\le n$, then the endemic equilibrium $E^{\ast }$ of system (1(1) $\begin{array}{rl}\frac{\mathrm{d}{S}_{if}}{\mathrm{d}t}& ={\epsilon }_i{\mathrm{\Lambda }}_i-{\beta }_i\frac{S_{if}{I}_{im}}{N_{im}}-{\mu }_{if}{S}_{if}+\sum _{j=1,j\ne i}^n(m_{ij}{S}_{jf}-m_{ji}{S}_{if}),\\ \frac{\mathrm{d}{I}_{if}}{\mathrm{d}t}& ={\beta }_i\frac{S_{if}{I}_{im}}{N_{im}}-({\mu }_{if}+{\gamma }_{if})I_{if}+\sum _{j=1,j\ne i}^n(m_{ij}{I}_{jf}-m_{ji}{I}_{if}),\\ \frac{\mathrm{d}{A}_{if}}{\mathrm{d}t}& ={\gamma }_{if}{I}_{if}-({\mu }_{if}+d_{if})A_{if}+\sum _{j=1,j\ne i}^n(m_{ij}{A}_{jf}-m_{ji}{A}_{if}),\\ \frac{\mathrm{d}{S}_{im}}{\mathrm{d}t}& =(1-{\epsilon }_i){\mathrm{\Lambda }}_i-{\beta }_i\frac{S_{im}{I}_{if}}{N_{if}}-{\mu }_{im}{S}_{im}+\sum _{j=1,j\ne i}^n(m_{ij}{S}_{jm}-m_{ji}{S}_{im}),\\ \frac{\mathrm{d}{I}_{im}}{\mathrm{d}t}& ={\beta }_i\frac{S_{im}{I}_{if}}{N_{if}}-({\mu }_{im}+{\gamma }_{im})I_{im}+\sum _{j=1,j\ne i}^n(m_{ij}{I}_{jm}-m_{ji}{I}_{im}),\\ \frac{\mathrm{d}{A}_{im}}{\mathrm{d}t}& ={\gamma }_{im}{I}_{im}-({\mu }_{im}+d_{im})A_{im}+\sum _{j=1,j\ne i}^n(m_{ij}{A}_{jm}-m_{ji}{A}_{im}),\end{array}$(1) ) is unique and globally asymptotically stable in int(X) when $R_0>1$.

The proof process of Theorem 3.4, see Appendix 1.

4. Application to two patches

In the case where the patch number n is 2, (1(1) $\begin{array}{rl}\frac{\mathrm{d}{S}_{if}}{\mathrm{d}t}& ={\epsilon }_i{\mathrm{\Lambda }}_i-{\beta }_i\frac{S_{if}{I}_{im}}{N_{im}}-{\mu }_{if}{S}_{if}+\sum _{j=1,j\ne i}^n(m_{ij}{S}_{jf}-m_{ji}{S}_{if}),\\ \frac{\mathrm{d}{I}_{if}}{\mathrm{d}t}& ={\beta }_i\frac{S_{if}{I}_{im}}{N_{im}}-({\mu }_{if}+{\gamma }_{if})I_{if}+\sum _{j=1,j\ne i}^n(m_{ij}{I}_{jf}-m_{ji}{I}_{if}),\\ \frac{\mathrm{d}{A}_{if}}{\mathrm{d}t}& ={\gamma }_{if}{I}_{if}-({\mu }_{if}+d_{if})A_{if}+\sum _{j=1,j\ne i}^n(m_{ij}{A}_{jf}-m_{ji}{A}_{if}),\\ \frac{\mathrm{d}{S}_{im}}{\mathrm{d}t}& =(1-{\epsilon }_i){\mathrm{\Lambda }}_i-{\beta }_i\frac{S_{im}{I}_{if}}{N_{if}}-{\mu }_{im}{S}_{im}+\sum _{j=1,j\ne i}^n(m_{ij}{S}_{jm}-m_{ji}{S}_{im}),\\ \frac{\mathrm{d}{I}_{im}}{\mathrm{d}t}& ={\beta }_i\frac{S_{im}{I}_{if}}{N_{if}}-({\mu }_{im}+{\gamma }_{im})I_{im}+\sum _{j=1,j\ne i}^n(m_{ij}{I}_{jm}-m_{ji}{I}_{im}),\\ \frac{\mathrm{d}{A}_{im}}{\mathrm{d}t}& ={\gamma }_{im}{I}_{im}-({\mu }_{im}+d_{im})A_{im}+\sum _{j=1,j\ne i}^n(m_{ij}{A}_{jm}-m_{ji}{A}_{im}),\end{array}$(1) ) reduces to (4) $\begin{array}{rl}\frac{\mathrm{d}{S}_{1f}}{\mathrm{d}t}& ={\epsilon }_1{\mathrm{\Lambda }}_1-{\beta }_1\frac{S_{1f}{I}_{1m}}{N_{1m}}-{\mu }_{1f}{S}_{1f}+m_{12}{S}_{2f}-m_{21}{S}_{1f},\\ \frac{\mathrm{d}{I}_{1f}}{\mathrm{d}t}& ={\beta }_1\frac{S_{1f}{I}_{1m}}{N_{1m}}-({\mu }_{1f}+{\gamma }_{1f})I_{1f}+m_{12}{I}_{2f}-m_{21}{I}_{1f},\\ \frac{\mathrm{d}{A}_{1f}}{\mathrm{d}t}& ={\gamma }_{1f}{I}_{1f}-({\mu }_{1f}+d_{1f})A_{1f}+m_{12}{A}_{2f}-m_{21}{A}_{1f},\\ \frac{\mathrm{d}{S}_{1m}}{\mathrm{d}t}& =(1-{\epsilon }_1){\mathrm{\Lambda }}_1-{\beta }_1\frac{S_{1m}{I}_{1f}}{N_{1f}}-{\mu }_{1m}{S}_{1m}+m_{12}{S}_{2m}-m_{21}{S}_{1m},\\ \frac{\mathrm{d}{I}_{1m}}{\mathrm{d}t}& ={\beta }_1\frac{S_{1m}{I}_{1f}}{N_{1f}}-({\mu }_{1m}+{\gamma }_{1m})I_{1m}+m_{12}{I}_{2m}-m_{21}{I}_{1m},\\ \frac{\mathrm{d}{A}_{1m}}{\mathrm{d}t}& ={\gamma }_{1m}{I}_{1m}-({\mu }_{1m}+d_{1m})A_{1m}+m_{12}{A}_{2m}-m_{21}{A}_{1m},\\ \frac{\mathrm{d}{S}_{2f}}{\mathrm{d}t}& ={\epsilon }_2{\mathrm{\Lambda }}_2-{\beta }_2\frac{S_{2f}{I}_{2m}}{N_{2m}}-{\mu }_{2f}{S}_{2f}+m_{21}{S}_{1f}-m_{12}{S}_{2f},\\ \frac{\mathrm{d}{I}_{2f}}{\mathrm{d}t}& ={\beta }_2\frac{S_{2f}{I}_{2m}}{N_{2m}}-({\mu }_{2f}+{\gamma }_{2f})I_{2f}+m_{21}{I}_{1f}-m_{12}{I}_{2f},\\ \frac{\mathrm{d}{A}_{2f}}{\mathrm{d}t}& ={\gamma }_{2f}{I}_{2f}-({\mu }_{2f}+d_{2f})A_{2f}+m_{21}{A}_{1f}-m_{12}{A}_{2f},\\ \frac{\mathrm{d}{S}_{2m}}{\mathrm{d}t}& =(1-{\epsilon }_2){\mathrm{\Lambda }}_2-{\beta }_2\frac{S_{2m}{I}_{2f}}{N_{2f}}-{\mu }_{2m}{S}_{2m}+m_{21}{S}_{1m}-m_{12}{S}_{2m},\\ \frac{\mathrm{d}{I}_{2m}}{\mathrm{d}t}& ={\beta }_2\frac{S_{2m}{I}_{2f}}{N_{2f}}-({\mu }_{2m}+{\gamma }_{2m})I_{2m}+m_{21}{I}_{1m}-m_{12}{I}_{2m},\\ \frac{\mathrm{d}{A}_{2m}}{\mathrm{d}t}& ={\gamma }_{2m}{I}_{2m}-({\mu }_{2m}+d_{2m})A_{2m}+m_{21}{A}_{1m}-m_{12}{A}_{2m},\end{array}$(4) where the parameters $m_{21}$ and $m_{12}$ denote, respectively, the migration rates from patches 1 to 2 and from patches 2 to 1.

If $R_0<1$, the unique disease-free equilibrium of (4(4) $\begin{array}{rl}\frac{\mathrm{d}{S}_{1f}}{\mathrm{d}t}& ={\epsilon }_1{\mathrm{\Lambda }}_1-{\beta }_1\frac{S_{1f}{I}_{1m}}{N_{1m}}-{\mu }_{1f}{S}_{1f}+m_{12}{S}_{2f}-m_{21}{S}_{1f},\\ \frac{\mathrm{d}{I}_{1f}}{\mathrm{d}t}& ={\beta }_1\frac{S_{1f}{I}_{1m}}{N_{1m}}-({\mu }_{1f}+{\gamma }_{1f})I_{1f}+m_{12}{I}_{2f}-m_{21}{I}_{1f},\\ \frac{\mathrm{d}{A}_{1f}}{\mathrm{d}t}& ={\gamma }_{1f}{I}_{1f}-({\mu }_{1f}+d_{1f})A_{1f}+m_{12}{A}_{2f}-m_{21}{A}_{1f},\\ \frac{\mathrm{d}{S}_{1m}}{\mathrm{d}t}& =(1-{\epsilon }_1){\mathrm{\Lambda }}_1-{\beta }_1\frac{S_{1m}{I}_{1f}}{N_{1f}}-{\mu }_{1m}{S}_{1m}+m_{12}{S}_{2m}-m_{21}{S}_{1m},\\ \frac{\mathrm{d}{I}_{1m}}{\mathrm{d}t}& ={\beta }_1\frac{S_{1m}{I}_{1f}}{N_{1f}}-({\mu }_{1m}+{\gamma }_{1m})I_{1m}+m_{12}{I}_{2m}-m_{21}{I}_{1m},\\ \frac{\mathrm{d}{A}_{1m}}{\mathrm{d}t}& ={\gamma }_{1m}{I}_{1m}-({\mu }_{1m}+d_{1m})A_{1m}+m_{12}{A}_{2m}-m_{21}{A}_{1m},\\ \frac{\mathrm{d}{S}_{2f}}{\mathrm{d}t}& ={\epsilon }_2{\mathrm{\Lambda }}_2-{\beta }_2\frac{S_{2f}{I}_{2m}}{N_{2m}}-{\mu }_{2f}{S}_{2f}+m_{21}{S}_{1f}-m_{12}{S}_{2f},\\ \frac{\mathrm{d}{I}_{2f}}{\mathrm{d}t}& ={\beta }_2\frac{S_{2f}{I}_{2m}}{N_{2m}}-({\mu }_{2f}+{\gamma }_{2f})I_{2f}+m_{21}{I}_{1f}-m_{12}{I}_{2f},\\ \frac{\mathrm{d}{A}_{2f}}{\mathrm{d}t}& ={\gamma }_{2f}{I}_{2f}-({\mu }_{2f}+d_{2f})A_{2f}+m_{21}{A}_{1f}-m_{12}{A}_{2f},\\ \frac{\mathrm{d}{S}_{2m}}{\mathrm{d}t}& =(1-{\epsilon }_2){\mathrm{\Lambda }}_2-{\beta }_2\frac{S_{2m}{I}_{2f}}{N_{2f}}-{\mu }_{2m}{S}_{2m}+m_{21}{S}_{1m}-m_{12}{S}_{2m},\\ \frac{\mathrm{d}{I}_{2m}}{\mathrm{d}t}& ={\beta }_2\frac{S_{2m}{I}_{2f}}{N_{2f}}-({\mu }_{2m}+{\gamma }_{2m})I_{2m}+m_{21}{I}_{1m}-m_{12}{I}_{2m},\\ \frac{\mathrm{d}{A}_{2m}}{\mathrm{d}t}& ={\gamma }_{2m}{I}_{2m}-({\mu }_{2m}+d_{2m})A_{2m}+m_{21}{A}_{1m}-m_{12}{A}_{2m},\end{array}$(4) ) is $E_0=(S_{1f}^0,0,0,S_{1m}^0,0,0,S_{2f}^0,0,0,S_{2m}^0,0,0),$where $\begin{array}{rl}{S}_{1f}^0& =\frac{({\mu }_{2f}+m_{12}){\epsilon }_1{\mathrm{\Lambda }}_1+m_{12}{\epsilon }_2{\mathrm{\Lambda }}_2}{({\mu }_{1f}+m_{21})({\mu }_{2f}+m_{12})-m_{12}{m}_{21}},\\ S_{2f}^0& =\frac{({\mu }_{1f}+m_{21}){\epsilon }_2{\mathrm{\Lambda }}_2+m_{21}{\epsilon }_1{\mathrm{\Lambda }}_1}{({\mu }_{1f}+m_{21})({\mu }_{2f}+m_{12})-m_{12}{m}_{21}},\\ S_{1m}^0& =\frac{({\mu }_{2m}+m_{12})(1-{\epsilon }_1){\mathrm{\Lambda }}_1+m_{12}(1-{\epsilon }_2){\mathrm{\Lambda }}_2}{({\mu }_{1m}+m_{21})({\mu }_{2m}+m_{12})-m_{12}{m}_{21}},\\ S_{2m}^0& =\frac{({\mu }_{1m}+m_{21})(1-{\epsilon }_2){\mathrm{\Lambda }}_2+m_{21}(1-{\epsilon }_1){\mathrm{\Lambda }}_1}{({\mu }_{1m}+m_{21})({\mu }_{2m}+m_{12})-m_{12}{m}_{21}}.\end{array}$Now $F_{12}{V}_{22}^{-1}{F}_{21}{V}_{11}^{-1}$ becomes $\left[\begin{array}{cc}a& b\\ c& d\end{array}\right]$, where $\begin{array}{rl}a& =\frac{1}{M_1{M}_2}[{\beta }_1^2({\mu }_{2f}+{\gamma }_{2f}+m_{12})({\mu }_{2m}+{\gamma }_{2m}+m_{12})+{\beta }_1{\beta }_2{m}_{12}{m}_{21}\frac{S_{1f}^0{S}_{2m}^0}{S_{1m}^0{S}_{2f}^0}],\\ b& =\frac{1}{M_1{M}_2}[{\beta }_1^2{m}_{12}({\mu }_{2m}+{\gamma }_{2m}+m_{12})+{\beta }_1{\beta }_2{m}_{12}({\mu }_{1f}+{\gamma }_{1f}+m_{21})\frac{S_{1f}^0{S}_{2m}^0}{S_{1m}^0{S}_{2f}^0}],\\ c& =\frac{1}{M_1{M}_2}[{\beta }_1{\beta }_2{m}_{21}({\mu }_{2f}+{\gamma }_{2f}+m_{12})\frac{S_{1m}^0{S}_{2f}^0}{S_{1f}^0{S}_{2m}^0}+{\beta }_2^2{m}_{21}({\mu }_{1m}+{\gamma }_{1m}+m_{21})],\\ d& =\frac{1}{M_1{M}_2}[{\beta }_1{\beta }_2{m}_{12}{m}_{21}\frac{S_{1m}^0{S}_{2f}^0}{S_{1f}^0{S}_{2m}^0}+{\beta }_2^2({\mu }_{1f}+{\gamma }_{1f}+m_{21})({\mu }_{1m}+{\gamma }_{1m}+m_{21})],\end{array}$where $\begin{array}{rl}{M}_1& =({\mu }_{1f}+{\gamma }_{1f}+m_{21})({\mu }_{2f}+{\gamma }_{2f}+m_{12})-m_{12}{m}_{21},\\ M_2& =({\mu }_{1m}+{\gamma }_{1m}+m_{21})({\mu }_{2m}+{\gamma }_{2m}+m_{12})-m_{12}{m}_{21}.\end{array}$Its characteristic equation is ${\lambda }^2-(a+d)\lambda +ad-bc=0,$ we have $\begin{array}{rl}{\lambda }_1& =\frac{a+d+\sqrt{(a+d{)}^2-4(ad-bc)}}{2},\\ {\lambda }_2& =\frac{a+d-\sqrt{(a+d{)}^2-4(ad-bc)}}{2}.\end{array}$It is easy to see that $R_0=\sqrt{\rho (F_{12}{V}_{22}^{-1}{F}_{21}{V}_{11}^{-1})}=\sqrt{{\lambda }_1}.$

When the relationship between $m_{21}$ and $m_{12}$ satisfies bc−ad + a + d−1 = 0, $R_0=1$.

The basic reproduction number for the ith patch in isolation (i.e. there is no migration between patch i and other patches) is given by $R_0^{(1)}=\frac{{\beta }_1}{\sqrt{({\mu }_{1f}+{\gamma }_{1f})({\mu }_{1m}+{\gamma }_{1m})}},R_0^{(2)}=\frac{{\beta }_2}{\sqrt{({\mu }_{2f}+{\gamma }_{2f})({\mu }_{2m}+{\gamma }_{2m})}}.$It is well known that $R_0^{(1)}$ is a reproduction number of the disease in the first patch and $R_0^{(2)}$ is a reproduction number of the disease in the second patch.

If $R_0>1$, and there exist ${\eta }_{11},{\eta }_{21},{\eta }_{12},{\eta }_{22}>0$ such that $S_{1f}^{\ast }={\eta }_{11}{I}_{1f}^{\ast },S_{2f}^{\ast }={\eta }_{21}{I}_{2f}^{\ast },S_{1m}^{\ast }={\eta }_{12}{I}_{1m}^{\ast }$ and $S_{2m}^{\ast }={\eta }_{22}{I}_{2m}^{\ast }$, the unique endemic equilibrium of (4(4) $\begin{array}{rl}\frac{\mathrm{d}{S}_{1f}}{\mathrm{d}t}& ={\epsilon }_1{\mathrm{\Lambda }}_1-{\beta }_1\frac{S_{1f}{I}_{1m}}{N_{1m}}-{\mu }_{1f}{S}_{1f}+m_{12}{S}_{2f}-m_{21}{S}_{1f},\\ \frac{\mathrm{d}{I}_{1f}}{\mathrm{d}t}& ={\beta }_1\frac{S_{1f}{I}_{1m}}{N_{1m}}-({\mu }_{1f}+{\gamma }_{1f})I_{1f}+m_{12}{I}_{2f}-m_{21}{I}_{1f},\\ \frac{\mathrm{d}{A}_{1f}}{\mathrm{d}t}& ={\gamma }_{1f}{I}_{1f}-({\mu }_{1f}+d_{1f})A_{1f}+m_{12}{A}_{2f}-m_{21}{A}_{1f},\\ \frac{\mathrm{d}{S}_{1m}}{\mathrm{d}t}& =(1-{\epsilon }_1){\mathrm{\Lambda }}_1-{\beta }_1\frac{S_{1m}{I}_{1f}}{N_{1f}}-{\mu }_{1m}{S}_{1m}+m_{12}{S}_{2m}-m_{21}{S}_{1m},\\ \frac{\mathrm{d}{I}_{1m}}{\mathrm{d}t}& ={\beta }_1\frac{S_{1m}{I}_{1f}}{N_{1f}}-({\mu }_{1m}+{\gamma }_{1m})I_{1m}+m_{12}{I}_{2m}-m_{21}{I}_{1m},\\ \frac{\mathrm{d}{A}_{1m}}{\mathrm{d}t}& ={\gamma }_{1m}{I}_{1m}-({\mu }_{1m}+d_{1m})A_{1m}+m_{12}{A}_{2m}-m_{21}{A}_{1m},\\ \frac{\mathrm{d}{S}_{2f}}{\mathrm{d}t}& ={\epsilon }_2{\mathrm{\Lambda }}_2-{\beta }_2\frac{S_{2f}{I}_{2m}}{N_{2m}}-{\mu }_{2f}{S}_{2f}+m_{21}{S}_{1f}-m_{12}{S}_{2f},\\ \frac{\mathrm{d}{I}_{2f}}{\mathrm{d}t}& ={\beta }_2\frac{S_{2f}{I}_{2m}}{N_{2m}}-({\mu }_{2f}+{\gamma }_{2f})I_{2f}+m_{21}{I}_{1f}-m_{12}{I}_{2f},\\ \frac{\mathrm{d}{A}_{2f}}{\mathrm{d}t}& ={\gamma }_{2f}{I}_{2f}-({\mu }_{2f}+d_{2f})A_{2f}+m_{21}{A}_{1f}-m_{12}{A}_{2f},\\ \frac{\mathrm{d}{S}_{2m}}{\mathrm{d}t}& =(1-{\epsilon }_2){\mathrm{\Lambda }}_2-{\beta }_2\frac{S_{2m}{I}_{2f}}{N_{2f}}-{\mu }_{2m}{S}_{2m}+m_{21}{S}_{1m}-m_{12}{S}_{2m},\\ \frac{\mathrm{d}{I}_{2m}}{\mathrm{d}t}& ={\beta }_2\frac{S_{2m}{I}_{2f}}{N_{2f}}-({\mu }_{2m}+{\gamma }_{2m})I_{2m}+m_{21}{I}_{1m}-m_{12}{I}_{2m},\\ \frac{\mathrm{d}{A}_{2m}}{\mathrm{d}t}& ={\gamma }_{2m}{I}_{2m}-({\mu }_{2m}+d_{2m})A_{2m}+m_{21}{A}_{1m}-m_{12}{A}_{2m},\end{array}$(4) ) is $E^{\ast }=(S_{1f}^{\ast },I_{1f}^{\ast },A_{1f}^{\ast },S_{1m}^{\ast },I_{1m}^{\ast },A_{1m}^{\ast },S_{2f}^{\ast },I_{2f}^{\ast },A_{2f}^{\ast },S_{2m}^{\ast },I_{2m}^{\ast },A_{2m}^{\ast }),$where $\begin{array}{rl}{I}_{1f}^{\ast }& =\frac{q_f[({\mu }_{2f}+m_{12}){\epsilon }_1{\mathrm{\Lambda }}_1+m_{12}{\epsilon }_2{\mathrm{\Lambda }}_2]+m_{12}{\gamma }_{2f}[m_{21}{\epsilon }_1{\mathrm{\Lambda }}_1+({\mu }_{1f}+m_{21}){\epsilon }_2{\mathrm{\Lambda }}_2]}{m_{21}{m}_{12}{\gamma }_{1f}{\gamma }_{2f}-p_f{q}_f},\\ I_{1m}^{\ast }& =\frac{\begin{array}{c}{q}_m[({\mu }_{2m}+m_{12})(1-{\epsilon }_1){\mathrm{\Lambda }}_1+m_{12}(1-{\epsilon }_2){\mathrm{\Lambda }}_2]+m_{12}{\gamma }_{2m}\\ [m_{21}(1-{\epsilon }_1){\mathrm{\Lambda }}_1+({\mu }_{1m}+m_{21})(1-{\epsilon }_2){\mathrm{\Lambda }}_2]\end{array}}{m_{21}{m}_{12}{\gamma }_{1m}{\gamma }_{2m}-p_m{q}_m},\\ I_{2f}^{\ast }& =\frac{m_{21}{\gamma }_{1f}[({\mu }_{2f}+m_{12}){\epsilon }_1{\mathrm{\Lambda }}_1+m_{12}{\epsilon }_2{\mathrm{\Lambda }}_2]+p_f[m_{21}{\epsilon }_1{\mathrm{\Lambda }}_1+({\mu }_{1f}+m_{21}){\epsilon }_2{\mathrm{\Lambda }}_2]}{m_{21}{m}_{12}{\gamma }_{1f}{\gamma }_{2f}-p_f{q}_f},\\ I_{2m}^{\ast }& =\frac{\begin{array}{c}{m}_{21}{\gamma }_{1m}[({\mu }_{2m}+m_{12})(1-{\epsilon }_1){\mathrm{\Lambda }}_1+m_{12}(1-{\epsilon }_2){\mathrm{\Lambda }}_2]\\ +p_m[m_{21}(1-{\epsilon }_1){\mathrm{\Lambda }}_1+({\mu }_{1m}+m_{21})(1-{\epsilon }_2){\mathrm{\Lambda }}_2]\end{array}}{m_{21}{m}_{12}{\gamma }_{1m}{\gamma }_{2m}-p_m{q}_m},\\ S_{1f}^{\ast }& =\frac{({\mu }_{2f}+m_{12}){\epsilon }_1{\mathrm{\Lambda }}_1+m_{12}{\epsilon }_2{\mathrm{\Lambda }}_2+[m_{21}{m}_{12}-({\mu }_{2f}+m_{12})({\mu }_{1f}+{\gamma }_{1f}+m_{21})]I_{1f}^{\ast }-m_{12}{\gamma }_{2f}{I}_{2f}^{\ast }}{({\mu }_{1f}+m_{21})({\mu }_{2f}+m_{12})-m_{21}{m}_{12}},\\ S_{1m}^{\ast }& =\frac{\begin{array}{c}({\mu }_{2m}+m_{12})(1-{\epsilon }_1){\mathrm{\Lambda }}_1+m_{12}(1-{\epsilon }_2){\mathrm{\Lambda }}_2\\ +[m_{21}{m}_{12}-({\mu }_{2m}+m_{12})({\mu }_{1m}+{\gamma }_{1m}+m_{21})]I_{1m}^{\ast }-m_{12}{\gamma }_{2m}{I}_{2m}^{\ast }\end{array}}{({\mu }_{1m}+m_{21})({\mu }_{2m}+m_{12})-m_{21}{m}_{12}},\\ S_{2f}^{\ast }& =\frac{m_{21}{\epsilon }_1{\mathrm{\Lambda }}_1+({\mu }_{1f}+m_{21}){\epsilon }_2{\mathrm{\Lambda }}_2-m_{21}{\gamma }_{1f}{I}_{1f}^{\ast }+[m_{21}{m}_{12}-({\mu }_{1f}+m_{21})({\mu }_{2f}+{\gamma }_{2f}+m_{12})]I_{2f}^{\ast }}{({\mu }_{1f}+m_{21})({\mu }_{2f}+m_{12})-m_{21}{m}_{12}},\\ S_{2m}^{\ast }& =\frac{\begin{array}{c}{m}_{21}(1-{\epsilon }_1){\mathrm{\Lambda }}_1+({\mu }_{1m}+m_{21})(1-{\epsilon }_2){\mathrm{\Lambda }}_2-m_{21}{\gamma }_{1m}{I}_{1m}^{\ast }\\ +[m_{21}{m}_{12}-({\mu }_{1m}+m_{21})({\mu }_{2m}+{\gamma }_{2m}+m_{12})]I_{2m}^{\ast }\end{array}}{({\mu }_{1m}+m_{21})({\mu }_{2m}+m_{12})-m_{21}{m}_{12}},\\ A_{1f}^{\ast }& =\frac{({\mu }_{2f}+d_{2f}+m_{12}){\gamma }_{1f}{I}_{1f}^{\ast }+m_{12}{\gamma }_{2f}{I}_{2f}^{\ast }}{({\mu }_{1f}+d_{1f}+m_{21})({\mu }_{2f}+d_{2f}+m_{12})-m_{21}{m}_{12}},\\ A_{1m}^{\ast }& =\frac{({\mu }_{2m}+d_{2m}+m_{12}){\gamma }_{1m}{I}_{1m}^{\ast }+m_{12}{\gamma }_{2m}{I}_{2m}^{\ast }}{({\mu }_{1m}+d_{1m}+m_{21})({\mu }_{2m}+d_{2m}+m_{12})-m_{21}{m}_{12}},\\ A_{2f}^{\ast }& =\frac{m_{21}{\gamma }_{1f}{I}_{1f}^{\ast }+({\mu }_{1f}+d_{1f}+m_{21}){\gamma }_{2f}{I}_{2f}^{\ast }}{({\mu }_{1f}+d_{1f}+m_{21})({\mu }_{2f}+d_{2f}+m_{12})-m_{21}{m}_{12}},\\ A_{2m}^{\ast }& =\frac{m_{21}{\gamma }_{1m}{I}_{1m}^{\ast }+({\mu }_{1m}+d_{1m}+m_{21}){\gamma }_{2m}{I}_{2m}^{\ast }}{({\mu }_{1m}+d_{1m}+m_{21})({\mu }_{2m}+d_{2m}+m_{12})-m_{21}{m}_{12}},\end{array}$

where $\begin{array}{rl}{p}_f& =({\eta }_{11}+1)m_{21}{m}_{12}-({\mu }_{2f}+m_{12})({\mu }_{1f}+{\gamma }_{1f}+m_{21})\\ & -{\eta }_{11}({\mu }_{1f}+m_{21})({\mu }_{2f}+m_{12}),\\ q_f& =({\eta }_{21}+1)m_{21}{m}_{12}-({\mu }_{1f}+m_{21})({\mu }_{2f}+{\gamma }_{2f}+m_{12})\\ & -{\eta }_{21}({\mu }_{1f}+m_{21})({\mu }_{2f}+m_{12}),\\ p_m& =({\eta }_{12}+1)m_{21}{m}_{12}-({\mu }_{2m}+m_{12})({\mu }_{1m}+{\gamma }_{1m}+m_{21})\\ & -{\eta }_{12}({\mu }_{1m}+m_{21})({\mu }_{2m}+m_{12}),\\ q_m& =({\eta }_{22}+1)m_{21}{m}_{12}-({\mu }_{1m}+m_{21})({\mu }_{2m}+{\gamma }_{2m}+m_{12})\\ & -{\eta }_{22}({\mu }_{1m}+m_{21})({\mu }_{2m}+m_{12}).\end{array}$

5. Numerical simulations

To have a more comprehensive understanding of the spread of HIV/AIDS between patches, we conduct numerical simulations for different situations to present various possible results. We introduce it from three aspects: in Section 5.1, we present time series diagrams under different conditions; in Section 5.2, we use contour plots to show the relationship between $R_0$ and migration rates; in Section 5.3, we make a comparison between no migration and migration.

We consider a particular case of system (1(1) $\begin{array}{rl}\frac{\mathrm{d}{S}_{if}}{\mathrm{d}t}& ={\epsilon }_i{\mathrm{\Lambda }}_i-{\beta }_i\frac{S_{if}{I}_{im}}{N_{im}}-{\mu }_{if}{S}_{if}+\sum _{j=1,j\ne i}^n(m_{ij}{S}_{jf}-m_{ji}{S}_{if}),\\ \frac{\mathrm{d}{I}_{if}}{\mathrm{d}t}& ={\beta }_i\frac{S_{if}{I}_{im}}{N_{im}}-({\mu }_{if}+{\gamma }_{if})I_{if}+\sum _{j=1,j\ne i}^n(m_{ij}{I}_{jf}-m_{ji}{I}_{if}),\\ \frac{\mathrm{d}{A}_{if}}{\mathrm{d}t}& ={\gamma }_{if}{I}_{if}-({\mu }_{if}+d_{if})A_{if}+\sum _{j=1,j\ne i}^n(m_{ij}{A}_{jf}-m_{ji}{A}_{if}),\\ \frac{\mathrm{d}{S}_{im}}{\mathrm{d}t}& =(1-{\epsilon }_i){\mathrm{\Lambda }}_i-{\beta }_i\frac{S_{im}{I}_{if}}{N_{if}}-{\mu }_{im}{S}_{im}+\sum _{j=1,j\ne i}^n(m_{ij}{S}_{jm}-m_{ji}{S}_{im}),\\ \frac{\mathrm{d}{I}_{im}}{\mathrm{d}t}& ={\beta }_i\frac{S_{im}{I}_{if}}{N_{if}}-({\mu }_{im}+{\gamma }_{im})I_{im}+\sum _{j=1,j\ne i}^n(m_{ij}{I}_{jm}-m_{ji}{I}_{im}),\\ \frac{\mathrm{d}{A}_{im}}{\mathrm{d}t}& ={\gamma }_{im}{I}_{im}-({\mu }_{im}+d_{im})A_{im}+\sum _{j=1,j\ne i}^n(m_{ij}{A}_{jm}-m_{ji}{A}_{im}),\end{array}$(1) ) with two patches. We study the impact of population migration on the spread of HIV/AIDS epidemic with heterosexual transmission. At the same time, we also verify the theoretical analysis. We choose the data from the Beijing City and Yunnan Province, Beijing City is regarded as the first patch, Yunnan Province is regarded as the second patch. We choose 2017 as the initial time. The initial population size is given as follows: $S_{1f}^0=10704298,I_{1f}^0=538A_{1f}^0=164S_{1m}^0=11236585I_{1m}^0=1833A_{1m}^0=582S_{2f}^0=23092948I_{2f}^0=1804A_{2f}^0=1248S_{2m}^0=24898433I_{2m}^0=6142A_{2m}^0=4425$ respectively. The data is from Chinese Center for Disease Control and Prevention [4], Beijing Municipal Bureau of Statistics [2] and Yunnan Provincial Bureau of Statistics [3]. Given that women generally live longer than men and men have higher death rates due to AIDS, we choose ${\mu }_{1f}=0.0120,{\mu }_{1m}=0.0127,{\mu }_{2f}=0.0122,{\mu }_{2m}=0.0130,d_{1f}=0.045,d_{1m}=0.049,d_{2f}=0.333,d_{2m}=0.369$, and ${\beta }_i$ is selected from the range of 0.011−0.95 [26]. For the transfer rates from stage $I_{if}$ to $A_{if}$ and from stage $I_{im}$ to $A_{im}$ in the two regions, we select the populations of I (the infected) and A (the AIDS) in the two regions from 2015 to 2017 from the China Center for Disease Control and prevention to calculate their average, consider that men have higher transfer rates from stage I to A than women, we choose ${\gamma }_{1f}=0.30,{\gamma }_{1m}=0.32,{\gamma }_{2f}=0.71,{\gamma }_{2m}=0.73$. For convenience, these model parameters values in numerical simulations are given in Table 2, with a time scale of a year.

Table 2. Model parameter values in numerical simulations.

Parameters	Description	Value
${\mathrm{\Lambda }}_i$	The recruitment rate of the population in patch i per unit time, i = 1, 2	$5\times {10}^5$
${\epsilon }_1$	The fraction of recruited susceptible $S_{1f}$	0.49
${\epsilon }_2$	The fraction of recruited susceptible $S_{2f}$	0.48
${\mu }_{1f}$	The natural death rate of women in patch 1	0.0120
${\mu }_{1m}$	The natural death rate of men in patch 1	0.0127
${\mu }_{2f}$	The natural death rate of women in patch 2	0.0122
${\mu }_{2m}$	The natural death rate of men in patch 2	0.0130
$d_{1f}$	The death rate due to AIDS of women in patch 1	0.045
$d_{1m}$	The death rate due to AIDS of men in patch 1	0.049
$d_{2f}$	The death rate due to AIDS of women in patch 2	0.333
$d_{2m}$	The death rate due to AIDS of men in patch 2	0.369
${\gamma }_{1f}$	The transfer rate from stage $I_{1f}$ to $A_{1f}$	0.30
${\gamma }_{1m}$	The transfer rate from stage $I_{1m}$ to $A_{1m}$	0.32
${\gamma }_{2f}$	The transfer rate from stage $I_{2f}$ to $A_{2f}$	0.71
${\gamma }_{2m}$	The transfer rate from stage $I_{2m}$ to $A_{2m}$	0.73

5.1. Time series diagrams

To fully understand the spread of disease between patches, we will present time series diagrams under different conditions. For system (1(1) $\begin{array}{rl}\frac{\mathrm{d}{S}_{if}}{\mathrm{d}t}& ={\epsilon }_i{\mathrm{\Lambda }}_i-{\beta }_i\frac{S_{if}{I}_{im}}{N_{im}}-{\mu }_{if}{S}_{if}+\sum _{j=1,j\ne i}^n(m_{ij}{S}_{jf}-m_{ji}{S}_{if}),\\ \frac{\mathrm{d}{I}_{if}}{\mathrm{d}t}& ={\beta }_i\frac{S_{if}{I}_{im}}{N_{im}}-({\mu }_{if}+{\gamma }_{if})I_{if}+\sum _{j=1,j\ne i}^n(m_{ij}{I}_{jf}-m_{ji}{I}_{if}),\\ \frac{\mathrm{d}{A}_{if}}{\mathrm{d}t}& ={\gamma }_{if}{I}_{if}-({\mu }_{if}+d_{if})A_{if}+\sum _{j=1,j\ne i}^n(m_{ij}{A}_{jf}-m_{ji}{A}_{if}),\\ \frac{\mathrm{d}{S}_{im}}{\mathrm{d}t}& =(1-{\epsilon }_i){\mathrm{\Lambda }}_i-{\beta }_i\frac{S_{im}{I}_{if}}{N_{if}}-{\mu }_{im}{S}_{im}+\sum _{j=1,j\ne i}^n(m_{ij}{S}_{jm}-m_{ji}{S}_{im}),\\ \frac{\mathrm{d}{I}_{im}}{\mathrm{d}t}& ={\beta }_i\frac{S_{im}{I}_{if}}{N_{if}}-({\mu }_{im}+{\gamma }_{im})I_{im}+\sum _{j=1,j\ne i}^n(m_{ij}{I}_{jm}-m_{ji}{I}_{im}),\\ \frac{\mathrm{d}{A}_{im}}{\mathrm{d}t}& ={\gamma }_{im}{I}_{im}-({\mu }_{im}+d_{im})A_{im}+\sum _{j=1,j\ne i}^n(m_{ij}{A}_{jm}-m_{ji}{A}_{im}),\end{array}$(1) ), we present some cases where the disease dies out or persists in each isolated patch but becomes endemic or extinct, respectively, when there are suitable migration rates between them. We simulate the changes of susceptible, infected and AIDS individuals over time in two patches in each case. $R_0^{(1)}$ and $R_0^{(2)}$ are the basic reproduction numbers for both patches in isolation. If $R_0^{(i)}<1$, then we call patch i a low-risk patch in isolation, while if $R_0^{(i)}>1$, then we call it a high-risk patch in isolation. $R_0$ is the basic reproduction number when two patches are connected.

5.1.1. $R_0^{(1)}<1,R_0^{(2)}<1$

Assuming ${\beta }_1=0.30$, ${\beta }_2=0.69$. We fix the migration rates by letting $m_{21}=0,m_{12}=0,$ thus the basic reproduction numbers for both patches in isolation are $R_0^{(1)}=0.9311<1$, $R_0^{(2)}=0.9419<1$. This shows when there is no migration between the two patches, they are low-risk patches and the disease dies out in both patches.

If we fix the migration rates by letting $m_{21}=0.03,m_{12}=0.01,$ and get $R_0=0.9407<1$. This shows that the disease still dies out in both patches. Time series diagrams of different individuals in two patches are presented in Figure 2.

5.1.2. $R_0^{(1)}>1,R_0^{(2)}>1$

Assuming ${\beta }_1=0.37$, ${\beta }_2=0.82$. We fix the migration rates by letting $m_{21}=0,m_{12}=0,$ thus the basic reproduction numbers for both patches in isolation are $R_0^{(1)}=1.1484>1$, $R_0^{(2)}=1.1194>1$. This shows when there is no migration between the two patches, they are high-risk patches and the disease becomes endemic in both patches. If we fix the migration rates by letting $m_{21}=0.03,m_{12}=0.01,$ and get $R_0=1.1241>1$. This shows that the disease still becomes endemic in both patches and there exists an endemic equilibrium. Time series diagrams of different individuals in two patches are presented in Figure 3.

Figure 2. Time series diagrams of different individuals in two patches with ${\beta }_1=0.30,{\beta }_2=0.69,m_{21}=0.03,m_{12}=0.01$, other parameters are shown in Table 2, and the basic reproduction number $R_0=0.9407<1$: (a) susceptible individuals, (b) infective individuals and (c) AIDS individuals.

5.1.3. $R_0^{(1)}>1,R_0^{(2)}<1$

Assuming ${\beta }_1=0.37$, ${\beta }_2=0.69$. We fix the migration rates by letting $m_{21}=0,m_{12}=0,$ that is, there is no migration between the two patches. Thus the basic reproduction numbers for both patches in isolation are $R_0^{(1)}=1.1484>1$, $R_0^{(2)}=0.9419<1$. This shows that the disease becomes endemic in high-risk patch 1 and dies out in low-risk patch 2. When migration occurs and migration rates are set differently, the spread of disease may be different. If we fix the migration rates by letting $m_{21}=0.48,m_{12}=0.10,$ thus $R_0=0.9616<1$. Therefore, the disease-free equilibrium of (1(1) $\begin{array}{rl}\frac{\mathrm{d}{S}_{if}}{\mathrm{d}t}& ={\epsilon }_i{\mathrm{\Lambda }}_i-{\beta }_i\frac{S_{if}{I}_{im}}{N_{im}}-{\mu }_{if}{S}_{if}+\sum _{j=1,j\ne i}^n(m_{ij}{S}_{jf}-m_{ji}{S}_{if}),\\ \frac{\mathrm{d}{I}_{if}}{\mathrm{d}t}& ={\beta }_i\frac{S_{if}{I}_{im}}{N_{im}}-({\mu }_{if}+{\gamma }_{if})I_{if}+\sum _{j=1,j\ne i}^n(m_{ij}{I}_{jf}-m_{ji}{I}_{if}),\\ \frac{\mathrm{d}{A}_{if}}{\mathrm{d}t}& ={\gamma }_{if}{I}_{if}-({\mu }_{if}+d_{if})A_{if}+\sum _{j=1,j\ne i}^n(m_{ij}{A}_{jf}-m_{ji}{A}_{if}),\\ \frac{\mathrm{d}{S}_{im}}{\mathrm{d}t}& =(1-{\epsilon }_i){\mathrm{\Lambda }}_i-{\beta }_i\frac{S_{im}{I}_{if}}{N_{if}}-{\mu }_{im}{S}_{im}+\sum _{j=1,j\ne i}^n(m_{ij}{S}_{jm}-m_{ji}{S}_{im}),\\ \frac{\mathrm{d}{I}_{im}}{\mathrm{d}t}& ={\beta }_i\frac{S_{im}{I}_{if}}{N_{if}}-({\mu }_{im}+{\gamma }_{im})I_{im}+\sum _{j=1,j\ne i}^n(m_{ij}{I}_{jm}-m_{ji}{I}_{im}),\\ \frac{\mathrm{d}{A}_{im}}{\mathrm{d}t}& ={\gamma }_{im}{I}_{im}-({\mu }_{im}+d_{im})A_{im}+\sum _{j=1,j\ne i}^n(m_{ij}{A}_{jm}-m_{ji}{A}_{im}),\end{array}$(1) ) $E_0$ is globally asymptotically stable, the disease dies out in two patches, which means that patch 1 changes from a high-risk patch to a low-risk patch. In this way, migration can decrease the spread of disease. Time series diagrams of different individuals in two patches are presented in Figure 4.

Figure 3. Time series diagrams of different individuals in two patches with ${\beta }_1=0.37,{\beta }_2=0.82,m_{21}=0.03,m_{12}=0.01$, other parameters are shown in Table 2, and the basic reproduction number $R_0=1.1241>1$: (a) susceptible individuals, (b) infective individuals and (c) AIDS individuals.

If we fix the migration rates by letting $m_{21}=0.02,m_{12}=0.35,$ thus $R_0=1.1303>1$. Hence, the endemic equilibrium of (1(1) $\begin{array}{rl}\frac{\mathrm{d}{S}_{if}}{\mathrm{d}t}& ={\epsilon }_i{\mathrm{\Lambda }}_i-{\beta }_i\frac{S_{if}{I}_{im}}{N_{im}}-{\mu }_{if}{S}_{if}+\sum _{j=1,j\ne i}^n(m_{ij}{S}_{jf}-m_{ji}{S}_{if}),\\ \frac{\mathrm{d}{I}_{if}}{\mathrm{d}t}& ={\beta }_i\frac{S_{if}{I}_{im}}{N_{im}}-({\mu }_{if}+{\gamma }_{if})I_{if}+\sum _{j=1,j\ne i}^n(m_{ij}{I}_{jf}-m_{ji}{I}_{if}),\\ \frac{\mathrm{d}{A}_{if}}{\mathrm{d}t}& ={\gamma }_{if}{I}_{if}-({\mu }_{if}+d_{if})A_{if}+\sum _{j=1,j\ne i}^n(m_{ij}{A}_{jf}-m_{ji}{A}_{if}),\\ \frac{\mathrm{d}{S}_{im}}{\mathrm{d}t}& =(1-{\epsilon }_i){\mathrm{\Lambda }}_i-{\beta }_i\frac{S_{im}{I}_{if}}{N_{if}}-{\mu }_{im}{S}_{im}+\sum _{j=1,j\ne i}^n(m_{ij}{S}_{jm}-m_{ji}{S}_{im}),\\ \frac{\mathrm{d}{I}_{im}}{\mathrm{d}t}& ={\beta }_i\frac{S_{im}{I}_{if}}{N_{if}}-({\mu }_{im}+{\gamma }_{im})I_{im}+\sum _{j=1,j\ne i}^n(m_{ij}{I}_{jm}-m_{ji}{I}_{im}),\\ \frac{\mathrm{d}{A}_{im}}{\mathrm{d}t}& ={\gamma }_{im}{I}_{im}-({\mu }_{im}+d_{im})A_{im}+\sum _{j=1,j\ne i}^n(m_{ij}{A}_{jm}-m_{ji}{A}_{im}),\end{array}$(1) ) $E^{\ast }$ is globally asymptotically stable, the disease becomes endemic in both patches and there exists an endemic equilibrium, which means that patch 2 changes from a low-risk patch to a high-risk patch. In this way, migration can increase the spread of disease. Time series diagrams of different individuals in two patches are presented in Figure 5.

5.1.4. $R_0^{(1)}<1,R_0^{(2)}>1$

Assuming ${\beta }_1=0.30$, ${\beta }_2=0.82$. We fix the migration rates by letting $m_{21}=0,m_{12}=0,$ that is, there is no migration between the two patches. Thus the basic reproduction numbers for both patches in isolation are $R_0^{(1)}=0.9311<1$, $R_0^{(2)}=1.1194>1$. This shows that the disease dies out in low-risk patch 1 and becomes endemic in high-risk patch 2. When migration occurs and migration rates are set differently, the spread of disease may be different.

If we fix the migration rates by letting $m_{21}=0.05,m_{12}=0.70,$ thus $R_0=0.9620<1$. Therefore, the disease-free equilibrium of (1(1) $\begin{array}{rl}\frac{\mathrm{d}{S}_{if}}{\mathrm{d}t}& ={\epsilon }_i{\mathrm{\Lambda }}_i-{\beta }_i\frac{S_{if}{I}_{im}}{N_{im}}-{\mu }_{if}{S}_{if}+\sum _{j=1,j\ne i}^n(m_{ij}{S}_{jf}-m_{ji}{S}_{if}),\\ \frac{\mathrm{d}{I}_{if}}{\mathrm{d}t}& ={\beta }_i\frac{S_{if}{I}_{im}}{N_{im}}-({\mu }_{if}+{\gamma }_{if})I_{if}+\sum _{j=1,j\ne i}^n(m_{ij}{I}_{jf}-m_{ji}{I}_{if}),\\ \frac{\mathrm{d}{A}_{if}}{\mathrm{d}t}& ={\gamma }_{if}{I}_{if}-({\mu }_{if}+d_{if})A_{if}+\sum _{j=1,j\ne i}^n(m_{ij}{A}_{jf}-m_{ji}{A}_{if}),\\ \frac{\mathrm{d}{S}_{im}}{\mathrm{d}t}& =(1-{\epsilon }_i){\mathrm{\Lambda }}_i-{\beta }_i\frac{S_{im}{I}_{if}}{N_{if}}-{\mu }_{im}{S}_{im}+\sum _{j=1,j\ne i}^n(m_{ij}{S}_{jm}-m_{ji}{S}_{im}),\\ \frac{\mathrm{d}{I}_{im}}{\mathrm{d}t}& ={\beta }_i\frac{S_{im}{I}_{if}}{N_{if}}-({\mu }_{im}+{\gamma }_{im})I_{im}+\sum _{j=1,j\ne i}^n(m_{ij}{I}_{jm}-m_{ji}{I}_{im}),\\ \frac{\mathrm{d}{A}_{im}}{\mathrm{d}t}& ={\gamma }_{im}{I}_{im}-({\mu }_{im}+d_{im})A_{im}+\sum _{j=1,j\ne i}^n(m_{ij}{A}_{jm}-m_{ji}{A}_{im}),\end{array}$(1) ) $E_0$ is globally asymptotically stable, the disease dies out in two patches, which means that patch 2 changes from a high-risk patch to a low-risk patch. In this way, migration can decrease the spread of disease. Time series diagrams of different individuals in two patches are presented in Figure 6. If we fix the migration rates by letting $m_{21}=0.01,m_{12}=0.02,$ thus $R_0=1.0947>1$. Therefore, the disease becomes endemic in both patches and there exists an endemic equilibrium, which means that patch 1 changes from a low-risk patch to a high-risk patch. In this way, migration can increase the spread of disease. Time series diagrams of different individuals in two patches are presented in Figure 7.

These cases suggest that the effect of population migration on disease transmission is related to the status of the two patches in isolation. If HIV/AIDS becomes extinct in each patch when two patches are isolated, the disease remains extinct in two patches when the population migration occurs; if HIV/AIDS spreads in each patch when two patches are isolated, the disease remains persistent in two patches when the population migration occurs; if the disease disappears in one patch and spreads in the other patch when they are isolated, the disease can spread or disappear in two patches if migration rates of individuals are suitably chosen.

Figure 4. Time series diagrams of different individuals in two patches with ${\beta }_1=0.37,{\beta }_2=0.69,m_{21}=0.48,m_{12}=0.10$, other parameters are shown in Table 2, and the basic reproduction number $R_0=0.9616<1$: (a) susceptible individuals, (b) infective individuals and (c) AIDS individuals.

Figure 5. Time series diagrams of different individuals in two patches with ${\beta }_1=0.37,{\beta }_2=0.69,m_{21}=0.02,m_{12}=0.35$, other parameters are shown in Table 2, and the basic reproduction number $R_0=1.1303>1$: (a) susceptible individuals, (b) infective individuals and (c) AIDS individuals.

Figure 6. Time series diagrams of different individuals in two patches with ${\beta }_1=0.30,{\beta }_2=0.82,m_{21}=0.05,m_{12}=0.70$, other parameters are shown in Table 2, and the basic reproduction number $R_0=0.9620<1$: (a) susceptible individuals, (b) infective individuals and (c) AIDS individuals.

Figure 7. Time series diagrams of different individuals in two patches with ${\beta }_1=0.30,{\beta }_2=0.82,m_{21}=0.01,m_{12}=0.02$, other parameters are shown in Table 2, and the basic reproduction number $R_0=1.0947>1$: (a) susceptible individuals, (b) infective individuals and (c) AIDS individuals.

5.2. $R_0$ versus migration rates $m_{21}$ and $m_{12}$

Sections 5.1.1 –5.1.4 show that the different migration rates between two patches may have different impacts on the spread of HIV/AIDS. To further understand the impacts of migration rates on the spread of HIV/AIDS, we use contour plots to show the relationship between $R_0$ and migration rates.

5.2.1. $R_0^{(1)}<1,R_0^{(2)}<1$

In Figure 8, when ${\beta }_1=0.30,{\beta }_2=0.69,$ the respective basic reproduction numbers of the disease for both patches in isolation are $R_0^{(1)}=0.9311<1,R_0^{(2)}=0.9419<1$, and no matter how migration rates $m_{21}$ and $m_{12}$ are set, $R_0<1$. Thus if HIV/AIDS becomes extinct in each patch when two patches are isolated, the disease remains extinct in two patches when population migration between two patches occurs.

5.2.2. $R_0^{(1)}>1,R_0^{(2)}>1$

In Figure 9, when ${\beta }_1=0.37,{\beta }_2=0.82,$ the respective basic reproduction numbers of the disease for both patches in isolation are $R_0^{(1)}=1.1484>1,R_0^{(2)}=1.1194>1$, and no matter how migration rates $m_{21}$ and $m_{12}$ are set, $R_0>1$. Thus if HIV/AIDS spreads in each patch when two patches are isolated, the disease remains persistent in two patches when population migration between two patches occurs.

5.2.3. $R_0^{(1)}>1,R_0^{(2)}<1$

In Figure 10, according to Section 5.1.3, when ${\beta }_1=0.37,{\beta }_2=0.69,$ $R_0^{(1)}=1.1484>1,R_0^{(2)}=0.9419<1$. Thus HIV/AIDS spreads in the first patch (Beijing City) and disappears in the second patch (Yunnan Province) when two patches are isolated. The red line denotes ‘$R_0=1$’, that is the relationship between $m_{21}$ and $m_{12}$ satisfies bc−ad + a + d−1 = 0. When the relationship between $m_{21}$ and $m_{12}$ is below the red line, that is bc−ad + a + d−1<0, then $R_0<1$, which means the disease will disappear in the two patches, otherwise it will spread.

5.2.4. $R_0^{(1)}<1,R_0^{(2)}>1$

In Figure 11, according to Section 5.1.4, when ${\beta }_1=0.30,{\beta }_2=0.82,$ $R_0^{(1)}=0.9311<1,R_0^{(2)}=1.1194>1$. Thus HIV/AIDS disappears in the first patch (Beijing City) and spreads in the second patch (Yunnan Province) when two patches are isolated. The red line denotes ‘$R_0=1$’, that is the relationship between $m_{21}$ and $m_{12}$ satisfies bc−ad + a + d−1 = 0. When the relationship between $m_{21}$ and $m_{12}$ is above the red line, that is bc−ad + a + d−1<0, then $R_0<1$, which means the disease will disappear in the two patches, otherwise it will spread.

If two patches are low-risk when isolated, they are still low-risk regardless of migration. If two patches are high-risk when isolated, they are still high-risk regardless of migration. If one patch is low-risk and the other is high-risk in isolation, migration may make two patches low-risk and two patches high-risk. Therefore, when two patches are low-risk in isolation, we hardly need to control the migration rates. When one patch is low-risk and the other patch is high-risk in isolation, we can control the migration rates to make $R_0<1$ according to Figures 10 and 11, so that two patches are low-risk. Otherwise, if the migration rates are not properly controlled, the low-risk patch may become high-risk patch. After that, no matter how the migration rates are controlled, it is difficult to reverse the situation, and the disease will spread in two patches.

Figure 8. The contour plot of the basic reproduction number $R_0$ versus migration rates $m_{21}$ and $m_{12}$ under parameters setting: ${\beta }_1=0.30,{\beta }_2=0.69$. Other parameters are shown in Table 2.

Figure 9. The contour plot of the basic reproduction number $R_0$ versus migration rates $m_{21}$ and $m_{12}$ under parameter setting: ${\beta }_1=0.37,{\beta }_2=0.82$. Other parameters are shown in Table 2.

Figure 10. The contour plot of the basic reproduction number $R_0$ versus migration rates $m_{21}$ and $m_{12}$ under parameters setting: ${\beta }_1=0.37,{\beta }_2=0.69$. Other parameters are shown in Table 2.

Figure 11. The contour plot of the basic reproduction number $R_0$ versus migration rates $m_{21}$ and $m_{12}$ under parameters setting: ${\beta }_1=0.30,{\beta }_2=0.82$. Other parameters are shown in Table 2.

5.3. Comparison between no migration and migration in two patches

We consider ${\beta }_1=0.37,{\beta }_2=0.69$. When there is no migration between the two patches, i.e. $m_{21}=m_{12}=0$, the basic reproduction numbers for both patches in isolation are $R_0^{(1)}=1.1484>1,R_0^{(2)}=0.9419<1$. This shows that when there is no migration between the two patches, the disease becomes endemic in patch 1 and dies out in patch 2.

Figure 12. Comparison between no migration and migration in patch 1 with ${\beta }_1=0.37,R_0^{(1)}=1.1484>1,{\beta }_2=0.69,R_0^{(2)}=0.9419<1,R_0=1.1188>1.$ Other parameters are shown in Table 2: (a) infective individuals and (b) AIDS individuals.

Figure 13. Comparison between no migration and migration in patch 2 with ${\beta }_1=0.37,R_0^{(1)}=1.1484>1,{\beta }_2=0.69,R_0^{(2)}=0.9419<1,R_0=1.1188>1.$ Other parameters are shown in Table 2: (a) infective individuals and (b) AIDS individuals.

Figure 14. Comparison between no migration and migration in patch 1 with ${\beta }_1=0.37,R_0^{(1)}=1.1484>1,{\beta }_2=0.69,R_0^{(2)}=0.9419<1,R_0=0.9800<1.$ Other parameters are shown in Table 2: (a) infective individuals and (b) AIDS individuals.

Figure 15. Comparison between no migration and migration in patch 2 with ${\beta }_1=0.37,R_0^{(1)}=1.1484>1,{\beta }_2=0.69,R_0^{(2)}=0.9419<1,R_0=0.9800<1.$ Other parameters are shown in Table 2: (a) infective individuals and (b) AIDS individuals.

We fix the migration rates by letting $m_{21}=0.01,m_{12}=0.02$, thus $R_0=1.1188>1$. Therefore, the disease becomes endemic in both patches and there exists an endemic equilibrium. Comparison between no migration and migration in patch 1 is presented in Figure 12 and comparison between no migration and migration in patch 2 is presented in Figure 13. We fix the migration rates by letting $m_{21}=0.25,m_{12}=0.10$, thus $R_0=0.9800<1$. Therefore, the disease dies out in two patches and there exists a disease-free equilibrium. Comparison between no migration and migration in patch 1 is presented in Figure 14, comparison between no migration and migration in patch 2 is presented in Figure 15.

It can be seen from Figures 12 and 13 that when migration occurs, the disease is still prevalent in patch 1, and the number of infected individuals becomes more; the disease also begin to spread in patch 2. It can be seen from Figures 14 and 15 that when migration occurs, the disease disappears in patch 1; the disease still disappears in patch 2.

6. Sensitivity analysis of the basic reproduction number $R_0$

To analyse the effects of the parameter values on the basic reproduction number $R_0$ of system (1(1) $\begin{array}{rl}\frac{\mathrm{d}{S}_{if}}{\mathrm{d}t}& ={\epsilon }_i{\mathrm{\Lambda }}_i-{\beta }_i\frac{S_{if}{I}_{im}}{N_{im}}-{\mu }_{if}{S}_{if}+\sum _{j=1,j\ne i}^n(m_{ij}{S}_{jf}-m_{ji}{S}_{if}),\\ \frac{\mathrm{d}{I}_{if}}{\mathrm{d}t}& ={\beta }_i\frac{S_{if}{I}_{im}}{N_{im}}-({\mu }_{if}+{\gamma }_{if})I_{if}+\sum _{j=1,j\ne i}^n(m_{ij}{I}_{jf}-m_{ji}{I}_{if}),\\ \frac{\mathrm{d}{A}_{if}}{\mathrm{d}t}& ={\gamma }_{if}{I}_{if}-({\mu }_{if}+d_{if})A_{if}+\sum _{j=1,j\ne i}^n(m_{ij}{A}_{jf}-m_{ji}{A}_{if}),\\ \frac{\mathrm{d}{S}_{im}}{\mathrm{d}t}& =(1-{\epsilon }_i){\mathrm{\Lambda }}_i-{\beta }_i\frac{S_{im}{I}_{if}}{N_{if}}-{\mu }_{im}{S}_{im}+\sum _{j=1,j\ne i}^n(m_{ij}{S}_{jm}-m_{ji}{S}_{im}),\\ \frac{\mathrm{d}{I}_{im}}{\mathrm{d}t}& ={\beta }_i\frac{S_{im}{I}_{if}}{N_{if}}-({\mu }_{im}+{\gamma }_{im})I_{im}+\sum _{j=1,j\ne i}^n(m_{ij}{I}_{jm}-m_{ji}{I}_{im}),\\ \frac{\mathrm{d}{A}_{im}}{\mathrm{d}t}& ={\gamma }_{im}{I}_{im}-({\mu }_{im}+d_{im})A_{im}+\sum _{j=1,j\ne i}^n(m_{ij}{A}_{jm}-m_{ji}{A}_{im}),\end{array}$(1) ), we perform sensitivity analysis by Latin square sampling and partial rank correlation coefficient (PRCC) methods [24]. We tested for significant PRCC for parameters of $R_0$. When $R_0^{(1)}<1$, PRCC is given in Table 3, and PRCC values of the parameters against the basic reproduction number $R_0$ are shown Figure 16. When $R_0^{(1)}>1$, PRCCs are given in Table 4, and PRCC values of the parameters against the basic reproduction number $R_0$ are shown in Figure 17. The parameters ${\beta }_1$ and ${\beta }_2$ of system (1(1) $\begin{array}{rl}\frac{\mathrm{d}{S}_{if}}{\mathrm{d}t}& ={\epsilon }_i{\mathrm{\Lambda }}_i-{\beta }_i\frac{S_{if}{I}_{im}}{N_{im}}-{\mu }_{if}{S}_{if}+\sum _{j=1,j\ne i}^n(m_{ij}{S}_{jf}-m_{ji}{S}_{if}),\\ \frac{\mathrm{d}{I}_{if}}{\mathrm{d}t}& ={\beta }_i\frac{S_{if}{I}_{im}}{N_{im}}-({\mu }_{if}+{\gamma }_{if})I_{if}+\sum _{j=1,j\ne i}^n(m_{ij}{I}_{jf}-m_{ji}{I}_{if}),\\ \frac{\mathrm{d}{A}_{if}}{\mathrm{d}t}& ={\gamma }_{if}{I}_{if}-({\mu }_{if}+d_{if})A_{if}+\sum _{j=1,j\ne i}^n(m_{ij}{A}_{jf}-m_{ji}{A}_{if}),\\ \frac{\mathrm{d}{S}_{im}}{\mathrm{d}t}& =(1-{\epsilon }_i){\mathrm{\Lambda }}_i-{\beta }_i\frac{S_{im}{I}_{if}}{N_{if}}-{\mu }_{im}{S}_{im}+\sum _{j=1,j\ne i}^n(m_{ij}{S}_{jm}-m_{ji}{S}_{im}),\\ \frac{\mathrm{d}{I}_{im}}{\mathrm{d}t}& ={\beta }_i\frac{S_{im}{I}_{if}}{N_{if}}-({\mu }_{im}+{\gamma }_{im})I_{im}+\sum _{j=1,j\ne i}^n(m_{ij}{I}_{jm}-m_{ji}{I}_{im}),\\ \frac{\mathrm{d}{A}_{im}}{\mathrm{d}t}& ={\gamma }_{im}{I}_{im}-({\mu }_{im}+d_{im})A_{im}+\sum _{j=1,j\ne i}^n(m_{ij}{A}_{jm}-m_{ji}{A}_{im}),\end{array}$(1) ) have positive and significant impact on the basic reproduction number $R_0$, which indicate that decreasing the transmission rates of sexual ${\beta }_1$ and ${\beta }_2$ in two patches is two effective measures to reduce the basic reproduction number $R_0$. The parameters ${\gamma }_{1f},{\gamma }_{1m},{\gamma }_{2f},{\gamma }_{2m},{\mu }_{1f},{\mu }_{1m},{\mu }_{2f},{\mu }_{2m}$ have negative impacts on the basic reproduction number $R_0$, of which ${\gamma }_{1f},{\gamma }_{1m},{\gamma }_{2f},{\gamma }_{2m}$ have more significant negative impacts on $R_0$. It means that if these transfer rates can be strictly controlled, $R_0$ will decrease, that is, if the number of infected individuals can be controlled, the disease will gradually disappear. We have shown that the impact of migration rates $m_{21}$ and $m_{12}$ on the basic reproduction number $R_0$ is always one positive and one negative. In Figure 16, when $R_0^{(1)}<1$, $m_{21}$ has a positive impact and $m_{12}$ has a negative impact. This is also shown in Figures 8 and 11 of the contour plots. However, in Figure 17, when $R_0^{(1)}>1$, $m_{21}$ has a negative impact and $m_{12}$ has a positive impact. This is also shown in Figures 9 and 10 of the contour plots. Therefore, we should effectively control the migration rates according to the actual situation of the two patches, so as to reduce $R_0$ and make the disease disappear gradually. From the perspective of the research process, the most effective measure to control and prevent HIV/AIDS spread is to reduce contact with infected individuals, thereby reducing the parameter values of ${\beta }_1$ and ${\beta }_2$. From the perspective of control, $R_0$ can be reduced by increasing publicity, education and using condoms, etc. In Figure 18, it is shown that decreasing transmission rates can reduce the spread of HIV/AIDS in two patches. Hence, to control the spread of HIV/AIDS, the intervention measure must be taken.

7. Conclusion

Cosner et al. [12] and Sattenspiel et al. [27] introduced two models of human movement, including Euler movement model and Lagrange movement model. Citron et al. [11] further introduced the two models and examples were given to illustrate how to combine the two models with the research problems. Therefore, Euler movement and Lagrange movement are usually considered for individual movement.

Figure 16. Value of correlation coefficient for outcome $R_0$ when $R_0^{(1)}<1.$

Figure 17. Value of correlation coefficient for outcome $R_0$ when $R_0^{(1)}>1.$

Figure 18. Influence of different transmission rate parameters on HIV/AIDS transmission in two patches with $m_{21}=0.01,m_{12}=0.02$. Other parameters are shown in Table 2: (a) infective individuals in patch 1, (b) AIDS individuals in patch 1, (c) infective individuals in patch 2 and (d) AIDS individuals in patch 2.

Table 3. The PRCC of the parameters in system (1(1) $\begin{array}{rl}\frac{\mathrm{d}{S}_{if}}{\mathrm{d}t}& ={\epsilon }_i{\mathrm{\Lambda }}_i-{\beta }_i\frac{S_{if}{I}_{im}}{N_{im}}-{\mu }_{if}{S}_{if}+\sum _{j=1,j\ne i}^n(m_{ij}{S}_{jf}-m_{ji}{S}_{if}),\\ \frac{\mathrm{d}{I}_{if}}{\mathrm{d}t}& ={\beta }_i\frac{S_{if}{I}_{im}}{N_{im}}-({\mu }_{if}+{\gamma }_{if})I_{if}+\sum _{j=1,j\ne i}^n(m_{ij}{I}_{jf}-m_{ji}{I}_{if}),\\ \frac{\mathrm{d}{A}_{if}}{\mathrm{d}t}& ={\gamma }_{if}{I}_{if}-({\mu }_{if}+d_{if})A_{if}+\sum _{j=1,j\ne i}^n(m_{ij}{A}_{jf}-m_{ji}{A}_{if}),\\ \frac{\mathrm{d}{S}_{im}}{\mathrm{d}t}& =(1-{\epsilon }_i){\mathrm{\Lambda }}_i-{\beta }_i\frac{S_{im}{I}_{if}}{N_{if}}-{\mu }_{im}{S}_{im}+\sum _{j=1,j\ne i}^n(m_{ij}{S}_{jm}-m_{ji}{S}_{im}),\\ \frac{\mathrm{d}{I}_{im}}{\mathrm{d}t}& ={\beta }_i\frac{S_{im}{I}_{if}}{N_{if}}-({\mu }_{im}+{\gamma }_{im})I_{im}+\sum _{j=1,j\ne i}^n(m_{ij}{I}_{jm}-m_{ji}{I}_{im}),\\ \frac{\mathrm{d}{A}_{im}}{\mathrm{d}t}& ={\gamma }_{im}{I}_{im}-({\mu }_{im}+d_{im})A_{im}+\sum _{j=1,j\ne i}^n(m_{ij}{A}_{jm}-m_{ji}{A}_{im}),\end{array}$(1) ) when $R_0^{(1)}<1$.

Parameters	p values	PRCC	Parameters	p values	PRCC
${\beta }_1$	0	0.8710	${\beta }_2$	0	0.9795
${\gamma }_{1f}$	0	$-0.6849$	${\gamma }_{1m}$	0	$-0.6967$
${\gamma }_{2f}$	0	$-0.9048$	${\gamma }_{2m}$	0	$-0.9126$
$m_{21}$	$1.0\ast {10}^{-12}\ast 0.3353$	0.2287	$m_{12}$	$1.0\ast {10}^{-13}\ast 0.2481$	$-0.2392$
${\mu }_{1f}$	0.2829	$-0.0342$	${\mu }_{1m}$	0.0698	$-0.0577$
${\mu }_{2f}$	0.7015	$-0.0122$	${\mu }_{2m}$	0.0374	$-0.0662$

Table 4. The PRCC of the parameters in system (1(1) $\begin{array}{rl}\frac{\mathrm{d}{S}_{if}}{\mathrm{d}t}& ={\epsilon }_i{\mathrm{\Lambda }}_i-{\beta }_i\frac{S_{if}{I}_{im}}{N_{im}}-{\mu }_{if}{S}_{if}+\sum _{j=1,j\ne i}^n(m_{ij}{S}_{jf}-m_{ji}{S}_{if}),\\ \frac{\mathrm{d}{I}_{if}}{\mathrm{d}t}& ={\beta }_i\frac{S_{if}{I}_{im}}{N_{im}}-({\mu }_{if}+{\gamma }_{if})I_{if}+\sum _{j=1,j\ne i}^n(m_{ij}{I}_{jf}-m_{ji}{I}_{if}),\\ \frac{\mathrm{d}{A}_{if}}{\mathrm{d}t}& ={\gamma }_{if}{I}_{if}-({\mu }_{if}+d_{if})A_{if}+\sum _{j=1,j\ne i}^n(m_{ij}{A}_{jf}-m_{ji}{A}_{if}),\\ \frac{\mathrm{d}{S}_{im}}{\mathrm{d}t}& =(1-{\epsilon }_i){\mathrm{\Lambda }}_i-{\beta }_i\frac{S_{im}{I}_{if}}{N_{if}}-{\mu }_{im}{S}_{im}+\sum _{j=1,j\ne i}^n(m_{ij}{S}_{jm}-m_{ji}{S}_{im}),\\ \frac{\mathrm{d}{I}_{im}}{\mathrm{d}t}& ={\beta }_i\frac{S_{im}{I}_{if}}{N_{if}}-({\mu }_{im}+{\gamma }_{im})I_{im}+\sum _{j=1,j\ne i}^n(m_{ij}{I}_{jm}-m_{ji}{I}_{im}),\\ \frac{\mathrm{d}{A}_{im}}{\mathrm{d}t}& ={\gamma }_{im}{I}_{im}-({\mu }_{im}+d_{im})A_{im}+\sum _{j=1,j\ne i}^n(m_{ij}{A}_{jm}-m_{ji}{A}_{im}),\end{array}$(1) ) when $R_0^{(1)}>1$.

Parameters	p values	PRCC	Parameters	p values	PRCC
${\beta }_1$	0	0.9664	${\beta }_2$	0	0.9276
${\gamma }_{1f}$	0	$-0.8467$	${\gamma }_{1m}$	0	$-0.8491$
${\gamma }_{2f}$	0	$-0.7972$	${\gamma }_{2m}$	0	$-0.8115$
$m_{21}$	$1.0\ast {10}^{-21}\ast 0.2110$	$-0.3027$	$m_{12}$	$1.0\ast {10}^{-14}\ast 0.7444$	0.2439
${\mu }_{1f}$	0.0250	$-0.0713$	${\mu }_{1m}$	$1.0\ast {10}^{-3}\ast 0.6598$	$-0.1081$
${\mu }_{2f}$	0.3532	$-0.0296$	${\mu }_{2m}$	0.1569	$-0.0450$

In this paper, we considered the movement of individuals as the form of Euler movement. And we have proposed an HIV/AIDS epidemic model with heterosexual transmission in a patchy environment to study the dynamics of disease transmission under the influence of population migration among patches. A patch can be interpreted as a city, a province or a country, etc. We have shown that $R_0<1$ implies that the disease-free equilibrium is unique and globally asymptotically stable and $R_0>1$ implies that the disease is uniformly persistent if migration matrix $M=(m_{ij})$ is irreducible. If $M=(m_{ij})$ is irreducible and there exist ${\eta }_{j1},{\eta }_{j2}>0$ such that $S_{jf}^{\ast }={\eta }_{j1}{I}_{jf}^{\ast },S_{jm}^{\ast }={\eta }_{j2}{I}_{jm}^{\ast }$ for all $1\le i,j\le n$, the endemic equilibrium $E^{\ast }$ is unique and globally asymptotically stable. We have selected two patches for numerical simulations. We have chosen the data from the Beijing City and Yunnan Province, Beijing City is regarded as the first patch, Yunnan Province is regarded as the second patch. Population migration may increase the spread of disease, turning a low-risk patch into a high-risk patch, or it may reduce the spread of disease, turning a high-risk patch into a low-risk patch. The effect of population migration on disease transmission is related to the status of the two patches in isolation. If HIV/AIDS becomes extinct in each patch when two patches are isolated, the disease remains extinct in two patches when the population migration occurs; if HIV/AIDS spreads in each patch when two patches are isolated, the disease remains persistent in two patches when the population migration occurs; if the disease disappears in one patch and spreads in the other patch when they are isolated, the disease can spread or disappear in two patches if migration rates of individuals are suitably chosen. We have shown that the impacts of migration rates $m_{21}$ and $m_{12}$ on the basic reproduction number $R_0$ is always one positive and one negative. When $R_0^{(1)}<1$, $m_{21}$ has a positive impact and $m_{12}$ has a negative impact; when $R_0^{(1)}>1$, $m_{21}$ has a negative impact and $m_{12}$ has a positive impact. When two patches are low-risk in isolation, we hardly need to control the migration rates. When one patch is low-risk and the other patch is high-risk in isolation, we can control the migration rates to make $R_0<1$ according to Figures 10–11, so that two patches are low-risk. Otherwise, if the migration rates are not properly controlled, the low-risk patch may become high-risk patch. After that, no matter how the migration rates are controlled, it is difficult to reverse the situation, and the disease will spread in two patches. Therefore, we should effectively control the migration rates according to the actual situation of the two patches, so as to reduce $R_0$ and make the disease reduce gradually. In addition, Figure 18 has shown that reducing contact with infected individuals, i.e. reducing the parameter values of ${\beta }_1$ and ${\beta }_2$ is also a key measure to reduce the spread of the disease. The model established in this paper is helpful to understand the spread of HIV/AIDS among regions when population migration occurs. In addition, the model can also be used in other regions. These studies enrich the research on patches.

Acknowledgments

The authors would like to thank the referees for their helpful suggestions that improved the presentation of this manuscript.

Disclosure statement

No potential conflict of interest was reported by the author(s).

Additional information

Funding

This work is supported by the National Sciences Foundation of China (No. 11971278), Shanxi Provincial Research and Development Program (No. 202102130501002) and Health Commission of Shanxi Province (2020XM18).

Appendix

Appendix 1.

Proof of Theorem 3.4

Proof.

By Theorem 3.3, the existence of an endemic equilibrium $E^{\ast }$ is ensured. We prove that $E^{\ast }$ is globally asymptotically stable in int(X).

First, we prove that if there exist ${\eta }_{j1},{\eta }_{j2}>0$ such that $S_{jf}^{\ast }={\eta }_{j1}{I}_{jf}^{\ast },S_{jm}^{\ast }={\eta }_{j2}{I}_{jm}^{\ast }$ for $1\le j\le n$, $E^{\ast }$ is the only endemic equilibrium. We consider the following system: (A1) $\begin{array}{rl}& {\epsilon }_i{\mathrm{\Lambda }}_i-{\beta }_i\frac{S_{if}{I}_{im}}{N_{im}}-{\mu }_{if}{S}_{if}+\sum _{j=1,j\ne i}^n(m_{ij}{S}_{jf}-m_{ji}{S}_{if})=0,\\ & {\beta }_i\frac{S_{if}{I}_{im}}{N_{im}}-({\mu }_{if}+{\gamma }_{if})I_{if}+\sum _{j=1,j\ne i}^n(m_{ij}{I}_{jf}-m_{ji}{I}_{if})=0,\\ & {\gamma }_{if}{I}_{if}-({\mu }_{if}+d_{if})A_{if}+\sum _{j=1,j\ne i}^n(m_{ij}{A}_{jf}-m_{ji}{A}_{if})=0,\\ & (1-{\epsilon }_i){\mathrm{\Lambda }}_i-{\beta }_i\frac{S_{im}{I}_{if}}{N_{if}}-{\mu }_{im}{S}_{im}+\sum _{j=1,j\ne i}^n(m_{ij}{S}_{jm}-m_{ji}{S}_{im})=0,\\ & {\beta }_i\frac{S_{im}{I}_{if}}{N_{if}}-({\mu }_{im}+{\gamma }_{im})I_{im}+\sum _{j=1,j\ne i}^n(m_{ij}{I}_{jm}-m_{ji}{I}_{im})=0,\\ & {\gamma }_{im}{I}_{im}-({\mu }_{im}+d_{im})A_{im}+\sum _{j=1,j\ne i}^n(m_{ij}{A}_{jm}-m_{ji}{A}_{im})=0.\end{array}$(A1) From the first two equations of (A1(A1) $\begin{array}{rl}& {\epsilon }_i{\mathrm{\Lambda }}_i-{\beta }_i\frac{S_{if}{I}_{im}}{N_{im}}-{\mu }_{if}{S}_{if}+\sum _{j=1,j\ne i}^n(m_{ij}{S}_{jf}-m_{ji}{S}_{if})=0,\\ & {\beta }_i\frac{S_{if}{I}_{im}}{N_{im}}-({\mu }_{if}+{\gamma }_{if})I_{if}+\sum _{j=1,j\ne i}^n(m_{ij}{I}_{jf}-m_{ji}{I}_{if})=0,\\ & {\gamma }_{if}{I}_{if}-({\mu }_{if}+d_{if})A_{if}+\sum _{j=1,j\ne i}^n(m_{ij}{A}_{jf}-m_{ji}{A}_{if})=0,\\ & (1-{\epsilon }_i){\mathrm{\Lambda }}_i-{\beta }_i\frac{S_{im}{I}_{if}}{N_{if}}-{\mu }_{im}{S}_{im}+\sum _{j=1,j\ne i}^n(m_{ij}{S}_{jm}-m_{ji}{S}_{im})=0,\\ & {\beta }_i\frac{S_{im}{I}_{if}}{N_{if}}-({\mu }_{im}+{\gamma }_{im})I_{im}+\sum _{j=1,j\ne i}^n(m_{ij}{I}_{jm}-m_{ji}{I}_{im})=0,\\ & {\gamma }_{im}{I}_{im}-({\mu }_{im}+d_{im})A_{im}+\sum _{j=1,j\ne i}^n(m_{ij}{A}_{jm}-m_{ji}{A}_{im})=0.\end{array}$(A1) ), we have $\begin{array}{rl}& {\epsilon }_i{\mathrm{\Lambda }}_i-({\mu }_{if}+{\gamma }_{if}+\sum _{j=1,j\ne i}^n{m}_{ji})I_{if}+\sum _{j=1,j\ne i}^n{m}_{ij}{I}_{jf}-({\mu }_{if}+\sum _{j=1,j\ne i}^n{m}_{ji})S_{if}+\sum _{j=1,j\ne i}^n{m}_{ij}{S}_{jf}=0,\\ & i=1,2,\dots ,n,\end{array}$or in the form of matrix $A{S}_f=\epsilon \mathrm{\Lambda }-(A+{\gamma }_f)I_f,$where $S_f=(S_{1f},S_{2f},\dots ,S_{nf}{)}^{\mathrm{T}},I_f=(I_{1f},I_{2f},\dots ,I_{nf}{)}^{\mathrm{T}},\epsilon =\mathrm{diag}({\epsilon }_1,{\epsilon }_2,\dots ,{\epsilon }_n),\mathrm{\Lambda }=({\mathrm{\Lambda }}_1,{\mathrm{\Lambda }}_2,\dots$, ${\mathrm{\Lambda }}_n{)}^{\mathrm{T}},$ and ${\gamma }_f=\mathrm{diag}({\gamma }_{1f},{\gamma }_{2f},\dots ,{\gamma }_{nf})$. Therefore, we have $A{S}_f^{\ast }=\epsilon \mathrm{\Lambda }-(A+{\gamma }_f)I_f^{\ast }$, where $S_f^{\ast }=(S_{1f}^{\ast },S_{2f}^{\ast },\dots ,S_{nf}^{\ast }{)}^{\mathrm{T}},I_f^{\ast }=(I_{1f}^{\ast },I_{2f}^{\ast },\dots ,I_{nf}^{\ast }{)}^{\mathrm{T}}$. Since $S_{jf}^{\ast }={\eta }_{j1}{I}_{jf}^{\ast },j=1,2,\dots ,n$, hence $\epsilon \mathrm{\Lambda }-(A+{\gamma }_f)I_f^{\ast }=AU{I}_f^{\ast }$, where $U=\mathrm{diag}({\eta }_{11},{\eta }_{21},\dots ,{\eta }_{n1})$. Since all off-diagonal entries of $A+{\gamma }_f+AU$ is non-positive and the sum of the entries in each column of $A+{\gamma }_f+AU$ is positive, so $A+{\gamma }_f+AU$ is a nonsingular M-matrix and $(A+{\gamma }_f+AU{)}^{-1}\ge 0$, hence $I_f^{\ast }=(I_{1f}^{\ast },I_{2f}^{\ast },\dots ,I_{nf}^{\ast }{)}^{\mathrm{T}}=(A+{\gamma }_f+AU{)}^{-1}\epsilon \mathrm{\Lambda }>0$. In the same way, we have $B{S}_m=(1-\epsilon )\mathrm{\Lambda }-(B+{\gamma }_m)I_m$, where $S_m=(S_{1m},S_{2m},\dots ,S_{nm}{)}^{\mathrm{T}},I_m=(I_{1m},I_{2m},\dots ,I_{nm}{)}^{\mathrm{T}}$ and ${\gamma }_m=\mathrm{diag}({\gamma }_{1m},{\gamma }_{2m},\dots ,{\gamma }_{nm})$. Therefore, we have $B{S}_m^{\ast }=(1-\epsilon )\mathrm{\Lambda }-(B+{\gamma }_m)I_m^{\ast }$, where $S_m^{\ast }=(S_{1m}^{\ast },S_{2m}^{\ast },\dots ,S_{nm}^{\ast }{)}^{\mathrm{T}},I_m^{\ast }=(I_{1m}^{\ast },I_{2m}^{\ast },\dots ,I_{nm}^{\ast }{)}^{\mathrm{T}}$. Since $S_{jm}^{\ast }={\eta }_{j2}{I}_{jm}^{\ast },j=1,2,\dots ,n$, hence $(1-\epsilon )\mathrm{\Lambda }-(B+{\gamma }_m)I_m^{\ast }=BW{I}_m^{\ast }$, where $W=\mathrm{diag}({\eta }_{12}{\eta }_{22},\dots ,{\eta }_{n2})$. We obtain $I_m^{\ast }=(I_{1m}^{\ast },I_{2m}^{\ast },\dots ,I_{nm}^{\ast }{)}^{\mathrm{T}}=(B+{\gamma }_m+BW{)}^{-1}(1-\epsilon )\mathrm{\Lambda }>0$.

We write the third and sixth equations of (A1(A1) $\begin{array}{rl}& {\epsilon }_i{\mathrm{\Lambda }}_i-{\beta }_i\frac{S_{if}{I}_{im}}{N_{im}}-{\mu }_{if}{S}_{if}+\sum _{j=1,j\ne i}^n(m_{ij}{S}_{jf}-m_{ji}{S}_{if})=0,\\ & {\beta }_i\frac{S_{if}{I}_{im}}{N_{im}}-({\mu }_{if}+{\gamma }_{if})I_{if}+\sum _{j=1,j\ne i}^n(m_{ij}{I}_{jf}-m_{ji}{I}_{if})=0,\\ & {\gamma }_{if}{I}_{if}-({\mu }_{if}+d_{if})A_{if}+\sum _{j=1,j\ne i}^n(m_{ij}{A}_{jf}-m_{ji}{A}_{if})=0,\\ & (1-{\epsilon }_i){\mathrm{\Lambda }}_i-{\beta }_i\frac{S_{im}{I}_{if}}{N_{if}}-{\mu }_{im}{S}_{im}+\sum _{j=1,j\ne i}^n(m_{ij}{S}_{jm}-m_{ji}{S}_{im})=0,\\ & {\beta }_i\frac{S_{im}{I}_{if}}{N_{if}}-({\mu }_{im}+{\gamma }_{im})I_{im}+\sum _{j=1,j\ne i}^n(m_{ij}{I}_{jm}-m_{ji}{I}_{im})=0,\\ & {\gamma }_{im}{I}_{im}-({\mu }_{im}+d_{im})A_{im}+\sum _{j=1,j\ne i}^n(m_{ij}{A}_{jm}-m_{ji}{A}_{im})=0.\end{array}$(A1) ) in matrix form respectively, i.e. $(A+D_f)A_f={\gamma }_f{I}_f$ and $(B+D_m)A_m={\gamma }_m{I}_m$, where $D_f=\mathrm{diag}(d_{1f},d_{2f},\dots ,d_{nf}),D_m=\mathrm{diag}(d_{1m},d_{2m},\dots ,d_{nm})$. Since all off-diagonal entries of $A+D_f$ and $B+D_m$ are non-positive and the sum of the entries in each column of $A+D_f$ and $B+D_m$ is positive, so they are two nonsingular M-matrices, $(A+D_f{)}^{-1}\ge 0,(B+D_m{)}^{-1}\ge 0$. Therefore, linear system $(A+D_f)A_f={\gamma }_f{I}_f$ has a unique positive solution $A_f^{\ast }=(A_{1f}^{\ast },A_{2f}^{\ast },\dots ,A_{nf}^{\ast }{)}^{\mathrm{T}}={\gamma }_f(A+D_f{)}^{-1}{I}_f^{\ast }>0$, linear system $(B+D_m)A_m={\gamma }_m{I}_m$ has a unique positive solution $A_m^{\ast }=(A_{1m}^{\ast },A_{2m}^{\ast },\dots ,A_{nm}^{\ast }{)}^{\mathrm{T}}={\gamma }_m(B+D_m{)}^{-1}{I}_m^{\ast }>0$, $S_f^{\ast }=(S_{1f}^{\ast },S_{2f}^{\ast },\dots ,S_{nf}^{\ast }{)}^{\mathrm{T}}=A^{-1}\epsilon \mathrm{\Lambda }-A^{-1}(A+{\gamma }_f)I_f^{\ast }=U{I}_f^{\ast }>0$, and $S_m^{\ast }=(S_{1m}^{\ast },S_{2m}^{\ast },\dots ,S_{nm}^{\ast }{)}^{\mathrm{T}}=B^{-1}(1-\epsilon )\mathrm{\Lambda }-B^{-1}(B+{\gamma }_m)I_m^{\ast }=W{I}_m^{\ast }>0$.

Then we prove that the endemic equilibrium $E^{\ast }$ of system (1(1) $\begin{array}{rl}\frac{\mathrm{d}{S}_{if}}{\mathrm{d}t}& ={\epsilon }_i{\mathrm{\Lambda }}_i-{\beta }_i\frac{S_{if}{I}_{im}}{N_{im}}-{\mu }_{if}{S}_{if}+\sum _{j=1,j\ne i}^n(m_{ij}{S}_{jf}-m_{ji}{S}_{if}),\\ \frac{\mathrm{d}{I}_{if}}{\mathrm{d}t}& ={\beta }_i\frac{S_{if}{I}_{im}}{N_{im}}-({\mu }_{if}+{\gamma }_{if})I_{if}+\sum _{j=1,j\ne i}^n(m_{ij}{I}_{jf}-m_{ji}{I}_{if}),\\ \frac{\mathrm{d}{A}_{if}}{\mathrm{d}t}& ={\gamma }_{if}{I}_{if}-({\mu }_{if}+d_{if})A_{if}+\sum _{j=1,j\ne i}^n(m_{ij}{A}_{jf}-m_{ji}{A}_{if}),\\ \frac{\mathrm{d}{S}_{im}}{\mathrm{d}t}& =(1-{\epsilon }_i){\mathrm{\Lambda }}_i-{\beta }_i\frac{S_{im}{I}_{if}}{N_{if}}-{\mu }_{im}{S}_{im}+\sum _{j=1,j\ne i}^n(m_{ij}{S}_{jm}-m_{ji}{S}_{im}),\\ \frac{\mathrm{d}{I}_{im}}{\mathrm{d}t}& ={\beta }_i\frac{S_{im}{I}_{if}}{N_{if}}-({\mu }_{im}+{\gamma }_{im})I_{im}+\sum _{j=1,j\ne i}^n(m_{ij}{I}_{jm}-m_{ji}{I}_{im}),\\ \frac{\mathrm{d}{A}_{im}}{\mathrm{d}t}& ={\gamma }_{im}{I}_{im}-({\mu }_{im}+d_{im})A_{im}+\sum _{j=1,j\ne i}^n(m_{ij}{A}_{jm}-m_{ji}{A}_{im}),\end{array}$(1) ) is globally asymptotically stable in int(X) when $R_0>1$.

Define $V_i:\{(S_{1f},S_{1m},I_{1f},I_{1m},\dots ,S_{nf},S_{nm},I_{nf},I_{nm})\in X|S_{if}>0,S_{im}>0,I_{if}>0,I_{im}>0i=1,2,\dots ,n\}⟶R$ by $\begin{array}{rl}{V}_i(t)& =S_{if}-S_{if}^{\ast }-S_{if}^{\ast }\ln \frac{S_{if}}{S_{if}^{\ast }}+\frac{S_{if}^{\ast }{I}_{im}^{\ast }{N}_{if}^{\ast }}{S_{im}^{\ast }{I}_{if}^{\ast }{N}_{im}^{\ast }}(S_{im}-S_{im}^{\ast }-S_{im}^{\ast }\ln \frac{S_{im}}{S_{im}^{\ast }})+I_{if}-I_{if}^{\ast }\\ & -I_{if}^{\ast }\ln \frac{I_{if}}{I_{if}^{\ast }}+\frac{S_{if}^{\ast }{I}_{im}^{\ast }{N}_{if}^{\ast }}{S_{im}^{\ast }{I}_{if}^{\ast }{N}_{im}^{\ast }}(I_{im}-I_{im}^{\ast }-I_{im}^{\ast }\ln \frac{I_{im}}{I_{im}^{\ast }}).\end{array}$From equilibrium equations of (1(1) $\begin{array}{rl}\frac{\mathrm{d}{S}_{if}}{\mathrm{d}t}& ={\epsilon }_i{\mathrm{\Lambda }}_i-{\beta }_i\frac{S_{if}{I}_{im}}{N_{im}}-{\mu }_{if}{S}_{if}+\sum _{j=1,j\ne i}^n(m_{ij}{S}_{jf}-m_{ji}{S}_{if}),\\ \frac{\mathrm{d}{I}_{if}}{\mathrm{d}t}& ={\beta }_i\frac{S_{if}{I}_{im}}{N_{im}}-({\mu }_{if}+{\gamma }_{if})I_{if}+\sum _{j=1,j\ne i}^n(m_{ij}{I}_{jf}-m_{ji}{I}_{if}),\\ \frac{\mathrm{d}{A}_{if}}{\mathrm{d}t}& ={\gamma }_{if}{I}_{if}-({\mu }_{if}+d_{if})A_{if}+\sum _{j=1,j\ne i}^n(m_{ij}{A}_{jf}-m_{ji}{A}_{if}),\\ \frac{\mathrm{d}{S}_{im}}{\mathrm{d}t}& =(1-{\epsilon }_i){\mathrm{\Lambda }}_i-{\beta }_i\frac{S_{im}{I}_{if}}{N_{if}}-{\mu }_{im}{S}_{im}+\sum _{j=1,j\ne i}^n(m_{ij}{S}_{jm}-m_{ji}{S}_{im}),\\ \frac{\mathrm{d}{I}_{im}}{\mathrm{d}t}& ={\beta }_i\frac{S_{im}{I}_{if}}{N_{if}}-({\mu }_{im}+{\gamma }_{im})I_{im}+\sum _{j=1,j\ne i}^n(m_{ij}{I}_{jm}-m_{ji}{I}_{im}),\\ \frac{\mathrm{d}{A}_{im}}{\mathrm{d}t}& ={\gamma }_{im}{I}_{im}-({\mu }_{im}+d_{im})A_{im}+\sum _{j=1,j\ne i}^n(m_{ij}{A}_{jm}-m_{ji}{A}_{im}),\end{array}$(1) ), we obtain $\begin{array}{rl}{\mu }_{if}{S}_{if}^{\ast }& ={\epsilon }_i{\mathrm{\Lambda }}_i-{\beta }_i\frac{S_{if}^{\ast }{I}_{im}^{\ast }}{N_{im}^{\ast }}+\sum _{j=1}^n(m_{ij}{S}_{jf}^{\ast }-m_{ji}{S}_{if}^{\ast }),\\ ({\mu }_{if}+{\gamma }_{if})I_{if}^{\ast }& ={\beta }_i\frac{S_{if}^{\ast }{I}_{im}^{\ast }}{N_{im}^{\ast }}+\sum _{j=1}^n(m_{ij}{I}_{jf}^{\ast }-m_{ji}{I}_{if}^{\ast }),\\ {\mu }_{im}{S}_{im}^{\ast }& =(1-{\epsilon }_i){\mathrm{\Lambda }}_i-{\beta }_i\frac{S_{im}^{\ast }{I}_{if}^{\ast }}{N_{if}^{\ast }}+\sum _{j=1}^n(m_{ij}{S}_{jm}^{\ast }-m_{ji}{S}_{im}^{\ast }),\\ ({\mu }_{im}+{\gamma }_{im})I_{im}^{\ast }& ={\beta }_i\frac{S_{im}^{\ast }{I}_{if}^{\ast }}{N_{if}^{\ast }}+\sum _{j=1}^n(m_{ij}{I}_{jm}^{\ast }-m_{ji}{I}_{im}^{\ast }).\end{array}$Note that $1-x+\ln x\le 0$ for x>0 and equality holds if and only if x = 1. Calculating the time derivative of $V_i(t)$ along solutions of system (1(1) $\begin{array}{rl}\frac{\mathrm{d}{S}_{if}}{\mathrm{d}t}& ={\epsilon }_i{\mathrm{\Lambda }}_i-{\beta }_i\frac{S_{if}{I}_{im}}{N_{im}}-{\mu }_{if}{S}_{if}+\sum _{j=1,j\ne i}^n(m_{ij}{S}_{jf}-m_{ji}{S}_{if}),\\ \frac{\mathrm{d}{I}_{if}}{\mathrm{d}t}& ={\beta }_i\frac{S_{if}{I}_{im}}{N_{im}}-({\mu }_{if}+{\gamma }_{if})I_{if}+\sum _{j=1,j\ne i}^n(m_{ij}{I}_{jf}-m_{ji}{I}_{if}),\\ \frac{\mathrm{d}{A}_{if}}{\mathrm{d}t}& ={\gamma }_{if}{I}_{if}-({\mu }_{if}+d_{if})A_{if}+\sum _{j=1,j\ne i}^n(m_{ij}{A}_{jf}-m_{ji}{A}_{if}),\\ \frac{\mathrm{d}{S}_{im}}{\mathrm{d}t}& =(1-{\epsilon }_i){\mathrm{\Lambda }}_i-{\beta }_i\frac{S_{im}{I}_{if}}{N_{if}}-{\mu }_{im}{S}_{im}+\sum _{j=1,j\ne i}^n(m_{ij}{S}_{jm}-m_{ji}{S}_{im}),\\ \frac{\mathrm{d}{I}_{im}}{\mathrm{d}t}& ={\beta }_i\frac{S_{im}{I}_{if}}{N_{if}}-({\mu }_{im}+{\gamma }_{im})I_{im}+\sum _{j=1,j\ne i}^n(m_{ij}{I}_{jm}-m_{ji}{I}_{im}),\\ \frac{\mathrm{d}{A}_{im}}{\mathrm{d}t}& ={\gamma }_{im}{I}_{im}-({\mu }_{im}+d_{im})A_{im}+\sum _{j=1,j\ne i}^n(m_{ij}{A}_{jm}-m_{ji}{A}_{im}),\end{array}$(1) ), we have $\begin{array}{rl}{V}_i^{\prime }(t)& =(1-\frac{S_{if}^{\ast }}{S_{if}})[{\epsilon }_i{\mathrm{\Lambda }}_i+\sum _{j=1}^n(m_{ij}{S}_{jf}-m_{ji}{S}_{if})\\ & -{\beta }_i\frac{S_{if}{I}_{im}}{N_{im}}-\frac{{\epsilon }_i{\mathrm{\Lambda }}_i-{\beta }_i\frac{S_{if}^{\ast }{I}_{im}^{\ast }}{N_{im}^{\ast }}+\sum _{j=1}^n(m_{ij}{S}_{jf}^{\ast }-m_{ji}{S}_{if}^{\ast })}{S_{if}^{\ast }}{S}_{if}]\\ & +\frac{S_{if}^{\ast }{I}_{im}^{\ast }{N}_{if}^{\ast }}{S_{im}^{\ast }{I}_{if}^{\ast }{N}_{im}^{\ast }}(1-\frac{S_{im}^{\ast }}{S_{im}})[(1-{\epsilon }_i){\mathrm{\Lambda }}_i+\sum _{j=1}^n(m_{ij}{S}_{jm}-m_{ji}{S}_{im})-{\beta }_i\frac{S_{im}{I}_{if}}{N_{if}}\\ & -\frac{(1-{\epsilon }_i){\mathrm{\Lambda }}_i-{\beta }_i\frac{S_{im}^{\ast }{I}_{if}^{\ast }}{N_{if}^{\ast }}+\sum _{j=1}^n(m_{ij}{S}_{jm}^{\ast }-m_{ji}{S}_{im}^{\ast })}{S_{im}^{\ast }}{S}_{im}]\\ & +(1-\frac{I_{if}^{\ast }}{I_{if}})[{\beta }_i\frac{S_{if}{I}_{im}}{N_{im}}+\sum _{j=1}^n(m_{ij}{I}_{jf}-m_{ji}{I}_{if})-\frac{{\beta }_i\frac{S_{if}^{\ast }{I}_{im}^{\ast }}{N_{im}^{\ast }}+\sum _{j=1}^n(m_{ij}{I}_{jf}^{\ast }-m_{ji}{I}_{if}^{\ast })}{I_{if}^{\ast }}{I}_{if}]\\ & +\frac{S_{if}^{\ast }{I}_{im}^{\ast }{N}_{if}^{\ast }}{S_{im}^{\ast }{I}_{if}^{\ast }{N}_{im}^{\ast }}(1-\frac{I_{im}^{\ast }}{I_{im}})[{\beta }_i\frac{S_{im}{I}_{if}}{N_{if}}+\sum _{j=1}^n(m_{ij}{I}_{jm}-m_{ji}{I}_{im})\\ & -\frac{{\beta }_i\frac{S_{im}^{\ast }{I}_{if}^{\ast }}{N_{if}^{\ast }}+\sum _{j=1}^n(m_{ij}{I}_{jm}^{\ast }-m_{ji}{I}_{im}^{\ast })}{I_{im}^{\ast }}{I}_{im}]\\ & ={\epsilon }_i{\mathrm{\Lambda }}_i(2-\frac{S_{if}}{S_{if}^{\ast }}-\frac{S_{if}^{\ast }}{S_{if}})+(1-{\epsilon }_i){\mathrm{\Lambda }}_i(2-\frac{S_{im}}{S_{im}^{\ast }}-\frac{S_{im}^{\ast }}{S_{im}})\\ & +{\beta }_i\frac{S_{if}^{\ast }{I}_{im}^{\ast }}{N_{im}^{\ast }}(2-\frac{S_{if}^{\ast }}{S_{if}}-\frac{I_{if}}{I_{if}^{\ast }}+\frac{I_{im}{N}_{im}^{\ast }}{I_{im}^{\ast }{N}_{im}}-\frac{S_{if}{I}_{im}{N}_{im}^{\ast }{I}_{if}^{\ast }}{S_{if}^{\ast }{I}_{im}^{\ast }{N}_{im}{I}_{if}})\\ & +{\beta }_i\frac{S_{if}^{\ast }{I}_{im}^{\ast }}{N_{im}^{\ast }}(2-\frac{S_{im}^{\ast }}{S_{im}}-\frac{I_{im}}{I_{im}^{\ast }}+\frac{I_{if}{N}_{if}^{\ast }}{I_{if}^{\ast }{N}_{if}}-\frac{S_{im}{I}_{if}{N}_{if}^{\ast }{I}_{im}^{\ast }}{S_{im}^{\ast }{I}_{if}^{\ast }{N}_{if}{I}_{im}})\\ & +\sum _{j=1}^n{m}_{ij}{S}_{jf}^{\ast }(1+\frac{S_{jf}}{S_{jf}^{\ast }}-\frac{S_{if}}{S_{if}^{\ast }}-\frac{S_{if}^{\ast }{S}_{jf}}{S_{if}{S}_{jf}^{\ast }})+\sum _{j=1}^n\frac{S_{if}^{\ast }{I}_{im}^{\ast }{N}_{if}^{\ast }}{S_{im}^{\ast }{I}_{if}^{\ast }{N}_{im}^{\ast }}{m}_{ij}{S}_{jm}^{\ast }(1+\frac{S_{jm}}{S_{jm}^{\ast }}-\frac{S_{im}}{S_{im}^{\ast }}-\frac{S_{im}^{\ast }{S}_{jm}}{S_{im}{S}_{jm}^{\ast }})\\ & +\sum _{j=1}^n{m}_{ij}{I}_{jf}^{\ast }(1+\frac{I_{jf}}{I_{jf}^{\ast }}-\frac{I_{if}}{I_{if}^{\ast }}-\frac{I_{if}^{\ast }{I}_{jf}}{I_{if}{I}_{jf}^{\ast }})+\sum _{j=1}^n\frac{S_{if}^{\ast }{I}_{im}^{\ast }{N}_{if}^{\ast }}{S_{im}^{\ast }{I}_{if}^{\ast }{N}_{im}^{\ast }}{m}_{ij}{I}_{jm}^{\ast }(1+\frac{I_{jm}}{I_{jm}^{\ast }}-\frac{I_{im}}{I_{im}^{\ast }}-\frac{I_{im}^{\ast }{I}_{jm}}{I_{im}{I}_{jm}^{\ast }})\\ & \le {\epsilon }_i{\mathrm{\Lambda }}_i(2-\frac{S_{if}}{S_{if}^{\ast }}-\frac{S_{if}^{\ast }}{S_{if}})+\frac{S_{if}^{\ast }{I}_{im}^{\ast }{N}_{if}^{\ast }}{S_{im}^{\ast }{I}_{if}^{\ast }{N}_{im}^{\ast }}(1-{\epsilon }_i){\mathrm{\Lambda }}_i(2-\frac{S_{im}}{S_{im}^{\ast }}-\frac{S_{im}^{\ast }}{S_{im}})\\ & +{\beta }_i\frac{S_{if}^{\ast }{I}_{im}^{\ast }}{N_{im}^{\ast }}[(\frac{I_{im}{N}_{im}^{\ast }}{I_{im}^{\ast }{N}_{im}}-\frac{I_{if}}{I_{if}^{\ast }}-\ln \frac{S_{if}^{\ast }}{S_{if}}-\ln \frac{S_{if}{I}_{im}{N}_{im}^{\ast }{I}_{if}^{\ast }}{S_{if}^{\ast }{I}_{im}^{\ast }{N}_{im}{I}_{if}})\\ & +(\frac{I_{if}{N}_{if}^{\ast }}{I_{if}^{\ast }{N}_{if}}-\frac{I_{im}}{I_{im}^{\ast }}-\ln \frac{S_{im}^{\ast }}{S_{im}}-\ln \frac{S_{im}{I}_{if}{N}_{if}^{\ast }{I}_{im}^{\ast }}{S_{im}^{\ast }{I}_{if}^{\ast }{N}_{if}{I}_{im}})]\\ & +\sum _{j=1}^n{m}_{ij}{S}_{jf}^{\ast }(1-\frac{S_{if}^{\ast }{S}_{jf}}{S_{if}{S}_{jf}^{\ast }}+\ln \frac{S_{if}^{\ast }{S}_{jf}}{S_{if}{S}_{jf}^{\ast }})+\sum _{j=1}^n{m}_{ij}{S}_{jf}^{\ast }(\frac{S_{jf}}{S_{jf}^{\ast }}-\frac{S_{if}}{S_{if}^{\ast }}+\ln \frac{S_{if}}{S_{if}^{\ast }}-\ln \frac{S_{jf}}{S_{jf}^{\ast }})\\ & +\sum _{j=1}^n\frac{S_{if}^{\ast }{I}_{im}^{\ast }{N}_{if}^{\ast }}{S_{im}^{\ast }{I}_{if}^{\ast }{N}_{im}^{\ast }}{m}_{ij}{S}_{jm}^{\ast }(1-\frac{S_{im}^{\ast }{S}_{jm}}{S_{im}{S}_{jm}^{\ast }}+\ln \frac{S_{im}^{\ast }{S}_{jm}}{S_{im}{S}_{jm}^{\ast }})\\ & +\sum _{j=1}^n\frac{S_{if}^{\ast }{I}_{im}^{\ast }{N}_{if}^{\ast }}{S_{im}^{\ast }{I}_{if}^{\ast }{N}_{im}^{\ast }}{m}_{ij}{S}_{jm}^{\ast }(\frac{S_{jm}}{S_{jm}^{\ast }}-\frac{S_{im}}{S_{im}^{\ast }}+\ln \frac{S_{im}}{S_{im}^{\ast }}-\ln \frac{S_{jm}}{S_{jm}^{\ast }})\\ & +\sum _{j=1}^n{m}_{ij}{I}_{jf}^{\ast }(1-\frac{I_{if}^{\ast }{I}_{jf}}{I_{if}{I}_{jf}^{\ast }}+\ln \frac{I_{if}^{\ast }{I}_{jf}}{I_{if}{I}_{jf}^{\ast }})+\sum _{j=1}^n{m}_{ij}{I}_{jf}^{\ast }(\frac{I_{jf}}{I_{jf}^{\ast }}-\frac{I_{if}}{I_{if}^{\ast }}+\ln \frac{I_{if}}{I_{if}^{\ast }}-\ln \frac{I_{jf}}{I_{jf}^{\ast }})\\ & +\sum _{j=1}^n\frac{S_{if}^{\ast }{I}_{im}^{\ast }{N}_{if}^{\ast }}{S_{im}^{\ast }{I}_{if}^{\ast }{N}_{im}^{\ast }}{m}_{ij}{I}_{jm}^{\ast }(1-\frac{I_{im}^{\ast }{I}_{jm}}{I_{im}{I}_{jm}^{\ast }}+\ln \frac{I_{im}^{\ast }{I}_{jm}}{I_{im}{I}_{jm}^{\ast }})\\ & +\sum _{j=1}^n\frac{S_{if}^{\ast }{I}_{im}^{\ast }{N}_{if}^{\ast }}{S_{im}^{\ast }{I}_{if}^{\ast }{N}_{im}^{\ast }}{m}_{ij}{I}_{jm}^{\ast }(\frac{I_{jm}}{I_{jm}^{\ast }}-\frac{I_{im}}{I_{im}^{\ast }}+\ln \frac{I_{im}}{I_{im}^{\ast }}-\ln \frac{I_{jm}}{I_{jm}^{\ast }})\\ & \le {\beta }_i\frac{S_{if}^{\ast }{I}_{im}^{\ast }}{N_{im}^{\ast }}[(\frac{I_{im}{N}_{im}^{\ast }}{I_{im}^{\ast }{N}_{im}}-\ln \frac{I_{im}{N}_{im}^{\ast }}{I_{im}^{\ast }{N}_{im}}+\ln \frac{I_{if}}{I_{if}^{\ast }}-\frac{I_{if}}{I_{if}^{\ast }})+(\frac{I_{if}{N}_{if}^{\ast }}{I_{if}^{\ast }{N}_{if}}-\ln \frac{I_{if}{N}_{if}^{\ast }}{I_{if}^{\ast }{N}_{if}}+\ln \frac{I_{im}}{I_{im}^{\ast }}-\frac{I_{im}}{I_{im}^{\ast }})]\\ & +\sum _{j=1}^n{m}_{ij}{S}_{jf}^{\ast }(\frac{S_{jf}}{S_{jf}^{\ast }}-\frac{S_{if}}{S_{if}^{\ast }}+\ln \frac{S_{if}}{S_{if}^{\ast }}-\ln \frac{S_{jf}}{S_{jf}^{\ast }})\\ & +\sum _{j=1}^n\frac{S_{if}^{\ast }{I}_{im}^{\ast }{N}_{if}^{\ast }}{S_{im}^{\ast }{I}_{if}^{\ast }{N}_{im}^{\ast }}{m}_{ij}{S}_{jm}^{\ast }(\frac{S_{jm}}{S_{jm}^{\ast }}-\frac{S_{im}}{S_{im}^{\ast }}+\ln \frac{S_{im}}{S_{im}^{\ast }}-\ln \frac{S_{jm}}{S_{jm}^{\ast }})\\ & +\sum _{j=1}^n{m}_{ij}{I}_{jf}^{\ast }(\frac{I_{jf}}{I_{jf}^{\ast }}-\frac{I_{if}}{I_{if}^{\ast }}+\ln \frac{I_{if}}{I_{if}^{\ast }}-\ln \frac{I_{jf}}{I_{jf}^{\ast }})\\ & +\sum _{j=1}^n\frac{S_{if}^{\ast }{I}_{im}^{\ast }{N}_{if}^{\ast }}{S_{im}^{\ast }{I}_{if}^{\ast }{N}_{im}^{\ast }}{m}_{ij}{I}_{jm}^{\ast }(\frac{I_{jm}}{I_{jm}^{\ast }}-\frac{I_{im}}{I_{im}^{\ast }}+\ln \frac{I_{im}}{I_{im}^{\ast }}-\ln \frac{I_{jm}}{I_{jm}^{\ast }})\\ & ={\beta }_i\frac{S_{if}^{\ast }{I}_{im}^{\ast }}{N_{im}^{\ast }}[(\frac{I_{im}{N}_{im}^{\ast }}{I_{im}^{\ast }{N}_{im}}-\ln \frac{I_{im}{N}_{im}^{\ast }}{I_{im}^{\ast }{N}_{im}}+\ln \frac{I_{if}}{I_{if}^{\ast }}-\frac{I_{if}}{I_{if}^{\ast }})\\ & +(\frac{I_{if}{N}_{if}^{\ast }}{I_{if}^{\ast }{N}_{if}}-\ln \frac{I_{if}{N}_{if}^{\ast }}{I_{if}^{\ast }{N}_{if}}+\ln \frac{I_{im}}{I_{im}^{\ast }}-\frac{I_{im}}{I_{im}^{\ast }})]\\ & +\sum _{j=1}^n{m}_{ij}{I}_{jf}^{\ast }[({\eta }_{j1}\frac{S_{jf}}{S_{jf}^{\ast }}+{\eta }_{j1}\ln \frac{S_{if}}{S_{if}^{\ast }}-{\eta }_{j1}\frac{S_{if}}{S_{if}^{\ast }}-{\eta }_{j1}\ln \frac{S_{jf}}{S_{jf}^{\ast }})+(\frac{I_{jf}}{I_{jf}^{\ast }}-\frac{I_{if}}{I_{if}^{\ast }}+\ln \frac{I_{if}}{I_{if}^{\ast }}-\ln \frac{I_{jf}}{I_{jf}^{\ast }})]\\ & +\sum _{j=1}^n\frac{S_{if}^{\ast }{I}_{im}^{\ast }{N}_{if}^{\ast }}{S_{im}^{\ast }{I}_{if}^{\ast }{N}_{im}^{\ast }}{m}_{ij}{I}_{jm}^{\ast }[({\eta }_{j2}\frac{S_{jm}}{S_{jm}^{\ast }}+{\eta }_{j2}\ln \frac{S_{im}}{S_{im}^{\ast }}-{\eta }_{j2}\frac{S_{im}}{S_{im}^{\ast }}-{\eta }_{j2}\ln \frac{S_{jm}}{S_{jm}^{\ast }})\\ & +(\frac{I_{jm}}{I_{jm}^{\ast }}-\frac{I_{im}}{I_{im}^{\ast }}+\ln \frac{I_{im}}{I_{im}^{\ast }}-\ln \frac{I_{jm}}{I_{jm}^{\ast }})]\\ & ={\beta }_i\frac{S_{if}^{\ast }{I}_{im}^{\ast }}{N_{im}^{\ast }}[(\frac{I_{if}{N}_{if}^{\ast }}{I_{if}^{\ast }{N}_{if}}-\ln \frac{I_{if}{N}_{if}^{\ast }}{I_{if}^{\ast }{N}_{if}}+\frac{I_{im}{N}_{im}^{\ast }}{I_{im}^{\ast }{N}_{im}}-\ln \frac{I_{im}{N}_{im}^{\ast }}{I_{im}^{\ast }{N}_{im}})-(\frac{I_{if}}{I_{if}^{\ast }}-\ln \frac{I_{if}}{I_{if}^{\ast }}+\frac{I_{im}}{I_{im}^{\ast }}-\ln \frac{I_{im}}{I_{im}^{\ast }})]\\ & +\sum _{j=1}^n{m}_{ij}{I}_{jf}^{\ast }[({\eta }_{j1}\frac{S_{jf}}{S_{jf}^{\ast }}-{\eta }_{j1}\ln \frac{S_{jf}}{S_{jf}^{\ast }}+\frac{I_{jf}}{I_{jf}^{\ast }}-\ln \frac{I_{jf}}{I_{jf}^{\ast }})-({\eta }_{j1}\frac{S_{if}}{S_{if}^{\ast }}-{\eta }_{j1}\ln \frac{S_{if}}{S_{if}^{\ast }}+\frac{I_{if}}{I_{if}^{\ast }}-\ln \frac{I_{if}}{I_{if}^{\ast }})]\\ & +\sum _{j=1}^n\frac{S_{if}^{\ast }{I}_{im}^{\ast }{N}_{if}^{\ast }}{S_{im}^{\ast }{I}_{if}^{\ast }{N}_{im}^{\ast }}{m}_{ij}{I}_{jm}^{\ast }[({\eta }_{j2}\frac{S_{jm}}{S_{jm}^{\ast }}-{\eta }_{j2}\ln \frac{S_{jm}}{S_{jm}^{\ast }}+\frac{I_{jm}}{I_{jm}^{\ast }}-\ln \frac{I_{jm}}{I_{jm}^{\ast }})\\ & -({\eta }_{j2}\frac{S_{im}}{S_{im}^{\ast }}-{\eta }_{j2}\ln \frac{S_{im}}{S_{im}^{\ast }}+\frac{I_{im}}{I_{im}^{\ast }}-\ln \frac{I_{im}}{I_{im}^{\ast }})]\\ & ={\beta }_i\frac{S_{if}^{\ast }{I}_{im}^{\ast }}{N_{im}^{\ast }}[G_i(\frac{I_{if}}{N_{if}},\frac{I_{im}}{N_{im}})-G_i(I_{if},I_{im})]+\sum _{j=1}^n{m}_{ij}{I}_{jf}^{\ast }[G_{jf}(S_{jf},I_{jf})-G_{if}(S_{if},I_{if})]\\ & +\sum _{j=1}^n\frac{S_{if}^{\ast }{I}_{im}^{\ast }{N}_{if}^{\ast }}{S_{im}^{\ast }{I}_{if}^{\ast }{N}_{im}^{\ast }}{m}_{ij}{I}_{jm}^{\ast }[G_{jm}(S_{jm},I_{jm})-G_{im}(S_{im},I_{im})]\\ & \le {\beta }_i\frac{S_{if}^{\ast }{I}_{im}^{\ast }}{N_{im}^{\ast }}[G_i(I_{if},I_{im})-G_i(I_{if},I_{im})]+\sum _{j=1}^n{m}_{ij}{I}_{jf}^{\ast }[G_{jf}(S_{jf},I_{jf})-G_{if}(S_{if},I_{if})]\\ & +\sum _{j=1}^n\frac{S_{if}^{\ast }{I}_{im}^{\ast }{N}_{if}^{\ast }}{S_{im}^{\ast }{I}_{if}^{\ast }{N}_{im}^{\ast }}{m}_{ij}{I}_{jm}^{\ast }[G_{jm}(S_{jm},I_{jm})-G_{im}(S_{im},I_{im})],\end{array}$

where $\begin{array}{rl}{G}_i(I_{if},I_{im})& =\frac{I_{if}}{I_{if}^{\ast }}-\ln \frac{I_{if}}{I_{if}^{\ast }}+\frac{I_{im}}{I_{im}^{\ast }}-\ln \frac{I_{im}}{I_{im}^{\ast }},\\ G_{if}(S_{if},I_{if})& ={\eta }_{j1}\frac{S_{if}}{S_{if}^{\ast }}-{\eta }_{j1}\ln \frac{S_{if}}{S_{if}^{\ast }}+\frac{I_{if}}{I_{if}^{\ast }}-\ln \frac{I_{if}}{I_{if}^{\ast }},\\ G_{im}(S_{im},I_{im})& ={\eta }_{j2}\frac{S_{im}}{S_{im}^{\ast }}-{\eta }_{j2}\ln \frac{S_{im}}{S_{im}^{\ast }}+\frac{I_{im}}{I_{im}^{\ast }}-\ln \frac{I_{im}}{I_{im}^{\ast }}.\end{array}$Consider two weight matrices $W_f=(w_{ij}^f)$ with entry $w_{ij}^f=m_{ij}{I}_{jf}^{\ast }$ and $W_m=(w_{ij}^m)$ with entry $w_{ij}^m=m_{ij}{I}_{jm}^{\ast }$, and denote the corresponding weighted digraphs as $(G_f,W_f)$ and $(G_m,W_m)$. Let $C_i^f=\sum _{\tau \in T_i}{\omega }_f(\tau )$, $C_i^m=\sum _{\tau \in T_i}{\omega }_m(\tau )$. Then $\begin{array}{rl}& \sum _{i=1}^n{C}_i^f\sum _{j=1}^n{m}_{ij}{I}_{jf}^{\ast }[G_{jf}(S_{jf},I_{jf})-G_{if}(S_{if},I_{if})]=0,\\ & \sum _{i=1}^n{C}_i^m\sum _{j=1}^n\frac{S_{if}^{\ast }{I}_{im}^{\ast }{N}_{if}^{\ast }}{S_{im}^{\ast }{I}_{if}^{\ast }{N}_{im}^{\ast }}{m}_{ij}{I}_{jm}^{\ast }[G_{jm}(S_{jm},I_{jm})-G_{im}(S_{im},I_{im})]=0.\end{array}$Set $\begin{array}{rl}V& =\sum _{i=1}^n[C_i^f(S_{if}-S_{if}^{\ast }-S_{if}^{\ast }\ln \frac{S_{if}}{S_{if}^{\ast }})+C_i^m\frac{S_{if}^{\ast }{I}_{im}^{\ast }{N}_{if}^{\ast }}{S_{im}^{\ast }{I}_{if}^{\ast }{N}_{im}^{\ast }}(S_{im}-S_{im}^{\ast }-S_{im}^{\ast }\ln \frac{S_{im}}{S_{im}^{\ast }})\\ & +C_i^f(I_{if}-I_{if}^{\ast }-I_{if}^{\ast }\ln \frac{I_{if}}{I_{if}^{\ast }})+C_i^m\frac{S_{if}^{\ast }{I}_{im}^{\ast }{N}_{if}^{\ast }}{S_{im}^{\ast }{I}_{if}^{\ast }{N}_{im}^{\ast }}(I_{im}-I_{im}^{\ast }-I_{im}^{\ast }\ln \frac{I_{im}}{I_{im}^{\ast }})],\end{array}$then $\begin{array}{rl}\frac{\mathrm{d}V}{\mathrm{d}t}& \le \sum _{i=1}^n{C}_i^f\sum _{j=1}^n{m}_{ij}{I}_{jf}^{\ast }[G_{jf}(S_{jf},I_{jf})-G_{if}(S_{if},I_{if})]\\ & +\sum _{i=1}^n{C}_i^m\sum _{j=1}^n\frac{S_{if}^{\ast }{I}_{im}^{\ast }{N}_{if}^{\ast }}{S_{im}^{\ast }{I}_{if}^{\ast }{N}_{im}^{\ast }}{m}_{ij}{I}_{jm}^{\ast }[G_{jm}(S_{jm},I_{jm})-G_{im}(S_{im},I_{im})]=0.\end{array}$Therefore, V is a Lyapunov function for system (1(1) $\begin{array}{rl}\frac{\mathrm{d}{S}_{if}}{\mathrm{d}t}& ={\epsilon }_i{\mathrm{\Lambda }}_i-{\beta }_i\frac{S_{if}{I}_{im}}{N_{im}}-{\mu }_{if}{S}_{if}+\sum _{j=1,j\ne i}^n(m_{ij}{S}_{jf}-m_{ji}{S}_{if}),\\ \frac{\mathrm{d}{I}_{if}}{\mathrm{d}t}& ={\beta }_i\frac{S_{if}{I}_{im}}{N_{im}}-({\mu }_{if}+{\gamma }_{if})I_{if}+\sum _{j=1,j\ne i}^n(m_{ij}{I}_{jf}-m_{ji}{I}_{if}),\\ \frac{\mathrm{d}{A}_{if}}{\mathrm{d}t}& ={\gamma }_{if}{I}_{if}-({\mu }_{if}+d_{if})A_{if}+\sum _{j=1,j\ne i}^n(m_{ij}{A}_{jf}-m_{ji}{A}_{if}),\\ \frac{\mathrm{d}{S}_{im}}{\mathrm{d}t}& =(1-{\epsilon }_i){\mathrm{\Lambda }}_i-{\beta }_i\frac{S_{im}{I}_{if}}{N_{if}}-{\mu }_{im}{S}_{im}+\sum _{j=1,j\ne i}^n(m_{ij}{S}_{jm}-m_{ji}{S}_{im}),\\ \frac{\mathrm{d}{I}_{im}}{\mathrm{d}t}& ={\beta }_i\frac{S_{im}{I}_{if}}{N_{if}}-({\mu }_{im}+{\gamma }_{im})I_{im}+\sum _{j=1,j\ne i}^n(m_{ij}{I}_{jm}-m_{ji}{I}_{im}),\\ \frac{\mathrm{d}{A}_{im}}{\mathrm{d}t}& ={\gamma }_{im}{I}_{im}-({\mu }_{im}+d_{im})A_{im}+\sum _{j=1,j\ne i}^n(m_{ij}{A}_{jm}-m_{ji}{A}_{im}),\end{array}$(1) ). Since $M=(m_{ij})$ is irreducible, from Theorem 2.1 in [23], we know that $C_i^f>0$ and $C_i^m>0$ for all i, and thus $V^{\prime }=0$ implies that $S_{if}=S_{if}^{\ast },S_{im}=S_{im}^{\ast }$ for all i. From the third and sixth equations of (1(1) $\begin{array}{rl}\frac{\mathrm{d}{S}_{if}}{\mathrm{d}t}& ={\epsilon }_i{\mathrm{\Lambda }}_i-{\beta }_i\frac{S_{if}{I}_{im}}{N_{im}}-{\mu }_{if}{S}_{if}+\sum _{j=1,j\ne i}^n(m_{ij}{S}_{jf}-m_{ji}{S}_{if}),\\ \frac{\mathrm{d}{I}_{if}}{\mathrm{d}t}& ={\beta }_i\frac{S_{if}{I}_{im}}{N_{im}}-({\mu }_{if}+{\gamma }_{if})I_{if}+\sum _{j=1,j\ne i}^n(m_{ij}{I}_{jf}-m_{ji}{I}_{if}),\\ \frac{\mathrm{d}{A}_{if}}{\mathrm{d}t}& ={\gamma }_{if}{I}_{if}-({\mu }_{if}+d_{if})A_{if}+\sum _{j=1,j\ne i}^n(m_{ij}{A}_{jf}-m_{ji}{A}_{if}),\\ \frac{\mathrm{d}{S}_{im}}{\mathrm{d}t}& =(1-{\epsilon }_i){\mathrm{\Lambda }}_i-{\beta }_i\frac{S_{im}{I}_{if}}{N_{if}}-{\mu }_{im}{S}_{im}+\sum _{j=1,j\ne i}^n(m_{ij}{S}_{jm}-m_{ji}{S}_{im}),\\ \frac{\mathrm{d}{I}_{im}}{\mathrm{d}t}& ={\beta }_i\frac{S_{im}{I}_{if}}{N_{if}}-({\mu }_{im}+{\gamma }_{im})I_{im}+\sum _{j=1,j\ne i}^n(m_{ij}{I}_{jm}-m_{ji}{I}_{im}),\\ \frac{\mathrm{d}{A}_{im}}{\mathrm{d}t}& ={\gamma }_{im}{I}_{im}-({\mu }_{im}+d_{im})A_{im}+\sum _{j=1,j\ne i}^n(m_{ij}{A}_{jm}-m_{ji}{A}_{im}),\end{array}$(1) ), we obtain $\begin{array}{rl}0=\frac{\mathrm{d}{A}_{if}^{\ast }}{\mathrm{d}t}& ={\gamma }_{if}{I}_{if}-({\mu }_{if}+d_{if})A_{if}^{\ast }+\sum _{j=1}^n(m_{ij}{A}_{jf}^{\ast }-m_{ji}{A}_{if}^{\ast })\\ & ={\gamma }_{if}{I}_{if}-({\mu }_{if}+d_{if})A_{if}^{\ast }+({\mu }_{if}+d_{if})A_{if}^{\ast }-{\gamma }_{if}{I}_{if}^{\ast }={\gamma }_{if}(I_{if}-I_{if}^{\ast }),\\ 0=\frac{\mathrm{d}{A}_{im}^{\ast }}{\mathrm{d}t}& ={\gamma }_{im}{I}_{im}-({\mu }_{im}+d_{im})A_{im}^{\ast }+\sum _{j=1}^n(m_{ij}{A}_{jm}^{\ast }-m_{ji}{A}_{im}^{\ast })\\ & ={\gamma }_{im}{I}_{im}-({\mu }_{im}+d_{im})A_{im}^{\ast }+({\mu }_{im}+d_{im})A_{im}^{\ast }-{\gamma }_{im}{I}_{im}^{\ast }={\gamma }_{im}(I_{im}-I_{im}^{\ast }),\end{array}$which implies that $I_{if}=I_{if}^{\ast },I_{im}=I_{im}^{\ast }$ for all i. From the first and fourth equations of (1(1) $\begin{array}{rl}\frac{\mathrm{d}{S}_{if}}{\mathrm{d}t}& ={\epsilon }_i{\mathrm{\Lambda }}_i-{\beta }_i\frac{S_{if}{I}_{im}}{N_{im}}-{\mu }_{if}{S}_{if}+\sum _{j=1,j\ne i}^n(m_{ij}{S}_{jf}-m_{ji}{S}_{if}),\\ \frac{\mathrm{d}{I}_{if}}{\mathrm{d}t}& ={\beta }_i\frac{S_{if}{I}_{im}}{N_{im}}-({\mu }_{if}+{\gamma }_{if})I_{if}+\sum _{j=1,j\ne i}^n(m_{ij}{I}_{jf}-m_{ji}{I}_{if}),\\ \frac{\mathrm{d}{A}_{if}}{\mathrm{d}t}& ={\gamma }_{if}{I}_{if}-({\mu }_{if}+d_{if})A_{if}+\sum _{j=1,j\ne i}^n(m_{ij}{A}_{jf}-m_{ji}{A}_{if}),\\ \frac{\mathrm{d}{S}_{im}}{\mathrm{d}t}& =(1-{\epsilon }_i){\mathrm{\Lambda }}_i-{\beta }_i\frac{S_{im}{I}_{if}}{N_{if}}-{\mu }_{im}{S}_{im}+\sum _{j=1,j\ne i}^n(m_{ij}{S}_{jm}-m_{ji}{S}_{im}),\\ \frac{\mathrm{d}{I}_{im}}{\mathrm{d}t}& ={\beta }_i\frac{S_{im}{I}_{if}}{N_{if}}-({\mu }_{im}+{\gamma }_{im})I_{im}+\sum _{j=1,j\ne i}^n(m_{ij}{I}_{jm}-m_{ji}{I}_{im}),\\ \frac{\mathrm{d}{A}_{im}}{\mathrm{d}t}& ={\gamma }_{im}{I}_{im}-({\mu }_{im}+d_{im})A_{im}+\sum _{j=1,j\ne i}^n(m_{ij}{A}_{jm}-m_{ji}{A}_{im}),\end{array}$(1) ), we obtain $\begin{array}{rl}0=\frac{\mathrm{d}{S}_{if}^{\ast }}{\mathrm{d}t}& ={\epsilon }_i{\mathrm{\Lambda }}_i-{\beta }_i\frac{S_{if}^{\ast }{I}_{im}}{N_{im}}-{\mu }_{if}{S}_{if}^{\ast }+\sum _{j=1}^n(m_{ij}{S}_{jf}^{\ast }-m_{ji}{S}_{if}^{\ast })\\ & ={\beta }_i\frac{S_{if}^{\ast }{I}_{im}^{\ast }}{N_{im}^{\ast }}+{\mu }_{if}{S}_{if}^{\ast }-\sum _{j=1}^n(m_{ij}{S}_{jf}^{\ast }-m_{ji}{S}_{if}^{\ast })-{\beta }_i\frac{S_{if}^{\ast }{I}_{im}}{N_{im}}-{\mu }_{if}{S}_{if}^{\ast }+\sum _{j=1}^n(m_{ij}{S}_{jf}^{\ast }-m_{ji}{S}_{if}^{\ast })\\ & ={\beta }_i{S}_{if}^{\ast }(\frac{I_{im}^{\ast }}{N_{im}^{\ast }}-\frac{I_{im}}{N_{im}}),\\ 0=\frac{\mathrm{d}{S}_{im}^{\ast }}{\mathrm{d}t}& ={\epsilon }_i{\mathrm{\Lambda }}_i-{\beta }_i\frac{S_{im}^{\ast }{I}_{if}}{N_{if}}-{\mu }_{im}{S}_{im}^{\ast }+\sum _{j=1}^n(m_{ij}{S}_{jm}^{\ast }-m_{ji}{S}_{im}^{\ast })\\ & ={\beta }_i\frac{S_{im}^{\ast }{I}_{if}^{\ast }}{N_{if}^{\ast }}+{\mu }_{im}{S}_{im}^{\ast }-\sum _{j=1}^n(m_{ij}{S}_{jm}^{\ast }-m_{ji}{S}_{im}^{\ast })\\ & -{\beta }_i\frac{S_{im}^{\ast }{I}_{if}}{N_{if}}-{\mu }_{im}{S}_{im}^{\ast }+\sum _{j=1}^n(m_{ij}{S}_{jm}^{\ast }-m_{ji}{S}_{im}^{\ast })\\ & ={\beta }_i{S}_{im}^{\ast }(\frac{I_{if}^{\ast }}{N_{if}^{\ast }}-\frac{I_{if}}{N_{if}}),\end{array}$which implies that $N_{if}=N_{if}^{\ast },N_{im}=N_{im}^{\ast }$ for all i, thus $A_{if}=A_{if}^{\ast },A_{im}=A_{im}^{\ast }$. The largest (and only) invariant set on which $V^{\prime }=0$ is the singleton $\{E^{\ast }\}$. Therefore, by LaSalle Invariance Principle [21], $E^{\ast }$ is globally asymptotically stable when $R_0>1$.