Reconsideration of the plague transmission in perspective of multi-host zoonotic disease model with interspecific interaction

1.   Introduction

Plague is an exceedingly virulent infectious disease that has a high mortality rate without treatment ^[1,2,3]. The most serious symptom of plague is its typical bubonic form, which will lead to 40%–70% mortality ^[2]. It has caused at least three significant pandemics with millions of deaths ^[1,2,4,5]. During the Second European Pandemic, plague was called "Black Death" because of its special clinical symptoms, namely blackening and death of tissue in human extremities ^[6]. Although the "Black Death", as one of the turning points of Western Civilization ^[1], should be listed as the most important biological environmental event in European history, people know little about the dynamic mechanism of plague transmission. More importantly, it seems not only of historical significance, but also plague is still threatening human health today ^[2,7,8]. According to the latest report, local health officials confirmed the two cases of pneumonic plague in Beijing (China) on November 13, 2019 ^[9]. Therefore, reconsideration of the spread of plague, especially in China, a country with a large natural plague foci, must be re-emphasized ^[10].

Plague is caused by the bacterium Yersinia pestis ^[11] and is transmitted primarily by the bite of adult fleas from an infected rodent to other mammals ^[1]. It is reported that more than 200 mammalian species are naturally infected with Yersinia pestis ^[1]. In which, rodents are the main enzootic (maintenance) host and epizootic (amplification) host of plague ^[1,12], and they play important roles in plague long-term survival. Global infectious disease eradication schemes have achieved two remarkable successes in smallpox and rinderpest, which spread only from person to person, but the plague eradication program has proven elusive in more complex ecosystems ^[13,14]. Most wild rodents are the main enzootic host ^[15], while most commensal rodents are the main epizootic host ^[1,12]. For humans, they are infected by wild rodents in rural settings, or by commensal rodents that move freely between villages, forests and urban area ^[2]. Although the transmission route of plague is very complex and may contain more than 200 kinds of mammals and their flea ectoparasites, wild rodents-commensal rodents-humans is the common transmission route for human beings ^[1,2,12], as is shown in Figure 1.

In the absence of dramatic climate change, the surviving species should maintain a stable population size in their respective communities, which is called environmental capacity or $K$ ^[16]. However, $K$ will be affected by interspecific competition between rodents ^[17,18], which determines the decrease of stable coexisting population size (Figure 2A) ^[16,19]. After the human colonization, this balance between rodents will be broken again(Figure 2B). If there was no plague-death, we assume that the presence of rodents has little effect on the changes in the size of human population. After all, rodents do not eat enough food to reduce the birth rate of humans in times other than famine. Therefore, the size of the red circle representing the human population in Figure 2B has not changed.

The development of agriculture and the accumulation of food waste in cities provide more food for commensal rodents ^[20,21]. The process of human urbanization is eating away the natural habitat of wild rodents, thus increasing the opportunities for rodents to come into contact with humans ^[22]. The human-animal interface promotes the cross-species spread of plague, leading to the epidemic of plague in humans. Therefore, the interaction between different rodents should be considered in plague transmission, which has been neglected in previous studies ^{[4,5,24,25,26,27]}. This paper establish a new multi-host mathematical model to describe the complex plague ecosystem composed of three host populations (wild rodent, commensal rodent and human). Then we redefine the environmental capacity $K$ with interspecific interaction and reconstruct the zoonotic basic reproductive number $R^0_Z$.

In addition, we also analyze the biological implication of our plague model and test some biological hypotheses, such as the reason for introducing and receding of the European plague pandemic ^[22]. Burrowing rodent(wild rodent)-black rat(commensal rodent)-human interface may result the plague spread into Europe. Biologists speculate that the black rat was identified as the culprit during the Second Plague Pandemic, which colonized western regions along trade routes(Silk Road) and flourished in the great late Roman cities ^[22,23]. However, the receding plague epidemics from the eighteenth century may be due to the colonization of Europe by the brown rat with stronger commensal proclivity, which does not harbor anthropophilic fleas and their chances of transmitting plague to humans through fleas are much smaller ^[22]. The competition between black rat and brown rat, and their commensalism with human may lead to the disappearance of plague in Europe. The absence of data may prevent us from finding out the truth of plague history, but the simulation of multi-host plague model may allow us to find out the possible plague transmission process and prove the rationality of speculation in plague historical research.

2.   Methods

A new multi-SIRs model is established to reflect the spread of plague among wild rodents, commensal rodents, and humans, as Figure 3. The notation $N_1(t)$ (wild rodents), $N_2(t)$ (commensal rodents) and $N_3(t)$ (humans) represent the population size of these populations at time $t$. Each population are assumed to be divided into 3 epidemiological compartments: susceptibles ($S_i(t)$), infectives ($I_i(t)$), and recovered individuals ($R_i(t)$) at time $t$, $i = 1, 2, 3$. Competition between wild rodents and commensal rodents is represented by the notation $B_{12}$ and $B_{21}$. Commensalism between rodents and humans is represented by the notation $B_{13}$ and $B_{23}$, where $B_{13} < B_{23}$. Wild rodents, commensal rodents and humans have the maximum environmental capacity at $K_1$, $K_2$ and $K_3$ respectively, and they satisfy the logistic growth curve, $dN_i(t)/dt = r_iN_i(t)(1-N_i(t)/K_i)$, where $r_i$ represents the intrinsic growth rate of the populations.

with the initial conditions $S_i(0) = S_0^i > 0, I_i(0) = I_0^i > 0, R_i(0) = R_0^i > 0. \ \ i = 1, 2, 3$.

The parameter $r_i = b_i-d_i$, in which $b_i$ represents the birth rate while $d_i$ represents the natural mortality rate. The notation $\gamma_i$ denotes the recovery rate, and $e_i$ is the disease-induced mortality rate. $\phi _1 (t) = - \frac{{a_1 N_1 (t)}}{{K_1 }} - \frac{{a'_1 B_{12} N_2 (t)}}{{K_1 }} + \frac{{a''_1 B_{13} N_3 (t)}}{{K_1 }}$, $\phi _2 (t) = - \frac{{a_2 N_2 (t)}}{{K_2 }} - \frac{{a'_2 B_{21} N_1 (t)}}{{K_2 }} + \frac{{a''_2 B_{23} N_3 (t)}}{{K_2 }}$ and $\phi _3 (t) = - \frac{{a_3 N_3 (t)}}{{K_3 }}$ represent the contribution of inter- and intraspecific interactions into the births change. $\psi _1 (t) = \frac{{(1 - a_1)N_1 (t)}}{{K_1 }} + \frac{{(1 - a'_1)B_{12} N_2 (t)}}{{K_1 }} - \frac{{(1 - a''_1)B_{13} N_3 (t)}}{{K_1 }}$, $\psi _2 (t) = \frac{{(1 - a_2)N_2 (t)}}{{K_2 }} + \frac{{(1 - a'_2)B_{21} N_1 (t)}}{{K_2 }} - \frac{{(1 - a''_2)B_{23} N_3 (t)}}{{K_2 }}$ and $\psi _3 (t) = \frac{{(1 - a_3)N_3 (t)}}{{K_3 }}$ represent the contribution of inter- and intraspecific interactions into the deaths change. The notation $a_i$ is the parameter subdividing the contribution of intraspecific competition into the births decrease ($a_i$) and deaths increase ($1-a_i$), with $0\leq a_i\leq 1$. $a^\prime_i$ is the parameter subdividing the contribution of interspecific competition into the births decrease (${a^\prime}_i$) and deaths increase ($1-{a^\prime}_i$), with $0\leq {a^\prime}_i\leq 1$. While $a^{\prime\prime}$ is the parameter subdividing the contribution of commensalism from human to rodents into the births increase (${a^{\prime\prime}}_i$) and deaths decrease ($1-{a^{\prime\prime}}_i$), with $0\leq {a^{\prime\prime}}_i\leq 1$.

