Analyzing an Epidemic of Human Infections with Two Strains of Zoonotic Virus

Abstract

Due to the existence and variation of various viruses, an epidemic in which different strains spread at the same time will occur. here, an avian–human epidemic model with two strain viruses are established and analyzed. Both theoretical and simulation results reveal that the mixed infections intensify the epidemic and the dynamics become more complex and sensitive. There are six equilibria. The trivial equilibrium point is a high-order singular point and will undergo the transcritical bifurcations to bifurcate three equilibria. The existence and stability of equilibria mainly depend on five thresholds. A bifurcation portrait for the existence and stability of equilibria is presented. Simulations suggest that the key control measure is to develop the identification technology to eliminate the poultry infected with a high pathogenic virus preferentially, then the infected poultry with a low pathogenic virus in the recruitment and on farms. Controlling contact between human and poultry can effectively restrain the epidemic and controlling contagions in poultry can avoid great infection in humans.

1. Introduction

In recent years, the number of new and recurrent high-threat pathogens has increased, such as SARS-CoV, MERS-CoV, Ebola virus, Nipa virus, avian influenza virus and the latest, SARS-CoV-2. Among them, avian influenza and other zoonotic influenza have always been the great threat to human beings. Take human infection with H5N1 and H7N9 avian influenza as examples. In 1997, Hong Kong Special Administrative Region of China reported cases of human infection with high pathogenic H5N1 avian influenza virus. Since 2003, this avian influenza virus has spread from Asia to Europe and Africa, and is deeply rooted in poultry in some countries. Other subtypes of A(H5) avian influenza viruses also cause poultry epidemics and human infection. Regarding H7N9, the A(H7N9) virus has been detected in birds for many years and had not attracted people’s attention because of its low pathogenicity for birds and no previous record of human infections. However, the first human case was confirmed in March 2013 in China, before the epidemic of human infections with the A(H7N9) avian influenza virus broke out [1,2,3]. Successively, the same epidemic happened every autumn–winter, and so far there have been five outbreaks. The wave of 2016–2017 was an epidemic of human infections with both low and high pathogenic A(H7N9) avian influenza virus; that is, the virus had mutated and both the original virus and the mutated virus infected humans and poultry simultaneously. In fact, the gene sequencing analysis of the viruses isolated from two human A(H7N9) cases in January 2017 showed that the amino acid insertion mutations happened and the A(H7N9) virus had mutated into a high pathogenic virus for poultry [4]. Consequently, the number of cases and the fatality rate increased alarmingly. There were 1565 cases from 2013 to 2017 in total; however, the infection number was 766 in the wave of 2016–2017 [3]. It contributed almost half of all cases in those years. Not only that, in December 2016, there were 106 cases of human infections in China, and there were 192 human cases, 79 of which died in January 2017 [5]. This shows that the number of cases in a single month exceeds that of every past year and there were obviously more deaths. Infections during these five waves are presented in Figure 1.

As well as the number of cases, the geographical distribution of cases in the wave of 2016–2017 was also far more than any previous wave. There were 12, 12, 13, 14 and 27 provinces, cities and municipalities with infections, respectively. Details are listed in

Table A1. Provinces, cities and municipalities where human cases occurred.

Wave	Place
1st	Shanghai, Anhui, Jiangsu, Zhejiang, Henan, Beijing, (Taiwan), Fujian, Hunan, Jiangxi, Shangdong, Hebei (2013.7), Guangdong (2013.7)
2nd	New added: Guangxi, Jilin, Guizhou(back from Zhejiang), Hongkong, Xinjiang (2014.7); Not reappearing: Jiangxi, Henan, Hebei
3rd	New added: Hubei, Guizhou; Not reappearing: Guangxi, Jixin, Henan, Hebei; Resumed: Jiangxi
4th	New added: Tianjin, Liaoning; Not appearing: Xinjiang, Guizhou; Resumed: Hebei
5th	New added: Sichuan, Yunnan, Tibet, Shaanxi, Shanxi, Chongqing, Gansu, Inner Mongolia, Macao

in Appendix A.1 and in Figure 2. We can see from Figure 2 that infections began in the coastal cities of the southeast China and had a tendency to expand vertically and to the inland regions. It also can be seen that the epidemic scope of the first four waves is relatively limited and fixed, while in the fifth wave there were human cases in almost all of the country, which suggests that the epidemic had expanded extensively and unusually.

Multiple strains are circulating in poultry, which makes it possible for human beings to be infected by different strains at the same time. Moreover, the increased number and the expanded geographical distribution in the fifth wave of human infections with the A(H7N9) avian influenza emphasize that it is well worth studying the spread and control of an epidemic of human infections with two strains of zoonotic influenza. Avian influenza viruses are classified as low pathogenic avian influenza (LPAI) viruses and high pathogenic avian influenza (HPAI) viruses according to their ability to cause disease in poultry. In view of their great threat to public health, there has been much research concerning the corresponding issues of genetic analysis, epidemiology and disease control. Arunachalam [6] found out the selection pressure of each amino acid site with the purpose of helping to develop effective vaccines and drugs for the A(H7N9) virus. Zhang et al. [7] discussed the occurrence of human infections with the A(H7N9) virus. Xing et al. [8] investigated the recurrent factors of the outbreak. Guo et al. [9] implemented the global dynamics of an avian–human influenza model. Li et al. [10] also obtained global stabilities of an A(H7N9) transmission model. Authors of this article, in their previous research [11,12] paid attention to the transmission of the virus in the environment. These studies were mainly finished before 2017 and focused on human infection with the low pathogenic avian influenza (LPAI). As we emphasize the simultaneous transmission of multiple viruses, we will take human infections with both the low and high A(H7N9) avian influenza viruses as an example to carry out this investigation. There are several studies concerning the transmission of two strains. For example, Tuncer et al. [13] introduced the dynamics of low and high pathogenic avian influenza in wild and domestic bird populations but emphasized that LPAI and HPAI can coexist in both populations. Kuddus et al. [14] investigated a two-strain disease model to simulate the prevalence of drug-susceptible and drug-resistant disease strains. In the present investigation, more complex and sensitive dynamics will be given out and the key control measures in the case of human infection with two strains of zoonotic virus also will be proposed.

The organization of this paper is as follows. In Section 2, we establish a mathematical model of avian–human infections in the fifth wave in mainland China. In Section 3, the existence of equilibria are investigated. Thresholds and the bifurcation portrait are given. In Section 4, local and global stabilities of equilibria are discussed dynamically. In Section 5, numerical simulations are implemented to clarify the dynamics of the model and the effect of control measures, and to carry out some comparative studies on human infections with low pathogenic, high pathogenic and mixed A(H7N9) viruses. In Section 6, we end this paper with conclusions and discussions.

2. The Model

In the present paper, the density of susceptible poultry population is denoted by $S_a\left(t\right)$, the density of infected poultry with low pathogenic A(H7N9) virus is denoted by $E_a\left(t\right)$ because of very mild or no disease symptoms and the density of infected poultry with high pathogenic A(H7N9) virus is denoted by $I_a\left(t\right)$. $N_a\left(t\right)=S_a\left(t\right)+E_a\left(t\right)+I_a\left(t\right)$. Densities of susceptible, infected and recovered human population are denoted by $S_h\left(t\right)$, $I_h\left(t\right)$ and $R_h\left(t\right)$, respectively.

In view of biosafety, the implemented practical measures are to close the live poultry markets (LPMs) and to kill all the poultry within 3 kilometers of the outbreak, which would severely knock down the poultry industry and cause huge economic loss. In fact, the positive rates of H5, H7 and H9 in farms are far less than that in the live poultry markets. Yan et al. in [15], Zhu et al. in [16] and Cao et al. in [17] were all convinced that the positive rate in live poultry markets was much higher than that in poultry farms. Chen et al. in [18], noted that 44.4% of retail LPMs and 50.0% of wholesale LPMs were confirmed to be contaminated while no positive samples were detected from poultry farms. It seems that maintaining normal production can be considered in the outbreak. Therefore, we carry out our study with the premise of maintaining normal production. The following assumptions are made to establish the two-strain avian–human infection model:

Hypothesis 1 (H1).

Since we are considering normal production, the constant recruitment in poultry should be considered in the model. Moreover, in modern society, due to the expansion of people’s activities, migration is certain. Let $A>0$ and $B>0$ be the constant recruitment rate of poultry population and human population, respectively. Because the infected poultry with the high pathogenic A(H7N9) virus is easy to distinguish and the infected poultry with the low pathogenic A(H7N9) virus continues to show no disease symptoms, we consider that the susceptible and the infected poultry with the low pathogenic A(H7N9) virus are all included in the recruitment of poultry. $0\le a\le 1$ represents the proportion of the latter in the recruitment.

Hypothesis 2 (H2).

The infected poultry with the low pathogenic A(H7N9) virus and the infected poultry with the high pathogenic A(H7N9) virus are all infectious, because they are carrying a virus.

Hypothesis 3 (H3).

Let $d>0$ be the output rate in the poultry industry including natural deaths and sales. $\rho >0$ is the natural death rate of the human population. $\xi >0$ is the additional death rate of $I_a\left(t\right)$ due to infection with the high pathogenic A(H7N9) virus. $\delta >0$ is the extra death rate of human cases caused by disease. $\gamma >0$ represents the recovery rate of human cases.

Hypothesis 4 (H4).

Suppose infections in poultry are due to contact. A saturated incidence function $x/(1+\alpha x)$ is used to describe the contact behaviors in poultry for the reason that the poultry population is large and densely populated. Because

$${\left(\frac{x}{1+\alpha x}\right)}^{\prime }=\frac{1}{1+\alpha x}-\frac{\alpha x}{{(1+\alpha x)}^2}=\frac{1}{{(1+\alpha x)}^2}>0,{\left(\frac{x}{1+\alpha x}\right)}^{″}=-\frac{2\alpha }{{(1+\alpha x)}^3}<0$$

and $x/(1+\alpha x)\to 1/\alpha$ for the larger x, it is a monotonically increasing convex function and will reach saturation for the larger x. Therefore, it can express the saturated contact effect in poultry. x only represents the argument of the function. The corresponding infection rate functions are ${\beta }_1{S}_a{E}_a/(1+{\alpha }_1{E}_a)$ and ${\beta }_2{S}_a{I}_a/(1+{\alpha }_2{I}_a)$.

Hypothesis 5 (H5).

As with human infections, many research results have confirmed that human behaviors, such as personal hygienic practice, reduced exposure to live poultry and disinfection of live poultry markets, will directly reduce the incidence of disease [19,20,21,22], because human behaviors will be affected by propaganda and psychosocial effects. At the beginning, human infections will expand due to people’s response lags. With the increase of human cases, the authorities concerned and mass media would respond and inform the public. Benefiting from protective awareness and behavioral changes, the incidence will decrease. Therefore, the incidence function $g\left(x\right)$ from poultry to humans should satisfy $g\left(0\right)=0$, $g^{\prime }\left(0\right)>0$ and $g\left(x\right)$ will decrease when x is relatively large. Let $g\left(x\right)=x/(1+\nu x^2)$, there is

$$g^{\prime }\left(x\right)=\frac{1-\nu x^2}{{(1+\nu x^2)}^2},g^{″}\left(x\right)=-\frac{4\nu x}{{(1+\nu x^2)}^3}<0.$$

