Transmission dynamics of brucellosis with patch model: Shanxi and Hebei Provinces as cases

1.   Introduction

Brucellosis is one of the common zoonotic diseases, and it mainly infects livestock and has been reported in various countries ^[1]. The main reason people suffer from brucellosis is eating food contaminated with Brucella or contacting the secretions and excreta of sick animals; people don't infect each other ^[2]. There are many kinds of Brucella, including Brucella of cattle, sheep, pigs, dogs and mice, of which the main infectious sources are sick sheep and cattle ^[3,4]. With the large-scale management of the sheep industry and the frequent trading of sheep products among regions, the liquidity of sick sheep has increased. In the past 30 years, brucellosis epidemic areas have gradually shifted from pastoral areas (i.e., Inner Mongolia, Xinjiang, Tibet, Qinghai and Ningxia) to grassland and agricultural areas (i.e., Shanxi, Liaoning, Hebei, Shandong and Jilin). Especially since 2004, the affected areas have expanded from the north to the south of China ^[5,6]. There is evidence that the epidemic of brucellosis in the south is caused by infected animals imported from other regions ^[7]. Therefore, the development of animal husbandry and the immigration of sheep are the most commonly accepted explanations for the prevalence of human brucellosis.

Mathematical modeling is a good tool to study diseases and give measures for disease control. For example, mathematical models have played a key role in predicting and controlling the disease for the novel coronavirus currently present globally, Ma et al. ^[8] studied the effect of mask use on controlling the spread of COVID-19 by constructing a model, and Asamoah et al. ^[9] studied the optimal control and cost effectiveness of COVID-19 through modeling. Some mathematical models have been used in the study of brucellosis. For example, Li et al. ^[2] studied the transmission mechanism of brucellosis in the Hinggan League of Inner Mongolia, and the results showed that banning the mixed feeding of basic ewes and other sheep, vaccination, detection and culling were the effective strategies. Hou and Sun ^[10] established a multi-stage dynamic model of sheep brucellosis transmission, and it was concluded that the birth rate, vaccination rate and culling rate of sheep play an important role in the transmission of brucellosis. From the investigation and comparison of the effects of vaccination and culling strategies, the latter is better than the former. Chen et al. ^[11] studied the spatial distribution of human brucellosis in Shanxi Province and found that the proportion of affected cities and towns increased from 31.5% in 2005 to 82.5% in 2014. These papers only studied the temporal and spatial characteristics of human brucellosis and established a dynamic model according to the pathogenesis and propagative law of brucellosis in different epidemic areas. The immigration of sheep among patches is an important factor causing the spread of brucellosis, but there is little work.

With regard to the research on the spatial transmission of brucellosis, Zhang et al. ^[12] established a multi-patch dynamic model of cattle brucellosis, and they obtained that the dispersal of susceptible populations of each patch and the centralization of infected cattle to patches with large breeding scale are conducive to the control of the disease. However, the basic reproduction number, the uniqueness of the positive equilibrium and the global asymptotic stability have not been further analyzed, and there is no numerical simulation of practical problems. Sick sheep are a major source of human brucellosis. Based on the above research, this paper established a patch model composed of people, sheep and Brucella. The transmission of sheep brucellosis among patches was analyzed qualitatively and quantitatively, and the prevention and control measures were put forward in combination with practical problems.

The article is organized as follows. In Section 2, an n-patch dynamics model was proposed. In Section 3, we gave the mathematical analysis of the model, including the basic reproduction number, disease-free equilibrium and positive equilibrium. In Section 4, taking two patches in Shanxi Province and Hebei Province as examples, the numerical simulation was given. In Section 5, a brief conclusion was made.

2.   Model formulation

The spatial distribution of human brucellosis shows that there is inter-provincial transmission of human brucellosis in China. Therefore, we regarded each province as an isolated patch. Assuming there are n patches, a multi-patch model composed of humans, sheep and Brucella was established. We made some assumptions about the model: that human brucellosis is mainly transmitted by sick sheep, and that environmental brucellosis also causes infection in susceptible people and sheep. There are few reports of human-to-human transmission of brucellosis, so human-to-human and human-to-animal transmission are ignored. Since sheep are the main source of infection, we considered only sheep immigration among patches and not human immigration among patches. In patch i, the population is divided into three classes: $S_{i}^h$, $I_{i}^h$ and $Y_{i}^h$, which are susceptible, acute and chronic brucellosis patients at time t, respectively. Sheep are divided into two classes: $S_i$ and $I_i$, which are the number of susceptible sheep and infected sheep at time t, respectively. During infection, the infected sheep release brucella into the environment, and the number of brucella in the environment is denoted by $W_i$. The flow chart of brucellosis is shown in Figure 1.

Based on the model flow chart, we established a multi-patch dynamics model including sheep, humans and bacteria in the environment:

where $A_{i}^h$ and $d_{ih}$ are the numbers of birth and natural mortality of people per unit time, respectively; $A_{i}^v$ and $d_{iv}$ are the birth number and sale rate of sheep per unit time, respectively; $p_{i}$ is the transfer rate from acute infections to susceptible individuals; $m_{i}$ is the transfer rate from acute infections to chronic infections; $\alpha_{i}$ is the transmission rate of Brucella to susceptible humans; $\beta_{i}$ is the transmission rate of infectious sheep to susceptible humans; $\delta_{i}$ is the transmission rate of infectious sheep to susceptible sheep; $\phi_{i}$ is the transmission rate of Brucella to susceptible sheep; $a_{iv}$ is the disease-related culling rate of infectious sheep; $a_{ij}$ and $b_{ij}$ are the immigration rates of the susceptible sheep and infectious sheep from the jth patch to the ith patch for $i \neq j$; $-a_{ii} = \sum_{i\neq j} a_{ji}$ and $-b_{ii} = \sum_{i\neq j} b_{ji}$ are the emigration rates of the susceptible sheep and infectious sheep; $k_i$ is the amount of Brucella released from infected sheep; and $\sigma$ is the decay rate of Brucella.