The parameters $\beta_{ij}$ represents the per capita incidence rate from population $j$ to population $i$, where $i, j = 1, 2, 3$. The transmission of the pathogen from rodents to humans is assumed to be unidirectional, with a low probability of occurrence for transmission in the other direction. Therefore, we assume that $\beta_{31} = \beta_{32} = 0$. We assume that all parameters are positive above. In this model, we will no longer set the fleas as the separate compartment, and the effect of them is attributed to the change of $\beta_{ij}$, which will make a final decision to the propagation speed of the pathogen and the value of plague basic reproductive number.

2.1.   Equilibrium points and zoonotic reproduction number $R^0_Z$

After calculation, we summarize that there are at most 15 equilibrium points of the plague model, $E_0$, $E_1$, $E_2$, $E_3$, $E_4$, $E_5$, $E_6$, $E_7$, $E_8$, $BE_1$, $BE_2$, $BE_3$, $BE_4$, $BE_5$, $BE_6$. Their existence and stability conditions are shown in Table 1, Appendix A and Appendix B.

Then, we define the zoonotic reproduction number $R^0_Z$ for the plague transmission, including all elements considered in the introduction of this article ^[28,29,30].

here $\rho(\cdot)$ is the spectral radius.

By the persistence theory ^[31,32,33], we can get the existence conditions of $E_8$. However, the exact expression of $E_8$ is unknown. Similarly, we did not get the exact expression of $BE_1$, $BE_2$, $BE_3$, $BE_4$, $BE_5$, and $BE_6$, either. Even if their exact expressions are unknown, we can obtain the uniqueness condition for their existence. Fortunately, if this condition is met, their stability is obvious when they already exist in the system. The unique variation trend of some fixed initial values and parameters can be obtained according to the continuity of the solution and the continuous dependence on the initial values. This is enough to reflect the spread of plague, which is the main issue to be discussed later.

3.   Results

3.1.   Wildlife extinction caused by human colonization without hunting or predation

Under stable and limited environmental conditions, the population size always tends to a constant value, which is usually called the environmental capacity $K$. However, if we consider more than two populations ^[16], the constant $K$ is no longer applicable due to interspecific interactions between them, such as predation($+, -$), competition($-, -$), mutualism($+, +$), commensalism($+, 0$), amensalism($-, 0$) ^[16,34]. Therefore, the equilibrium point mentioned above is proposed to redefine the capacity $K$. The equilibrium point is the mathematical result of ODE model representing the final population size. $E_1$, $E_2$, and $E_3$ in Table 1 reflect the environmental capacity $K_1$, $K_2$, and $K_3$ corresponding to the three host populations, respectively. $E_4$, $E_5$, and $E_6$ describe the environmental capacity of coexistence of ($N_1$, $N_2$), ($N_1$, $N_3$), and ($N_2$, $N_3$). And $E_7$ describe the environmental capacity of coexistence of the three host populations without plague.

According to the calculation of our plague model, not all of the equilibrium points could exist all the time. Their existence should be satisfied with the conditions shown in Table 1. The change of these conditions would be the quantitative and qualitative evidence to answer the questions shown in Figure 2. For example, if we do not consider the infectious disease, $E_7$ is globally asymptotically stable in its domain of definition when $B_{12} B_{21} < 1$. $E_7$ is a saddle point when $B_{12} B_{21} > 1$. Compared with the classic Lotka-Volterra competitive mathematical model ^[34], the participation of the third population will not change the stability condition of the internal equilibrium point, but will change its existence condition. The threshold values turn from $B_{12} = \frac{K_1}{ K_2}$ and $B_{21} = \frac{K_2}{K_1}$ in ^[34] (the notations $B_{12}$ and $B_{21}$ instead of $b_{12}$ and $b_{21}$ in ^[34]) to $B_{12} = \frac{K_1+B_{13} K_3}{K_2+B_{23} K_3}$ and $B_{21} = \frac{K_2+B_{23} K_3}{K_1+B_{13} K_3}$ in our plague model. When taking appropriate experimental parameters, Figure 4(a) and (b) show that the existence of humans would break the original balance between the first two competitive rodent species. The wild rodents ($N_1$) have a more exceptional ability to survive in wild field. However, after human colonization, $N_1$ declines from dominant species to inferior species, even to extinction, as shown from point $A$ to point $B$ in Figure 4(a). However, the commensal rodents $N_2$ develop to a bigger population size. Therefore, we can infer that the original ecological balance will be changed only due to human colonization through commensalism of animals in different degrees, which has not been discussed in previous quantitative studies.

Then, considering the plague into its host populations, the incidence of plague deaths leads to the change of equilibrium points, shown by the point C in Figure 4(a) and the green curve in Figure 4(b). Disease-free boundary equilibrium points $E_1$, $E_2$, $E_3$, $E_4$, $E_5$, $E_6$, $E_7$ change to $BE_1$, $BE_2$, $BE_3$, $BE_4$, $BE_5$, $BE_6$, $E_8$, respectively. The threshold values then turn to $B_{12} = \frac{K_1+B_{13} \overline{N}_3}{K_2+B_{23} \overline{N}_3}$ and $B_{21} = \frac{K_2+B_{23} \overline{N}_3}{K_1+B_{13} \overline{N}_3}$, in which $\overline{N}_3$ represents the equilibrium point of $N_3$ without plague deaths. It is easy to know that $\overline{N}_3 < K_3$. However, because of the complexity of the plague model, we cannot give the exact mathematical expression of $\overline{N}_3$. More simulation results will be discussed in the following part.

3.2.   Effect of zoonotic infectious disease on host population size change