Obviously, $g\left(0\right)=0,g^{\prime }\left(0\right)>0$, and $g^{\prime }\left(x\right)>0$ for the small x while $g^{\prime }\left(x\right)<0$ for the larger x. In fact, $g^{\prime }\left(x\right)$ will obtain the maximum value at $x=1/\sqrt{\nu }$ which means the infections will increase at the beginning and then decrease after a certain time. The corresponding infection rate functions are ${\eta }_1{E}_a{S}_h/(1+{\nu }_1{E}_a^2)$ and ${\eta }_2{I}_a{S}_h/(1+{\nu }_2{I}_a^2)$.

Based on the above analysis, the disease transmission mechanism is presented in Figure 3 and we establish the two-strain avian–human infection model as follows.

$$\left\{\begin{array}{c}\begin{array}{ccc}{\displaystyle \frac{\mathrm{d}{S}_a}{\mathrm{d}t}}\hfill & =\hfill & (1-a)A-{\beta }_1{S}_a{E}_a/(1+{\alpha }_1{E}_a)-{\beta }_2{S}_a{I}_a/(1+{\alpha }_2{I}_a)-d{S}_a,\hfill \\ {\displaystyle \frac{\mathrm{d}{E}_a}{\mathrm{d}t}}\hfill & =\hfill & aA+{\beta }_1{S}_a{E}_a/(1+{\alpha }_1{E}_a)-d{E}_a,\hfill \\ {\displaystyle \frac{\mathrm{d}{I}_a}{\mathrm{d}t}}\hfill & =\hfill & {\beta }_2{S}_a{I}_a/(1+{\alpha }_2{I}_a)-d{I}_a-\xi I_a,\hfill \\ {\displaystyle \frac{\mathrm{d}{S}_h}{\mathrm{d}t}}\hfill & =\hfill & B-{\eta }_1{E}_a{S}_h/(1+{\nu }_1{E}_a^2)-{\eta }_2{I}_a{S}_h/(1+{\nu }_2{I}_a^2)-\rho S_h,\hfill \\ {\displaystyle \frac{\mathrm{d}{I}_h}{\mathrm{d}t}}\hfill & =\hfill & {\eta }_1{E}_a{S}_h/(1+{\nu }_1{E}_a^2)+{\eta }_2{I}_a{S}_h/(1+{\nu }_2{I}_a^2)-\delta I_h-\rho I_h-\gamma I_h,\hfill \\ {\displaystyle \frac{\mathrm{d}{R}_h}{\mathrm{d}t}}\hfill & =\hfill & \gamma I_h-\rho R_h,\hfill \end{array}\hfill \end{array}\right.$$

(1)

with the initial conditions $S_a\left(0\right)=S_{a_0}>0,E_a\left(0\right)=E_{a_0}\ge 0,I_a\left(0\right)=I_{a_0}>0,S_h\left(0\right)=S_{h_0}>0,I_h\left(0\right)=I_{h_0}\ge 0,R_h\left(0\right)=R_{h_0}\ge 0$.

It is obvious that all the solutions initiating in ${\mathbb{R}}_{+}^6$ exist continuously for all $t\ge 0$ and are unique. Where ${\mathbb{R}}_{+}^6=\{(x,y,z,u,v,w)\in {\mathbb{R}}^6:x\ge 0,y\ge 0,z\ge 0,u\ge 0,v\ge 0,w\ge 0\}.$

Firstly, the biological validity should be guaranteed. Let $N_a\left(t\right)=S_a\left(t\right)+E_a\left(t\right)+I_a\left(t\right),N_h\left(t\right)=S_h\left(t\right)+I_h\left(t\right)+R_h\left(t\right).$ Then, we have

$$\frac{\mathrm{d}X}{\mathrm{d}t}=A-d{N}_a-\xi I_a\le A-d{N}_a.$$

$$\frac{\mathrm{d}Y}{\mathrm{d}t}=B-\rho N_h-\delta I_h\le B-\rho N_h.$$

By the differential inequality, we can obtain

$$0\le N_a\left(t\right)\le \frac{A}{d}+N_a(S_{a_0},E_{a_0},I_{a_0})e^{-dt},0\le N_h\left(t\right)\le \frac{B}{\rho }+N_h(S_{h_0},I_{h_0},R_{h_0})e^{-\rho t}.$$

Thus, $0\le N_a\le A/d$ and $0\le N_h\le B/\rho$ as $t\to +\infty .$

Therefore, all the solutions of System (1) that initiate in $R_{+}^6$ are confined in the region

$$\Omega =\{(S_a,E_a,I_a,S_h,I_h,R_h)\in {\mathbb{R}}_{+}^6:0<S_a+E_a+I_a<A/d+ϵ,0<S_h+I_h+R_h<B/\rho +ϵ\}$$

for any $ϵ>0$ and $t\to +\infty .$ Hence, all the solutions of System (1) are ultimately uniform-bounded and the system is dissipative.

3. Existence of Equilibria

For a dynamic system, an equilibrium point corresponds to its constant solution. In epidemiology, equilibria can be classified as the disease-free equilibrium point and the disease-existence equilibrium point. Let us discuss the existence of equilibria firstly. As the constant solution of a differential equation $\dot{x}=f\left(x\right)$, it can be obtained by solving the equation $f\left(x\right)=0$. Therefore, the Equilibria of System (1) depend on that of the system’s avian-subsystem which is as follows.

$$\left\{\begin{array}{c}\begin{array}{ccc}{\displaystyle \frac{\mathrm{d}{S}_a}{\mathrm{d}t}}\hfill & =\hfill & (1-a)A-{\beta }_1{S}_a{E}_a/(1+{\alpha }_1{E}_a)-{\beta }_2{S}_a{I}_a/(1+{\alpha }_2{I}_a)-d{S}_a,\hfill \\ {\displaystyle \frac{\mathrm{d}{E}_a}{\mathrm{d}t}}\hfill & =\hfill & aA+{\beta }_1{S}_a{E}_a/(1+{\alpha }_1{E}_a)-d{E}_a,\hfill \\ {\displaystyle \frac{\mathrm{d}{I}_a}{\mathrm{d}t}}\hfill & =\hfill & {\beta }_2{S}_a{I}_a/(1+{\alpha }_2{I}_a)-d{I}_a-\xi I_a,\hfill \end{array}\hfill \end{array}\right.$$

(2)

Natually, all solutions of System (2) that initiate in $R_{+}^3$ are confined in the region

$${\Omega }_a=\{(S_a,E_a,I_a)\in {\mathbb{R}}_{+}^3:0<S_a+E_a+I_a<A/d+ϵ\}$$

for any $ϵ>0$ and $t\to +\infty ,$ and System (2) is dissipative.

In view of solving equilibria, the third question of System (2) hints there is $I_a=0$. Therefore, in the case of $a=0$ there is an axis equilibrium point $M_0(A/d,0,0)$ and a boundary equilibrium point $M_1(A/d-E_a^{\left(1\right)},E_a^{\left(1\right)},0)$ if $A{\beta }_1>d^2$, where $E_a^{\left(1\right)}=(A{\beta }_1-d^2)/d({\alpha }_{1}d+{\beta }_1)$. It is obvious that there is $M_1\to M_0$ as $A{\beta }_1\to d^2$. While in the case of $a\ne 0$, by the direct calculation, we can obtain

$$({\alpha }_{1}d+{\beta }_1)E_a^2+(-\frac{A{\beta }_1}{d}+d-aA{\alpha }_1)E_a-aA=0.$$

(3)

Denote

$$h\left(x\right)=({\alpha }_{1}d+{\beta }_1)x^2+(-\frac{A{\beta }_1}{d}+d-aA{\alpha }_1)x-aA,$$

then $h\left(0^{+}\right)=-aA<0$ and

$$h\left({\left(\frac{A}{d}\right)}^{-}\right)=\frac{{\alpha }_1{A}^2}{d}(1-a)+A(1-a)>0.$$

According to the properties of quadratic function, Equation (3) has an unique positive root $\widetilde{E_a^{\left(1\right)}}\in (0,A/d)$. Therefore, there is another boundary equilibrium point $\widetilde{M_1}(A/d-\widetilde{E_a^{\left(1\right)}},\widetilde{E_a^{\left(1\right)}},0)$. Equation (3) also hints $\widetilde{M_1}\to M_1$ as $a\to 0$. However, the existence of $M_1$ is conditional. $\widetilde{M_1}$’s unconditional existence and $M_1$’s conditional existence indicate that it is possible to avoid the outbreak in poultry if we can put an end to the infected poultry in recruitment.

Suppose $I_a\ne 0$, then in the case of $a=0$ there is a boundary equilibrium point

$$M_2(\frac{A{\alpha }_2+d+\xi }{d{\alpha }_2+{\beta }_2},0,\frac{A{\beta }_2-d(d+\xi )}{(d+\xi )(d{\alpha }_2+{\beta }_2)})$$

under the condition of $A{\beta }_2>d(d+\xi )$. It is obvious that $M_2\to M_0$ as $A{\beta }_2\to d(d+\xi )$. If also $E_a\ne 0$, by the detailed calculation, we can obtain

$$I_a^{\left(3\right)}=\frac{(d+A{\alpha }_1){\beta }_2-(d+\xi )({\beta }_1+d{\alpha }_1)}{(d+\xi )({\alpha }_2{\beta }_1+{\alpha }_1{\beta }_2+d{\alpha }_1{\alpha }_2)},E_a^{\left(3\right)}=\frac{(d+\xi +A{\alpha }_2){\beta }_1-d({\beta }_2+d{\alpha }_2)}{d({\alpha }_2{\beta }_1+{\alpha }_1{\beta }_2+d{\alpha }_1{\alpha }_2)},$$

then

$$E_a^{\left(3\right)}>0,I_a^{\left(3\right)}>0if\left\{\begin{array}{c}(d+A{\alpha }_1){\beta }_2>(d+\xi )({\beta }_1+d{\alpha }_1),\hfill \\ (d+\xi +A{\alpha }_2){\beta }_1>d({\beta }_2+d{\alpha }_2).\hfill \end{array}\right.$$

From the expression of Equation (2), it can be seen that

$$\left\{\begin{array}{l}{E}_a^{\left(3\right)}=({\beta }_1{S}_a^{\left(3\right)}/d-1)/{\alpha }_1,\\ I_a^{\left(3\right)}=[{\beta }_2{S}_a^{\left(3\right)}/(d+\xi )-1]/{\alpha }_2,\end{array}\right.$$

then it should be $S_a^{\left(3\right)}>d/{\beta }_1$ and $S_a^{\left(3\right)}>(d+\xi )/{\beta }_2$. Additionally, $S_a^{\left(3\right)}<A/d$ also should be satisfied. Therefore,

$$\frac{A}{d}>S_a^{\left(3\right)}>\frac{d}{{\beta }_1},\frac{A}{d}>S_a^{\left(3\right)}>\frac{d+\xi }{{\beta }_2}.$$

That is, there is a positive equilibrium point

$$M_3(\frac{(d+\xi )(1+{\alpha }_2{I}_a^{\left(3\right)})}{{\beta }_2},E_a^{\left(3\right)},I_a^{\left(3\right)})if\left\{\begin{array}{c}A{\beta }_1>d^2,A{\beta }_2>d(d+\xi ),\hfill \\ (d+A{\alpha }_1){\beta }_2>(d+\xi )({\beta }_1+d{\alpha }_1),\hfill \\ (d+\xi +A{\alpha }_2){\beta }_1>d({\beta }_2+d{\alpha }_2).\hfill \end{array}\right.$$

In the case of $a=0$, the equilibrium point also can be obtained from the second and third equations of System (2) by eliminating $S_a$. Then there is

$$\left\{\begin{array}{c}{E}_a^{\left(3\right)}=[(1+{\alpha }_2{I}_a^{\left(3\right)})·{\beta }_1(d+\xi )/{\beta }_{2}d-1]/{\alpha }_1,\hfill \\ I_a^{\left(3\right)}=[(1+{\alpha }_1{E}_a^{\left(3\right)})·{\beta }_{2}d/{\beta }_1(d+\xi )-1]/{\alpha }_2.\hfill \end{array}\right.$$

And

$$I_a^{\left(3\right)}>0if(d+A{\alpha }_1){\beta }_2>(d+\xi )({\beta }_1+d{\alpha }_1).$$

Therefore,

$$E_a^{\left(3\right)}>\frac{1}{{\alpha }_1}[\frac{{\beta }_1(d+\xi )}{{\beta }_{2}d}-1],$$

thus $E_a^{\left(3\right)}>0$ if ${\beta }_1(d+\xi )/{\beta }_{2}d\ge 1$. Hence $M_3$ exists if

$$\frac{{\beta }_1+d{\alpha }_1}{A{\alpha }_1+d}<\frac{{\beta }_2}{d+\xi }\le \frac{{\beta }_1}{d}.$$

Similarly,

$$E_a^{\left(3\right)}>0if(d+\xi +A{\alpha }_2){\beta }_1>d({\beta }_2+d{\alpha }_2)$$

and

$$I_a^{\left(3\right)}>\frac{1}{{\alpha }_2}[\frac{{\beta }_{2}d}{{\beta }_1(d+\xi )}-1].$$

Therefore, $M_3$ exists if

$$\frac{{\beta }_2+d{\alpha }_2}{A{\alpha }_2+d+\xi }<\frac{{\beta }_1}{d}\le \frac{{\beta }_2}{d+\xi }.$$

In the case of $I_a\ne 0$ and $a\ne 0$, a direct calculation yields

$$d[{\alpha }_1({\beta }_2+{\alpha }_{2}d)+{\alpha }_2{\beta }_1]E_a^2+[({\beta }_2+{\alpha }_{2}d)(d-{\alpha }_{1}aA)-{\beta }_1(A{\alpha }_2+d+\xi )]E_a-aA({\beta }_2+{\alpha }_{2}d)=0$$

(4)

Denote

$$g\left(x\right)=d[{\alpha }_1({\beta }_2+{\alpha }_{2}d)+{\alpha }_2{\beta }_1]x^2+[({\beta }_2+{\alpha }_{2}d)(d-{\alpha }_{1}aA)-{\beta }_1(A{\alpha }_2+d+\xi )]x-aA({\beta }_2+{\alpha }_{2}d),$$

then $g\left(0^{+}\right)=-aA({\beta }_2+{\alpha }_{2}d)<0,$

$$\begin{array}{ccc}g\left({\left(\frac{A}{d}\right)}^{-}\right)\hfill & =\hfill & (1-a)A^2/d·{\alpha }_1({\beta }_2+{\alpha }_{2}d)+(1-a)A({\beta }_2+{\alpha }_{2}d)-A{\beta }_1(d+\xi )/d\hfill \\ & =\hfill & (1-a)A/d·({\beta }_2+{\alpha }_{2}d)(A{\alpha }_1+d)-A{\beta }_1(d+\xi )/d\hfill \\ & =\hfill & A/d·[(1-a)({\beta }_2+{\alpha }_{2}d)(A{\alpha }_1+d)-{\beta }_1(d+\xi )].\hfill \end{array}$$

It is obvious that $g\left({\left(\frac{A}{d}\right)}^{-}\right)>0$ if

$$(1-a)(A{\alpha }_1+d)({\beta }_2+{\alpha }_{2}d)>{\beta }_1(d+\xi ).$$

Of course, $g\left({\left(\frac{A}{d}\right)}^{-}\right)>0$ if

$$(1-a)d({\beta }_2+{\alpha }_{2}d)>{\beta }_1(d+\xi )or(1-a)A{\alpha }_1({\beta }_2+{\alpha }_{2}d)>{\beta }_1(d+\xi ).$$

Therefore, Equation (4) has an unique positive root $E_a^{*}\in (0,A/d)$. Furthermore

$$g\left(\frac{aA}{d}\right)=\frac{aA}{d}[(a-1)A{\alpha }_2{\beta }_1-{\beta }_1(d+\xi )]<0.$$

That is $E_a^{*}\in (aA/d,A/d)$. Then, there is a positive equilibrium point $M^{*}(S_a^{*},E_a^{*},I_a^{*})$ with $S_a^{*}=(d{E}_a^{*}-aA)(1+{\alpha }_1{E}_a^{*})/{\beta }_1{E}_a^{*}>0$ and $I_a^{*}={\beta }_2{S}_a^{*}/{\alpha }_2(d+\xi )-1/{\alpha }_2$. Therefore, it still requires $A/d>S_a^{*}>(d+\xi )/{\beta }_2$. Thus, $M^{*}$ exists if

$$A{\beta }_2>d(d+\xi )and(1-a)(A{\alpha }_1+d)({\beta }_2+{\alpha }_{2}d)>{\beta }_1(d+\xi ).$$

Equation (4) also hints $M^{*}\to M_3$ as $a\to 0$.

Denote

$$R_1=\frac{A{\beta }_1}{d^2},R_2=\frac{A{\beta }_2}{d(d+\xi )},R_3=\frac{(A{\alpha }_1+d){\beta }_2}{({\beta }_1+d{\alpha }_1)(d+\xi )},$$

$$\widetilde{R_3}=\frac{(A{\alpha }_2+d+\xi ){\beta }_1}{d({\beta }_2+d{\alpha }_2)},R_4=\frac{(1-a)(A{\alpha }_1+d)({\beta }_2+{\alpha }_{2}d)}{{\beta }_1(d+\xi )}.$$

Noticing that if $R_1=1$ and $R_2=1$, there are

$$R_3=\frac{(A{\alpha }_1+d){\beta }_2}{({\beta }_1+d{\alpha }_1)(d+\xi )}=\frac{(A{\alpha }_1+d){\beta }_2·A}{({\beta }_1+d{\alpha }_1)(d+\xi )·A}=\frac{(A{\alpha }_1+d)d}{A({\beta }_1+d{\alpha }_1)}=\frac{A{\alpha }_{1}d+d^2}{A{\alpha }_{1}d+A{\beta }_1}=1,$$

$$\widetilde{R_3}=\frac{(A{\alpha }_2+d+\xi ){\beta }_1}{d({\beta }_2+d{\alpha }_2)}=\frac{(A{\alpha }_2+d+\xi ){\beta }_{1}A}{({\beta }_2+d{\alpha }_2)dA}=\frac{(A{\alpha }_2+d+\xi )d}{A({\beta }_2+d{\alpha }_2)}=\frac{A{\alpha }_{2}d+d(d+\xi )}{A{\alpha }_{2}d+A{\beta }_2}=1.$$

Then $E_a^{\left(3\right)}\to 0,I_a^{\left(3\right)}\to 0,M_3\to M_0$ as $R_1\to 1$ and $R_2\to 1$.

Summarizing up, the equilibria and their existence conditions are listed in

Table 1. Existence of equilibria.

Case	Equilibria	Threshold	Condition Characteristic	Remarks
$a=0$	$M_0(\frac{A}{d},0,0)$	None	Unconditional existence
	$M_1(S_a^{\left(1\right)},E_a^{\left(1\right)},0)$	$R_1=\frac{A{\beta }_1}{d^2}>1$	$\frac{d}{A}<\frac{{\beta }_1}{d}$	For larger ${\beta }_1$, $M_1$ exists.
	$M_2(S_a^{\left(2\right)},0,I_a^{\left(2\right)})$	$R_2=\frac{A{\beta }_2}{d(d+\xi )}>1$	$\frac{d}{A}<\frac{{\beta }_2}{d+\xi }$	For larger ${\beta }_2$, $M_2$ exists.
	$M_3(S_a^{\left(3\right)},E_a^{\left(3\right)},I_a^{\left(3\right)})$	$R_1>1$, $R_2>1$	$\frac{d}{A}<\frac{{\beta }_1}{d}$, $\frac{{\beta }_2+d{\alpha }_2}{A{\alpha }_2+d+\xi }<\frac{{\beta }_1}{d}$	For larger ${\beta }_1$ and ${\beta }_2$,
		$R_3>1$, $\widetilde{R_3}>1$	$\frac{d}{A}<\frac{{\beta }_2}{d+\xi }$, $\frac{{\beta }_1+d{\alpha }_1}{A{\alpha }_1+d}<\frac{{\beta }_2}{d+\xi }$	$M_3$ may exist.
		or $R_3>1,R_1\ge R_2$	$\frac{{\beta }_1+d{\alpha }_1}{A{\alpha }_1+d}<\frac{{\beta }_2}{d+\xi }<\frac{{\beta }_1}{d}$
		or $\widetilde{R_3}>1,R_2\ge R_1$	$\frac{{\beta }_2+d{\alpha }_2}{A{\alpha }_2+d+\xi }<\frac{{\beta }_1}{d}<\frac{{\beta }_2}{d+\xi }$
		or $R_1>1,R_3>1,\widetilde{R_3}>1$	$R_3<R_2$ if $R_1>1$
		or $R_2>1,\widetilde{R_3}>1,R_3>1$	$\widetilde{R_3}<R_1$ if $R_2>1$
$a\ne 0$	$\widetilde{M_1}(\widetilde{S_a^{\left(1\right)}},\widetilde{E_a^{\left(1\right)}},0)$	None	Unconditional existence
	$M^{*}(S_a^{*},E_a^{*},I_a^{*})$	$R_2>1,R_4>1$	$\frac{{\beta }_1}{A{\alpha }_1+d}<\frac{(1-a)({\beta }_2+{\alpha }_{2}d)}{d+\xi }$	$R_4=\frac{(1-a)(A{\alpha }_1+d)({\beta }_2+{\alpha }_{2}d)}{{\beta }_1(d+\xi )}$
			$\frac{d}{A}<\frac{{\beta }_2}{d+\xi }$	For larger ${\beta }_2$, $M^{*}$ exists.

Where $\left\{\right({\beta }_1,{\beta }_2\left)\right|R_1>1,R_2>1,R_3>1,\widetilde{R_3}>1\}=\{({\beta }_1,{\beta }_2)|R_1>1,R_3>1,\widetilde{R_3}>1\}=\left\{\right({\beta }_1,{\beta }_2\left)\right|R_2>1,\widetilde{R_3}>1,R_3>1\}=\{({\beta }_1,{\beta }_2)|R_3>1,R_1\ge R_2\}\cup \left\{\right({\beta }_1,{\beta }_2\left)\right|\widetilde{R_3}>1,R_2\ge R_1\}$.

, and the corresponding conclusions are as follows:

Theorem 1.

(i) 

Suppose $a=0$, System (2) always has an axis equilibrium point $M_0(A/d,0,0)$;

And, there also exist two boundary equilibria $M_1(A/d-E_a^{\left(1\right)},E_a^{\left(1\right)},0)$ if $R_1>1$ and $M_2((A{\alpha }_2+d+\xi )/(d{\alpha }_2+{\beta }_2),0,(A{\beta }_2-d(d+\xi ))/(d+\xi )(d{\alpha }_2+{\beta }_2))$ if $R_2>1$, where $E_a^{\left(1\right)}=(A{\beta }_1-d^2)/d({\alpha }_{1}d+{\beta }_1)$;

And, there is a positive equilibrium point $M_3((d+\xi )(1+{\alpha }_2{I}_a^{\left(3\right)})/{\beta }_2,E_a^{\left(3\right)},I_a^{\left(3\right)})$ if $R_1>1,R_2>1,R_3>1$ and $\widetilde{R_3}>1$, where $I_a^{\left(3\right)}=[(d+A{\alpha }_1){\beta }_2-(d+\xi )({\beta }_1+d{\alpha }_1)]/(d+\xi )({\alpha }_2{\beta }_1+{\alpha }_1{\beta }_2+d{\alpha }_1{\alpha }_2),E_a^{\left(3\right)}=[(d+\xi +A{\alpha }_2){\beta }_1-d({\beta }_2+d{\alpha }_2)]/d({\alpha }_2{\beta }_1+{\alpha }_1{\beta }_2+d{\alpha }_1{\alpha }_2)$.

(ii) 

Suppose $a\ne 0$, System (2) has a boundary equilibrium point $\widetilde{M_1}(\frac{A}{d}-\widetilde{E_a^{\left(1\right)}},\widetilde{E_a^{\left(1\right)}},0)$, where $\widetilde{E_a^{\left(1\right)}}$ is the unique positive root of Equation (3), and have a positive equilibrium point $M^{*}$ if $R_2>1$ and $R_4>1$.

(iii) 

$\widetilde{M_1}\to M_1$ and $M^{*}\to M_3$ as $a\to 0$.

(iv) 

$M_1\to M_0$ if $R_1\to 1$; $M_2\to M_0$ if $R_2\to 1$; $M_3\to M_0$ if $R_1\to 1$ and $R_2\to 1$.

Corollary 1.

As $a=0$, System (2) has a positive equilibrium point $M_3$ if $R_3>1,R_1\ge R_2$ or $\widetilde{R_3}>1,R_2\ge R_1$.

Corollary 2.

$M_0$ undergoes the transcritical bifurcations to bifurcate equilibria $M_1,M_2$ and $M_3$ as $R_{1,2}$ increase from $R_{1,2}=1$ to $R_{1,2}>1$.

By directly calculation, we also find that $R_3=R_2$ if $R_1=1$; $R_3<R_2$ if $R_1>1$; $\widetilde{R_3}=R_1$ if $R_2=1$; $\widetilde{R_3}<R_1$ if $R_2>1$. The calculations are as follows:

$$if{R}_1\ge 1,R_3=\frac{(A{\alpha }_1+d){\beta }_2}{({\beta }_1+d{\alpha }_1)(d+\xi )}=\frac{(A{\alpha }_1+d){\beta }_{2}A}{(A{\beta }_1+Ad{\alpha }_1)(d+\xi )}\le \frac{(A{\alpha }_1+d)A{\beta }_2}{d(d+A{\alpha }_1)(d+\xi )}=R_2,$$

$$if{R}_2\ge 1,\widetilde{R_3}=\frac{(A{\alpha }_2+d+\xi ){\beta }_1}{d({\beta }_2+d{\alpha }_2)}=\frac{(A{\alpha }_2+d+\xi )A{\beta }_1}{d(A{\beta }_2+Ad{\alpha }_2)}\le \frac{(A{\alpha }_2+d+\xi )A{\beta }_1}{d^2(d+\xi +A{\alpha }_2)}=R_1.$$

Thus, there is

Corollary 3.

Suppose $a=0$, System (2) has a positive equilibrium point $M_3$ if $R_1>1,R_3>1,\widetilde{R_3}>1$ or $R_2>1,\widetilde{R_3}>1,R_3>1$.

As for these thresholds, there are the following relationships.

$$R_1>1\iff \frac{d}{A}<\frac{{\beta }_1}{d},R_2>1\iff \frac{d}{A}<\frac{{\beta }_2}{d+\xi },R_3>1\iff \frac{{\beta }_1+d{\alpha }_1}{A{\alpha }_1+d}<\frac{{\beta }_2}{d+\xi },$$

$$\widetilde{R_3}>1\iff \frac{{\beta }_2+d{\alpha }_2}{A{\alpha }_2+d+\xi }<\frac{{\beta }_1}{d},R_4>1\iff \frac{{\beta }_1}{A{\alpha }_1+d}<\frac{(1-a)({\beta }_2+{\alpha }_{2}d)}{d+\xi }.$$

Details of the existence of equilibria are presented in Table 1.

The above results indicate that $M_1$ and $M_2$ exist with the larger ${\beta }_1$ and ${\beta }_2$, respectively. $M^{*}$ will exist for larger ${\beta }_2$, but the necessary existence conditions of $M_3$ may be stronger. In other words, the existence of $M_3$ is sensitive to parameters. Additionally, although $\widetilde{M_1}\to M_1$ and $M^{*}\to M_3$ as $a\to 0$, which we can detect from Equations (3) and (4), both the satisfied conditions are stronger. The reason is

$$\frac{{\beta }_1}{A{\alpha }_1+d}<\frac{{\beta }_1+d{\alpha }_1}{A{\alpha }_1+d}<\frac{{\beta }_2}{d+\xi }<\frac{{\beta }_2+{\alpha }_{2}d}{d+\xi }.$$

(5)

The unequal relationship between the two ends of (5) corresponds to $R_4>1$ at $a=0$, while the intermediate unequal relationship of (5) corresponds to $R_3>1$. Moreover, $\widetilde{M_1}$ exists unconditionally while $M_1$ exists if $R_1>1$. In other words, these equilibria more likely to not exist if $a=0$. It also indicates that it will be beneficial to disease control if we can eliminate the infected poultry in the recruitment.

Let $M(S_a,E_a,I_a)$ be an equilibrium point of System (2), a direct calculation yields a corresponding equilibrium point $\overline{M}(S_a,E_a,I_a,S_h,I_h,R_h)$ of System (1), where

$$S_h=B/[\frac{{\eta }_1{E}_a}{1+{\nu }_1{E}_a^2}+\frac{{\eta }_2{I}_a}{1+{\nu }_2{I}_a^2}+\rho ],I_h=[\frac{{\eta }_1{E}_a{S}_h}{1+{\nu }_1{E}_a^2}+\frac{{\eta }_2{I}_a{S}_h}{1+{\nu }_2{I}_a^2}]/(\gamma +\delta +\rho ),R_h=\frac{\gamma I_h}{\rho }.$$

In the following, we will hold names of equilibria unchanged. Correspondingly, the existence of equilibria in System (1) is as follows.

Theorem 2.

(i) 

Suppose $a=0$, System (1) always has an axis equilibrium point $M_0$; and, there also exists two boundary equilibria $M_1$ if $R_1>1$ and $M_2$ if $R_2>1$; and, there is a positive equilibrium point $M_3$ if $R_1>1,R_2>1,R_3>1$ and $\widetilde{R_3}>1$.

(ii) 

Suppose $a\ne 0$, System (1) has a boundary equilibrium point $\widetilde{M_1}$ and may have a positive equilibrium point $M^{*}$ if $R_4>1$ and $R_2>1$.

Only the axis equilibrium point $M_0(A/d,0,0,B/\rho ,0,0)$ is the disease-free equilibrium point. All the boundary equilibria $M_1$, $\widetilde{M_1}$ and $M_2$ and the positive equilibria $M_3$ and $M^{*}$ are the disease-existence equilibria. Epidemiologically, the existence of a disease-free equilibrium point indicates that it is possible to eliminate the disease, the existence of a disease-existence equilibrium point means the disease may exist. Therefore, all $R_1,R_2,R_3,\widetilde{R_3}$ and $R_4$ are thresholds of disease existence and the disease will exist with great probability. Results also indicate that we will do our best to put an end to the infected poultry in the recruitment and control the transmission in poultry so as to keep ${\beta }_i,i=1,2$ lower.

Dynamically, $M_0$ will undergo a series of transcritical bifurcations. $M_0$ will bifurcate an equilibrium point $M_1$ as $R_1$ increases from $R_1=1$ to $R_1>1$ and will bifurcate an equilibrium point $M_2$ as $R_2$ increases from $R_2=1$ to $R_2>1$. $M_0$ also may bifurcate an equilibrium point $M_3$ as $R_{1,2}$ increase from $R_{1,2}=1$ to $R_{1,2}>1$ when $R_3>1$ and $\widetilde{R_3}>1$ are satisfied. The bifurcation portrait of equilibria in the case of $a=0$ is depicted in Figure 4.

4. Stability of Equilibria

In terms of the dynamics of infectious disease, the stability of a disease-existence equilibrium point means the pandemic will occur, and the stability of a disease-free equilibrium point means the extinction of the disease. Now, we will discuss the stability of equilibria. Firstly, we investigate System (1). The Jacobian matrix of System (1) is given as

$$J=\left(\begin{array}{cccccc}-\frac{{\beta }_1{E}_a}{1+{\alpha }_1{E}_a}-\frac{{\beta }_2{I}_a}{1+{\alpha }_2{I}_a}-d& -\frac{{\beta }_1{S}_a}{{(1+{\alpha }_1{E}_a)}^2}& -\frac{{\beta }_2{S}_a}{{(1+{\alpha }_2{I}_a)}^2}& 0& 0& 0\\ \frac{{\beta }_1{E}_a}{1+{\alpha }_1{E}_a}& \frac{{\beta }_1{S}_a}{{(1+{\alpha }_1{E}_a)}^2}-d& 0& 0& 0& 0\\ \frac{{\beta }_2{I}_a}{1+{\alpha }_2{I}_a}& 0& \frac{{\beta }_2{S}_a}{{(1+{\alpha }_2{I}_a)}^2}-d-\xi & 0& 0& 0\\ 0& \frac{{\eta }_1{\nu }_1{S}_h{E}_a^2-{\eta }_1{S}_h}{{(1+{\nu }_1{E}_a^2)}^2}& \frac{{\eta }_2{\nu }_2{S}_h{I}_a^2-{\eta }_2{S}_h}{{(1+{\nu }_2{I}_a^2)}^2}& -\frac{{\eta }_1{E}_a}{1+{\nu }_1{E}_a^2}-\frac{{\eta }_2{I}_a}{1+{\nu }_2{I}_a^2}-\rho & 0& 0\\ 0& \frac{{\eta }_1{S}_h-{\eta }_1{\nu }_1{S}_h{E}_a^2}{{(1+{\nu }_1{E}_a^2)}^2}& \frac{{\eta }_2{S}_h-{\eta }_2{\nu }_2{S}_h{I}_a^2}{{(1+{\nu }_2{I}_a^2)}^2}& \frac{{\eta }_1{E}_a}{1+{\nu }_1{E}_a^2}+\frac{{\eta }_2{I}_a}{1+{\nu }_2{I}_a^2}& -\omega & 0\\ 0& 0& 0& 0& \gamma & -\rho \end{array}\right).$$

Denote

$$J=\left(\begin{array}{cc}H& O\\ L& Q\end{array}\right),$$

where

$$H=\left(\begin{array}{ccc}-\frac{{\beta }_1{E}_a}{1+{\alpha }_1{E}_a}-\frac{{\beta }_2{I}_a}{1+{\alpha }_2{I}_a}-d& -\frac{{\beta }_1{S}_a}{{(1+{\alpha }_1{E}_a)}^2}& -\frac{{\beta }_2{S}_a}{{(1+{\alpha }_2{I}_a)}^2}\\ \frac{{\beta }_1{E}_a}{1+{\alpha }_1{E}_a}& \frac{{\beta }_1{S}_a}{{(1+{\alpha }_1{E}_a)}^2}-d& 0\\ \frac{{\beta }_2{I}_a}{1+{\alpha }_2{I}_a}& 0& \frac{{\beta }_2{S}_a}{{(1+{\alpha }_2{I}_a)}^2}-d-\xi \end{array}\right),$$

$$Q=\left(\begin{array}{ccc}-\frac{{\eta }_1{E}_a}{1+{\nu }_1{E}_a^2}-\frac{{\eta }_2{I}_a}{1+{\nu }_2{I}_a^2}-\rho & 0& 0\\ \frac{{\eta }_1{E}_a}{1+{\nu }_1{E}_a^2}+\frac{{\eta }_2{I}_a}{1+{\nu }_2{I}_a^2}& -\omega & 0\\ 0& \gamma & -\rho \end{array}\right),\omega =\gamma +\delta +\rho .$$

Therefore, J evaluated at every equilibrium point is stable if and only if so are H and Q as J is a block triangular matrix. Obviously, all the eigenvalues of the submatrix Q at every equilibrium point have the negative real parts, then the local stability of every equilibrium point depends on the evaluation of submatrix H. The Jacobian submatrix H corresponds to the poultry subsystem (2) which is independent of System (1). H evaluated at every equilibrium point $M_0,M_1,\widetilde{M_1},M_2,M_3,M^{*}$ is as follows, respectively.

$${J;H|}_{M_0}=\left(\begin{array}{ccc}-d& -\frac{A{\beta }_1}{d}& -\frac{A{\beta }_2}{d}\\ 0& \frac{A{\beta }_1}{d}-d& 0\\ 0& 0& \frac{A{\beta }_2}{d}-d-\xi \end{array}\right),$$

$${J;H|}_{M_1,\widetilde{M_1}}={\left(\begin{array}{ccc}-\frac{{\beta }_1{E}_a}{1+{\alpha }_1{E}_a}-d& -\frac{{\beta }_1{S}_a}{{(1+{\alpha }_1{E}_a)}^2}& -{\beta }_2{S}_a\\ \frac{{\beta }_1{E}_a}{1+{\alpha }_1{E}_a}& \frac{{\beta }_1{S}_a}{{(1+{\alpha }_1{E}_a)}^2}-d& 0\\ 0& 0& {\beta }_2{S}_a-d-\xi \end{array}\right)}_{M_1,\widetilde{M_1}},$$

$${J;H|}_{M_2}=\left(\begin{array}{ccc}-\frac{{\beta }_2{I}_a^{\left(2\right)}}{1+{\alpha }_2{I}_a^{\left(2\right)}}-d& -{\beta }_1{S}_a^{\left(2\right)}& -\frac{{\beta }_2{S}_a^{\left(2\right)}}{{(1+{\alpha }_2{I}_a^{\left(2\right)})}^2}\\ 0& {\beta }_1{S}_a^{\left(2\right)}-d& 0\\ \frac{{\beta }_2{I}_a^{\left(2\right)}}{1+{\alpha }_2{I}_a^{\left(2\right)}}& 0& \frac{{\beta }_2{S}_a^{\left(2\right)}}{{(1+{\alpha }_2{I}_a^{\left(2\right)})}^2}-d-\xi \end{array}\right),$$

$${J;H|}_{M_3,M^{*}}={\left(\begin{array}{ccc}-\frac{{\beta }_1{E}_a}{1+{\alpha }_1{E}_a}-\frac{{\beta }_2{I}_a}{1+{\alpha }_2{I}_a}-d& -\frac{{\beta }_1{S}_a}{{(1+{\alpha }_1{E}_a)}^2}& -\frac{{\beta }_2{S}_a}{{(1+{\alpha }_2{I}_a)}^2}\\ \frac{{\beta }_1{E}_a}{1+{\alpha }_1{E}_a}& \frac{{\beta }_1{S}_a}{{(1+{\alpha }_1{E}_a)}^2}-d& 0\\ \frac{{\beta }_2{I}_a}{1+{\alpha }_2{I}_a}& 0& \frac{{\beta }_2{S}_a}{{(1+{\alpha }_2{I}_a)}^2}-d-\xi \end{array}\right)}_{M_3,M^{*}}.$$

The corresponding eigenvalues of $M_0$ at submatrix H are $-d,A{\beta }_1/d-d$ and $A{\beta }_2/d-d-\xi$. Then the axis equilibrium point $M_0$ is locally asymptotically stable (LAS) if thresholds $R_1<1$ and $R_2<1$; $M_0$ is an unstable saddle if the threshold $R_1>1$ or $R_2>1$; $M_0$ is locally stable if $R_1=1$ and $R_2<1$ or $R_1<1$ and $R_2=1$ because the unique zero eigenvalue is a simple eigenvalue; $M_0$ is a high-order singular point if $R_1=1$ and $R_2=1$.

Because coordinates of $\widetilde{M_1}$ satisfy the condition

$$aA+\frac{{\beta }_1\widetilde{S_a^{\left(1\right)}}\widetilde{E_a^{\left(1\right)}}}{1+{\alpha }_1\widetilde{E_a^{\left(1\right)}}}-d\widetilde{E_a^{\left(1\right)}}=0,$$

then

$$\frac{{\beta }_1\widetilde{S_a^{\left(1\right)}}}{{(1+{\alpha }_1\widetilde{E_a^{\left(1\right)}})}^2}-d<\frac{{\beta }_1\widetilde{S_a^{\left(1\right)}}}{1+{\alpha }_1\widetilde{E_a^{\left(1\right)}}}-d=-\frac{aA}{\widetilde{E_a^{\left(1\right)}}}<0.$$

For $M_1$, there is ${\beta }_1{S}_a^{\left(1\right)}/{(1+{\alpha }_1{E}_a^{\left(1\right)})}^2-d<{\beta }_1{S}_a^{\left(1\right)}/(1+{\alpha }_1{E}_a^{\left(1\right)})-d=0$. So that ${J;H|}_{M_1,\widetilde{M_1}}$ can characterize as

$${J;H|}_{M_1,\widetilde{M_1}}:\left(\begin{array}{ccc}-& -& -\\ +& -& 0\\ 0& 0& {\beta }_2{S}_a-d-\xi \end{array}\right)$$

For the first second-order block matrix $M_{12}=\left(\begin{array}{cc}{a}_{11}& a_{12}\\ a_{21}& a_{22}\end{array}\right):\left(\begin{array}{cc}-& -\\ +& -\end{array}\right)$, it is easy to notice that there are two negative real part eigenvalues. Therefore, the local stability of either $M_1$ or $\widetilde{M_1}$ depends on the sign of ${\beta }_2{S}_a-d-\xi$. Additionally ${\beta }_2{S}_a-d-\xi <A{\beta }_2/d-d-\xi .$ Thus, $M_1$ and $\widetilde{M_1}$ are all locally asymptotically stable (LAS) if the threshold $R_2<1$.

Because $M_2$ satisfies ${\beta }_2{S}_a^{\left(2\right)}{I}_a^{\left(2\right)}/(1+{\alpha }_2{I}_a^{\left(2\right)})=d{I}_a^{\left(2\right)}+\xi I_a^{\left(2\right)}$, so that ${\beta }_2{S}_a^{\left(2\right)}/{(1+{\alpha }_2{I}_a^{\left(2\right)})}^2-d-\xi =(d+\xi )/(1+{\alpha }_2{I}_a^{\left(2\right)})-d-\xi <d+\xi -d-\xi =0$. Similarly, ${J;H|}_{M_2}$ can characterize as

$${J;H|}_{M_2}:\left(\begin{array}{ccc}-& -& -\\ 0& {\beta }_1{S}_a^{\left(2\right)}-d& 0\\ +& 0& -\end{array}\right)$$

Moreover, the second-order block matrix $M_{13}=\left(\begin{array}{cc}{a}_{11}& a_{13}\\ a_{31}& a_{33}\end{array}\right):\left(\begin{array}{cc}-& -\\ +& -\end{array}\right)$, as the same reason that the local stability of $M_2$ depends on the sign of ${\beta }_1{S}_a^{\left(2\right)}-d$. And ${\beta }_1{S}_a^{\left(2\right)}-d<A{\beta }_1/d-d.$ Thus, $M_2$ is locally asymptotically stable (LAS) if the threshold $R_1<1$.

Similar to the above analysis, for ${J;H|}_{M_3,M^{*}}$, there are $a_{22}<0$ and $a_{33}<0$. The Jacobian matrix ${J;H|}_{M_3,M^{*}}$ can be characterized as

$${J;H|}_{M_3,M^{*}}\triangleq \left(\begin{array}{ccc}{a}_{11}& a_{12}& a_{13}\\ a_{21}& a_{22}& 0\\ a_{31}& 0& a_{33}\end{array}\right):\left(\begin{array}{ccc}-& -& -\\ +& -& 0\\ +& 0& -\end{array}\right)$$

The corresponding characteristic polynomial is

$$\begin{array}{ccc}\left|\begin{array}{ccc}\lambda -a_{11}& -a_{12}& -a_{13}\\ -a_{21}& \lambda -a_{22}& 0\\ -a_{13}& 0& \lambda -a_{33}\end{array}\right|\hfill & =\hfill & {\lambda }^3-(a_{11}+a_{22}+a_{33}){\lambda }^2+(a_{11}{a}_{22}+a_{11}{a}_{33}+a_{22}{a}_{33}-a_{13}{a}_{31}-a_{12}{a}_{21})\lambda \hfill \\ & +\hfill & a_{13}{a}_{22}{a}_{31}+a_{12}{a}_{21}{a}_{33}-a_{11}{a}_{22}{a}_{33},\hfill \end{array}$$

where $a_{11}<0,a_{22}<0,a_{33}<0,a_{12}<0,a_{13}<0,a_{21}>0,a_{31}>0$. Then

$$-(a_{11}+a_{22}+a_{33})>0,a_{11}{a}_{22}+a_{11}{a}_{33}+a_{22}{a}_{33}-a_{13}{a}_{31}-a_{12}{a}_{21}>0,a_{13}{a}_{22}{a}_{31}+a_{12}{a}_{21}{a}_{33}-a_{11}{a}_{22}{a}_{33}>0.$$

And the Routh–Hurwitz table is as follows.

${\lambda }^3$	1	$a_{11}{a}_{22}+a_{11}{a}_{33}+a_{22}{a}_{33}-a_{13}{a}_{31}-a_{12}{a}_{21}$
${\lambda }^2$	$-(a_{11}+a_{22}+a_{33})$	$a_{13}{a}_{22}{a}_{31}+a_{12}{a}_{21}{a}_{33}-a_{11}{a}_{22}{a}_{33}$
${\lambda }^1$	$\Delta /(a_{11}+a_{22}+a_{33})$	0
${\lambda }^0$	$a_{13}{a}_{22}{a}_{31}+a_{12}{a}_{21}{a}_{33}-a_{11}{a}_{22}{a}_{33}$	0

where

$$\begin{array}{ccc}\Delta \hfill & =\hfill & \left|\begin{array}{cc}1& a_{11}{a}_{22}+a_{11}{a}_{33}+a_{22}{a}_{33}-a_{13}{a}_{31}-a_{12}{a}_{21}\\ -(a_{11}+a_{22}+a_{33})& a_{13}{a}_{22}{a}_{31}+a_{12}{a}_{21}{a}_{33}-a_{11}{a}_{22}{a}_{33}\end{array}\right|\hfill \\ & =\hfill & \left|\begin{array}{cc}1& a_{11}{a}_{22}+a_{11}{a}_{33}+a_{22}{a}_{33}-a_{13}{a}_{31}-a_{12}{a}_{21}\\ -a_{11}& a_{22}(a_{11}{a}_{22}+a_{11}{a}_{33}+a_{22}{a}_{33}-a_{12}{a}_{21})+a_{33}(a_{11}{a}_{33}+a_{22}{a}_{33}-a_{13}{a}_{31})\end{array}\right|\hfill \\ & :\hfill & \left|\begin{array}{cc}1& +\\ +& -\end{array}\right|<0,\hfill \end{array}$$

Then

$$\Delta <0,\frac{\Delta }{a_{11}+a_{22}+a_{33}}>0.$$

Therefore, items in the first column of the Routh–Hurwitz table are all positive. Thus by the Routh–Hurwitz criterion, all eigenvalues of the submatrix ${J;H|}_{M_3,M^{*}}$ have negative real parts. That is, positive equilibra $M_3$ and $M^{*}$ are locally asymptotically stable if they exist.

Summarizing up the above analysis, we have following conclusions. Results are also represented in

Table 2. Global stabilities of equilibria.

Case	Equilibria	Stabilities and Conditions	Remarks
$a=0$	$M_0(\frac{A}{d},0,0)$	$R_1\le 1,R_2<1$ or $R_2\le 1,R_1<1$
	$M_1(S_a^{\left(1\right)},E_a^{\left(1\right)},0)$	$R_1>1,R_2\le 1$	Exists if $R_1>1$.
	$M_2(S_a^{\left(2\right)},0,I_a^{\left(2\right)})$	$R_2>1,R_1\le 1$	Exists if $R_2>1$.
	$M_3(S_a^{\left(3\right)},E_a^{\left(3\right)},I_a^{\left(3\right)})$	For larger ${\beta }_1$ and ${\beta }_2$	$\frac{2{\beta }_1{\alpha }_1{k}_e{k}_s}{{(1+{\alpha }_1{K}_e)}^2}+\frac{{\beta }_2{\alpha }_2{k}_i{k}_s}{{(1+{\alpha }_2{K}_i)}^2}>2d+\xi$
		and $R_1>1,R_2>1,R_3>1,\widetilde{R_3}>1$	Exists if $R_1>1,R_2>1,R_3>1$ and $\widetilde{R_3}>1$
$a\ne 0$	$\widetilde{M_1}(\widetilde{S_a^{\left(1\right)}},\widetilde{E_a^{\left(1\right)}},0)$	$R_2\le 1$
	$M^{*}(S_a^{*},E_a^{*},I_a^{*})$	For larger ${\beta }_1$ and ${\beta }_2$	$\frac{2{\beta }_1{\alpha }_1{k}_e{k}_s}{{(1+{\alpha }_1{K}_e)}^2}+\frac{{\beta }_2{\alpha }_2{k}_i{k}_s}{{(1+{\alpha }_2{K}_i)}^2}+\frac{2aA}{K_e}>2d+\xi$
		and $R_2>1,R_4>1$	Exists if $R_2>1,R_4>1$

.

Theorem 3.

(i) 

The axis equilibrium point $M_0$ is locally asymptotically stable if $R_1<1$ and $R_2<1$;

(ii) 

$M_1$ is locally asymptotically stable if $R_1>1$ and $R_2<1$;

(iii) 

$\widetilde{M_1}$ is locally asymptotically stable if $R_2<1$;

(iv) 

$M_2$ is locally asymptotically stable if $R_2>1$ and $R_1<1$;

(v) 

$M_3$ is locally asymptotically stable if $R_1>1,R_2>1,R_3>1$ and $\widetilde{R_3}>1$;

(vi) 

$M^{*}$ is locally asymptotically stable if $R_2>1$ and $R_4>1$.

The locally asymptotical stability of an equilibrium point means trajectories initiating in its neighborhood will converge to it eventually, which reflects the local dynamics of the system. Next, we will investigate the global stabilities of equilibria which will represent the global dynamics of the system. We investigate the avian subsystem (2) firstly. Choosing the Liapunov function $V=I_a$, then

$$\dot{V}{|}_{\left(2\right)}=(\frac{{\beta }_2{S}_a}{1+{\alpha }_2{I}_a}-d-\xi )I_a\le ({\beta }_2{S}_a-d-\xi )I_a\le (\frac{A{\beta }_2}{d}-d-\xi )I_a.$$

Therefore $\dot{V}{|}_{\left(2\right)}\le 0$ as $R_2\le 1$. Let $D=\{(S_a,E_a,I_a)\in {\mathbb{R}}_{+}^3:\dot{V}=0\}=\{I_a=0\}$. By the LaSalle Invariance Principle, all solutions of the avian subsystem (2) will approach the $S_a-E_a$ plane.

Now, considering the limit system of the avian subsystem (2) about $I_a$, we obtain

$$\left\{\begin{array}{ccc}{\displaystyle \frac{\mathrm{d}{S}_a}{\mathrm{d}t}}\hfill & =\hfill & (1-a)A-{\beta }_1{S}_a{E}_a/(1+{\alpha }_1{E}_a)-d{S}_a\triangleq F(S_a,E_a),\hfill \\ {\displaystyle \frac{\mathrm{d}{E}_a}{\mathrm{d}t}}\hfill & =\hfill & aA+{\beta }_1{S}_a{E}_a/(1+{\alpha }_1{E}_a)-d{E}_a\triangleq G(S_a,E_a).\hfill \end{array}\right.$$

By theories of limit system, the stabilities of equilibria $M_0,M_1$ and $\widetilde{M_1}$ convert to that of the corresponding limit system’s equilibria (still marked as $M_0,M_1$ and $\widetilde{M_1}$, respectively.). It is easy to obtain the corresponding equilibria $M_0=(A/d,0)$, $M_1(S_a^{\left(1\right)},E_a^{\left(1\right)})$(as $R_1>1$ ) and $\widetilde{M_1}(\widetilde{S_a^{\left(1\right)}},\widetilde{E_a^{\left(1\right)}})$ of the limit system. Furthermore, $M_0$ is locally asymptotically stable as $R_1<1$ and $R_2<1$, $M_1$ and $\widetilde{M_1}$ are locally asymptotically stable as $R_2<1$, respectively.

Next, we select the Dulac function $B(S_a,E_a)=E_a^{-1}$, then

$$\frac{\partial \left(BF\right)}{\partial S_a}+\frac{\partial \left(BG\right)}{\partial E_a}=-\frac{{\beta }_1}{1+{\alpha }_1{E}_a}-\frac{d}{E_a}-\frac{aA}{{\left(E_a\right)}^2}-\frac{{\alpha }_1{\beta }_1{S}_a}{{(1+{\alpha }_1{E}_a)}^2}<0.$$

The Dulac criterion implies that the existence of periodic trajectories are excluded in the region $\tilde{D}=\{(S_a,E_a)\in {\mathbb{R}}_{+}^2:0\le S_a+E_a\le A/d+ϵ\}$. Therefore, by Poincaré-Bendixson theorem, equilibria $M_0,M_1$ and $\widetilde{M_1}$ are globally asymptotically stable in $\tilde{D}$, respectively. That is, $M_0(A/d,0,0)$ is globally asymptotically stable as $R_2\le 1$ and $R_1<1$; $M_1(S_a^{\left(1\right)},E_a^{\left(1\right)},0)$ is globally asymptotically stable as $R_1>1$ and $R_2\le 1$; $\widetilde{M_1}(\widetilde{S_a^{\left(1\right)}},\widetilde{E_a^{\left(1\right)}},0)$ is globally asymptotically stable as $R_2\le 1$. Similarly, $M_0(A/d,0,0)$ is also globally asymptotically stable as $R_1\le 1$ and $R_2<1$.

Using the same discussion for the limiting human-subsystem, we can obtain the globally asymptotical stability of $M_2$ as $R_1\le 1$ and $R_2>1$. Details can refer to the Appendix A.2. Therefore, by the LaSalle Invariance Principle and the theory of limit system, we can conclude that

Theorem 4.

(i) 

$M_0$ is globally asymptotically stable as $R_2\le 1$ and $R_1<1$ or $R_1\le 1$ and $R_2<1$;

(ii) 

$M_1$ is globally asymptotically stable as $R_1>1$ and $R_2\le 1$;

(iii) 

$\widetilde{M_1}$ is globally asymptotically stable as $R_2\le 1$;

(iv) 

$M_2$ is globally asymptotically stable as $R_1\le 1$ and $R_2>1$.

The globally asymptotical stability of a disease-free equilibrium point means the disease will be eradicated finally. This theorem tells us that for the lower ${\beta }_1$ and ${\beta }_2$ it is possible to eliminate the epidemic. Therefore, we should put an end to the infected individuals in the recruitment firstly, then try to control the spread of disease in poultry. Of course, properly reducing the circulation will help control the epidemic.

In the following, we will investigate global stabilities of positive equilibria $M_3$ and $M^{*}$. Moreover, it suffices to discuss System (2). $M_3$ is the unique equilibrium point in ${\Omega }_a$ as $a=0$ and there is an unique equilibrium point $M^{*}$ in ${\Omega }_a$ as $a\ne 0$. Based on the above analysis, all solutions of System (2) are ultimately uniform bounded and System (2) is dissipative. Therefore, System (2) is uniformly persistent and domain ${\Omega }_a$ is a bounded cone. That is, there is a compact subset $K\subset {\Omega }_a$.

Let $x\mapsto f\left(x\right)\in {\mathbf{R}}^n$ be a $C^1$ function for x in an open set ${\Omega }_a\subset {\mathbf{R}}^n$. Consider the differential equation

$$x^{\prime }=f\left(x\right).$$

(6)

Denote $x(t,x_0)$ as the solution of (6) with respect to $x(0,x_0)=x_0$. We make the following two basic assumptions of Equation (6):

Hypothesis 6 (H6).

There exists a compact absorbing set $K\subset {\Omega }_a$.

Hypothesis 7 (H7).

There exists an unique equilibrium point $\overline{x}$ in ${\Omega }_a$.

Therefore, the System (2) satisfies conditions Hypotheses 6 and 7. In order to investigate the global stability of the positive equilibria, we introduce the Bendixson criterion [23]:

Let $x\mapsto P\left(x\right)$ be an $\left(\begin{array}{c}n\\ 2\end{array}\right)$×$\left(\begin{array}{c}n\\ 2\end{array}\right)$ matrix-valued function that is $C^1$ for $x\in \Omega$.

Assume that $P^{-1}\left(x\right)$ exists and is continuous for $x\in K$. Denote

$$B=P_f{P}^{-1}+P\frac{\partial f^{\left[2\right]}}{\partial x}{P}^{-1},$$

where the matrix $P_f$ is obtained by replacing each entry $p_{ij}$ of P by its derivative in the direction of f, $\frac{\partial P_{ij}^{*}}{\partial x}f=\frac{d{P}_{ij}}{dt}{|}_{\left(6\right)}$. $\frac{\partial f^{\left[2\right]}}{\partial x}$ is the second additive compound matrix [24,25] of the Jacobian matrix $\frac{\partial f}{\partial x}$ of f, $\mu \left(B\right)$ is the Lozinskiľ measure of B with respect to a vector norm $\left|\right|·\left|\right|$ in ${\mathbf{R}}^N$, $N=\left(\begin{array}{c}n\\ 2\end{array}\right)$, defined by

$$\mu \left(B\right)=\underset{h\to 0^{+}}{lim}\frac{\left|\right|I+hB\left|\right|-1}{h}.$$

A quantity q is defined as

$$q=\underset{t\to +\infty }{lim\; sup}\underset{x_0\in K}{sup}\frac{1}{t}{\int }_0^t\mu \left(B\left(x(s,x_0)\right)\right)ds.$$

Then,

Lemma 1

([23]). if ${\Omega }_a$ is simply connected and the assumptions $\left(H_1\right),\left(H_2\right)$ hold, the unique equilibrium point $\overline{x}$ of (6) is globally stable in ${\Omega }_a$ if $q<0$.

Here, we denote the right hand of System (2) by $f\left(x\right)$ with $x=(S_a,E_a,I_a)$, the Jacobian matrix $J=\frac{\partial f}{\partial x}{=J;H|}_{M_3,M^{*}}$ of System (2) has been listed above, its second additive compound matrix $J^{\left[2\right]}$ is

$$J^{\left[2\right]}=\left(\begin{array}{ccc}\frac{{\beta }_1{S}_a}{{(1+{\alpha }_1{E}_a)}^2}-\frac{{\beta }_1{E}_a}{1+{\alpha }_1{E}_a}-\frac{{\beta }_2{I}_a}{1+{\alpha }_2{I}_a}-2d& 0& \frac{{\beta }_2{S}_a}{{(1+{\alpha }_2{I}_a)}^2}\\ 0& \frac{{\beta }_2{S}_a}{{(1+{\alpha }_2{I}_a)}^2}-\frac{{\beta }_1{E}_a}{1+{\alpha }_1{E}_a}-\frac{{\beta }_2{I}_a}{1+{\alpha }_2{I}_a}-2d-\xi & -\frac{{\beta }_1{S}_a}{{(1+{\alpha }_1{E}_a)}^2}\\ -\frac{{\beta }_2{I}_a}{1+{\alpha }_2{I}_a}& \frac{{\beta }_1{E}_a}{1+{\alpha }_1{E}_a}& \frac{{\beta }_1{S}_a}{{(1+{\alpha }_1{E}_a)}^2}+\frac{{\beta }_2{S}_a}{{(1+{\alpha }_2{I}_a)}^2}-2d-\xi \end{array}\right).$$

Define the function $P\left(x\right)=P(S_a,E_a,I_a)=diag(1,S_a/I_a,S_a/I_a)$, then $P_f{P}^{-1}=diag(0,S_a^{\prime }/S_a-I_a^{\prime }/I_a,S_a^{\prime }/S_a-I_a^{\prime }/I_a).$ The matrix

$$B=P_f{P}^{-1}+P\frac{\partial f^{\left[2\right]}}{\partial x}{P}^{-1}=\left(\begin{array}{cc}{B}_{11}& B_{12}\\ B_{21}& B_{22}\end{array}\right),$$

where

$$B_{11}=\frac{{\beta }_1{S}_a}{{(1+{\alpha }_1{E}_a)}^2}-\frac{{\beta }_1{E}_a}{1+{\alpha }_1{E}_a}-\frac{{\beta }_2{I}_a}{1+{\alpha }_2{I}_a}-2d,B_{12}=(0,\frac{{\beta }_2{I}_a}{{(1+{\alpha }_2{I}_a)}^2}),B_{21}={(0,-\frac{{\beta }_2{S}_a}{1+{\alpha }_2{I}_a})}^T,$$

$$B_{22}=\left(\begin{array}{cc}{S}_a^{\prime }/S_a-I_a^{\prime }/I_a+\frac{{\beta }_2{S}_a}{{(1+{\alpha }_2{I}_a)}^2}-\frac{{\beta }_1{E}_a}{1+{\alpha }_1{E}_a}-\frac{{\beta }_2{I}_a}{1+{\alpha }_2{I}_a}-2d-\xi & -\frac{{\beta }_1{S}_a}{{(1+{\alpha }_1{E}_a)}^2}\\ \frac{{\beta }_1{E}_a}{1+{\alpha }_1{E}_a}& S_a^{\prime }/S_a-I_a^{\prime }/I_a+\frac{{\beta }_1{S}_a}{{(1+{\alpha }_1{E}_a)}^2}+\frac{{\beta }_2{S}_a}{{(1+{\alpha }_2{I}_a)}^2}-2d-\xi \end{array}\right).$$

Let $(u,v,w)$ be the vector in ${\mathbf{R}}^3\cong {\mathbf{R}}^{\left(\begin{array}{c}3\\ 2\end{array}\right)}$, we take a norm in ${\mathbf{R}}^3$ as $\left|\right|(u,v,w)\left|\right|=max\left\{\right|u|,|v|+|w\left|\right\}$ and $\mu$ denotes the Lozinskiľ measure with respect to this norm. By the method in Fiedler [24], we have $\mu \left(B\right)\le sup\{g_1,g_2\}$, where

$$g_1={\mu }_1\left(B_{11}\right)+|B_{12}|,g_2=\left|B_{21}\right|+{\mu }_1\left(B_{22}\right).$$

$|B_{12}|$ and $|B_{21}|$ are the matrix norms with respect to the $l_1$ vector norm, ${\mu }_1$ is the Lozinskiľ measure with respect to the $l_1$ norm.

There are

$$|B_{12}|=\frac{{\beta }_2{I}_a}{{(1+{\alpha }_2{I}_a)}^2},\left|B_{21}\right|=\frac{{\beta }_2{S}_a}{1+{\alpha }_2{I}_a}$$

and

$$g_1=\frac{{\beta }_1{S}_a}{{(1+{\alpha }_1{E}_a)}^2}-\frac{{\beta }_1{E}_a}{1+{\alpha }_1{E}_a}-\frac{{\alpha }_2{\beta }_2{I}_a^2}{{(1+{\alpha }_2{I}_a)}^2}-2d,$$

$$\begin{array}{ccc}{g}_2\hfill & =\hfill & {\beta }_2{S}_a/(1+{\alpha }_2{I}_a)+S_a^{\prime }/S_a-I_a^{\prime }/I_a\hfill \\ & & +max\{{\beta }_2{S}_a/{(1+{\alpha }_2{I}_a)}^2-{\beta }_2{I}_a/(1+{\alpha }_2{I}_a)-2d-\xi ,2{\beta }_1{S}_a/{(1+{\alpha }_1{E}_a)}^2+{\beta }_2{S}_a/{(1+{\alpha }_2{I}_a)}^2-2d-\xi \}\hfill \\ & =\hfill & {\beta }_2{S}_a/(1+{\alpha }_2{I}_a)+S_a^{\prime }/S_a-I_a^{\prime }/I_a+2{\beta }_1{S}_a/{(1+{\alpha }_1{E}_a)}^2+{\beta }_2{S}_a/{(1+{\alpha }_2{I}_a)}^2-2d-\xi .\hfill \end{array}$$

Because

$$\left\{\begin{array}{ccc}{\displaystyle \frac{{\mathrm{S}}_a^{\prime }}{{\mathrm{S}}_a}}\hfill & =\hfill & (1-a)A/S_a-{\beta }_1{E}_a/(1+{\alpha }_1{E}_a)-{\beta }_2{I}_a/(1+{\alpha }_2{I}_a)-d,\hfill \\ {\displaystyle \frac{{\mathrm{E}}_a^{\prime }}{{\mathrm{E}}_a}}\hfill & =\hfill & aA/E_a+{\beta }_1{S}_a/(1+{\alpha }_1{E}_a)-d,\hfill \\ {\displaystyle \frac{{\mathrm{I}}_a^{\prime }}{{\mathrm{I}}_a}}\hfill & =\hfill & {\beta }_2{S}_a/(1+{\alpha }_2{I}_a)-d-\xi ,\hfill \end{array}\right.$$

Then

$$g_1=\frac{E_a^{\prime }}{E_a}-\frac{{\beta }_1{E}_a}{1+{\alpha }_1{E}_a}-\frac{{\beta }_1{\alpha }_1{E}_a{S}_a}{{(1+{\alpha }_1{E}_a)}^2}-\frac{{\alpha }_2{\beta }_2{I}_a^2}{{(1+{\alpha }_2{I}_a)}^2}-d-\frac{aA}{E_a},$$

$$g_2=\frac{S_a^{\prime }}{S_a}+2\frac{E_a^{\prime }}{E_a}+\frac{I_a^{\prime }}{I_a}-\frac{2{\beta }_1{\alpha }_1{E}_a{S}_a}{{(1+{\alpha }_1{E}_a)}^2}-\frac{{\beta }_2{\alpha }_2{I}_a{S}_a}{{(1+{\alpha }_2{I}_a)}^2}-\frac{2aA}{E_a}+2d+\xi .$$

It is obvious that $g_1\le E_a^{\prime }/E_a-d$, and there are constants $k_s,K_s,k_e,K_e,k_i,K_i>0$ and $T>0$ independent of $(S_a\left(0\right),E_a\left(0\right),I_a\left(0\right))\in K$ such that $k_s\le S_a\left(t\right)\le K_s,k_e\le E_a\left(t\right)\le K_e,k_i\le I_a\left(t\right)\le K_i$ for $t>T$ because of the uniform persistence of the system. Then

$$g_2\le \frac{S_a^{\prime }}{S_a}+2\frac{E_a^{\prime }}{E_a}+\frac{I_a^{\prime }}{I_a}-(\frac{2{\beta }_1{\alpha }_1{k}_e{k}_s}{{(1+{\alpha }_1{K}_e)}^2}+\frac{{\beta }_2{\alpha }_2{k}_i{k}_s}{{(1+{\alpha }_2{K}_i)}^2}+\frac{2aA}{K_e}-2d-\xi ).$$

Denote

$$\tilde{d}=\frac{2{\beta }_1{\alpha }_1{k}_e{k}_s}{{(1+{\alpha }_1{K}_e)}^2}+\frac{{\beta }_2{\alpha }_2{k}_i{k}_s}{{(1+{\alpha }_2{K}_i)}^2}+\frac{2aA}{K_e}-2d-\xi .$$

Therefore, along each solution $(S_a\left(t\right),E_a\left(t\right),I_a\left(t\right))$ of System (2) such that $(S_a\left(0\right),E_a\left(0\right),I_a\left(0\right))\in K$ and $t>T$, we have

$$\frac{1}{t}{\int }_0^t{g}_{1}ds\le \frac{1}{t}{\int }_0^T{g}_{1}ds+\frac{1}{t}ln\frac{E_a\left(t\right)}{E_a\left(T\right)}-d\frac{t-T}{t},$$

$$\frac{1}{t}{\int }_0^t{g}_{2}ds\le \frac{1}{t}{\int }_0^T{g}_{2}ds+\frac{1}{t}ln\frac{S_a\left(t\right)}{S_a\left(T\right)}+\frac{2}{t}ln\frac{E_a\left(t\right)}{E_a\left(T\right)}+\frac{1}{t}ln\frac{I_a\left(t\right)}{I_a\left(T\right)}-\tilde{d}\frac{t-T}{t}if\tilde{d}>0.$$

Let $b=max\{d,\tilde{d}\}$, then $q\le -b/2<0$.

Hence, by Lemma 1, we have $(S_a\left(t\right),E_a\left(t\right),I_a\left(t\right))\to (S_a^{*},E_a^{*},I_a^{*})$ or $M_3(S_a^{\left(3\right)},E_a^{\left(3\right)},I_a^{\left(3\right)})(a=0)$ as $t\to +\infty$. Then, we get the conclusion:

Theorem 5.

Either the positive equilibrium point $M_3$ (in the case of $a=0$) or $M^{*}$ (in the case of $a\ne 0$) exists and will be globally asymptotically stable if $2{\beta }_1{\alpha }_1{k}_e{k}_s/{(1+{\alpha }_1{K}_e)}^2+{\beta }_2{\alpha }_2{k}_i{k}_s/{(1+{\alpha }_2{K}_i)}^2+2aA/K_e>2d+\xi$.

The global stability of the positive equilibrium point means the pandemic will occur. Therefore, it is possible for the global stability of the disease-existence equilibrium point $M_3$ or $M^{*}$ for the larger ${\beta }_1$ and ${\beta }_2$. Results also indicate the pandemic risk is increased if there is infected poultry in the recruitment.

Summing up the above analysis, the global stabilities of equilibria are presented in Table 2.

5. Simulations

In simulations, for mathematical convenience, we suppose ${\alpha }_1={\alpha }_2=\alpha$, ${\eta }_1={\eta }_2=\eta$ and ${\nu }_1={\nu }_2=\nu$. There are reports from the Chinese CDC that in April 2017, 8500 poultry in a farm in Hebei were infected with the new H7N9 and 5000 died; in May, 7500 poultry in a farm in Henan were infected and 5770 died; in a farm in Tianjin, 10,000 were infected and 6000 died; in a farm in Shangxi, were infected and 22,000 died; in a farm in Hohhot, Inner Mongolia, 59,556 were infected and 35,526 died; in June, in a farm in Batou, Inner Mongolia, 3850 were infected and 2056 died; in a farm in Heilongjiang, 20,150 were infected and 19,500 died, and so on. These data indicate that poultry can be infected with high pathogenic viruses in large numbers. Therefore, without loss of generality, we suppose ${\beta }_2>{\beta }_1$. For an outbreak region, let us refer to A = 10,000 and $B=500$. For poultry infection with the high pathogenic A(H7N9) virus, the extra death rate caused by the disease might be $\xi =0.6$ as well based on the above data. Poultry on farms usually are kept for 50 days in the farm industry in China, so the output rate is $d=1/50=0.02$. The natural death rate of people can be given as $\rho =1/(75\times 365)$ as far as life expectancy is 75 years. The death rate of confirmed cases is about 40% and the treatment time of confirmed cases would be 15 days to 20 days according to documents from the Chinese CDC. Because the severe critical pneumonia mostly occurs about 3 to 7 days after the onset of the disease for a confirmed case, we may take an average of 5 days as reference, then the death rate due to disease is $\delta =1/5\times 0.4=0.08$ and the recovery rate is $\gamma =1/15\times 0.6=0.04$. Let $p=0.2,{\beta }_1=0.003,\alpha =0.0001,\eta =0.0001,\nu =0.0005$ as well.

Now we carry out simulations with the purpose of investigating system dynamics, explaining control measures and doing some comparative studies for human infections with the low, high and mixed A(H7N9) viruses. Firstly, the dynamic behaviors of System (2) are illustrated in Figure 5. For the above parameters, all five thresholds are greater than 1 and it can be seen from Figure 5a that the positive equilibrium point $M^{*}$ is globally asymptotically stable. If we let $A=1000$, there is $R_3<1$. Therefore, the positive equilibrium point $M^{*}$ does not exist and we can see from Figure 5b that the trajectories tend to the boundary equilibrium point $\widetilde{M_1}$.

Furthermore, the dynamics of System (1) is depicted in Figure 6 which presents curves of $I_h\left(t\right)$ for human infections with the low, high and mixed A(H7N9) viruses. We can see that the mixed A(H7N9) viruses will cause more serious infections than the other two, so the unusually serious epidemic in the fifth wave is explicated. The curve of $I_h\left(t\right)$ corresponding to human infections with the low pathogenic A(H7N9) virus shows a different trend, that the curve will reach a small peak in a short time, then increase further after falling. This hints that people should respond quickly to control the epidemic before reinflating for the low pathogenic A(H7N9) infections. Moreover, it also suggests that human infection with the mixed A(H7N9) viruses is not the simple superposition of human infection with the low and high pathogenic A(H7N9) viruses.

In view of the practical control measures, we can see from Figure 7 and Figure 8 that there are serious infections if $I_a\left(0\right)\ne 0$, which can be detected from the representations of the three up curves and the two down curves which correspond to $I_a\left(0\right)=0$. This suggests that the more important thing is to identify and eliminate the infected poultry with the high pathogenic A(H7N9) virus. It is possible for us to eliminate the infected poultry with symptoms by killing and burying all sick poultry. As for the identification of infected poultry with low pathogenic A(H7N9) viruses on farms and in the recruitment, which one is more important? Figure 7 and Figure 8 tell us that it is different in the short term and in the long term, which can be extracted from the two down curves in Figure 7 and Figure 8. It seems that in the short term, the identification and elimination of the infected poultry on farms is prioritized. However, if we cannot eliminate the infected poultry with the low pathogenic A(H7N9) virus in the recruitment, there is no difference for $E_a\left(0\right)=0$ or not which can be seen from the overlap curves. However, without the identification technology, it is also important to eliminate the infected poultry with the low pathogenic A(H7N9) virus in the recruitment which can be captured from the two up curves in Figure 8. In short, the identification technology is very important for the control of disease. It is more important for the identification of infected poultry with the high pathogenic A(H7N9) virus, than the infected poultry with the low pathogenic A(H7N9) virus on farms and in the recruitment.

From Figure 9 and Figure 10, we can see that the rapid infection of all poultry will occur whether by the low pathogenic virus or the high pathogenic virus. However, the infection caused by the low pathogenic virus of poultry to human beings far exceeds that caused by the high pathogenic virus in terms of time and quantity. There are two waves of human infections in the case of the mixed transmission. This may be the reason that the Chinese government decisively implemented strict vaccine measures against all poultry in the country when the mixed A(H7N9) infections seriously threatened people in 2017.

Figure 11 suggests that controlling contact both in poultry and between human and poultry can effectively control the epidemic in a month and losing controlling contacts in poultry will cause a serious pandemic. It seems it is very sensitive of parameters $\eta$ and $\nu$ embodying human response and the psychosocial effect, compared to parameters $\beta$ and $\alpha$ which correspond to the transmission in poultry. Losing control of the contact between poultry will lead to a surge of human cases, which is depicted in Figure 12. Therefore, controlling the contact between poultry is as significant as controlling the contact for people with poultry. In short, we should respond immediately and take measures to reduce the contagions in poultry, and at the same time publicize the information so as to avoid the dangerous levels of exposure once a human case is confirmed.

6. Conclusions and Discussion

In 2013, the avian influenza A(H7N9) virus, which previously circulated only in poultry, began to infect humans. Hereafter, a new wave of epidemics emerged in every autumn and winter, severely discouraging the poultry industry and threatening public health. Frightened by the surge in number of human cases and the extension of the distribution of infected areas in the wave of 2016–2017, considering the fact that human infections with both the original low pathogenic and the newly emerging high pathogenic A(H7N9) viruses, we established and analyzed a two-strain avian–human infection model in order to investigate the dynamics of the epidemic, to explain and explore disease control measures and to make some comparative studies for human infections with the low or high pathogenic viruses and the mixed infections.

Because the low and high pathogenic A(H7N9) viruses are susceptible to both poultry and humans, the dynamics of the epidemic are more complex and sensitive. Theoretical results suggest there are four equilibria in the case of $a=0$ and two equilibria in the case of $a\ne 0$. In addition to the trivial equilibrium point $M_0$, the others are all the disease-existence equilibria, including the boundary equilibria. It greatly increases the risk of pandemic. We also obtain the bifurcation portrait of the existence and stability of equilibria at $a=0$. From the perspective of the existence conditions, basically equilibria $M_1$ and $M_2$ will exist for the larger ${\beta }_i,i=1,2$. However, for the positive equilibrium point $M_3$, the situation becomes more complicated, which also hints at the complex and sensitive nature of the dynamics. In fact, the existence of $M_3$ is not inevitable, which can be detected from the expressions of its existence thresholds. We also know the existing equilibria $\widetilde{M_1}$ and $M^{*}$ in the case of $a\ne 0$ will converge to the corresponding $M_1$ and $M_3$ at $a=0$, although the required existence conditions are stronger. In any case, theoretically it is still the prior to eliminate the infected poultry in the recruitment so that $a=0$, which may be beneficial to the disease eradication. The trivital equilibrium point $M_0$ is a high-order singular point as $R_1=1$ and $R_2=1$ and will bifurcates equilibria $M_1$, $M_2$ and $M_3$ as $R_i,i=1,2$ increase from $R_i=1$ to $R_i>1,i=1,2$. Of course, $R_i,i=1,2$ are the bifurcation thresholds. $M_0$ undergoes the transcritical bifurcations. Furthermore, if we can put an end to the infected poultry in recruitment, the stabilities of equilibria tell us that we have the possibility to eradicate the epidemic by controlling the transmission in poultry such that $R_i<1,i=1,2$. The related results are presented in Theorems 1–3. The global stabilities of the positive equilibria $M_3$ and $M^{*}$ indicate again the importance of eliminating the infected poultry in recruitment which can be extracted from Theorem 4.

The implemented simulations verify the theoretical results and reveal that the dynamics of human infections with the low pathogenic A(H7N9) virus are different from that of the high pathogenic A(H7N9) virus and the mixed infections. The mixed low and high pathogenic A(H7N9) viruses would cause more serious infections in humans. Moreover, human infections with mixed low and high pathogenic A(H7N9) viruses is not the simple superposition of human infections with low and high pathogenic A(H7N9) viruses. For disease control, simulations suggest that the more important thing is to develop the identification technology to eliminate the infected poultry. Again, controlling the contact between human and poultry can effectively control the epidemic, and controlling the contagions in poultry can avoid a terrible rise in infections in humans. Evacuation stocking, isolated chicken houses and public notification mechanisms should be strengthened. This can give reference to the case of the simultaneous spreading of two strains between humans and poultry.

Because of the diversity of disease transmission, the widely used SEIR model is extended to include more categories of population, such as the asymptomatic, the hospitalized patients and the isolated individuals, etc. Additionally, more and more factors are also considered in modeling, such as medical resources, traffic power, contact behaviors, social and economic factors, and so on. Among these newly emerging models, it is worth mentioning those researches based on complex networks and networked populations. For example, Wang et al. in [26] gave a detailed and valuable introduction on epidemics in networked populations, which contributed the understanding of the complex networked disease transmission. Li et al. in [27] proposed a modified signed-susceptible-infectious-susceptible epidemiological model with positive and negative transmission rates and structural balance. Research on epidemic spreading via complex systems would be a successful way for us to face the challenge of disease transmission.

Obviously, this study is mainly focused on theoretical dynamics. In view of the realities of a complex epidemic, there are some factors which should be taken into account. For example, factor of the transmissions caused by the latent poultry infected with the high pathogenic virus should be considered, and the transmission via environment and objects should also be incorporated into models. From the point of view of disease control, the isolation of poultry production and trading/marketplaces are also worth considering. Therefore, further research will be carried out. There are also issues of the cotransmission and alternate prevalence of multivirulent strains in human seasonal influenza epidemic and the cross-transmission of sudden infectious diseases in the influenza season. In view of the interpersonal transmission of diseases, a networked approach will help to understand the spread of the disease. We will consider these in our subsequent research.