Since the last three equations of system (2.1) are independent of the first three equations, we only studied the following system for the dynamic analysis of (2.1):

3.   Dynamic analysis of model

Firstly, we considered the existence and uniqueness of the disease-free equilibrium of (2.2).

Let $E^0 = (S_1^0, S_2^0, ..., S_n^0, 0, ..., 0, 0, ..., 0)$ be the disease-free equilibrium of (2.2); then $S^0 = (S_1^0, S_2^0, ..., S_n^0)$ is the positive equilibrium of the following subsystem:

According to $A_i^v-d_{iv} S_{i}+\sum_{j = 1}^{n} a_{ij} S_{j} = 0$, we define an auxiliary matrix

Then, $M_1S = A^v$, where $S = (S_1, S_2, ..., S_n)^T, A^v = (A_1^v, A_2^v, .., A_n^v)^T$. Note that $M_1$ is an irreducible $M$-matrix (Appendix A ^[13]), and $M_1S = A^v$ has a unique solution $S$. According to corollary 4.3.2 ^[14], $M_1^{-1}$is a positive matrix; then, $S = M_1^{-1}A^v > 0$, so $S^0 = (S_1^0, S_2^0, ..., S_n^0)$ is the only positive equilibrium of (3.1).

Define $s (M)$ as the spectral bound of matrix $M$.

$s(M)$ = max{Re$\lambda$:$\lambda$ is an eigenvalue of $M$}.

$M_1S = A^v$, so $-M_1S = -A^v$. Since $-M_1$ is an irreducible Metzler matrix, by the Perron-Frobenius Theorem ^[15], $s(-M_1)$ is a zero solution of the characteristic polynomial, there is no other eigenvalue $\lambda$ so that $R(\lambda) = s(-M_1)$, and there is only one positive eigenvector for the eigenvalue $s(-M_1)$. Let $V = (v_1, v_2, …, v_n)$ be a positive eigenvector associated with $s(-M_1)$; then, $-M_1V = s(-M_1)V = -A^v$, so $s(-M_1) < 0$, and hence all eigenvalues of $-M_1$ have negative real parts. Since (3.1) is a linear system, $S^0 = (S_1^0, S_2^0, …, S_n^0)$ is globally asymptotically stable on $S\in R_+^n \setminus {0}$. Thus, $E^0 = (S_1^0, S_2^0, ..., S_n^0, 0, 0, ..., 0, 0, 0, ..., 0)$ is the disease-free equilibrium of (2.2).

3.1.   Basic reproduction number

To derive the basic reproduction number $R_0$ for (2.2), we ordered the infected variables first by disease state, then by patch, i.e.,

Applying the method of the next generation matrix ^[13], we can obtain the expression of the basic reproduction number; define

$F = \begin{bmatrix}{} F_1 & F_2 \\ 0 & 0 \end{bmatrix}$, $V = \begin{bmatrix}{} V_1 & 0 \\ V_3&V_4 \end{bmatrix}$, where $F_1 = (\delta_{ij}(\delta_iS_i^0))_{n \times n}$, $F_2 = (\delta_{ij}(\phi_iS_i^0))_{n \times n}$, $V_1 = (\delta_{ij}(d_{iv}+a_{iv})-B)_{n \times n}$, $B = (b_{ij})_{n \times n}$, $V_3 = (\delta_{ij}(-k_i))_{n \times n}$, $V_4 = (\delta_{ij}(\sigma))_{n \times n}$.

$\delta_{ij}$ denotes the Kronecker delta (i.e., 1 when i = j and 0 otherwise). Define $M = F-V$.

then, $R_0 = \rho (FV^{-1}) = \rho (F_1 V_1^{-1}-F_2 V_4^{-1} V_3 V_1^{-1})$, where $\rho (FV^{-1})$ represents the spectral radius of the matrix $FV^{-1}$. In particular, the basic reproduction number of the single patch model (2.2) is $R_{0i} = \frac{\delta_iS_i^0\sigma+k_i\phi_iS_i^0}{(d_{iv}+a_{iv})\sigma}.$

Corollary 3.1. $\min_{1\leq i \leq n}R_{0i}\leq R_0 \leq \max_{1\leq i \leq n} R_{0i}.$

Proof. Through the expression of $R_0$, we gained $F_1 V_1^{-1}-F_2 V_4^{-1} V_3 V_1^{-1} = (F_1 -F_2 V_4^{-1} V_3)V_1^{-1}$.

Let

then, $R_0 = \rho(HV_1^{-1})$.

Since $V_1$ is the $M$-matrix, $V_1^{-1}$ is a nonnegative matrix. We apply Fischer's inequality (Theorems 2.5.4(e) ^[16]) and 3.4 ^[17] to estimate the diagonal entries of matrix $V_1^{-1}$.