Wild animals cause most zoonotic infectious diseases, but they are not toxic to them. Some wild plague infections can only cause inapparent to mild illness ^[1]. However, in humans, the plague occurs in bubonic form and pneumonic form, which has a high mortality rate without treatment. The simulations about the effect of plague on population size change are shown in Figure 5(A1) and (A2). The two simulations reflect the population change of rodents and humans with $B_{12} < \frac{K_1+B_{13} K_3}{K_2+B_{23} K_3}$ and $B_{21} < \frac{K_2+B_{23} K_3}{K_1+B_{13} K_3}$, which guarantees the existence of $E_7$. If $R^0_Z < 1$, $E_8$ does not exist and $E_7$ is stable. However, the process of population size change could be different with disease-induced mortality rate $e_i = 0$ or not. The time to the equilibrium points would get longer with bigger $e_i$ as the direction of the arrows in Figure 5(A1). If $R^0_Z > 1$, $E_8$ is existent and stable. The value of equilibrium points depend on $e_i = 0$ or not. If $e_i = 0$, $(N_1(t), N_2(t), N_3(t)) \rightarrow E_7$ with $t \rightarrow \infty$. However, if $e_i > 0$, the equilibrium points would change depending on $I_i^{**}$, $e_i$ and $r_i$. Choosing appropriate experimental parameters, Figure 5(A2) shows that if $e_i > 0$, inferior species $N_1$ once again turn to be the dominant one. Plague is the protection of $N_1$'s existence. It proves the importance of One Health, which emphasizes the coordinated health of human beings, animals and ecosystems as the best way to face a new zoonosis ^[22]. Furthermore, when humans invade the original habitat of animals, they will not only encounter physical defense of animals, but also attack by poisonous pathogens.

3.3.   Effect of interspecific interaction on the zoonotic infectious disease

The colonization of Europe by the brown rats might be regarded as one of the biological hypotheses for the receding plague epidemics from the eighteenth century ^[22]. The brown rats would depend more on human beings than the black rats, which flourished in the great late Roman cities in Europe ^[22]. Because of their strong commensal proclivity, they drove the black rats out of the city. However, brown mice do not harbor anthropophilic fleas, so their chances of transmitting plague to humans through fleas are much smaller. In our model, we take $N_1$ as the black rats and $N_2$ as the brown rats. Then we suppose that $B_{23} > B_{13} > 0$ and $\beta_{31} > \beta_{32} > 0$. Because we do not get the exact data of all animals from the Second Plague Pandemic in the Middle Ages, we only choose some special experimental parameters in ^[24] in order to verify the different results with or without interspecific interaction. Figure 5(B1) and (B2) show that the existence of $N_2$ in our model would lead to the extinction of the plague epidemic. If $N_2 = 0$, the colonization of new areas by humans would continued to be troubled by plague without any prevention and control measures. However, if $N_2 > 0$, the competition and the commensalism would lead to the disappearance of plague in Europe. In Figure 5(B1), the human infections diminish gradually with $N_1$ replaced by $N_2$ in cities. And the deaths caused by plague would remain unchanged shown in Figure 5(B2). In the absence of historical data, through dynamic process simulation, it is proved that the mathematical model can prove the rationality of speculation in plague historical research.

3.4.   The advantage of zoonoses mathematical model with interspecific interaction

Comparing to the human infectious disease mathematical model ^[24], our plague model has more complicated zoonotic reproduction number $R^0_Z$.

In order to discuss the effect of zoonotic reproduction number $R^0_Z$, we suppose that $e_i = 0$ to eliminate the effect of population size variation by zoonoses and $I^3_0 = 0$ to intensify the infection from animal hosts. Furthermore, we suppose that at least one of $\beta_{31}$ and $\beta_{32}$ is bigger than 0, or else there is no transmission between animals and humans. Figure 5(C1) shows that the infected population size change with $R^0_Z > 1$, $R^0_{12} > 1$ and $0 < \beta_{31}, \beta_{32} \ll 1$. Even if $\beta_{31}$ and $\beta_{32}$ are small enough to transmit the pathogen from animals to humans and $R^0_3 = 0$, we still can not neglect the infected population with low incidence due to $I(\infty) > 0$ with $R^0_Z > 1$. Moreover, the human morbidity will increase with the increasing of $R^0_3$, which represents the basic production number of human plague model. In the previous studies, we have always believed that human epidemics will only occur when there is continually transmission in humans, that is $R^0_3 > 0$. The results in Figure 5(C1) seem to be evidence of low infection rates. However, if we choose $\beta_{31}$ and $\beta_{32}$ bigger enough as is shown in Figure 5(B2), the zoonoses would break out in humans similar to the results with $R^0_3 > 0$ in Figure 5(B1). Therefore, we can get the opposite conclusion to the results. Figure 5(C3) shows the coefficient of $\beta_{31}$, $\beta_{32}$ and $R^0_3$. Although the parameters $\beta_{31}$ and $\beta_{32}$ does not appear on the expression of $R^0_Z$, and they do not determine the outbreak of zoonoses in the humans, these two parameters could affect on the equilibrium of $I_3$. And then, they can determine the density of human infection, which is essential in the prevention and control of zoonotic diseases.

It is worth noting that the variation of $\beta_{31}$ and $\beta_{32}$ in Figure 5(C2) and (C3) also reflects the influence of human ectoparasitic fleas in plague transmission. Different from the method in ^[24], the effect of ectoparasitic fleas is explored by the addition of flea compartments and the selection of the models by fitting historical data, we judge the impact of fleas by simulating the sensitivity of parameters $\beta_{31}$ and $\beta_{32}$.

Meanwhile, with the continuous acceleration of urbanization, the contradiction between economic development and ecological environment is becoming increasingly prominent. In previous studies, the quantitative results of the impact of urbanization on zoonosis are rarely involved. By our plague model, we get that $R^0_Z$ is increasing with the increase of $K_3$ by the exact mathematical expression, which is the maximum environmental capacity of humans. The urban environmental capacity would increase as the efficient land utilization with increasing food production and convenient transportation ^[35]. $K_3$ value will increase with the process of urbanization, and then with a higher $R^0_Z$ value to increase the risk of infectious disease transmission. To sum up, many dynamic processes of disease ecology can be reflected by research on zoonosis models.

4.   Discussion

Our study establish a new plague mathematical model to reflect the dynamic process of plague transmission with interspecific action between wild rodents, commensal rodents and humans, in which the roles of different species will no longer be at the same level. Mathematical models in epidemiology can elucidate the mechanisms underlying plague transmission among its hosts and provide a way to more quickly and easily reflect the dynamic process of plague transmission ^[36]. By our new plague model, we find that the human-rodent interface has promoted the cross-species transmission of plague and resulted in the prevalence of plague in humans. In addition, the threshold values of population development and disease transmission are also discussed in order to provide scientific basis for future health decision makers in plague prevention and control.

The introduction of multiple hosts into zoonotic disease model leads to the failure of constant $K$ in describing the environmental capacity. Therefore, we refine it by the boundary disease-free equilibrium points calculated by our plague model. These equilibrium points represent the property of population size in a particular environment, including interspecific interaction. Correspondingly, we have to redefine the threshold conditions that affect their existence and stability conditions. And then we propose the zoonotic reproduction number $R^0_Z$, which is more applicable to the study of plague transmission. These results reflect the methodological value of plague model.

