1. Introduction

With the development of transportation and urbanisation, population migration across regions becomes more frequent, and more and more rural population crowded into cities. The increasing mobility among regions might lead to the spread of the infectious diseases regionally and globally much faster than ever before [19]. For example, SARS was first reported in Guangdong Province of China in November of 2002, and in late June of 2003, the emerging infectious disease had spread to 32 countries and regions due to the human mobility [21,25]. In February 2014 Ebola virus appeared in Guinea and then due to the human mobility the disease spread very quickly to other countries including the United States, Spain and the United Kingdom et al [12], and has caused about 6070 reported deaths and 17145 reported cases of Ebola virus disease up to December 3,2014 according to the report from the World Health Organization (WHO) [6]. All the above facts show that the population dispersal can affect transmission dynamics of the infectious diseases.

In the recent years, the impact of population dispersal has received increasing attention, and many mathematical patch models are formulated to investigate this hot issue (see [24,3,14] and the references cited therein). Here, the patches can be cities, towns, states, countries or other appropriate community divisions. Wang and Zhao [24] proposed an epidemic model with population dispersal to describe the dynamics of disease spread between two patches and n patches. Arino and van den Driessche [3] developed a multi-city epidemic model to analyze the spatial spread of infectious diseases. In 2011, Gao and Ruan [15] formulated an SIS patch model with non-constant transmission coefficients to investigate the effect of media coverage and human movement on the spread of infectious diseases among patches, and soon after, Gao and Ruan [13] proposed a multi-patch model to study the impacts of population dispersal on the spatial spread of malaria between patches. All the above mathematical models have provided useful information about the effect of host mobility on transmission dynamics of infectious diseases, but almost these models do not include the control measures, such as vaccination in it.

There is no doubt that the top priority of global public security is to prevent and contain the spread of infectious diseases. Thus it is important to study how to control the spread of infectious diseases in patchy environment and how the increasing mobility of hosts affects the current public health security. In this paper, we will use a mathematical model to explore this important issue. As we all know, vaccination is one of the most effective biological means of containing the outbreak of infectious diseases, which inoculates antigenic material into the individuals to stimulate immune system to develop adaptive immunity to a pathogen. Since Edward Jenner, the founder of vaccinology, inoculated a 13 year-old-boy with vaccinia virus (cowpox) and demonstrated immunity to smallpox [18] in 1796, vaccination has played an important role in controlling and preventing the outbreak of infectious diseases. The widespread immunity due to vaccination is largely responsible for the worldwide eradication of smallpox and the restriction of infectious diseases, such as polio, measles, and tetanus from much of the world [17]. Over the past two decades, many modeling studies have been conducted the effect of vaccination on transmission dynamics of infectious diseases (see [1,2] and reference therein). However, most of the epidemic models with vaccination are formulated in an isolated patch, ignoring spatial heterogeneity both for populations and disease transmissions.

The main purpose of the paper is to formulate an SIR epidemic model to study the impact of vaccination on transmission dynamics of infectious disease in patchy environment and the impact of the increasing mobility of hosts on the current immunization strategy. The paper is organized as follows. In Section 2, based on the SIR model with birth targeted vaccination we propose an SIR epidemic model with vaccination in patchy environment. In Section 3, we mainly present some preliminary results and derive the reproduction number. A classification of the equilibria of system on two patches and its the local dynamical behavior is provided in Section 4. We conclude with some numerical simulations in Section 5 and give a brief conclusion in the final section.

---

2. Model formulation

In this section, we employ an $n$-patch SIR epidemic model capable of informing vaccination policy to illustrate the impact of population migration between patches on the transmission dynamics of an infectious disease.

First, let us formulate a model for the spread of the disease in the $i$-th patch. Suppose there is no host migration among patches, i.e., the patches are isolated. We assume that the total host population $N_i(t)$ in the $i-$th patch is partitioned into three distinct epidemiological subclasses which are susceptible, infectious and removed either by recovery from infection or by vaccination, with sizes denoted by $S_i(t), I_i(t)$ and $R_i(t)$, respectively. Our assumptions on the dynamical transfer of the host population in the $i-$th patch are demonstrated in Figure 1.

We assume that the hosts are recruited at a rate $\mu_i N_i$ into the susceptible class and die of natural causes at a rate $\mu_i$. All individuals are born susceptible, and a fraction $p_i$ of newborns will receive vaccination. The newborns that successfully take the vaccine will then develop immunity to infection. The force of infection for susceptible is $\beta_iS_i\frac{I_i}{N_i}$, where $\beta_i$ denotes the transmission coefficient. An infected host recovers at the rate $\gamma_i$.

Based on the transfer diagram 1, the spread of an infectious disease in the $i$-th patch can be described by the following equations:

$\left\{ \begin{array}{ll} \displaystyle \frac{dS_i}{dt}=(1-p_i)\mu_i N_i-\beta_i\frac{I_i}{N_i} S_i-\mu_i S_i,\\[2ex] \displaystyle \frac{dI_i}{dt}=\beta_i\frac{I_i}{N_i} S_i-(\mu_i+\gamma_i) I_i,\\[2ex] \displaystyle \frac{dR_i}{dt}=p_i\mu_i N_i+\gamma_i I_i-\mu_i R_i.\\[2ex] \end{array}\label{model:0}\right.$

(1)

When $n$ patches are connected, we assume that only susceptible and recovered hosts can migrate among the patches and the infected hosts cannot migrate among patches due to health problems or strict border inspection. Let $m_{ij}\geq 0$ denotes the per capita rate that susceptible or recovered hosts of patch $i$ leave for patch $j$, where $i \neq j$; $d$ represents the migration period, and $l_{ij}$ denotes the fraction of the hosts in patch $i$ that move to patch $j$ and satisfies $\sum_{j=1}^nl_{ij}=1$, then the migration rate $m_{ij}$ can be calculated or estimated from available data. If we given the fraction $l_{ij}$ of the hosts migrate from patch $i$ to patch $j$ within $d$ days, the migration rate $m_{ij}$ can be determined from using the relationship $1-e^{-m_{ij}d}=l_{ij}$, that is

$m_{ij}=-\ln(1-l_{ij})^\frac{1}{d},\ \ i,j=1,2,\cdots,n,\ i\neq j.$

(2)

Then the dynamics of the hosts with migration is governed by the following model:

$\left\{ \begin{array}{ll} \displaystyle \frac{dS_i}{dt}=(1-p_i)\mu_i N_i-\beta_i\frac{I_i}{N_i} S_i-\mu_i S_i+\sum\limits_{j\neq i}^n(m_{ji}S_j-m_{ij}S_i),\\[2ex] \displaystyle \frac{dI_i}{dt}=\beta_i\frac{I_i}{N_i} S_i-(\mu_i+\gamma_i) I_i,\\[2ex] \displaystyle \frac{dR_i}{dt}=p_i\mu_i N_i+\gamma I_i-\mu_i R_i+\sum\limits_{j\neq i}^n(m_{ji}R_j-m_{ij}R_i),\\[2ex] \displaystyle N_i=S_i+I_i+R_i, i=1,2,\cdots,n. \end{array}\label{model:1}\right.$

(3)

In this paper, we will use the system (3) to investigate the effect of vaccination on transmission dynamics of infectious disease in patchy environment and the impact of the increasing mobility of hosts on the current immunization strategy.

---

3. Preliminary results and the reproduction number

We first introduce some notations which will be used throughout this paper. Let $x,y\in \mathbb{R}_+^n$, where $\mathbb{R}_+^n=\{x\in\mathbb{R}^n:x_i\ge0\ \mbox{for}\ 1\le i\le n\}$ be the positive orthant in $\mathbb{R}^n$. We write $x*y=(x_1y_1,x_2y_2,\cdots,x_ny_n)$, $y\ge x$ whenever $y-x\in \mathbb{R}_+^n$, $y>x$ whenever $y-x\in \mathbb{R}_+^n$ and $x\neq y$, and ${\rm Diag}(x)$ the $n\times n$ matrix whose diagonal is given by the components of $x$ and the other terms are zero.

Let $N_i(t)=S_i(t)+I_i(t)+R_i(t)$ be the total population in patch $i$ at time $t$ and let the new infection term in patch $i$ equal zero if $N_i=0$, see Adding all equations in (3) together leads to $(N_1(t)+N_2(t)+\cdots+N_n(t))'=0$, giving that the total population $N(t)=\sum_{i=1}^nN_i(t)$ is always a constant. If the total initial population given by $N(0)$, then the fesible region

$\Gamma=\{(S_1,I_1,R_1,\cdots,S_n,I_n,R_n)\in \mathbb{R}_+^{3n}:\sum\limits_{i=1}^n(S_i+I_i+R_i)\le N(0),i=1,2,\cdots,n \},$

is positively invariant with respect to system (3).

Define movement matrix

$M=\left(\begin{array}{cccc} \sum_{j\neq1}^nm_{1j}&-m_{21}&\cdots&-m_{n1} \\[2ex] -m_{12}&\sum_{j\neq 2}^nm_{2j}&\cdots&-m_{n2} \\[2ex] \vdots&\vdots&\ddots&\vdots \\[2ex] -m_{1n}& -m_{2n}& \cdots&\sum_{j\neq n}^nm_{nj} \end{array} \right).$

(4)

In this paper, we always assume that the movement matrix is irreducible, that is, the graph of the patches are strongly connected through the movement of hosts with respect to disease. If the movement matrix is reducible, the system may be decoupled into several samll systems (see [11] and reference therein).

To find the disease-free equilibrium with all $I_i=0$ of system (3), consider the following linear systems

$\left\{ \begin{array}{l} \displaystyle (1-p_i)\mu_i N_i -\mu_i S_i+\sum\limits_{j\neq i}^n(m_{ji}S_j-m_{ij}S_i)=0,\\[2ex] \displaystyle p_i\mu_i N_i -\mu_i R_i+\sum\limits_{j\neq i}^n(m_{ji}R_j-m_{ij}R_i)=0,\\[2ex] \displaystyle N_i=S_i+R_i, i=1,2,\cdots,n, \end{array}\right.$

(5)

or in the form of matrix systems

$\left\{\begin{array}{l} {\rm Diag}(({\bf{1}}-{\bf{p}})*\mu) {\bf{N}}-(M+{\rm Diag}(\mu)) {\bf{S}}=0,\\[2ex] {\rm Diag}({\bf{p}}*\mu) {\bf{N}}-(M+{\rm Diag}(\mu)) {\bf{R}}=0,\\[2ex] M {\bf{N}}=0,\\%[2ex] \end{array}\right.$

(6)

where $\mu=(\mu_1,\mu_2,\cdots,\mu_n)$, ${\bf{1}}=(1,1,\cdots,1)$ and ${\bf{X}}=(X_1,X_2,\cdots, X_n)^T$ (subscript $T$ denotes transpose), $X$ represents $S,I,R,N$ and $p$.

We first solve the third equation of (5) or (6) which independent of the first two equations. Applying the results presented in Lemma 2.1 [16], the general solution to the third equation of (6) can be given as ${\bf{N}}=k(c_{11},c_{22},\cdots,c_{nn})^T$, where $c_{kk}>0,\ (k=1,2,\cdots,n)$ is the co-factor of the $k$-th diagonal entry of $M$ and $k$ is a constant to be specified later. It follows from $\sum_{i=1}^n N_i = N(0)$ that $k=\frac{N(0)}{\sum_{i=1}^nc_{ii}}>0$, then the third equation of (6) has a unique positive solution

${\bf{N}}^0\triangleq(N_1^0,N_2^0,\cdots,N_n^0)^T=\frac{N(0)}{\sum_{i=1}^nc_{ii}}(c_{11},c_{22},\cdots,c_{nn})^T.$

Substituting ${\bf{N}}^0$ back into the first two equations of (6) yields that

$\begin{array}{l} {\bf{S}}^0\triangleq(S_1^0,S_2^0,\cdots,S_n^0)=(M+{\rm Diag}(\mathbf{\mu}))^{-1}{\rm Diag}(({\bf{1}}-{\bf{p}})*\mu) {\bf{N}}^0,\\[2ex] {\bf{R}}^0\triangleq(R_1^0,R_2^0,\cdots,R_n^0)=(M+{\rm Diag}(\mu))^{-1}{\rm Diag}({\bf{p}}*\mu){\bf{N}}^0. \end{array}$

(7)

Since all off-diagonal entries of matrix $M+{\rm Diag}(\mu)$ are negative (i.e., the Z-sign pattern) and the sum of the entries in each column is positive, it follows from Chapter 6 in [5] that $M+{\rm Diag}(\mu)$ is a nonsingular irreducible M-matrix and $(M+{\rm Diag}(\mu))^{-1}>0$. Therefore, the linear system (6) has a unique positive solution ${\bf{S}}^0, {\bf{R}}^0, {\bf{N}}^0$, that is, system (3) always has unique disease-free equilibrium $E^0=({\bf{S}}^0,0, {\bf{R}}^0)$.

In absence of infectious disease, adding the three equations of system (3) together leads to

$\frac{dN_i}{dt}=\sum\limits_{j\neq i}^n(m_{ji} N_j -m_{ij}N_i),\ \ i=1,2,\cdots,n,$

(8)

or in the form of matrix system

$\frac{d{\bf{N}}(t) }{dt}=-M {\bf{N}} (t).$

(9)

It follows from Theorem 2.1 in [4] that the positive equilibrium ${\bf{N}}^0= \frac{N(0)}{\sum_{i=1}^nc_{ii}} (c_{11}, c_{22},$ $\cdots, c_{nn})^T$ of system (8) is globally asymptotically stable on the affine hyperplane $\sum_{i=1}^n N_i^0 = N(0)$. Using the expressions of ${\bf{S}}^0,{\bf{I}}^0, {\bf{N}}^0$ and the condition $I_i=0$ for all $i$, we can transfer the stability of system (3) into the following limit system

$\left\{\begin{array}{l} ({\bf{S}}-{\bf{S}}^0)'=-(M+{\rm Diag}(\mu))({\bf{S}}-{\bf{S}}^0),\\[2ex] ({\bf{R}}-{\bf{R}}^0)'=-(M+{\rm Diag}(\mu))({\bf{R}}-{\bf{R}}^0). \end{array}\right.$

(10)

Since the Gerschgorin circular disc theorem implies that matrix $M+{\rm Diag}(\mu)$ is stable, i.e., all the eigenvalues of $M+{\rm Diag}(\mu)$ have negative real parts, then equilibrium ${\bf{S}}={\bf{S}}^0, {\bf{R}}={\bf{R}}^0$ is global stability of system (10). Following the Theorem 2.1 in [8], we can directly obtain that the unique positive equilibrium $E^0$ of linear system (3) is globally asymptotically stable on the affine hyperplane $\sum_{i=1}^n N_i^0 = N(0)$. Summarizing the above discussions, we can obtain the following result.

Theorem 3.1. System (3) always has a disease-free equilibrium $E^0({\bf{S}}^0,{\bf{0}},{\bf{R}}^0)$, where

$\begin{array}{l} \mathit{\boldsymbol{S}}^0=(M+{\rm Diag}(\mu))^{-1}{\rm Diag}(({\bf{1}}-\mathit{\boldsymbol{p}})\ast\mu) \mathit{\boldsymbol{N}}^0, \\[2ex] \mathit{\boldsymbol{R}}^0=(M+{\rm Diag}(\mu))^{-1}{\rm Diag}(\mathit{\boldsymbol{p}}*\mu) \mathit{\boldsymbol{N}}^0,\\[2ex] \mathit{\boldsymbol{N}}^0= \frac{N(0)}{\sum_{i=1}^nc_{ii}}(c_{11},c_{22},\cdots,c_{nn})^T, \end{array}$

and $N(0)$ is the initial total population and $c_{kk}>0,\ (k=1,2,\cdots,n)$ is the co-factor of the $k$-th diagonal entry of the movement matrix $M$. Moreover, the disease-free equilibrium $E^0({\bf{S}}^0,{\bf{0}},{\bf{R}}^0)$ is globally asymptotically stable in

$\Gamma^0\mspace{-3mu}=\mspace{-3mu}\{(S_1,\mspace{-3mu}\cdots,\mspace{-3mu}S_n,\mspace{-3mu}I_1,\mspace{-3mu}\cdots, \mspace{-3mu}I_n,\mspace{-3mu}R_1,\mspace{-3mu}\cdots,\mspace{-3mu} R_n):\mspace{-3mu}\sum\limits_{i=1}^n\mspace{-3mu} (S_i+R_i)\mspace{-3mu}=\mspace{-3mu}N(0),\mspace{-3mu} I_i=0,i=1,2,\cdots,n\}.$

Note that the system (3) has $n$ infected variables, namely, $I_1,I_2,\cdots, I_n$, it then follows that, using the notation of Driessche and Watmough [10], the matrices $F$ and $V$ (corresponding to the new infection terms and the remaining transfer terms, respectively) for entire population are given by

$F={\rm Diag}(\beta_1\frac{S_1^0}{N_1^0},\beta_2\frac{S_2^0}{N_2^0},\cdots,\beta_n\frac{S_n^0}{N_n^0})\ {\rm and}\ V={\rm Diag}(\mu+\gamma).$

From literature [10], the reproduction number $\mathfrak{R}_v$ is defined as the spectral radius of the next generation matrix $FV^{-1}$, that is,

$\mathfrak{R}_v= \rho( \mathcal{F}\mathcal{V}^{-1})=\max\{\mathfrak{R}_{v1},\mathfrak{R}_{v2},\cdots,\mathfrak{R}_{vn}\},$

(11)

where $\rho(M)$ represents the spectral radius of the nonnegative matrix $M$ and

$\mathfrak{R}_{vi}=\frac{\beta_i}{\mu_i+\gamma_i}\frac{S_i^0}{N_i^0},$

(12)

which represents the reproduction number in the $i$-th patch. It is clearly that the domain $\mathfrak{R}_v$ implicitly depend on the movement of susceptible individuals and vaccination coverage through $S_i^0$. For disease-free equilibrium, we then have the following result based on the Theorem 2 in [10].

Theorem 3.2. The disease-free equilibrium $E^0$ of system (3) is locally asymptotically stable if $\mathfrak{R}_v<1$, whereas it is unstable if $\mathfrak{R}_v>1$.

In the special case $n=1$, the control reproduction number has the explicit expression

$\mathfrak{R}_v=(1-p_1)\frac{\beta_1}{\mu_1+\gamma_1}.$

(13)

which represents the numbers of secondary cases directly produced by infectious diease during the period of infection in a susceptible population.

In the special case of no movement between patches (i.e., $M=0$), the control reproduction number $\mathfrak{R}_v$ defined in (11) is given by the maximum value of control reproduction numbers $\mathfrak{R}_{vi}$ in all patches. Namely,

$\mathfrak{R}_v= \max\{\mathfrak{R}_{v1},\mathfrak{R}_{v2},\cdots,\mathfrak{R}_{vn}\},$

(14)

with $\mathfrak{R}_{vi}=(1-p_i)\frac{\beta_i}{\mu_i+\gamma_i}.$

Theorem 3.3. If $\mathfrak{R}_{v}<1$, then system (3) exists only one equilibrium whose coordinates includes zero, that is, the disease-free equilibrium $E^0$.

Proof. For any equilibrium $E({\bf{S}},{\bf{I}},{\bf{R}})$ of system (3), it must satisfy the matrix system:

$\left\{ \begin{array}{ll} \displaystyle {\rm Diag}((1-p)*\mu) {\bf{N}}-{\rm Diag}(\mu+\gamma) {\bf{I}}-(M+{\rm Diag}(\mu)){\bf{S}}=0,\\[2ex] \displaystyle\displaystyle {\rm Diag}({\bf{I}})({\bf{S}}-B{\bf{N}})=0,\\[2ex] \displaystyle {\rm Diag}(p*\mu) {\bf{N}}+{\rm Diag}(\gamma) {\bf{I}}-(M+{\rm Diag}(\mu) {\bf{R}}=0,\\[2ex] \displaystyle {\bf{N}} = {\bf{S}}+{\bf{I}}+{\bf{R}}. \end{array}\right.$

(15)

where $B={\rm Diag}(\frac{\mu_1+\gamma_1}{\beta_1},\frac{\mu_2+\gamma_2}{\beta_2},\cdots, \frac{\mu_n+\gamma_n}{\beta_n})$.

Adding the first three equations of (15) together yields $M({\bf{N}}-{\bf{I}})=0$, whose general solution can be given as (see Lemma 2.1 [16])

${\bf{N}}-{\bf{I}}=k(c_{11},c_{22},\cdots,c_{nn})^T,$

where $c_{kk}>0,\ (k=1,2,\cdots,n)$ is the co-factor of the $k$-th diagonal entry of $M$ and $k$ is constant to be specified later. Since $N(0)=\sum_{i=1}^nN_i$, then direct calculation implies that $k=\frac{1}{\sum_{i=1}^nc_{ii}}N(0)-\frac{1}{\sum_{i=1}^nc_{ii}}{\bf{1}} \cdot {\bf{I}}.$ Therefore, we can rewrite ${\bf{N}}$ as

$\begin{array}{rl} {\bf{N}}={\bf{N}}^0+(\mathbb{E}-C){\bf{I}}, \end{array}$

(16)

where ${\bf{N}}^0= \frac{N(0)}{\sum_{i=1}^nc_{ii}}(c_{11},c_{22},\cdots,c_{nn})^T$, $C=\big(\frac{c_{11}}{\sum_{i=1}^nc_{ii}},\frac{c_{22}}{\sum_{i=1}^nc_{ii}}, \cdots,\frac{c_{nn}}{\sum_{i=1}^nc_{ii}}\big)^T\cdot {\bf{1}}$ and $\mathbb{E}$ is the identity matrix. Substituting (16) into the first equation of (15) leads to

$\begin{array}{rl} {\bf{S}}=&{\bf{S}}^0+(M+{\rm Diag}(\mu))^{-1}\Big({\rm Diag}(({\bf{1}}-{\bf{p}})*\mu)(\mathbb{E}-C)-{\rm Diag}(\mu+\gamma)\Big){\bf{I}}, \end{array}$

(17)

where ${\bf{S}}^0=(M+{\rm Diag}(\mu))^{-1}{\rm Diag}(({\bf{1}}-{\bf{p}})\mu) {\bf{N}}^0$.

Substituting (16), (17) into the second equation of (15), the system of equation (15) can be reduced to the following equation with one single equation of ${\bf{I}}$

$\begin{array}{l} \displaystyle {\rm Diag}({\bf{I}})\Big({\bf{S}}^0-B{\bf{N}}^0-((M+{\rm Diag}(\mu))^{-1}{\rm Diag}(({\bf{1}}-{\bf{p}})\mu)C-B C){\bf{I}}\\ -(B+(M+{\rm Diag}(\mu))^{-1}({\rm Diag}(\mu+\gamma)-{\rm Diag}(({\bf{1}}-{\bf{p}})\mu)){\bf{I}}\Big)=0. \end{array}$

(18)

In the following, we only need to solve (18) for ${\bf{I}}$ and then back-substitute into (15), the solutions for other variables will be found.

Since $CN(0)=N(0)C={\bf{N}}^0\cdot {\bf{1}}$, it follows from the relationship between ${\bf{S}}^0$ and ${\bf{N}}^0$ that

$\begin{array}{rl} &\displaystyle(M+{\rm Diag}(\mu))^{-1}{\rm Diag}(({\bf{1}}-{\bf{p}})*\mu)C-BC \\ &=\displaystyle (M+{\rm Diag}(\mu))^{-1}\left({\rm Diag}(({\bf{1}}-{\bf{p}})*\mu){\bf{N}}^0-(M+{\rm Diag}(\mu))B{\bf{N}}^0\right)\frac{1}{N(0)}{\bf{1}} \\ &=\displaystyle \left({\bf{S}}^0-B{\bf{N}}^0\right)\frac{1}{N(0)}{\bf{1}}, \end{array}$

(19)

and the expression for $\mathfrak{R}_{vi}$ given in (12) that

${\bf{S}}^0-B{\bf{N}}^0=\displaystyle \left( \begin{array}{c} \displaystyle\frac{\mu_1+\gamma_1}{\beta_1}N_1^0(\mathfrak{R}_{v1}-1)\\ \displaystyle\frac{\mu_2+\gamma_2}{\beta_2}N_2^0(\mathfrak{R}_{v2}-1)\\ \vdots\\ \displaystyle\frac{\mu_n+\gamma_n}{\beta_n}N_n^0(\mathfrak{R}_{vn}-1) \end{array} \right).$

(20)

Therefore, the equation (18) can be expressed as

${\rm Diag}({\bf{I}})(M+{\rm Diag}(\mu))^{-1}({\bf{b}}-A{\bf{I}})=0,$

(21)

where

$\label{2} {\bf{b}}\triangleq(b_1,b_2,\cdots,b_n)^T=(M+{\rm Diag}(\mu))({\bf{S}}^0-B{\bf{N}}^0),$

and

$\label{1} \begin{array}{rl} A\triangleq (a_{ij})_{n\times n}=&{\rm Diag}(\gamma+{\bf{p}}\mu)+(M+{\rm Diag}(\mu))B\\[2ex] &\displaystyle +(M+{\rm Diag}(\mu))\left({\bf{S}}^0-B{\bf{N}}^0\right)\frac{1}{N(0)}{\bf{1}}.\\[2ex] \end{array}$

Note that $M+{\rm Diag}(\mu)$ is a nonsingular M-matrix, it follows from Chapter 6 in [5] and equation (20) that ${\bf{b}}<{\bf{0}}$ if $\mathfrak{R}_v<1$. Since all off-diagonal entries are negative and every column sum of matrix $A$ is positive if $\mathfrak{R}_v<1$, then matrix $A$ is a nonsingular M-matrix. Denote sub-matrix $A(i_1,i_2,\cdots,i_k)$ which composed by the $i_1$-th, $i_2$-th, $\cdots,i_k$-th rows and columns of $A$, similarly, one can verify that sub-matrix $A(i_1,i_2,\cdots,i_k)$ is also a non-singular M-matrix.

It is easily to see that ${\bf{I}}={\bf{0}}$ is one solution of (21). Except ${\bf{I}}={\bf{0}}$, the other solution of equation (21) should be ${\bf{I}}\neq{\bf{0}}$. For ease of presentation, we discuss the cases $I_i\neq0$ for all $i$ and $I_i\neq0$ for some $i$ separately. If $I_i\neq 0$ for all $i$, then the solution of ${\bf{I}}$ must satisfy ${\bf{b}}-A{\bf{I}}=0$, it follows from ${\bf{b}}<{\bf{0}}$ and $A$ is a nonsingular M-matrix that ${\bf{I}}=A^{-1}{\bf{b}}<{\bf{0}}$. If $I_i\neq0$ for some $i$, without loss of generality, we assume that $I_{i1}\neq0,I_{i2}\neq0,\cdots,I_{ik}\neq0$ and $I_{i,k+1}=I_{i,k+2}=\cdots=I_{in}=0$, then solution of $I_{i1},I_{i2},\cdots,I_{ik}$ satisfy ${\bf{b}}(i_1,i_2,\cdots,i_k)-A(i_1,i_2,\cdots,i_k)(I_{i1},I_{i2},\cdots,I_{ik})^T=0$, where ${\bf{b}}(i_1,i_2,\cdots,i_k)$ composed with the $i_1$-th, $i_2$-th, $\cdots,i_k$-th rows of ${\bf{b}}$. It follows from ${\bf{b}}(i_1,i_2,\cdots,i_k)<{\bf{0}}$ and $A(i_1,i_2,\cdots,i_k)$ is a nonsingular M-matrix that $I_{ij}<0$ for $j=1,2,\cdots,k$. That is, the solutions of (21) when $\mathfrak{R}_v<1$ either equal to zero or less than zero. Considering the biological significance, we omit the solution including the element $I_i<0$ $(i=1,2,\cdots,n)$. Then back-substituting ${\bf{I}}=0$ into (21), we obtain the expression of $E({\bf{S}},{\bf{I}},{\bf{R}})$ which in fact is the disease-free equilibrium $E^0$. Therefore, the disease-free equilibrium $E^0$ is the unique equilibrium under the condition $\mathfrak{R}_v<1$ for system (3). This completes the proof.

---

4. The SIR model on two patches

In this section, we mainly consider the dynamic behaviors for system (3) with $n=2$ due to it is hard to find the explicit solutions of (21) when $n$ is very large as $\mathfrak{R}_v>1$. In this case, following the results presented in the previous section, we know that disease-free equilibrium $E^0=(S_1^0,S_2^0,0,0,R_1^0,R_2^0)$ always exists and the explicit expression can be computed as

$\label{dfe2} \begin{array}{l} \displaystyle S_1^0=\frac{((1-p_1)(\mu_1\mu_2+\mu_1m_{21})+(1-p_2)\mu_2m_{12})m_{21}N(0)} {(\mu_1\mu_2+\mu_1m_{21}+\mu_2m_{12})(m_{12}+m_{21})},\\[2ex] \displaystyle S_2^0=\frac{((1-p_1)\mu_1m_{21}+(1-p_2)(\mu_1\mu_2+\mu_2m_{12})m_{12}N(0)} {(\mu_1\mu_2+\mu_1m_{21}+\mu_2m_{12})(m_{12}+m_{21})},\\[2ex] \displaystyle R_1^0=\frac{(p_1(\mu_1\mu_2+\mu_1m_{21})+p_2\mu_2m_{12})m_{21}N(0)} {(\mu_1\mu_2+\mu_1m_{21}+\mu_2m_{12})(m_{12}+m_{21})},\\[2ex] \displaystyle R_2^0=\frac{(p_1\mu_1m_{21}+p_2(\mu_1\mu_2+\mu_2m_{12})m_{12}N(0)} {(\mu_1\mu_2+\mu_1m_{21}+\mu_2m_{12})(m_{12}+m_{21})}.\\[2ex] \end{array}$

From (11) and (12), the control reproduction number for this case can be given by

$\mathfrak{R}_v=\max\{\mathfrak{R}_{v1},\mathfrak{R}_{v2}\},$

(22)

where

$\begin{array}{l} \displaystyle\mathfrak{R}_{v1}=\frac{\beta_1}{\mu_1+\gamma_1}\frac{(1-p_1)(\mu_1\mu_2+\mu_1m_{21})+(1-p_2)\mu_2m_{12}}{\mu_1\mu_2+\mu_1m_{21}+\mu_2m_{12}},\\ \displaystyle \mathfrak{R}_{v2}=\frac{\beta_2}{\mu_2+\gamma_2}\frac{(1-p_1)\mu_1m_{21}+(1-p_2)(\mu_1\mu_2+\mu_2m_{12})}{\mu_1\mu_2+\mu_1m_{21}+\mu_2m_{12}}, \end{array}$

(23)

represent the control reproduction number correspond to the sub-patch 1 and 2, respectively.

Like in the single patch model (1) or many other epidemic models, we have the global stability of the disease-free equilibrium for system (3) with $n=2$ as $\mathfrak{R}_v<1$.

Theorem 4.1. If $\mathfrak{R}_v<1$, then $E^0$ is globally asymptotically stable for system (3) without vaccination, whereas if $\mathfrak{R}_v>1$, $E^0$ is unstable

The proof of Theorem (4.1) is analogous to those of Theorem 2.4 in Gao and Ruan [15] and Theorem 3.2 in Sun et al. [23]. We omit the details here.

Following Theorem 4.1 and the proof of Theorem 4.1, when $\mathfrak{R}_v<1$, there does not exist any endemic equilibrium. We now turn to the case where $\mathfrak{R}_v>1$. We first study the existence of other equilibria for system (3) with $n=2$ when $\mathfrak{R}_v>1$, and then establish its the stability.

For convenience of presentation, set

$\begin{array}{l} \xi_1=\mu_1\mu_2+\mu_1m_{21}+\mu_2m_{12}, \ \xi_2=(\gamma_2+p_2\mu_2)\beta_2, \ \xi_3=(\gamma_1+p_1\mu_1)\beta_1, \\ \xi_4=(\mu_1+\gamma_1)(\mu_2+\gamma_2)(\mu_2+m_{21})+(\mu_1+\gamma_1)\xi_2+(\mu_2+\gamma_2)^2m_{12},\\ \xi_5=(\mu_1+\gamma_1)(\mu_2+\gamma_2)(\mu_1+m_{12})+(\mu_2+\gamma_2)\xi_3+(\mu_1+\gamma_1)^2m_{21}, \end{array}$

and define

$\begin{array}{rl} \displaystyle \mathfrak{\bar{R}}_{v1}=\frac{((\mu_2+\gamma_2)(\mu_2+m_{21})+\xi_2)(1-p_1)\mu_1\beta_1+(\mu_2+\gamma_2)^2m_{12}\beta_1} {(\mu_1+\gamma_1)((\mu_2+\gamma_2)\xi_1+(\mu_1+m_{12})\xi_2)},\\ \displaystyle \mathfrak{\bar{R}}_{v2}=\frac{((\mu_1+\gamma_1)(\mu_1+m_{12})+\xi_3)(1-p_2)\mu_2\beta_2+(\mu_1+\gamma_1)^2m_{21}\beta_2} {(\mu_2+\gamma_2)((\mu_1+\gamma_1)\xi_1+(\mu_2+m_{21})\xi_3)}, \end{array}$

(24)

which can be considered as a second threshold for epidemic invasion of sub-populations 1 and 2, respectively.

Theorem 4.2. The system (3) can have other three equilibria, and we have the following results:

1. Boundary equilibria $\hat{E}=(\hat{S}_1,\hat{S}_2,\hat{I}_1,0,\hat{R}_1,\hat{R}_2)$ exists if and only if $\mathfrak{R}_{v1}>1$, and $\bar{E}=(\bar{S}_1,\bar{S}_2,0,\bar{I}_2,\bar{R}_1,\bar{R}_2)$ exists if and only if $\mathfrak{R}_{v2}>1$. Here,

$\label{en:boundaryI2}\begin{array}{l} \displaystyle \hat{S}_1\mspace{-3mu}=\mspace{-3mu}\frac{((\mu_1+\gamma_1)(\mu_2+m_{21})\mspace{-3mu}+\mspace{-3mu}m_{12}(1-p_2)\mu_2)(\mu_1+\gamma_1)m_{21}N(0)} {(\mu_1\mspace{-3mu}+\mspace{-3mu}\gamma_1)(m_{12}\xi_1\mspace{-3mu}+\mspace{-3mu} (\mu_2\mspace{-3mu}+\mspace{-3mu}m_{21})m_{21}\beta_1)\mspace{-3mu}+\mspace{-3mu}m_{12} (\mu_2\mspace{-3mu}+\mspace{-3mu}m_{21})\xi_3\mspace{-3mu}+\mspace{-3mu}m_{12}m_{21} (1-p_2)\mu_2\beta_1},\\ \displaystyle\hat{S}_2\mspace{-3mu}=\mspace{-3mu}\frac{((\mu_1\mspace{-3mu}+\mspace{-3mu}\gamma_1)^2m_{21}\mspace{-3mu} +\mspace{-3mu}((\mu_1\mspace{-3mu}+\mspace{-3mu}\gamma_1)(\mu_1\mspace{-3mu} +\mspace{-3mu}m_{12})\mspace{-3mu}+\mspace{-3mu}(\gamma_1\mspace{-3mu} +\mspace{-3mu}p_1\mu_1)\beta_1)(1\mspace{-3mu}-\mspace{-3mu}p_2)\mu_2)m_{12}N(0)} {(\mu_1\mspace{-3mu}+\mspace{-3mu}\gamma_1)(m_{12}\xi_1\mspace{-3mu}+\mspace{-3mu}(\mu_2\mspace{-3mu}+\mspace{-3mu}m_{21})m_{21}\beta_1)\mspace{-3mu}+\mspace{-3mu}m_{12}(\mu_2\mspace{-3mu}+\mspace{-3mu}m_{21})\xi_3\mspace{-3mu}+\mspace{-3mu}m_{12}m_{21} (1-p_2)\mu_2\beta_1},\\ \displaystyle\hat{I}_1\mspace{-3mu}=\mspace{-3mu}\frac{(\mu_1\mspace{-3mu}+\mspace{-3mu}\gamma_1)\xi_1m_{21}(\mathfrak{R}_{v1}-1)N(0)} {(\mu_1\mspace{-3mu}+\mspace{-3mu}\gamma_1)(m_{12}\xi_1\mspace{-3mu}+\mspace{-3mu} (\mu_2\mspace{-3mu}+\mspace{-3mu}m_{21})m_{21}\beta_1)\mspace{-3mu}+\mspace{-3mu}m_{12} (\mu_2\mspace{-3mu}+\mspace{-3mu}m_{21})\xi_3\mspace{-3mu}+m_{12}m_{21} (1-p_2)\mu_2\beta_1},\\ \displaystyle\hat{R}_1\mspace{-3mu}=\mspace{-3mu}\frac{((\mu_2\mspace{-3mu}+\mspace{-3mu}m_{21})(\xi_3-(\mu_1\mspace{-3mu}+\mspace{-3mu}\gamma_1)\gamma_1) \mspace{-3mu}+\mspace{-3mu}m_{12}(\mu_1\mspace{-3mu}+\mspace{-3mu}\gamma_1)p_2\mu_2)m_{21}N(0)} {(\mu_1\mspace{-3mu}+\mspace{-3mu}\gamma_1)(m_{12}\xi_1\mspace{-3mu}+\mspace{-3mu}(\mu_2\mspace{-3mu}+\mspace{-3mu}m_{21})m_{21}\beta_1)\mspace{-3mu}+\mspace{-3mu}m_{12}(\mu_2\mspace{-3mu}+\mspace{-3mu}m_{21})\xi_3\mspace{-3mu}+\mspace{-3mu}m_{12}m_{21} (1-p_2)\mu_2\beta_1},\\ \displaystyle\hat{R}_2\mspace{-3mu}=\mspace{-3mu}\frac{(m_{21}(\xi_3-(\mu_1\mspace{-3mu}+\mspace{-3mu}\gamma_1)\gamma_1)\mspace{-3mu}+\mspace{-3mu}((\mu_1\mspace{-3mu}+\mspace{-3mu}\gamma_1)(\mu_1\mspace{-3mu}+\mspace{-3mu}m_{12}) \mspace{-3mu}+\mspace{-3mu}(\gamma_1\mspace{-3mu}+\mspace{-3mu}p_1\mu_1)\beta_1)p_2\mu_2)m_{12}N(0)} {(\mu_1\mspace{-3mu}+\mspace{-3mu}\gamma_1)(m_{12}\xi_1\mspace{-3mu}+\mspace{-3mu}(\mu_2\mspace{-3mu}+\mspace{-3mu}m_{21})m_{21}\beta_1)\mspace{-3mu}+\mspace{-3mu}m_{12}(\mu_2\mspace{-3mu}+\mspace{-3mu}m_{21})\xi_3\mspace{-3mu}+\mspace{-3mu}m_{12}m_{21} (1-p_2)\mu_2\beta_1}, \end{array}$

and

$\label{en:boundaryI1} \begin{array}{l} \displaystyle \displaystyle \bar{S}_1\mspace{-3mu}=\mspace{-3mu}\frac{((\mu_2\mspace{-3mu}+\mspace{-3mu}\gamma_2)^2m_{12}\mspace{-3mu}+\mspace{-3mu}((\mu_2\mspace{-3mu}+\mspace{-3mu}\gamma_2)(\mu_2\mspace{-3mu}+\mspace{-3mu}m_{21})\mspace{-3mu}+\mspace{-3mu}(\gamma_2\mspace{-3mu}+\mspace{-3mu}p_2\mu_2)\beta_2)(1-p_1)\mu_1)m_{21}N(0)} {(\mu_2\mspace{-3mu}+\mspace{-3mu}\gamma_2)(m_{21}\xi_1\mspace{-3mu}+\mspace{-3mu}(\mu_1\mspace{-3mu}+\mspace{-3mu}m_{12})m_{12}\beta_2)\mspace{-3mu}+\mspace{-3mu}m_{21}(\mu_1\mspace{-3mu}+\mspace{-3mu}m_{12})\xi_2\mspace{-3mu}+\mspace{-3mu}m_{12}m_{21} (1-p_1)\mu_1\beta_2}, \\ \displaystyle \bar{S}_2\mspace{-3mu}=\mspace{-3mu}\frac{((\mu_2\mspace{-3mu}+\mspace{-3mu}\gamma_2)(\mu_1\mspace{-3mu}+\mspace{-3mu}m_{12})\mspace{-3mu}+\mspace{-3mu}m_{21}(1-p_1)\mu_1)(\mu_2\mspace{-3mu}+\mspace{-3mu}\gamma_2)m_{12}N(0)} {(\mu_2\mspace{-3mu}+\mspace{-3mu}\gamma_2)(m_{21}\xi_1\mspace{-3mu}+\mspace{-3mu}(\mu_1\mspace{-3mu}+\mspace{-3mu}m_{12})m_{12}\beta_2)\mspace{-3mu}+\mspace{-3mu}m_{21}(\mu_1\mspace{-3mu}+\mspace{-3mu}m_{12})\xi_2\mspace{-3mu}+\mspace{-3mu}m_{12}m_{21} (1-p_1)\mu_1\beta_2},\\ \displaystyle\bar{I}_2\mspace{-3mu}=\mspace{-3mu}\frac{(\mu_2\mspace{-3mu}+\mspace{-3mu}\gamma_2)\xi_1m_{12}(\mathfrak{R}_{v2}-1)N(0)} {(\mu_2\mspace{-3mu}+\mspace{-3mu}\gamma_2)(m_{21}\xi_1\mspace{-3mu}+\mspace{-3mu}(\mu_1\mspace{-3mu}+\mspace{-3mu}m_{12})m_{12}\beta_2)\mspace{-3mu}+\mspace{-3mu}m_{21}(\mu_1\mspace{-3mu}+\mspace{-3mu}m_{12})\xi_2\mspace{-3mu}+\mspace{-3mu}m_{12}m_{21} (1-p_1)\mu_1\beta_2},\\%[11pt] \end{array}$

$\begin{array}{l} \displaystyle \bar{R}_1\mspace{-3mu}=\mspace{-3mu}\frac{(m_{12}(\xi_2-(\mu_2\mspace{-3mu}+\mspace{-3mu}\gamma_2)\gamma2)\mspace{-3mu}+\mspace{-3mu}((\mu_2\mspace{-3mu}+\mspace{-3mu}\gamma_2)(\mu_2\mspace{-3mu}+\mspace{-3mu}m_{21}) \mspace{-3mu}+\mspace{-3mu}(\gamma_2\mspace{-3mu}+\mspace{-3mu}p_2\mu_2)\beta_2)p_1\mu_1)m_{21}N(0)} {(\mu_2\mspace{-3mu}+\mspace{-3mu}\gamma_2)(m_{21}\xi_1\mspace{-3mu}+\mspace{-3mu}(\mu_1\mspace{-3mu}+\mspace{-3mu}m_{12})m_{12}\beta_2)\mspace{-3mu}+\mspace{-3mu}m_{21}(\mu_1\mspace{-3mu}+\mspace{-3mu}m_{12})\xi_2\mspace{-3mu}+\mspace{-3mu}m_{12}m_{21} (1-p_1)\mu_1\beta_2},\\ \displaystyle\bar{R}_2\mspace{-3mu}=\mspace{-3mu}\frac{((\mu_1\mspace{-3mu}+\mspace{-3mu}m_{12})(\xi_2-(\mu_2\mspace{-3mu}+\mspace{-3mu}\gamma_2)\gamma_2) \mspace{-3mu}+\mspace{-3mu}m_{21}(\mu_2\mspace{-3mu}+\mspace{-3mu}\gamma_2)p_1\mu_1)m_{12}N(0)} {(\mu_2\mspace{-3mu}+\mspace{-3mu}\gamma_2)(m_{21}\xi_1\mspace{-3mu}+\mspace{-3mu}(\mu_1\mspace{-3mu}+\mspace{-3mu}m_{12})m_{12}\beta_2)\mspace{-3mu}+\mspace{-3mu}m_{21}(\mu_1\mspace{-3mu}+\mspace{-3mu}m_{12})\xi_2\mspace{-3mu}+\mspace{-3mu}m_{12}m_{21} (1-p_1)\mu_1\beta_2}.\\%[11pt] \end{array}$

2. System (3) has a unique endemic equilibrium $E^*=(S_1^*,S_2^*,I_1^*,I_2^*,R_1^*,R_2^*)$ if and only if $\mathfrak{\bar{R}}_{v1}>1$ and $\mathfrak{\bar{R}}_{v2}>1$. Here

$\begin{split} \displaystyle S_1^*=&\frac{((\mu_1+\gamma_1)((\mu_2+\gamma_2)(\mu_2+m_{21})+\xi_2)+(\mu_2+\gamma_2)^2m_{12})(\mu_1+\gamma_1)m_{21}N(0)} {m_{21}\beta_1\xi_4+m_{12}\beta_2\xi_5},\\ \displaystyle S_2^*=&\frac{((\mu_2+\gamma_2)((\mu_1+\gamma_1)(\mu_1+m_{12})+\xi_3)+(\mu_1+\gamma_1)^2m_{21})(\mu_2+\gamma_2)m_{12}N(0)} {m_{21}\beta_1\xi_4+m_{12}\beta_2\xi_5},\\ \displaystyle I_1^*=&\frac{(\mu_1+\gamma_1)((\mu_1+m_{12})\xi_2+(\mu_2+\gamma_2)\xi_1)(\mathfrak{\bar{R}}_{v1}-1)m_{21}N(0)} {m_{21}\beta_1\xi_4+m_{12}\beta_2\xi_5},\\ \displaystyle I_2^*=&\frac{(\mu_2+\gamma_2)((\mu_2+m_{21})\xi_3+(\mu_1+\gamma_1)\xi_1)(\mathfrak{\bar{R}}_{v2}-1)m_{12}N(0)} {m_{21}\beta_1\xi_4+m_{12}\beta_2\xi_5},\\ \displaystyle R_1^*=&\frac{(\xi_2 +(\mu_2+ \gamma_2) (\mu_2 + m_{21}) ) (\xi_3- \gamma_1(\mu_1 + \gamma_1) ) m_{21} N (0)} {m_{21}\beta_1\xi_4 + m_{12}\beta_2\xi_5}\\ &+\frac{m_{12}(\mu_1 + \gamma_1) (\xi_2 - \gamma_2(\mu_2 + \gamma_2) ) m_{21} N (0)}{m_{21}\beta_1\xi_4 + m_{12}\beta_2\xi_5},\\ \displaystyle R_2^*=&\frac{ (\xi_3 + (\mu_1 + \gamma_1) (\mu_1 + m_{12}) ) (\xi_2 - \gamma_2 (\mu_2 + \gamma_2) )m_{12} N (0)} {m_{21}\beta_1\xi_4 + m_{12}\beta_2\xi_5}\\ &+\frac{ m_{21}(\mu_2 + \gamma_2) (\xi_3 - \gamma_1(\mu_1 + \gamma_1) ) m_{12} N (0)} {m_{21}\beta_1\xi_4 + m_{12}\beta_2\xi_5}.\\%[9pt] \end{split}$

It follows from the expressions of $\mathfrak{\bar{R}}_{vi}$ and $\mathfrak{R}_{vi}$ ($i=1,2$) that

$\label{d} \begin{array}{l} \displaystyle \mathfrak{\bar{R}}_{v1}-\mathfrak{R}_{v1}=\frac{(\mu_2+\gamma_2)(\gamma_2+p_2\mu_2)m_{12}\beta_1} {(\mu_1+\gamma_1)((\mu_2+\gamma_2) \xi_1+(\mu_1+m_{12})\xi_2)}(1-\mathfrak{R}_{v2}),\\ \displaystyle \mathfrak{\bar{R}}_{v2}-\mathfrak{R}_{v2}=\frac{(\mu_1+\gamma_1)(\gamma_1+p_1\mu_1)m_{21}\beta_2} {(\mu_2+\gamma_2)((\mu_1+\gamma_1) \xi_1+(\mu_2+m_{21})\xi_3)}(1-\mathfrak{R}_{v1}). \end{array}$

Thus, if $\mathfrak{R}_v>1$, the inequalities $\mathfrak{R}_{v1}>\mathfrak{\bar{R}}_{v1}>1$ and $\mathfrak{R}_{v2}>\mathfrak{\bar{R}}_{v2}>1$ hold. Therefore, system (3) coexists at most four equilibriums if $\mathfrak{\bar{R}}_{v1}>1$ and $\mathfrak{\bar{R}}_{v2}>1$.

Using $N(0)=N_1+N_2$, we can reduce system (3) with $n=2$ as follows

$\left\{ \begin{array}{ll} \displaystyle \frac{dS_1}{dt}=(1-p_1)\mu_1 N_1-\beta_1\frac{I_1}{N_1} S_1-\mu_1 S_1+m_{21}S_2-m_{12} S_1,\\ \displaystyle \frac{dS_2}{dt}=(1-p_2)\mu_2 (N(0)-N_1)-\beta_2\frac{I_2}{N(0)-N_1} S_2-\mu_2 S_1+m_{12}S_1-m_{21} S_2,\\ \displaystyle \frac{dI_1}{dt}=\beta_1\frac{I_1}{N_1} S_1-(\mu_1+\gamma_1) I_1, \displaystyle \frac{dI_2}{dt}=\beta_2\frac{I_2}{N(0)-N_1} S_2-(\mu_2+\gamma_2) I_2,\\ \displaystyle \frac{dN_1}{dt}=m_{21}(N(0)-N_1-I_2)-m_{12}(N_1-I_1), \end{array}\label{model:n2}\right.$

(25)

which can be used to study the local behavior of system (3) near the boundary equilibria. By considering the linear system for (25), we have the following theorems.

Theorem 4.3. If $\mathfrak{R}_{v1}>1$ and $\mathfrak{\bar{R}}_{v2}<1$, then $\hat{E}$ is locally asymptotically stable, whereas if $\mathfrak{\bar{R}}_{v2}>1$, then $\hat{E}$ is unstable.

Proof. Evaluating system (25) at boundary equilibrium $\hat{E}$, we have the following Jacobian matrix

$\begin{split} J\mspace{-2mu}(\hat{E})\mspace{-5mu}=\mspace{-5mu}\left(\begin{array}{cccccc} -\mspace{-4mu}2\mu_1\mspace{-5mu} -\mspace{-5mu}\gamma_1\mspace{-5mu}-\mspace{-5mu}m_{12}\mspace{-3mu}&m_{21}& \mspace{-3mu}-\mspace{-3mu}\beta_1\mspace{-3mu}\frac{\hat{S}_1}{\hat{N}_1}&0& (\mspace{-3mu}1\mspace{-5mu}-\mspace{-5mu}p_1)\mu_1\mspace{-5mu}+\mspace{-5mu} (\mu_1\mspace{-5mu}+\mspace{-5mu}\gamma_1)\mspace{-3mu}\frac{\hat{I}_1}{\hat{N}_1}\\ m_{12}& \mspace{-4mu}-\mspace{-4mu}\mu_2\mspace{-5mu}-\mspace{-5mu}m_{21}&0&-\beta_2\frac{\hat{S}_2}{\hat{N}_2}&-(1-p_2)\mu_2\\ \beta_1\frac{\hat{I}_1}{\hat{N}_1}&0&0&0&-(\mu_1+\gamma_1)\frac{\hat{I}_1}{\hat{N}_1}\\ 0&0&0&\beta_2\mspace{-4mu}\frac{\hat{S}_2}{\hat{N}_2}\mspace{-5mu}-\mspace{-5mu} \mu_2\mspace{-5mu}-\mspace{-5mu}\gamma_2&0\\ 0&0&m_{12}&-m_{21}&-(m_{12}+m_{21}) \end{array}\right). \end{split}$

It is clearly that one of the eigenvalue of $J(\hat{E})$ is

$\lambda_1=\beta_2\frac{\hat{S}_2}{\hat{N}_2}-(\mu_2+\gamma_2)=(\mu_2+\gamma_2)(\mathfrak{\bar{R}}_{v2}-1),$

where $\mathfrak{\bar{R}}_{v2}$ given in (24). After complexity calculation, the other eigenvalues of $J(\hat{E})$ are the solutions of quartic equation

$\lambda^4+\alpha_3\lambda_3+\alpha_2\lambda_2+\alpha_1\lambda_1+\alpha_0=0,$

(26)

where

$\begin{array}{rl} \alpha_0=&m_{12}(\mu_1\mu_2+\mu_1m_{21}+\mu_2m_{12}(\mu_1+\gamma_1)\hat{N}_1^3\hat{I}_1+m_{12}m_{21}(1-p_2)\mu_2\hat{N}_1^4\\ &+(\mu_1\gamma_1)^2(\mu_2+m_{21})m_{21}\hat{N}_1^4+(\gamma_1+p_1\mu_1)(\mu_2+m_{21})m_{12}\hat{N}_1^4>0,\\ \alpha_1=&(\mu_1+\gamma_1)(\mu_1+\mu_2+m_{12}+m_{21})m_{12}\hat{N}_1^2\hat{I}_1+((2\mu_1+\gamma_1)m_{21}^2+\mu_2m_{12}^2)\hat{N}_1^2\\ &+(\mu_2(\mu_1+\gamma_1)^2+(m_{12}+2m_{21})\gamma_1^2+2\mu_1(\mu_1+\mu_2)+(4\mu_1+\mu_2)\gamma_1)\hat{N}_1^2\\ &+(m_{21}(2\mu_1+\mu_2+\gamma_1)+\mu_1(\gamma_1+p_1\mu_1+p_1\gamma_1)+\mu_2(\gamma_1+2\mu_1))\hat{N}_1^2>0,\\ \alpha_2=&m_{12}(\mu_1+\gamma_1)\hat{N}_1\hat{I}_1+(m_{12}^2+m_{21}^2+(\mu_1+\gamma_1)^2+\mu_2\gamma_1+2\mu_1\mu_2)\hat{N}_1^2\\ &+(m_{21}(2\gamma_1+4\mu_1+\mu_2) +m_{12}(2m_{21}+\gamma_1+2\mu_1+2\mu_2))\hat{N}_1^2>0,\\ \alpha_3=&(2\mu_1+\mu_2+\gamma_1+2m_{12}+2m_{21})\hat{N}_1>0. \end{array}$

By Routh-Hurwitz theorem, (26) has roots with negative real parts only requires that $\alpha_1>0,\alpha_3>0$ and $\alpha_1\alpha_2-\alpha_0\alpha_3>0$. Then we only need to prove the last inequality. Direct computation gives that

$\begin{split} \alpha_1\alpha_2\mspace{-4mu}-\mspace{-4mu}\alpha_0\alpha_3 &= m_{12}^2(\mu_1+\gamma_2)^2(\mu_1+\mu_2+m_{12}+m_{21})\hat{I}_1^2\hat{N}_1^3 +(\mu_1+\gamma_1)(m_{12}^3+m_{21}^3\\ &+\mu_1\gamma_1^2 +2\mu_1^2+\gamma_1+\mu_1^3+2\mu_2\gamma_1^2+4\mu_1\mu_2\gamma_1+2\mu_1^2\mu_2 +\mu_2^2\gamma_1+\mu_1\mu_2^2\\ &+\mspace{-3mu}(3\mu_1\mspace{-3mu}+\mspace{-3mu}2\mu_2\mspace{-3mu}+\mspace{-3mu}\gamma_1\mspace{-3mu}+\mspace{-3mu}2m_{21})m_{12}^2\mspace{-3mu}+\mspace{-3mu}(5\mu_1+2\mu_2+3\gamma_1)m_{21}^2 +(5\mu_1^2 \mspace{-3mu}+\mspace{-3mu}7\mu_1\gamma_1\\ &+\mspace{-3mu}2\gamma_1^2+4\mu_2\gamma_1+6\mu_1\mu_2+\mu_2^2)m_{21}^2+(3\mu_1^2+2\gamma_1^2+\mu_2^2+3m_{21}^2+p_1\mu_1^2\\ &+\mspace{-3mu}4\mu_1\mu_2 \mspace{-3mu}+\mspace{-3mu}4(2\mu_1\mspace{-3mu}+\mspace{-3mu}\mu_2\mspace{-3mu}+\mspace{-3mu}\gamma_1)m_{21}+(p_1+4)\mu_1\gamma_1+2\mu_2\gamma_1)m_{12})m_{12}\hat{N}_1^4\hat{I}_1\\ &+\mspace{-3mu}( m_{21}^4(2\mu_1\mspace{-3mu}+\mspace{-3mu}\gamma_1)\mspace{-3mu}+\mspace{-3mu}\mu_2m_{12}^4\mspace{-3mu} +\mspace{-3mu}2m_{21}^3(2\mu_1\mspace{-3mu}+\mspace{-3mu}\gamma_1)(2\mu_1\mspace{-3mu}+\mspace{-3mu} \gamma_1\mspace{-3mu}+\mspace{-3mu}\mu_2) m_{21}^2(2\mu_1\mspace{-3mu}+\mspace{-3mu}\gamma_1) \end{split}$

$\begin{split}&\times\mspace{-3mu}(4(\mu_1\mspace{-3mu}+\mspace{-3mu}\gamma_1)^2\mspace{-3mu}+\mspace{-3mu} 3\mu_2\gamma_1\mspace{-3mu}+\mspace{-3mu}6\mu_1\mu_2\mspace{-3mu}+\mspace{-3mu}\mu_2^2) \mspace{-3mu}+\mspace{-3mu}(\mu_1+\gamma_1)^2\mu_2(\gamma_1^2+(2\mu_1+\mu_2)\gamma_1\\ &+\mspace{-3mu}(\mu_1\mspace{-3mu}+\mspace{-3mu}2\mu_2)\mu_1) \mspace{-3mu}+\mspace{-3mu}m_{12}^3(\gamma_1^2\mspace{-3mu}+\mspace{-3mu}p_1\mu_1^2 \mspace{-3mu}+\mspace{-3mu}4\mu_1\mu_2+2\mu_2^2+(\mu_1+p_1\mu_1+2\mu_2)\gamma_1\\ &+\mspace{-3mu}(2\mu_1\mspace{-3mu}+\mspace{-3mu}\gamma_1\mspace{-3mu}+\mspace{-3mu}3\mu_2)m_{21}) \mspace{-3mu}+\mspace{-3mu}m_{21}(2\gamma_1^4\mspace{-3mu}+\mspace{-3mu}4\gamma_1^3(2\mu_1 \mspace{-3mu}+\mspace{-3mu}\mu_2)\mspace{-3mu}+\mspace{-3mu}4\mu_1\gamma_1(2\mu_1^2\mspace{-3mu}+\mspace{-3mu}5\mu_1\mu_2\\ &+\mu_2^2)+\gamma_1^2 (12\mu_1+16\mu_1\mu_2+\mu_2^2)+2\mu_1^2(\mu_1^2+4\mu_1\mu_2+2\mu_2^2)) m_{12}^2(\gamma_1^3\\ &+3m_{21}^2(2\mu_1+\mu_2+\gamma_1)+\gamma_1^2(3\mu_2+(p_1+3)\mu_1)+\gamma_1(8\mu_1\mu_2+2\mu_2^2\\ &+\mspace{-4mu}(2\mspace{-4mu}+\mspace{-4mu}3p_1)\mu_1^2) \mspace{-3mu}+\mspace{-3mu}2\mu_1(p_1\mu^2\mspace{-4mu}+\mspace{-4mu} 3\mu_2(\mu_1\mspace{-4mu}+\mspace{-4mu}\mspace{-3mu}\mu_2))\mspace{-3mu}+\mspace{-3mu}m_{21}(3\gamma_1^2 \mspace{-3mu}+\mspace{-3mu}6\mu_1^2\mspace{-3mu}+\mspace{-3mu}2(p_2\mspace{-3mu}+\mspace{-3mu}7)\mu_1\mu_2\\ &+\mspace{-3mu}3\mu_2^2\mspace{-3mu}+\mspace{-3mu}2\gamma_1(4\mu_1\mspace{-3mu}+\mspace{-3mu} (p_2\mspace{-3mu}+\mspace{-3mu}3)\mu_2))) (\gamma_1^4\mspace{-3mu}+\mspace{-3mu}m_{21}(3\gamma_1\mspace{-3mu}+\mspace{-3mu}6\mu_1 \mspace{-3mu}+\mspace{-3mu}\mu_2)\mspace{-3mu}+\mspace{-3mu}\gamma_1^3((p_1\mspace{-3mu}+\mspace{-3mu}3)\mu_1\\ &+\mspace{-3mu}2\mu_2)\mspace{-3mu}+\mspace{-3mu}\gamma_1^2(3(p_1\mspace{-3mu}+\mspace{-3mu}1)\mu_1^2+8\mu_1\mu_2+2\mu_2^2) +\mu_1\gamma_1((1+3p_1)\mu_1^2+10\mu_1\mu_2\\ &+\mspace{-3mu}(7\mspace{-3mu}-\mspace{-3mu}p_1)\mu_2^2)\mspace{-3mu}+\mspace{-3mu} \mu_1^2(p_1\mu_1^2\mspace{-3mu}+\mspace{-3mu}4\mu_1\mu_2\mspace{-3mu}+\mspace{-3mu}(6-p_2)\mu_2^2) +m_{21}^2(4\gamma_1^2+(14-p_1)\mu_1^2\\ &+\mspace{-3mu}2(p_2\mspace{-3mu}+\mspace{-3mu}7)\mu_1\mu_2\mspace{-3mu}+\mspace{-3mu}\mu_2^2\mspace{-3mu}+\mspace{-3mu} \gamma_1((15\mspace{-3mu}-\mspace{-3mu}p_1)\mu_1\mspace{-3mu}+\mspace{-3mu}2(p_2\mspace{-3mu}+\mspace{-3mu}3)\mu_2)) \mspace{-3mu}+\mspace{-3mu}m_{21}(\gamma_1^2((15\mspace{-3mu}+\mspace{-3mu}p_1)\mu_1\\ &+\mspace{-3mu}4\gamma_1^3\mspace{-3mu}+\mspace{-3mu}(6\mspace{-3mu}+\mspace{-3mu}p_2)\mu_2)\mspace{-3mu}+\mspace{-3mu}\gamma_1((17\mspace{-3mu}+\mspace{-3mu}3p_1)\mu_1^2\mspace{-3mu}+\mspace{-3mu}(21-2p_1\mspace{-3mu}+\mspace{-3mu}3p_2)\mu_1\mu_2 \mspace{-3mu}+\mspace{-3mu}(p_2\mspace{-3mu}+\mspace{-3mu}2)\mu_2^2)\\ &+\mspace{-3mu}\mu_1(2(p_1+3)\mu_1^2+(19-2p_1+2p_2)\mu_1\mu_2+(p_2+7)\mu_2^2))) )\hat{N}_1^5>0. \end{split}$

Then all solution of (26) have negative real parts. Therefore, based on the above discussion, we know that if $\mathfrak{R}_{v1}>1$, there exists a boundary equilibrium $\hat{E}$ and it is locally asymptotically stable if $\mathfrak{\bar{R}}_{v2}<1$. This completes the proof.

Similar results hold for boundary equilibrium $\bar{E}$ if $\mathfrak{R}_{v2}>1$.

Theorem 4.4. If $\mathfrak{R}_{v2}>1$ and $\mathfrak{\bar{R}}_{v1}<1$, then $\bar{E}$ is locally asymptotically stable, whereas if $\mathfrak{\bar{R}}_{v1}>1$, then $\bar{E}$ is unstable.

---

5. Numerical simulations

To complement the mathematical analysis carried out in the previous sections, we now investigate some of the numerical properties of system (3). We take the default parameter values as: $N(0)=10,000$, $\beta_1=0.5, \beta_2=0.3, \gamma=1/7, \mu=0.0006$, $d=360$ and the initial condition for system (3) considered as $(S_{1}(0),S_{2}(0),I_{1}(0),I_{2}(0),R_{1}(0),R_{2}(0))=(6489,3489,10,10,1,1)$. In this section, we mainly change migration rate and the vaccination coverage to simulate how migration and vaccination affect the outbreak of an infectious diseases. Since the annual (or quarter) floating population number for a region or a country (e.g. 2015 China Statistical Yearbook) determined the migration proportion $l_{ij}$, it follows from (2) that we can choose $d=360$ (or $d=90$) to calculate the migration rate $m_{ij}$. In this paper, we focus on the migration duration $d=360$, the corresponding results for other migration durations can be similarly obtained.

Time evolution of system (3) in the special case of $n=2$ with multi-initial conditions are presented in Fig. 2. One can observe from the first figure that the trajectories of the two-patch system converge to the disease-free equilibrium $E^0$ when $\mathfrak{R}_{v1}=0.792<1$ and $\mathfrak{R}_{v2}=0.421<1$. This means that the disease disappears in the whole population as proved in Theorem 4.1. If $\mathfrak{R}_{v1}=2.39>1,\ \mathfrak{R}_{v2}=0.85<1$ and $\bar{\mathfrak{R}}_{v1}=2.46>1,\ \bar{\mathfrak{R}}_{v2}=0.395<1$, the trajectories of this two-patch system converge to the boundary equilibrium $\hat{E}=(2046,542, 14, 0,7131,2869)$ as depicted in Fig. 2(b), while if $\mathfrak{R}_{v1}=0.938<1,\ \mathfrak{R}_{v2}=1.24>1$ and $\bar{\mathfrak{R}}_{v1}=0.805<1,\ \bar{\mathfrak{R}}_{v2}=1.265>1$, the trajectories converge to the boundary equilibrium $\bar{E}=(1645,1375, 0, 3,7125,2875)$ as shown in Fig. 2(c). This suggests that the disease persist in one patch and disappear in another patch when just one patch control reproduction number larger than one as proved in Theorem 4.3 and 4.4. If the both the control reproduction number and the invasion threshold are all larger than one, as shown in Fig. 2(d) $\mathfrak{R}_{v1}=2.99>\bar{\mathfrak{R}}_{v1}=2.54>1$ and $\mathfrak{R}_{v2}=1.84>\bar{\mathfrak{R}}_{v2}=1.81>1$, then the trajectories converge to the endemic equilibrium $E^*=(2046,1372, 20, 2,7131,2869)$. Similarly to the model considered in [15,23], the disease persist in both sub-populations if the control reproduction numbers in each patch are all larger than one, to which complement the mathematical analysis carried out in the previous section.

The theoretical and numerical results all show that $\mathfrak{R}_v=1$ acting as a sharp threshold between disease extinction and persist state, to study the influence of vaccination and migration, we can first study how the control reproduction number $\mathfrak{R}_v$ depends on the vaccination rate $p_1,p_2$ and the migration rate $m_{12},m_{21}$ (which determined by $l_{12},l_{21}$). Following Fig. 3, we can observe that if there is no migration between patches, then vaccines always helpful to control the spread of an infectious disease and the optimal vaccination strategy should be vaccine more individuals in patch 1 and relatively lower in patch 2. This is very reasonable, when more people in higher transmission patch (i.e., patch 1) get vaccinated, then the number of people to be infected is smaller, hence the threshold $\mathfrak{R}_v$ will be reduced more. For the specific vaccination rate $(p_1,p_2)=(0.8,0.6)$, if there is no migration between patches then $\mathfrak{R}_v=0.83$. If there exist migration between patches, then more people migrate from patch 2 to patch 1 could make the threshold $\mathfrak{R}_v=0.83$ marked in red curve in Fig. 3(b). Since migration rate $m_{12}$ ($m_{21}$) is an increasing function of variable $l_{12}$ ($l_{21}$), therefore, this figure also implies that the threshold $\mathfrak{R}_v$ increases with the variable $l_{12}$ and decreases with variable $l_{21}$. However, if we fix $(p_1,p_2)=(0.6,0.8)$, it follows from Fig. 3(c) that threshold $\mathfrak{R}_v$ decreases with the variable $l_{12}$ and increases with variable $l_{21}$. This suggests that if we want to control the threshold $\mathfrak{R}_v$ less than one, we should consider not only the effectiveness of vaccines but also the impact of migration.

To explore the effect of vaccination and migration, we also compare the second peak size and second peak time with various vaccination coverage and migration rate. For the case of $p_1=0.2, p_2=0.1$, although the endemic equilibrium $E^*$ is stable due to $\mathfrak{R}_{v1}=2.79>1$ and $\mathfrak{R}_{v1}=1.88>1$, we can obtain that the second peak size decrease form 330 at 1112 day to 294 at 1450 day for patch 1, decrease from 119 at 1645 day decrease to 92 at 1906 day for patch 2, and decrease from 342 at 1110 day to 291 at 1490 day compared to the case of $p_1=p_2=0$ as shown in Fig.4. This figure also depicts that the disease will disappear at the first 2000 days if the vaccination rate add up to $p_1=0.65,p_2=0.7$, at this time, the threshold $\mathfrak{R}_{v1}=1.22>1$ and $\mathfrak{R}_{v1}=0.627<1$ and the boundary equilibrium is stable. This tells us that vaccine plays a critical role in reducing the second peak size and delay the second peak time, since the control reproduction number is decreasing with vaccination rate $p_1,p_2$. If there exist migration between patches in absence vaccination, then the control reproduction number $\mathfrak{R}_{v1}=2.09>1$ and $\mathfrak{R}_{v2}=3.48>1$, that is, the disease persist in both patches. In this case, we can observe from Fig. 5 that higher migration rate in patch 1 is beneficial to individuals in both patches and higher migration rate in patch 2 will magnify the second peak size of patch 1 and the entire population. This implies that, if we want to control the outbreak size during the first 2000 days, the scale of migration between patches should be considered. These figures all show that the second peaks happened after three or four years later, this may be because the decay of the effectiveness of the vaccination or individuals migration change the spatial structure of sub-populations.

We also compared the residual values of the first peak size to investigate the impact of vaccination and migration, which shown in the histogram 6. The results show that migration can reduce the first peak size for each patches and the entire population as long as the migration rate $m_{12}$ is less than or equal to migration rate $m_{21}$. Such as, for the entire population, the residual of first peak sizes almost equal to -8 in the case of $l_{12}=l_{21}=0.2$ and equal to -29 in the case of $l_{12}=0.3, l_{21}=0.1$, while the residual of first peak sizes almost equal to 24 in the case of $l_{12}=0.04,l_{21}=0.36$ as shown in histogram 6(c).

---

6. Conclusion

In this paper, we proposed a multi-patch SIR model with vaccination to study the influence of vaccination coverage and human mobility on disease transmission. Our theoretical results show that the control reproduction number $\mathfrak{R}_v$ is a threshold parameter of the disease dynamics. It founds that system only exists a disease-free equilibrium $E^0$ if $\mathfrak{R}_v<1$ and is locally asymptotically stable. Particularly, in the case of $n=2$, boundary equilibrium $\hat{E}$ ($\bar{E}$) exists if $\mathfrak{R}_{v1}>1$ ($\mathfrak{R}_{v2}>1$) and it is globally stability if $\mathfrak{\bar{R}}_{v2}<1$ ($\mathfrak{\bar{R}}_{v1}<1$), and endemic equilibrium $E^*$ exists only when $\mathfrak{\bar{R}}_{v1}>1$ and $\mathfrak{\bar{R}}_{v2}>1$. The simulation results show that, for parameter values considered, vaccines always can shrink the outbreak of an infectious while mobility restriction (change the migration rate) dose not necessarily always have a positive impact on the overall spread of disease. An increase of migration rate from one patch to the other sometimes may prevent the outbreak of infectious in one patch while intensify the disease spread in other patch (see Figs. 5 and 6).

In our model, we assume that the infective do not move between patches, corresponding to either a very severe disease so that infective are not able to move or move is forbidden in order to control outbreak of disease. In the further, we can generalize the current model with infective move between patches.

---

Acknowledgments

This project has been partially supported by grants from National Natural Science Foundation of China (Nos. 11671206,11271190) and Scientific Research Innovation Project of Jiangsu Province (No. KYZZ15_0130). We also thank two anonymous referees for their helpful comments.

---