For example, let $V_1 = (a_{ij})_{n \times n}$ and $V_1^{-1} = (\alpha_{ij})_{n \times n}$; then, $\frac{1}{a_{ii}}\leq \alpha_{ii}$, $i = 1, 2, ..., n$. Therefore, $0 \leq diag\{\frac{1}{d_{1v}+a_{1v}}, ..., \frac{1}{d_{nv}+a_{nv}}\}\leq diag\{\frac{1}{d_{1v}+a_{1v}-b_{11}}, ..., \frac{1}{d_{nv}+a_{nv}-b_{nn}}\}$ $\leq diag\{\alpha_{11}, ..., \alpha_{nn}\} \leq V_1^{-1}$.

So, $\min_{1\leq i \leq n}R_{0i}\leq R_0$.

Let $G = V_1+B = diag\{d_{1v}+a_{1v}, ..., d_{nv}+a_{nv}\}$; then, $I(V_1G^{-1}) = I \Rightarrow I(GV_1^{-1}) = I$, where $I = (1, 1, ..., 1)_{1 \times n}$, which implies that the spectral radius of $GV_1^{-1}$ is 1, and hence $\rho(V_1^{-1}) = \rho(G^{-1}GV_1^{-1})\leq \rho(G^{-1})\rho(GV_1^{-1}) = \rho(G^{-1})$.

So, $\rho(HV_1^{-1})\leq \rho(H) \rho(V_1^{-1})\leq \rho(H) \rho(G^{-1})\leq \max_{1\leq i \leq n} R_{0i}$.

Based on the above, $\min_{1\leq i \leq n}R_{0i}\leq R_0 \leq \max_{1\leq i \leq n} R_{0i}$.

Lemma 3.1. There hold two equivalences ^[13]:

By Theorem 2 ^[13], the disease-free equilibrium $E^0$ is locally asymptotically stable if $R_0 < 1$ and is unstable if $R_0 > 1$.

Next, we studied the global dynamics of (2.2), the disease-free equilibrium $E^0$ is globally attractive if $R_0 < 1$, and (2.2) has a positive equilibrium if $R_0 > 1$.

Lemma 3.2. Let $k = \max\{k_i:1\leq i\leq n\}, d_v = \min\{d_{iv}:1\leq i\leq n\}, N^* = \frac{A}{d_v}$, where $A = \sum_{i = 1}^{n}A_i^v$. Every forward orbit of $(2.2)$ eventually enters into $G = \{(S, I, W)\in R_+^{3n}:\sum_{i = 1}^{n}(S_i+I_i)\leq N^*, \sum_{i = 1}^{n}W_i \leq \frac{kN^*}{\sigma}\}$, and G is positively invariant for $(2.2)$.

Theorem 3.1. When $R_0 < 1$, the disease-free equilibrium of system $(2.2)$ is globally asymptotically stable in G ^[12].

Theorem 3.2. When $R_{0} > 1$, then $(2.2)$ admits at least one positive equilibrium, and there is a positive constant $\kappa$ such that every solution $\phi_t(\chi_0) = (S(t), I(t), W(t))$ of $(2.2)$ with $\chi_0 = (S_1(0), ..., S_n(0), I_1(0), ..., I_n(0), W_1(0), ..., W_n(0))\in R_{+}^n \times R_{+}^n\setminus \{0\} \times R_{+}^n\setminus \{0\}$ satisfies

Proof. Let

Then, we showed that (2.2) is uniformly persistent with respect to $(X_0, \partial X_0)$.

Clearly, $X$ and $X_0$ are positively invariant, and $\partial X_0$ is relatively closed in $X$. Furthermore, system (2.2) is point dissipative ^[18]. Set $M_\partial = \{(S(0), I(0), W(0)):(S(t), I(t), W(t))$ satisfies (2.2) and $(S(t), I(t), W(t)) \in \partial X_0, \forall t \geq 0\}$.We next show that $M_\partial = \{(S, 0, 0):S\geq 0\}$.

Set $(S(0), I(0), W(0))\in M_\partial$; then, $I(t) = 0, W(t) = 0, \forall t \geq 0$. Suppose not, then there exist an $i_0, 1\leq i_0\leq n$, and a $t_0\geq0$ such that $I_{i0}(t_0) > 0$. Here we only analyzed $I(t)$, because the change of $W(t)$ depends on the change of $I(t)$. If $I_i(t) > 0$, then $W_i(t) > 0$; if $I_i(t) = 0$, then $W_i(t)$ eventually tends to 0.

We partition $\{1, 2, ..., n\}$ into two sets $Q_1$ and $Q_2$ such that

$I_i(t_0) = 0, \forall i\in Q_1$,

$I_i(t_0) > 0, \forall i\in Q_2$.

Defined by $M_\partial$, $Q_1$ is a non-empty set. Since $I_{i0} > 0$, $Q_2$ is non-empty. For any $j\in Q_1$, we have $I_j'(t_0) = \delta _{j} S_{j}(t_0) I_{j}(t_0)+\phi _{j} W_{j}(t_0) S_{j}(t_0)-d_{v} I_{j}(t_0)-a_{v} I_{j}(t_0)+$$\sum_{i = 1}^{n} b_{ji} I_{i}(t_0) = \phi _{j} W_{j}(t_0) S_{j}(t_0)+\sum_{i = 1}^{n} b_{ji} I_{i}(t_0)\geq b_{ji_0}I_{i0}(t_0) > 0$.