In the biological application of our plague model, we focus more on the interaction between interspecific relationships and zoonotic disease transmission. Our model uses a combination of quantitative and qualitative methods to study many biological hypotheses. For example, we prove the change of threshold conditions for rodent coexistence after human colonization. Without hunting or predation, human colonization may lead to the extinction of wild animals due to commensal proclivity. Furthermore, we simulate the plague epidemic with or without commensal rodents, which may answer the rationality of one of the hypotheses for disappearance of the Second Plague Pandemic.

Acknowledgments

This work is supported by the National Natural Science Foundation of China, 11771044.

Conflict of interest

The authors declare no conflicts of interest.

Appendix A

If $e_i = 0$, we integrate $N_i (t) = S_i (t) +I_i (t) +R_i (t), i = 1, 2, 3$ of system (2.1), then we get system

with the initial conditions $N_i(0) = N_0^i = S_0^i+I_0^i+R_0^i > 0, i = 1, 2, 3$.

System (4.1) always has the following equilibrium points:

If $\frac{B_{12}K_2-K_1}{B_{12}B_{21}-1} > 0, \frac{B_{21}K_1-K_2}{B_{12}B_{21}-1} > 0$, then there exists a unique positive boundary equilibrium:

If $\frac{B_{12}K_2-K_1+(B_{12}B_{23}-B_{13})K_3}{B_{12}B_{21}-1} > 0, \frac{B_{21}K_1-K_2+(B_{21}B_{13}-B_{23})K_3}{B_{12}B_{21}-1} > 0$, then there exists a unique positive internal equilibrium:

Lemma 1. The set $\; \Omega = \{(N_1, N_2, N_3) | N_i \geq 0, i = 1, 2, 3, N_j \leq K_j + B_{j3} K_3, j = 1, 2, N_3\leq K_3\}$ is a positively invariant region for Model (4.1). Moreover, every trajectory of Model (4.1) is eventually staying in a compact subset of $\Omega$.

Theorem 1. For Model (4.1), $E_0$, $E_1$, $E_2$, $E_3$ and $E_4$ are always unstable. As for the existence and the stability of $E_5$, $E_6$ and $E_7$, there are 4 cases:

1. If $B_{12} < \frac{K_1+B_{13} K_3}{K_2+B_{23} K_3}$ and $B_{21} > \frac{K_2+B_{23} K_3}{K_1+B_{13} K_3}$, there is no internal equilibrium point $E_7$, and $E_5$ is locally asymptotically stable in ${\Omega}$, $E_6$ is a saddle point.

2. If $B_{12} > \frac{K_1+B_{13} K_3}{K_2+B_{23} K_3}$ and $B_{21} < \frac{K_2+B_{23} K_3}{K_1+B_{13} K_3}$, there is no internal equilibrium point $E_7$, and $E_6$ is locally asymptotically stable in ${\Omega}$, $E_5$ is a saddle point.

3. If $B_{12} > \frac{K_1+B_{13} K_3}{K_2+B_{23} K_3}$ and $B_{21} > \frac{K_2+B_{23} K_3}{K_1+B_{13} K_3}$, there exists a unique internal equilibrium point $E_7$, which is a saddle point. Both $E_5$ and $E_6$ are locally asymptotically stable.

4. If $B_{12} < \frac{K_1+B_{13} K_3}{K_2+B_{23} K_3}$ and $B_{21} < \frac{K_2+B_{23} K_3}{K_1+B_{13} K_3}$, there exists also a unique internal equilibrium point $E_7$, which is locally asymptotically stable in ${\Omega}$, while either $E_5$ or $E_6$ is unstable.

Proof. The Jacobian of system (4.1) is

$J$ = $\left(\begin{array}{ccc} r_1 (1 - \frac{2N_1}{K_1} - \frac{B_{12} N_2}{K_1} + \frac{B_{13} N_3 }{K_1 }) & { - r_1 N_1 \frac{{B_{12} }}{{K_1 }}} & {r_1 N_1 \frac{{B_{13} }}{{K_1 }}} \\ { - r_2 N_2 \frac{{B_{21} }}{{K_2 }}} & {r_2 (1 - \frac{{2N_2 }}{{K_2 }} - \frac{{B_{21} N_1 }}{{K_2 }} + \frac{{B_{23} N_3 }}{{K_2 }})} & {r_2 N_2 \frac{{B_{23} }}{{K_2 }}} \\ 0 & 0 & {r_3 (1 - \frac{{2N_3 }}{{K_3 }})} \\ \end{array} \right)$

For $\; E_0 = (0, 0, 0)$, we have $J_{E_0}$ = $\left(\begin{matrix}r_1 & 0 & 0\\0&r_2 & 0\\0 & 0&r_3\\\end{matrix}\right)$. So $E_0$ is unstable.

For $E_1 = (K_1, 0, 0)$, we have $J_{E_1}$ = $\left(\begin{matrix}-r_1 & -r_1B_{12}&r_1B_{13}\\0&r_2(1-\frac{B_{21}K_1}{K_2}) & 0\\0 & 0&r_3\\\end{matrix}\right)$. So $E_1$ is unstable.

For $E_2 = (0, K_2, 0)$, we have $J_{E_2} = \left(\begin{matrix}r_1(1-\frac{B_{12}K_2}{K_1}) & 0 & 0\\-r_2B_{21} & -r_2&r_2B_{23}\\0 & 0&r_3\\\end{matrix}\right)$. So $E_2$ is unstable.

For $E_3 = (0, 0, K_3)$, we have $J_{E_3} = \left(\begin{matrix}r_1(1+\frac{B_{13}K_3}{K_1}) & 0 & 0\\0&r_2(1+\frac{B_{23}K_3}{K_2}) & 0\\0 & 0 & -r_3\\\end{matrix}\right)$. So $E_3$ is unstable.

For $\; E_4 = ({N^\prime}_1, {N^\prime}_2, 0) = (\frac{B_{12} K_2 - K_1 } { B_{12} B_{21} -1}, \frac{B_{21} K_1 - K_2} {B_{12} B_{21}-1}, 0)$, we have

$\; J_{E_4} = \left(\begin{matrix} - r_1 \frac{{N^\prime}_1}{K_1} & - r_1{N^\prime}_1 \frac{B_{12}}{K_1} & r_1{N^\prime}_1\frac{B_{13}}{K_1} \\ - r_2{N^\prime}_2\frac{B_{21}}{K_2} & -r_2\frac{{N^\prime}_2}{K_2} & r_2{N^\prime}_2\frac{B_{23}}{K_2} \\ 0 & 0 & r_3\\ \end{matrix}\right)$. So $E_4$ is unstable.

For $\; E_5 = \left(\overline{N_1}, 0, \overline{N_3}\right) = \left(K_1+B_{13}K_3, 0, K_3\right)$, we have

$J_{E_5} = \left(\begin{matrix}r_1(-1-\frac{B_{13}K_3}{K_1}) & -r_1(K_1+B_{13}K_3)\frac{B_{12}}{K_1}&r_1(K_1+B_{13}K_3)\frac{B_{13}}{K_1}\\0&r_2(1-\frac{B_{21}(K_1+B_{13}K_3)}{K_2}+\frac{B_{23}K_3}{K_2}) & 0\\0 & 0 & -r_3\\\end{matrix}\right)$.

If $B_{21}K_1-K_2+(B_{21}B_{13}-B_{23})K_3 > 0$, $E_5$ is locally asymptotically stable.

If $B_{21}K_1-K_2+(B_{21}B_{13}-B_{23})K_3 < 0$, $E_5$ is a saddle point, which is unstable.

For $\; E_6 = \left(0, \overline{N_2}, \overline{N_3}\right) = \left(0, K_2+B_{23}K_3, K_3\right)$, we have

$J_{E_6} = \left(\begin{matrix}r_1(1-\frac{B_{12}(K_2+B_{23}K_3)}{K_1}+\frac{B_{13}K_3}{K_1}) & 0 & 0\\-r_2(K_2+B_{23}K_3)\frac{B_{21}}{K_2}&r_2(-1-\frac{B_{23}K_3}{K_2})&r_2(K_2+B_{23}K_3)\frac{B_{23}}{K_2}\\0 & 0 & -r_3\\\end{matrix}\right)$.

If $B_{12}K_2-K_1+(B_{12}B_{23}-B_{13})K_3 > 0$, $E_6$ is locally asymptotically stable.

If $B_{12}K_2-K_1+(B_{12}B_{23}-B_{13})K_3 < 0$, $E_6$ is a saddle point, which is unstable.

For $E_7 = (N_1^\ast, N_2^\ast, N_3^\ast) = (\frac{B_{12}K_2-K_1+(B_{12}B_{23}-B_{13})K_3}{B_{12}B_{21}-1}, \frac{B_{21}K_1-K_2+(B_{21}B_{13}-B_{23})K_3}{B_{12}B_{21}-1}, K_3)$, we have

$J_{E_7} = \left(\begin{matrix}-r_1\frac{N_1^\ast}{K_1} & -r_1N_1^\ast\frac{B_{12}}{K_1}&r_1N_1^\ast\frac{B_{13}}{K_1}\\-r_2N_2^\ast\frac{B_{21}}{K_2} & -r_2\frac{N_2^\ast}{K_2}&r_2N_2^\ast\frac{B_{23}}{K_2}\\0 & 0 & -r_3\\\end{matrix}\right)$.

If $B_{12}B_{21} < 1$, $E_7$ is locally asymptotically stable.

If $B_{12}B_{21} > 1$, $E_7$ is a saddle point, which is unstable.

Appendix B

If $e_i > 0$, we focus on the equilibria of model (2.1) and study their stability.

At first, we only consider the third species $N_3$ (human species) of model (2.1). If $I_1(0) = I_2(0) = 0$, then we get model as follows:

Lemma 2. (^[42,43])

For model (4.2), if $\; R^0_3 = \frac{ \beta_{33} K_3} {b_3 + \gamma_3 + e_3 - a_3 r_3} < 1$, the free disease equilibrium $\; E_3^\ast = (K_3, 0, 0)$ is globally asymptotically stable in $\; \Omega_3 = \{(S_3, I_3, R_3)\mid S_3 \geq 0, I_3 \geq 0, R_3 \geq 0, 0 \leq S_3 + I_3 + R_3 \leq K_3 \}$.

If $\; R^0_3 = \frac{ \beta_{33} K_3} {b_3 + \gamma_3 + e_3 - a_3 r_3} > 1$, there exists a unique equilibrium $\; E_3^{\mathrm{**}} = (S_3^\ast, I_3^\ast, R_3^\ast)$, which is globally asymptotically stable in $\; \Omega_3 = \{(S_3, I_3, R_3)\mid S_3 \geq 0, I_3 \geq 0, R_3 \geq 0, 0 \leq S_3 + I_3 + R_3 \leq K_3 \}$.

According to the definition of basic reproduction number in ^[30]. If $J$ is the Jacobian matrix of infective compartments, then let $J = F-V$, $F$ be the rate of appearance of new infections in compartment $I$, $V$ be the rate of transfer of individuals out of compartment $I$. We call $FV^{-1}$ be the next generation matrix for the model (2.1) and set $\rho(FV^{-1})$ at epidemic equilibrium point be the basic reproduction number, where $\rho(A$) denotes the spectral radius of a matrix $A$.

For the single population, such as human population $N_3$, $F = \beta_{33} S_3$, $V = b_3+\gamma_3+e_3-a_3 r_3$, and its epidemic equilibrium point is $(K_3, 0, 0)$, then we can obtain that $R^0_3 = FV^{-1} = \frac{ \beta_{33} K_3} {b_3 + \gamma_3 + e_3 - a_3 r_3}$.

Similarly, for the basic reproduction number defined in two populations, such as $N_1$ and $N_2$, $R^0_{12}$ can be obtained by

Furthermore, let $R^0_1 = \frac{\beta_{11}K_1}{b_1+\gamma_1+e_1-a_1r_1}$; $R^0_2 = \frac{\beta_{22}K_2}{b_2+\gamma_2+e_2-a_2r_2}$; $R^0_3 = \frac{\beta_{33}K_3}{b_3+\gamma_3+e_3-a_3r_3}$;

$R^0_{12} = \rho \left (\begin{matrix} \frac{\beta_{11}{N^\prime}_1}{\alpha_3} & \frac{\beta_{12}{N^\prime}_2}{\alpha_4} \\ \frac{\beta_{21}{N^\prime}_1}{\alpha_3} & \frac{\beta_{22}{N^\prime}_2}{\alpha_4} \\ \end{matrix} \right)$;

$R^0_{13} = \max \left \{ \frac{ \beta_{33} K_3} {b_3 + \gamma_3 + e_3 -a_3 r_3}, \frac {\beta_{11}(K_1 + B_{13}K_3)}{d_1 + e_1 + \gamma_1 + r_1 \left (-\frac{(1 - {a''}_1) B_{13} K_3}{K_1} + \frac {(1 - a_1) (K_1 + B_{13} K_3)}{K_1} \right)} \right\}$;

$R^0_{23} = \max \left \{ \frac{\beta_{33} K_3}{b_3 + \gamma_3 + e_3 - a_3 r_3}, \frac{\beta_{22}(K_2 + B_{23}K_3)}{d_2 + e_2 + \gamma_2 + r_2 \left (-\frac{(1 - {a''}_2) B_{23} K_3}{K_2} + \frac{(1 - a_2) (K_2 + B_{23} K_3)} {K_2}\right)}\right \}$;

with $\alpha_3 = d_1 + e_1 + \gamma_1 + r_1 (\frac{(1-a_1){N^\prime}_1}{K_1} + \frac{(1-{a^\prime}_1)B_{12}{N^\prime}_2}{K_1})$, $\alpha_4 = d_2 + e_2 + \gamma_2 + r_2(\frac{(1-a_2){N^\prime}_2}{K_2} + \frac{(1-{a^\prime}_2)B_{21}{N^\prime}_1}{K_2})$.