Then exist an $\varepsilon_0 > 0$, when $t_0 < t < t_0+\varepsilon_0$, there is $I_j(t) > 0, j\in Q_1$, we can restrict $\varepsilon_0 > 0$ small enough such that $t_0 < t < t_0+\varepsilon_0, I_i(t) > 0, i\in Q_2$. This shows that when $t_0 < t < t_0+\varepsilon_0$, $(S(t), I(t), W(t))\notin \partial X_0$, which contradicts $(S(0), I(0), W(0))\in M_\partial$, so $M_\partial = \{(S, 0, 0):S\geq 0\}$.

Obviously, $M_\partial$ has only one equilibrium $E^0$. We chose $\eta > 0$ small enough such that $s(M-M_\eta) > 0$, where $M_\eta = \begin{bmatrix} A_\eta&B_\eta \\ 0 & 0 \end{bmatrix}$, $A_\eta = \delta_{ij}(\delta_i\eta)_{n\times n}, B_\eta = \delta_{ij}(\phi_i \eta)_{n\times n}$. Consider the perturbed system of (3.1)

First, with regard to our previous analysis of the system (3.1), restrict $\varepsilon_1 > 0$ small enough so that the system (3.3) has a unique equilibrium point $S^*(\varepsilon_1)$ and is globally asymptotically stable, and $S^*(\varepsilon_1)$ is continuous in $\varepsilon_1$. Therefore, we can further restrict $\varepsilon_1 > 0$ small enough such that $S^*(\varepsilon_1) > S^0-\eta$.

Let us consider an arbitrary positive solution $(S(t), I(t), W(t))$ of (2.2); then, $\underset{t \to \infty }{\lim} \sup \underset{i}{\max}\{I_i(t)\}$ $> \varepsilon_1$. Suppose there is a $T > 0$ such that $I_i(t) \leq\varepsilon_1, i = 1, 2, ..., n,$ for $t\geq T$; then for $t \geq T$, we have

Since the equilibrium $S^*(\varepsilon_1)$ of system (3.3) is globally asymptotically stable, and $S^*(\varepsilon_1) > S^0-\eta$, there is $T_1 > 0$ such that $S(t)\geq S^0-\eta$ for $t > T_1+T$. Therefore, when $t > T_1+T$,

Since the matrix $(M-M_\eta)$ has a positive eigenvalue $s(M-M_\eta)$ with a positive eigenvector, according to the comparison principle ^[19], $\underset{t \to \infty }{\lim} I_i(t) = \infty$; then, $\underset{t \to \infty }{\lim}W_i(t) = \infty, i = 1, 2, ..., n$, which leads to a contradiction.

For the system (3.1), we noted that $S^0$ is globally asymptotically stable. From the above, we can see that $E^0$ is an isolated invariant set in $X$, $W^s(E^0)\cap X_0 = \varnothing$. Clearly, each orbit in $M_\partial$ converges to $E^0$, and $E^0$ is acyclic in $M_\partial$. According to the theorem 4.6 ^[20], the system (2.2) is uniformly persistent with respect to $(X_0, \partial X_0)$. By the theorem 2.4 ^[21], (2.2) has an equilibrium $E^* = (S_1^*, ..., S_n^*, I_1^*, ..., I_n^*, W_1^*, ..., W_n^*)$ $\in X_0, S^*\in R_+^n, I^*\in int(R_+^n), W^*\in int(R_+^n)$. Where $S^*\in R_+^n \setminus \{0\}$, supposed that $S^* = 0$, from the second equation of system (2.2), we can get $0 = -\sum_{i = 1}^{n}(d_{iv}+a_{iv})I_i^*$, since $d_{iv}+a_{iv}\neq 0$; then, $I_i^* = 0, i = 1, 2, ..., n$, a contradiction. Through the first equation of system (2.2) and the irreducibility of matrix $(a_{ij})_{n\times n}$, $S^*\in int(R_+^n), \forall t > 0$; then, $(S^*, I^*, W^*)$ is the positive equilibrium of system (2.2).

3.2.   Uniqueness and uniform persistence of endemic equilibrium

We restricted the system (2.2) by assuming that the disease-related culling rate is 0, and the immigration rate of susceptible sheep and infected sheep is the same, that is, $a_{ij} = b_{ij}, i = 1, 2, ..., n$; then, the model of system (2.2) becomes

Add the first two equations of the system (3.5):

By the conclusion of system (3.1), system (3.6) has a unique positive equilibrium $N^* = (N_1^*, N_2^*, ..., N_n^*)$, and system (3.5) has the following limit system:

Lemma 3.3. For system $(3.7)$, the set $G_1 = \{(I, W)\in R_+^{2n}:I_i \leq N_i, W_i\leq \frac{k_iN_i}{\sigma}, i = 1, 2, ..., n \}$ is positively invariant.

Theorem 3.3. When $R_0 > 1$, the system $(3.7)$ admits a unique endemic equilibrium $\overline{E}^* = \{I_1^*, ..., I_n^*, W_1^*, ..., W_n^*\}$, which is globally asymptotically stable with respect to $(I(0), W(0))\in G_1$.

Proof. Through the definition on the right side of the system (3.7), $F\colon G_1\to G _1$. Obviously, $F$ is continuously differentiable and is cooperative on $G_1$, and $DF(I, W)$ is irreducible for every $(I, W)\in G_1$. $F(0) = 0$, and $F_i(I, W)\geq 0$ for all $(I, W)\in G_1$ with $I_i = 0, W_i = 0, i = 1, 2, ..., n$.

For $\forall \alpha \in(0, 1)$ and $(I, W) = (I_1, ..., I_n, W_1, ..., W_n)\in G_1$, we have