$R^0_i$, $i = 1, 2, 3$, is the reproduction number in single population $N_i$, while $R^0_{ij}$ is the basic reproduction number in two populations, $N_i$ and $N_j$, $i, j = 1, 2, 3$.

The choice of the basic reproductive numbers of different groups above references ^[29,36]

Lemma 3. The set $\; \Omega_0 = \{(S_1, I_1, R_1, S_2, I_2, R_2, S_3, I_3, R_3) | S_i, I_i, R_i \geq 0, i = 1, 2, 3, S_j + I_j + R_j \leq K_j + B_{j3} K_3, j = 1, 2, S_3 + I_3 + R_3\leq K_3\}$ is a positively invariant region for Model (2.1). Moreover, every trajectory of Model (2.1) is eventually staying in a compact subset of $\Omega_0$.

System (2.1) always has the equilibria $E_0^\prime = (0, 0, 0, 0, 0, 0, 0, 0, 0)$, $E_1^\prime = (K_1, 0, 0, 0, 0, 0, 0, 0, 0)$, $E_2^\prime = (0, 0, 0, 0, K_2, 0, 0, 0, 0)$, $E_3^\prime = (0, 0, 0, 0, 0, 0, 0, 0, K_3)$. $E_0^\prime$, $E_1^\prime$, $E_2^\prime$, and $E_3^\prime$ are always unstable.

If $R^0_1 > 1$, there exists the boundary equilibrium point $BE_1 = (S_1^\ast, I_1^\ast, R_1^\ast, 0, 0, 0, 0, 0, 0)$ with $N_2(0) = N_3(0) = 0$ by Lemma 3.

If $R^0_2 > 1$, there exists the boundary equilibrium point $BE_2 = (0, 0, 0, S_2^\ast, I_2^\ast, R_2^\ast, 0, 0, 0)$ with $N_1(0) = N_3(0) = 0$ by Lemma 3.

If $R^0_3 > 1$, there exists the boundary equilibrium point $BE_3 = (0, 0, 0, 0, 0, 0, S_3^\ast, I_3^\ast, R_3^\ast)$ with $N_1(0) = N_2(0) = 0$ by Lemma 3.

If $\frac{B_{12}K_2-K_1}{B_{12}B_{21}-1} > 0$ and $\frac{B_{21}K_1-K_2}{B_{12}B_{21}-1} > 0$, there exists the boundary equilibrium point $E_4^\prime = ({N^\prime}_1, 0, 0, {N^\prime}_2, 0, 0, 0, 0, 0) = (\frac{B_{12}K_2-K_1}{B_{12}B_{21}-1}, 0, 0, \frac{B_{21}K_1-K_2}{B_{12}B_{21}-1}, 0, 0, 0, 0, 0)$. $E_4^\prime$ is always unstable.

If $R^0_{12} > 1$ and $\frac{B_{12}K_2-K_1}{B_{12}B_{21}-1} > 0$, $\frac{B_{21}K_1-K_2}{B_{12}B_{21}-1} > 0$, there exists the boundary equilibrium point $BE_4 = ({S^\prime}_1, {I^\prime}_1, {R^\prime}_1, {S^\prime}_2, {I^\prime}_2, {R^\prime}_2, 0, 0, 0)$ ^[32]. Here we always assume that $I_i(0)\le N_i(0)$, $i = 1, 2, 3$. The infective compartment is just the part of human species. $BE_1$, $BE_2$, $BE_3$ and $BE_4$ are always unstable.

For $E_5$ and $E_6$ in model (4.1), we can also isolate the boundary equilibrium point $E_5^\prime$, $BE_5$, $E_6^\prime$ and $BE_6$ in model (2.1).

$E_5^\prime = \left(\overline{N_1}, 0, 0, 0, 0, 0, \overline{N_3}, 0, 0 \right) = \left(K_1 + B_{13} K_3, 0, 0, 0, 0, 0, K_3, 0, 0\right)$,

$BE_5 = \left(\overline{S_1}, \overline{\ I_1}, \overline{R_1}, 0, 0, 0, \overline{S_3}, \overline{I_3}, \overline{R_3}\right)$,

$E_6^\prime = \left(0, 0, 0, \overline{N_2}, 0, 0, \overline{N_3}, 0, 0\right) = \left(0, 0, 0, K_2+B_{23}K_3, 0, 0, K_3, 0, 0\right)$,

$BE_6 = \left(0, 0, 0, \overline{\overline{S_2}}, \overline{\overline{I_2}}, \overline{\overline{R_2}}, \overline{\overline{S_3}}, \overline{\overline{I_3}}, \overline{\overline{R_3}}\right)$.

The existence of boundary equilibrium $E_7^\prime$ and internal equilibrium $E_8^\prime$ depends on the conditions that $\frac{B_{12}K_2-K_1+(B_{12}B_{23}-B_{13})K_3}{B_{12}B_{21}-1} > 0$ and $\frac{B_{21}K_1-K_2+(B_{21}B_{13}-B_{23})K_3}{B_{12}B_{21}-1} > 0$.

If $B_{12}B_{21} < 1$, $B_{12}K_2-K_1+\left(B_{12}B_{23}-B_{13}\right)K_3 > 0$ and $B_{21}K_1-K_2+\left(B_{21}B_{13}-B_{23}\right)K_3 > 0$, there exists stable boundary equilibrium $E_7^\prime$ by Theorem 1.

However, the existence of internal equilibrium $E_8^\prime$ with infection is more complex. Inspired by the ^[32], we discuss the existence of $E_8^\prime$ by its persistence isolated from other equilibria.

The internal equilibrium $E_8^\prime$ of model (2.1) could satisfy the equations as follows:

Integrate $N_i (t) = S_i(t) + I_i(t) + R_i(t)$, $i = 1, 2, 3$ of equations in (4.3), then we get

From (4.4), we get $N_3 = K_3(1-\frac{e_3}{r_3}\frac{I_3}{N_3})$, $N_1 = \frac{B_{12}K_2-K_1+\left(B_{12}B_{23}-B_{13}\right)N_3+\frac{e_1}{r_1}\frac{I_1}{N_1}-B_{12}\frac{e_2}{r_2}\frac{I_2}{N_2}}{B_{12}B_{21}-1}$, and $N_1 = \frac{B_{21}K_1-K_2+\left(B_{21}B_{13}-B_{23}\right)N_3+\frac{e_2}{r_2}\frac{I_2}{N_2}-B_{21}\frac{e_1}{r_1}\frac{I_1}{N_1}}{B_{12}B_{21}-1}$. Set $Y_i: = \frac{I_i}{N_i}$, $i = 1, 2, 3$. Note that conditions in Lemma 3 imply that $0 \leq Y_i \leq 1$, since $0\leq I_i \leq N_i$. $r_3 > e_3$ implies the existence of steady state of $N_3$ and ${K_3(1-\frac{e_3}{r_3})\leq N}_3{\leq K}_3$. Next, we get the following results for the existence of internal equilibrium in the view of persistence theory.

Let

with

Theorem 2. For system (2.1), if $\frac{B_{12}K_2-K_1+(B_{12}B_{23}-B_{13})K_3}{B_{12}B_{21}-1} > 0$, $\frac{B_{21}K_1-K_2+(B_{21}B_{13}-B_{23})K_3}{B_{12}B_{21}-1} > 0$, $B_{12}B_{21} < 1$ and $R^0_Z > 1$, there exists an $\epsilon > 0$ such that $\lim{{inf}_{t\rightarrow\infty}}min\{I_1(t), I_2(t), I_3(t)\} > \epsilon$, for any solution with $N_1(0) > 0$, $N_2(0) > 0$, $N_3(0) > 0$ and $I_1(0) > 0$ or $I_2(0) > 0$ or $I_3(0) > 0$.

Proof. Define

$D_2$ and ${\widetilde{D}}_2$ are forward invariant.

Let $\Omega^\ast$ consist of equilibria $E_0^\prime$, $E_1^\prime$, $BE_1$, $E_2^\prime$, $BE_2$, $E_3^\prime$, $BE_3$, $E_4^\prime$, $BE_4$, $E_5^\prime$, $BE_5$, $E_6^\prime$, $BE_6$, and $E_7^\prime$. These equilibria cannot be chained to each other in $D_1$. By analyzing the flow in neighborhood of each equilibrium, it is easy to see that $\Omega^\ast$ is isolated in $D$ and $D_1$ is a uniform strong repeller for ${\widetilde{D}}_2$.

If the solution $x(t) = (S_1(t), I_1(t), S_2(t), I_2(t))$ of system (2.1) stays close to $E_0^\prime$, we have two cases:

● if $N_1(0) = N_2(0) = N_3(0) = 0$, then $N_1(t) = N_2(t) = N_3(t) = 0$;

● if $N_1(0) > 0\ or\ N_2(0) > 0\ or\ N_3(0) > 0$, then $N_1(t) > 0\ or\ N_2(t) > 0\ or\ N_3(t) > 0$. $E_0^\prime$ is isolated in D.

If $x(t)$ stays in a small neighborhood of $E_1^\prime$, we have three cases:

● if $I_1(0) = N_2(0) = N_3(0) = 0$, then $I_1(t) = N_2(t) = N_3(t) = 0$;

● if $N_2(0) > 0\ or\ N_3(0) > 0$, then $\ N_2(t) > 0\ or\ N_3(t) > 0$, since $\frac{B_{21}K_1-K_2+(B_{21}B_{13}-B_{23})K_3}{B_{12}B_{21}-1} > 0$, $B_{12}B_{21} < 1$ by Lemma (3);

● if $I_1(0) > 0$, then $I_1(t) > 0, \forall\ t > 0$. Since $(S_1(t), I_1(t), R_1(t))$ satisfying system (2.1) has no invariant subset other than $E_1^\prime$ in its neighborhood. $E_1^\prime$ is isolated in $D$.

Similarly, we can prove that $BE_1$, $E_2^\prime$, $BE_2$, $E_3^\prime$, $BE_3$, $E_4^\prime$, $BE_4$, $E_5^\prime$, $BE_5$, $E_6^\prime$, $BE_6$, and $E_7^\prime$ are isolated in $D$.

Using Proposition 4.3 in ^[31], we can prove that $D_1$ is a uniform weak repeller for ${\widetilde{D}}_2$; and using Theorem 4.5 in ^[31], we can prove that $D_1$ is a uniform strong repeller for ${\widetilde{D}}_2$.

Then we get that there exists an $\epsilon > 0$ such that

with $N_1(0) > 0$, $N_2(0) > 0$, $N_3(0) > 0$ and $I_1(0) > 0$ or $I_2(0) > 0$ or $I_3(0) > 0$.

More details have been shown in ^[32,33,31], we won't repeat the process in this paper.

Theorem 3. For system (2.1), if $\frac{B_{12}K_2-K_1+(B_{12}B_{23}-B_{13})K_3}{B_{12}B_{21}-1} > 0$, $\frac{B_{21}K_1-K_2+(B_{21}B_{13}-B_{23})K_3}{B_{12}B_{21}-1} > 0$, $B_{12}B_{21} < 1$ and $R^0_Z < 1$, the equilibrium $E_7^\prime$ is stable; and system (2.1) is not persistent with $N_1(0) > 0$, $N_2(0) > 0$, $N_3(0) > 0$ and $I_1(0) > 0$ or $I_2(0) > 0$ or $I_3(0) > 0$.

Corollary 1. (^[32,33]) If $\frac{B_{12}K_2-K_1+(B_{12}B_{23}-B_{13})K_3}{B_{12}B_{21}-1} > 0$, $\frac{B_{21}K_1-K_2+(B_{21}B_{13}-B_{23})K_3}{B_{12}B_{21}-1} > 0$, $B_{12}B_{21} < 1$ and $R^0_Z > 1$, there exists at least one internal equilibrium of system (2.1).

System (2.1) is subdivided from system (4.1). The local stability of $E_0^\prime$ is similar with $E_0$, which is unstable. Likewise, $E_1^\prime$, $BE_1$, $E_2^\prime$, $BE_2$, $E_3^\prime$, $BE_3$, $E_4^\prime$ and $BE_4$ are unstable. However, the local stability of $E_5^\prime$, $BE_5$, $E_6^\prime$, $BE_6$, $E_7^\prime$ and $E_8^\prime$ are more complex for the complexity of $E_5$, $E_6$ and $E_7$.

By Theorem 1 and the classical theory about basic reproductive number in ^[29], we can get the following results.

If $B_{21}K_1-K_2+(B_{21}B_{13}-B_{23})K_3 > 0$, $E_5$ is stable. So $BE_5$ is stable when $R^0_{13} < 1$, $E_5^\prime$ is when $R^0_{13} > 1$. If $B_{21}K_1-K_2+(B_{21}B_{13}-B_{23})K_3 < 0$, $E_5$ is unstable. So $E_5^\prime$, $BE_5$ are unstable.

If $B_{12}K_2-K_1+(B_{12}B_{23}-B_{13})K_3 > 0$, $E_6$ is stable. So $BE_6$ is stable when $R^0_{23} < 1$, $E_6^\prime$ is when $R^0_{23} > 1$. If $B_{12}K_2-K_1+(B_{12}B_{23}-B_{13})K_3 < 0$, $E_6$ is unstable. So $E_6^\prime$, $BE_6$ are unstable.

If $\frac{B_{12}K_2-K_1+(B_{12}B_{23}-B_{13})K_3}{B_{12}B_{21}-1} > 0$, $\frac{B_{21}K_1-K_2+(B_{21}B_{13}-B_{23})K_3}{B_{12}B_{21}-1} > 0$ and $B_{12}B_{21} > 1$, $E_7$ is a saddle point, which is unstable. So $E_7^\prime$, $E_8^\prime$ are unstable.