$\alpha \delta _{i} (N_i^*-\alpha I_i) I_{i}+\alpha \phi _{i} W_{i} (N_i^*-\alpha I_i)-\alpha d_{iv} I_{i}+$$\alpha \sum_{j = 1}^{n} b_{ij} I_{j} > \alpha[\delta _{i} (N_i^*-I_i) I_{i}+\phi _{i} W_{i} (N_i^*-I_i)-d_{iv} I_{i}+$$\sum_{j = 1}^{n} b_{ij} I_{j}], i = 1, 2, ..., n.$

$k_i \alpha I_i-\sigma \alpha W_i = \alpha (k_iI_i-\sigma W_i), i = 1, 2, ..., n.$

Therefore, $F$ is strictly sublinear on $G_1$ ^[22].

Let $s(DF(0)) = \overline{M} = \begin{bmatrix} F_1-V_1 & F_2 \\ -V_3 & -V_4 \end{bmatrix}$, where $F_1 = (\delta_{ij}(\delta_i N_i^*))_{n \times n}$, $F_2 = (\delta_{ij}(\phi_iN_i^*))_{n \times n}$, $V_1 = (\delta_{ij}(d_{iv}-B)_{n \times n}$, $B = (b_{ij})_{n \times n}$, $V_3 = (\delta_{ij}(-k_i))_{n \times n}$, $V_4 = (\delta_{ij}(\sigma))_{n \times n}$. Then, $\overline{M}$ is an irreducible Metzler matrix. According to the Perron-Frobenius theorem, $s(\overline{M})$ is an eigenvalue with a positive eigenvector, and let its positive eigenvector be $x = (x_1, ..., x_n, x_{n+1}, ..., x_{2n})$, so $\overline{M}x = s(\overline{M})x = (\delta_1 x_1^2+\phi_1 x_1x_{n+1}, ..., \delta_n x_n^2+\phi_n x_nx_{2n}, 0, ..., 0)^T$; then, $s(\overline{M}) > 0$. From lemma $3.1$ and corollary 3.2 ^[22], the system (3.7) has a unique positive equilibrium $\overline{E}^* = \{I_1^*, ..., I_n^*, W_1^*, ..., W_n^*\}$.

Theorem 3.4. When $R_0 > 1$, the system $(3.5)$ has a unique endemic equilibrium $E^* = (S_1^*, ..., S_n^*, I_1^*, ..., I_n^*, W_1^*, ..., W_n^*)$, which is globally asymptotically stable with respect to $(S(0), I(0),$ $W(0))\in G$.

Proof. Let $\psi (t)$ be the corresponding flow of system (3.7). According to the strong monotonicity of $\psi (t)$, $S_i^* = N_i^*-I_i^* > 0, i = 1, 2, ..., n$, so the system (3.5) has a unique positive equilibrium $E^* = (S_1^*, ..., S_n^*, I_1^*, ..., I_n^*, W_1^*, ..., W_n^*)$. Next, we prove the globally asymptotic stability of the positive equilibrium $E^*$.

Let $\phi(t)\colon\ R_+^{3n}\to R_+^{3n}\ $ be the solution semiflow of system (3.5), and let $\omega$ be the $\omega$ limit set of $\phi(S(0), I(0),$ $W(0))\in G$. By lemma $2.1'$ ^[21], $\omega$ is an internal chain transitive set for $\phi(t)$. Clearly, when $R_0 > 1$, the system (3.5) has only two equilibria, $E^0$ and $E^*$, by $S^0 = (S_1^0, S_2^0, ..., S_n^0)$ is globally asymptotically stable on $R_+^n\setminus \{0\}$ and theorem $3.1$, it is easy to know that $\phi(t)$ satisfies theorem 1.2.2 ^[23]. Thus, $\omega$ is $E^0$ or $E^*$. Next, we show that $\omega = \{E^*\}$.

Let's assume $\omega = \{E^0\}$; then, $\underset{t \to \infty }{\lim}S_i(t) = S_i^0, \underset{t \to \infty }{\lim}I_i(t) = 0, \underset{t \to \infty }{\lim}W_i(t) = 0, (i = 1, 2, ..., n)$.

Since $s(M) > 0$, we can choose a small enough $\epsilon > 0$ so that $s(M-M_\epsilon) > 0$, where $M_\epsilon = \begin{bmatrix} A_{01}&B_{01} \\ 0 & 0 \end{bmatrix}, A_{01} = (\delta_{ij}(\delta_i \epsilon))_{n \times n}, B_{01} = (\delta_{ij}(\phi_i \epsilon))_{n \times n}$. It follows that there exists a $T$ such that $S_i(t) > S_i^0-\epsilon$ for $t > T$; then, $\frac{dI_i}{dt} > \delta _{i} (S_i^0-\epsilon) I_{i}+\phi _{i} W_{i} (S_i^0-\epsilon)-d_{v} I_{i}+\sum_{j = 1}^{n} b_{ij} I_{j}$. Let $v = (v_1, v_2, ..., v_n, v_{n+1}, ..., v_{2n})$ be the positive eigenvector associated with $s(M-M_\epsilon)$, and choose a small enough $\alpha$ to satisfy $(I, W) = (I_1, ..., I_n, W_1, ..., W_n)\geq \alpha v$. By the comparison theorem, $(I, W)\geq \alpha v e^{s(M-M_\epsilon)(t-T)}$, and thus $\underset{t \to \infty }{\lim}I_i(t) = \infty, \underset{t \to \infty }{\lim}W_i(t) = \infty (i = 1, 2, ..., n)$, a contradiction. Therefore, $E^*$ is the only endemic equilibrium and is globally asymptotically stable.

4.   Numerical simulation

Assuming $n = 2$, the effect of sheep immigration on the transmission of brucellosis was studied by numerical simulation. Let $a_{21} = -a_{11}, a_{12} = -a_{22}, b_{21} = -b_{11}, b_{12} = -b_{22}$, and hence, system (2.1) reduces to

For system (4.1), the disease-free equilibrium is $P^0 = (S_{1h}^0, 0, 0, S_1^0, 0, 0, S_{2h}^0, 0, 0, S_2^0, 0, 0)$, where $S_{1h}^0 = \frac{A_1^h}{d_{1h}}$, $S_{2h}^0 = \frac{A_2^h}{d_{2h}}$, $S_1^0 = \frac{A_1^v a_{12}+A_2^v a_{12}+A_1^v d_{2v}}{a_{12}d_{1v}+a_{21}d_{2v}+d_{1v}d_{2v}}$, $S_2^0 = \frac{A_1^v a_{21}+A_2^v a_{21}+A_2^v d_{1v}}{a_{12}d_{1v}+a_{21}d_{2v}+d_{1v}d_{2v}}$. Then, the basic reproduction number of the two-patch submodels is calculated as

$R_0 = \frac{1}{2A}(C+D+\frac{\phi_{1} S_1^0B+E}{\sigma})+\frac{1}{2A}\sqrt{(C-D+\frac{\phi_{1} S_1^0B-E}{\sigma})^2+4b_{12}b_{21}\delta_{1}S_1^0\delta_{2}S_2^0-4b_{12}b_{21}k_2\phi_{1}S_1^0\phi_{2}S_2^0}$, where $A = (d_{1v}+a_{1v}+b_{21})(d_{2v}+a_{2v}+b_{12})-b_{12}b_{21},$ $B = (d_{2v}+a_{2v}+b_{12})k_1,$ $C = (d_{2v}+a_{2v}+b_{12})\delta_1S_1^0, D = (d_{1v}+a_{1v}+b_{21})\delta_2S_2^0, E = (d_{1v}+a_{1v}+b_{21})k_2\phi_2S_2^0$.

(A) According to The China Statistical Yearbook, $A_{1h} = 201,164, d_{1h} = 0.0056, A_{2h} = 467,718, d_{2h} = 0.0065$.

(B) According to The China Animal Husbandry Statistical Yearbook, in the past eight years, the average stock in Shanxi Province was $8,754,222$, the average sale rate was $d_{1v} = 0.51,$ and the birth number of sheep per unit time was $A_{1v} = 8,754,222\times0.51 = 4,464,653.$ The average stock in Hebei Province was 14,129,285, the average sale rate was $d_{2v} = 1.44,$ and the birth number of sheep per unit time was $A_{2v} = 14,129,285\times1.44 = 20,346,171.$

4.1.   Parameter estimation

Shanxi and Hebei provinces are adjacent provinces with high incidence rates and similar time series. We have chosen the two provinces as the two patches in model (4.1). The rationality of the model is verified by the least squares method using MATLAB. C1 (t) and C2 (t) are defined as the theoretical cumulative numbers of Shanxi Province and Hebei Province in the t year, that is, the solutions of models $\frac{dC1}{dt} = \beta_{1}S_{1}^h I_1+\alpha_{1} S_{1}^h W_{1}$ and $\frac{dC2}{dt} = \beta_{2}S_{2}^h I_2+\alpha_{2} S_{2}^h W_{2}$. By fitting the model solution with the cumulative number of human brucellosis reported from 2010 to 2018, the values of parameters $\alpha_{1}$, $\alpha_{2}$, $\beta_{1}$, $\beta_{2}$, $\delta_{1}$, $\delta_{2}$, $\phi_{1}$ and $\phi_{2}$ are estimated to obtain $\alpha_{1} = \alpha_{2} = 1.9e^{-13}$, $\beta_{1} = \beta_{2} = 6.64e^{-09}$, $\delta_{1} = 5.64e^{-08}$, $\delta_{2} = 5.648e^{-08}$, $\phi_{1} = 1e^{-08}$ and $\phi_{2} = 3e^{-09}$. The data used to simulate the model are from article ^[4], and the model parameter values are from Table 1. Figure 2(a), (b) shows the numerical simulation and 95% confidence interval of the cumulative number of human brucellosis cases in Shanxi and Hebei provinces from 2010 to 2018, respectively. As can be seen from Figure 2, the solution of the model is consistent with the reported data. In addition, we assumed that the number of newly infected people in the two regions obeys the Poisson distribution at time t, and the cumulative cases are C1 (t) and C2 (t), t = 2010, ..., 2018. If 1000 samples are taken from the Poisson distribution, we have 1000 groups of data samples and make least squares fitting for each group of data. At the same time, 1000 groups of corresponding estimated parameter values can be obtained. We can assume that the parameters obey the normal distribution, and then we can obtain the 95% confidence interval of the parameters (Table 1). According to the parameter values in the simulation and the basic reproduction number formula in the two patch models, $R_0 = 0.5457$, $R_ {01} = 0.7497$, and $R_ {02} = 0.5022$, which means that the disease will disappear in Shanxi and Hebei Provinces. This shows that the two regions are very effective in the prevention and control of brucellosis.

4.2.   Sensitivity analysis

As a zoonotic infection, the best way to control human brucellosis is to control the disease in animals. The control measures of livestock brucellosis include detection, vaccination and elimination of infected animals. We should also strengthen information dissemination and health education on brucellosis, and improve veterinary and public health supervision ^[25]. In this part, we used the PRCC (partial rank correlation coefficient) to analyze the sensitivity of $R_0$ on control parameters. The control parameters are disease-related culling rates $a_{1v}$ and $a_{2v}$, the infection rates of infected sheep to susceptible sheep $\delta_{1}$ and $\delta_{2}$, infection rates of Brucella to susceptible sheep $\phi_{1}$ and $\phi_{2}$, sheep sales rates $d_{1v}$ and $d_{2v}$, and sheep immigration rates $a_{12}$ and $a_{21}$. The sensitivity analysis of basic reproduction number $R_0$ to parameters is shown in Figure 3. We can observe that $R_0$ is more sensitive to $\delta_{1}$, $\delta_{2}$, $a_{12}$, $a_{1v}$, $d_{1v}$ and $d_{2v}$ (|PRCC| > 0.5), in which $\delta_{1}$, $\delta_{2}$ and $a_{12}$ are positively correlated (PRCC > 0.5), and $a_{1v}$, $a_{21}$, $d_{1v}$ and $d_{2v}$ are negatively correlated (PRCC < -0.5). There may be many factors interacting with human brucellosis. Therefore, when eradicating this disease, we can take a variety of measures, such as timely handling of the infected sheep in the two regions, reducing the infection rate from infected sheep to susceptible sheep in the two regions, decreasing the immigration of sheep from Hebei Province to Shanxi Province and improving the culling rate of infected sheep in Shanxi Province and increasing the sheep immigration rate from Shanxi Province to Hebei Province and the sales rate of sheep in the two regions. The relationships between the basic reproductive number $R_0$ and the disease-related culling rate $a_{2v}$ of sheep in Hebei Province and the infection rates $\phi_{1}$ and $\phi_{2}$ from Brucella to susceptible sheep of the two regions are weak, so we do not select these three parameters as control parameters. In particular, it can be observed that $R_0$ has a positive correlation with the sheep immigration rate ($a_{12}$) from Hebei Province to Shanxi Province and a negative correlation with the sheep immigration rate ($a_{21}$) from Shanxi Province to Hebei Province. Therefore, next, we studied the impact of the sheep immigration rate on human brucellosis.

4.3.   The influence of sheep immigration on human brucellosis infection

Due to $R_0$ having the opposite correlation with $a_{12}$ and $a_{21}$, in this part, we investigated the effects of sheep immigration in Shanxi and Hebei Provinces on human brucellosis infection. Figure 4(a), (b) shows that $R_0$ is more sensitive to $a_{12}$. We can reduce the immigration of sheep from Hebei Province to Shanxi Province to control the disease of these two patches. In addition, we chose different immigration rates $a_{12}$ and $a_{21}$ to study the changes of cumulative cases C1, C2 in the two provinces with time t (Figure 5). Figure 5(a) shows that as $a_{12}$ increases, the number of infected cases in the two regions is increasing, and when $a_{12} = 0$, the two regions have the least number of infected individuals. Figure 5(b) shows that as $a_{21}$ increases, the cases in Shanxi Province are decreasing, while the infected in Hebei Province are increasing. When $a_{21} = 0.44$, the least are infected in Shanxi Province, and the most are infected in Hebei Province. When $a_{21} = 0$, the most are infected in Shanxi Province, and the least are infected in Hebei Province. This indicates that in order to have fewer infected, the immigration of sheep cannot be completely absent, the immigration of sheep from Hebei Province to Shanxi Province should be reduced, and the immigration of sheep from Shanxi Province to Hebei Province should be controlled within a reasonable range.

The basic reproduction number $R_0$ is the expected number of secondary cases produced, in a completely susceptible population, by a typical infective individual ^[13]. For Shanxi Province and Hebei Province, $R_{01} = 0.7497 > R_{02} = 0.5022$, that is, in the population of all susceptible people, the average number of infected people in Shanxi Province is greater than that in Hebei Province. Here, we refer to Shanxi Province as the high-infected area and Hebei Province as the low-infected area. When there is only one-way immigration ($a_{12} = 0$ or $a_{21} = 0$), the change of cumulative cases C1, C2 with time t in both provinces is shown in Figure 6. When $a_{21} = 0$, only sheep immigrate from the low-infected area to the high-infected area, and Figure 6(a) shows that different values of $a_{12}$ increase human brucellosis in the high-infected area to different degrees. When $a_{12} = 0$, the increase is minimal. However, the situation of brucellosis in the low-infected zone shown in Figure 6(a) is exactly opposite to that shown in Figure 5(a), where increasing $a_{12}$ will reduce the occurrence of brucellosis in the low-infected zone but increase the occurrence of brucellosis in the high-infected zone. When $a_{12} = 0$ and only sheep immigrate from the high-infected area to the low-infected area, the situation shown in Figure 6(b) is similar to that shown in Figure 5(a) ($a_{12} = 0.07$) for C1, C2 with time t. This indicates that no sheep immigrate from the low-infected area to the high-infected area, or a small amount of immigration has no effect on brucellosis in these two areas. In conclusion, to better reduce human brucellosis in these two areas, there should be some degree of sheep immigration from high- to low-infection areas, and sheep immigration from low- to high-infection areas should be reduced.

From the above analysis, it can be seen that in order to better control human brucellosis in the two provinces, there needs to be some sheep immigration between the two provinces. So, what is the amount of sheep immigration under the best condition of disease control? We thought of the case that sheep immigration in these two provinces is a constant C and to minimize $R_0$ subject to the condition of $a_{21}A_{1v}+a_{12}A_{2v} = C$. We brought $a_{12} = \frac{C-a_{21}A_{1v}}{A_{2v}}$ into the $R_0$ expression, and then $R_0$ can be regarded as a function of $a_{21}$. Different immigration C values were selected to obtain different $\min R_0$. Figure 7 shows that $\min R_0$ first decreases and then increases with the increase of immigration. The maximum immigration under the condition of the best disease control effect is $C = 10,000,000$. The immigration rates of Shanxi Province and Hebei Province are $a_{21} = 0.88$ and $a_{12} = 0.3$, respectively, and the immigration rates are $C_1 = 3,928,894$ and $C_2 = 6,071,105$, respectively.

In fact, we collected the data of human brucellosis in Shanxi Province and Hebei Province from 2010 to 2020, and we observed that the number of cases in Shanxi Province and Hebei Province increased sharply in 2019 and 2020. To facilitate our study, we used the number of human brucellosis cases from 2010 to 2018. According to our research, human brucellosis will disappear in the two provinces. Based on the parameters in Table 1, MATLAB was applied to derive the numerical solutions of cumulative cases C1(t) and C2(t) from 2010 to 2020, and the predicted values of new infections in 2019 and 2020 were obtained according to $I_i^h (t) = Ci(t)-Ci(t-1), i = 1, 2, t = 2019, 2020.$ It is predicted that the numbers of human brucellosis cases in Shanxi Province may be 1912 and 1373, and the numbers of human brucellosis cases in Hebei Province may be 1732 and 1258, in 2019 and 2020 (see Figure 8). However, in 2019 and 2020, the actual numbers of cases in Shanxi Province were 3279 and 3365, and in Hebei Province they were 3236 and 2968. There is a huge difference between the predicted data and the actual data. According to our analysis, the main reasons are the occurrence of African swine fever in 2018, the closure of a large number of domestic pig farms and slaughterhouses, a large reduction in the source of pigs, a shortage of pork and a significant increase in prices. Most people chose to buy mutton, and the immigration of sheep in various provinces had also increased greatly, including Shanxi Province and Hebei Province. The frequent trading of sheep products between different regions and the increased mobility of infected animals led to a significant increase in the number of infected people in Shanxi and Hebei provinces in 2019 and 2020, indicating that sheep immigration has a great impact on human brucellosis infection.

5.   Conclusions

In the past 30 years, brucellosis epidemic areas have gradually shifted from pastoral areas to grassland and agricultural areas. Especially since 2004, the affected areas have expanded from the north to the south of China ^[5,6], and studies have confirmed that part of the epidemic in the southern region was caused by infected animals imported from other regions ^[7]. Therefore, the geographical transmission of brucellosis is caused by the immigration of sheep.

In this paper, a patch model was proposed to describe the spatial transmission dynamics of brucellosis and to study the impact of sheep immigration on the geographical transmission of brucellosis. Firstly, we analyzed the dynamics of the model, including the basic reproduction number and the existence, uniqueness and stability of the positive equilibrium.

Secondly, taking Shanxi Province and Hebei Province as examples, numerical simulation was carried out. The parameters were estimated by the least squares method, and $R_{0} = 0.5457$, $R_{01} = 0.7497$ and $R_{02} = 0.0.5022$ were obtained, which indicate that brucellosis will disappear in the two provinces. Sensitivity analysis of $R_0$ found that the infection rates $\delta_{1}$ and $\delta_{2}$ from infected sheep to susceptible sheep, the sheep sale rates $d_{1v}$ and $d_{2v}$, the sheep immigration rate $a_{12}$ from Hebei Province to Shanxi Province and the disease-related culling rate $a_{1v}$ of sheep in Shanxi Province had greater impacts on $R_0$. Therefore, when eradicating the disease, we can take a variety of measures, such as vaccinating sheep, timely dealing with infected sheep in these two areas, descreasing the infection rate of infected sheep to susceptible sheep in the two areas, reducing the immigration of sheep from Hebei Province to Shanxi Province, improving the culling rate of infected sheep in Shanxi Province and increasing the sheep immigration rate from Shanxi Province to Hebei Province and the sales rate of sheep in the two regions.

Finally, we studied the influence of sheep immigration rate on the occurrence of disease. In terms of sheep immigration between the two regions, we should minimize the sheep immigration from Hebei Province to Shanxi Province and control the sheep immigration from Shanxi Province to Hebei Province within a reasonable range. According to $R_{01} = 0.7497$ and $R_{02} = 0.5022$, when we only considered one-way immigration, we found that there should be a certain degree of immigration of sheep from high-infection area to low-infection area, and we should lessen the immigration of sheep from low-infection area to high-infection area.

We studied the geographic transmission of sheep brucellosis using a deterministic system and simulated case data from 2010 to 2018 in Shanxi and Hebei provinces. Our model does not consider stochasticity, periodicity and age structure. Next, the patch model combined with these factors can be studied. Meanwhile, only the two-patch model was used to simulate the data of two provinces. Complex transmission among three or more provinces needs to be researched.

Acknowledgments

This work is supported by the National Natural Science Foundation of China under grants (12101443, 11801398) and the Natural Science Foundation of Shanxi Province grants (20210302124260).

Conflict of interest

The authors declare there is no conflict of interest.