If $\frac{B_{12}K_2-K_1+(B_{12}B_{23}-B_{13})K_3}{B_{12}B_{21}-1} > 0$, $\frac{B_{21}K_1-K_2+(B_{21}B_{13}-B_{23})K_3}{B_{12}B_{21}-1} > 0$ and $B_{12}B_{21} < 1$, $E_7$ is globally asymptotically stable. So $E_7^\prime$ is stable when $R^0_Z < 1$, $E_8^\prime$ is stable when $R^0_Z > 1$.

Theorem 4. For system (2.1), the equilibria $E_0^\prime$, $E_1^\prime$, $E_2^\prime$ and $E_3^\prime$ always exist.

If $R^0_1 = \frac{\beta_{11}K_1}{b_1+\gamma_1+e_1-a_1r_1} > 1$ (or $R^0_2 = \frac{\beta_{22}K_2}{b_2+\gamma_2+e_2-a_2r_2} > 1$ or $R^0_3 = \frac{\beta_{33}K_3}{b_3+\gamma_3+e_3-a_3r_3} > 1)$, the equilibrium $BE_1 (or \ BE_2\ or \ BE_3)$ could exist.

If $\frac{B_{12}K_2-K_1}{B_{12}B_{21}-1} > 0$ and $\frac{B_{21}K_1-K_2}{B_{12}B_{21}-1} > 0$, the equilibrium $E_4^\prime$ exists.

If $R^0_{12} = \rho\left(\begin{matrix}\frac{\beta_{11}N_1}{b_1+\gamma_1+e_1-a_1r_1\frac{N_1}{K_1}} & \frac{\beta_{12}N_2}{b_2+\gamma_2+e_2-a_2r_2\frac{N_2}{K_2}}\\\frac{\beta_{21}N_1}{b_1+\gamma_1+e_1-a_1r_1\frac{N_1}{K_1}} & \frac{\beta_{22}N_2}{b_2+\gamma_2+e_2-a_2r_2\frac{N_2}{K_2}}\\\end{matrix}\right) > 1$, the equilibrium $BE_4$ could exist.

All these equilibria above are unstable.

The equilibrium $E_5^\prime$ and $E_6^\prime$ always exist.

As for the existence and the stability of $BE_5$, $BE_6$, $E_7^\prime$ and $E_8^\prime = (S_1^{\ast\ast}, I_1^{\ast\ast}, R_1^{\ast\ast}, S_2^{\ast\ast}, I_2^{\ast\ast}, R_2^{\ast\ast}, S_3^{\ast\ast}, I_3^{\ast\ast}, R_3^{\ast\ast})$, there are 4 cases:

● If $B_{12} > \frac{K_1+B_{13} K_3}{K_2+B_{23} K_3}$ and $B_{21} > \frac{K_2+B_{23} K_3}{K_1+B_{13} K_3}$, then $E_5^\prime$ is stable with $R^0_{13} < 1$, $E_6^\prime$ is stable with $R^0_{23} < 1$; $BE_5$ exists with $R^0_{13} > 1$, $BE_6$ exists with $R^0_{23} > 1$; $BE_5$ and $BE_6$ are stable; $E_7^\prime$ always exists; $E_8^\prime$ exists with $R^0_Z > 1$; $E_7^\prime$ and $E_8^\prime$ are unstable.

● If $B_{12} < \frac{K_1+B_{13} K_3}{K_2+B_{23} K_3}$ and $B_{21} > \frac{K_2+B_{23} K_3}{K_1+B_{13} K_3}$, then $E_5^\prime$ is stable with $R^0_{13} < 1$; $BE_5$ exists with $R^0_{13} > 1$ and $BE_5$ is stable; $E_6^\prime$ always exists; $BE_6$ exists with $R^0_{23} > 1$; $E_6^\prime$ and $BE_6$ are unstable; $E_7^\prime$ and $E_8^\prime$ do not exist.

● If $B_{12} > \frac{K_1+B_{13} K_3}{K_2+B_{23} K_3}$ and $B_{21} < \frac{K_2+B_{23} K_3}{K_1+B_{13} K_3}$, then $E_5^\prime$ always exists; $B_5$ exists with $R^0_{13} > 1$; $E_5^\prime$ and $BE_5$ are unstable; $E_6^\prime$ is stable with $R^0_{23} < 1$; $BE_6$ exists with $R^0_{23} > 1$ and $BE_6$ is stable; $E_7^\prime$ and $E_8^\prime$ do not exist.

● If $B_{12} < \frac{K_1+B_{13} K_3}{K_2+B_{23} K_3}$ and $B_{21} < \frac{K_2+B_{23} K_3}{K_1+B_{13} K_3}$, then $E_5^\prime$ and $E_6^\prime$ always exist; $BE_5$ exists with $R^0_{13} > 1$, $BE_6$ exists with $R^0_{23} > 1$; $E_5^\prime$, $E_6^\prime$, $BE_5$ and $BE_6$ are unstable; $E_7^\prime$ is stable with $R^0_Z < 1$; $E_8^\prime$ exists with $R^0_Z > 1$ and $E_8^\prime$ is stable.

Proof. The existence of equilibria in system (2.1) has been discussed in section 2.1. Next, we use Theorem 2 in ^[30] to discuss the reproduction numbers of system (2.1).

The Jacobian of $(I_1, I_2, I_3)$ is

with

Let $J = F-V$, $F$ be the rate of appearance of new infections in compartment $I$, $V$ be the rate of transfer of individuals out of compartment $I$. Then, we get

Set $R^0_Z = \rho(FV^{-1})$ at $E_7^\prime$, where $\rho(A$) denotes the spectral radius of a matrix $A$.

Then we get

with $\alpha_1 = d_1+e_1+\gamma_1+r_1(-\frac{(1-{a''}_1)B_{13}N_3^\ast}{K_1}+\frac{(1-a_1)N_1^\ast}{K_1}+\frac{(1-{a^\prime}_1)B_{12}N_2^\ast}{K_1})$, $\alpha_2 = d_2+e_2+\gamma_2+r_2(-\frac{(1-{a''}_2)B_{23}N_3^\ast}{K_2}+\frac{(1-a_2)N_2^\ast}{K_2}+\frac{(1-{a^\prime}_2)B_{21}N_1^\ast}{K_2})$.

Similarly, set

with $\alpha_1 = d_1 + e_1 + \gamma_1 + r_1(\frac {(1 - a_1) N'_1} {K_1} + \frac {(1-{a^\prime}_1) B_{12} {N'}_2} {K_1})$, $\alpha_2 = d_2+e_2+\gamma_2+r_2(\frac{(1-a_2){N'}_2}{K_2}+\frac{(1-{a^\prime}_2)B_{21}{N^\prime}_1}{K_2})$. Set

Using Theorem 2 in ^[30], we can prove Theorem 4.