# Dynamical Analysis of theSEIB Model for Brucellosis Transmission to the Dairy Cows with Immunological Threshold

Authors: Wei Zhang; Juan Zhang; Yong-Ping Wu; Li Li
Journal: Complexity
 (2019)

---

## Abstract

As we all know, bacteria is different from virus which with certain types can be killed by the immune cells in the body. The brucellosis, a bacterial disease, can invade the body by indirect transmission from environment, which has not been researched by combining with immune cells. Considering the effects of immune cells, we put a minimum infection dose of brucellosis invading into the dairy cows as an immunological threshold and get a switch model. In this paper, we accomplish a thorough dynamics analysis of aSEIB switch model. On the one hand, we can get a disease-free and bacteria-free steady state and up to three endemic steady states which may be thoroughly analyzed in different cases of a minimum infection dose in a switch model. On the other hand, we calculate the basic reproduction number R0 and know that the disease-free and bacteria-free steady state is a global stability when R0<1, and the one of the endemic steady state is a conditionally global stability when R0>1. We find that different amounts of R0 may lead to different steady states of brucellosis, and considering the effects of immunology is more serious in mathematics and biology.


---

## Body

## 1. Introduction

Brucellosis is the zoonotic sex contagion which is named as Mediterranean relax heat, wave heat, or wave form heat [1–3]. It is characterized by extensive host, strong infectivity, and difficulty in radical treatment after infection [4, 5]. Besides, it has a serious harm to economy, society, and public health. Therefore, brucellosis is listed as one of the communicable diseases that must be notified in the World Organization for Animal Health and is classified as second kind of animal diseases in China. Meanwhile, it is as the first zoonotic disease to be controlled in the National Medium and Long Term Plan for Animal Disease Prevention and Control 2011-2020 formulated by the ministry of agriculture. Brucellosis is acute or chronic infectious disease caused by Brucella which is a group of small bulbous gram-negative bacteria [6, 7]. In 1985, the World Health Organization (WHO) divided Brucella into six species and nineteen biological types [7]. They are more popular with Br. melitensis, Br. brovis, and Br. suis in China, among which the most popular Brucella is Br. melitensis and the second is Br. brovis. The most noticeable symptom of brucellosis is miscarriage in cows and orchitis in bulls. In addition, it can cause arthritis in the knee and wrist frequently [2, 3]. The route of spread is through direct contact with diseased animals and indirect infection with bacteria in the environment [8, 9]. And it can be divided into three ways. The first transmitted way is contacting with skin, such as direct contact with droppings, vaginal secretions, and vaginal delivery content of sick animals. It also can be indirectly exposed to the environment and objects contaminated by sick animals. The second infected way is via the digestive tract, like eating food, water, or milk contaminated with pathogens. And the last way of transmission is through the respiratory tract.When bacteria enter the body, mainly in the liver, spleen, bone marrow, lymphatic tissue, and other cells, those will grow, multiply, and produce endotoxin which can cause damaged tissues. In the meantime, it often releases bacteria and toxins into the blood causing systemic bacteremia and allergic reactions. But the Brucella does not take effect immediately when it enters the body. The bacteria with a small number will be killed by immune cells and only amounts of bacteria over a certain threshold called a minimum infection dose (MID) in environment can do harm to body [10, 11]. And different bacteria have variant thresholds which should be tested by experiment. Some bacteria threshold can be seen in [12]. Also, there were many authors researching the brucellosis, who had been unconscious of the impact on immunization [13–19]. By referring to [20, 21], we will put the pathogeny of brucellosis and MID into this paper. The authors [12, 20, 21] incorporated a MID, c, into the incidence term α(B), in a cholera transmission model, which was a piecewise continuous function which was zero under the MID and was a Holling-II response curve over that. (1)αB=0ifB<c,aB-cB-c+HifB≥c.where the parameter c represents MID, H is half-saturation pathogen density, and the description of parameter a is maximum rate of infection. One of typical mathematical models given in [22] is as follows:(2)S˙=A+qS+E+I-β1SI-β2SV-dS,E˙=β1SI+β2SV-dE-δE,I˙=δE-dI-cI,B˙=mE+I-w+ehB.where all parameters are positive and the description of parameters can be found in Table 1.Table 1
The description of parameters in original model.
Parameter Description Dimension A The annual introduction number of dairy cows year-1 q The annual birth rate of dairy cows year-1 d The annual natural elimination rate year-1 δ Clinical outcome rate year-1 c The annual elimination rate for the positive cows y e a r - 1 m The annual quantity of Brucella y e a r - 1 w The annual natural mortality of Brucella y e a r - 1 β 1 Cow-to-cow transmission rate none β 2 Brucella-to-cow transmission rate none h The sterilizing rate in a disinfection none e The number of disinfections every year onceIn [22], the authors gave a detailed demonstration of the basic reproduction number R0(R0=A(w+eh)β1δ+Am(d+c+δ)β2/(d+c)(d+δ)(w+eh)(d-q)), which can estimate the occurrence of the epidemic, and made brief dynamical analysis on the global stability.Basing on model (2) and considering the immunology of Brucella in the environment, we get a switch system with MID (M). If the number of Brucella in environment (B) is less than or equal to the MID, then we can get a globally asymptotically stable disease-free and bacteria-free steady state while the basic reproduction number is less than or equal to 1, an unstable disease-free and bacteria-free equilibrium and a conditionally globally stable endemic steady state when the basic reproduction number is more than 1. Else if B>M, it is going to be a little bit complicated. Combining with the relationship between parameter M and Mg (which is the formula received from the system), we can attain the three different switch systems. One of them can attain four equilibria and can produce a backward bifurcation [23, 24] and might be two steady states at most under certain parameter conditions.


## 2. Dynamical Analysis of the Model

### 2.1. The Dynamical Model

We classify the dairy cows into three compartments: the susceptible compartmentS, the exposed compartment E, and the infectious compartment I in Figure 1. And the Brucella in environment is denoted by B.Figure 1
A flow diagram of brucellosis transmission.One of the key differences of (2) is the incidence term α(B)S of indirect transmission of Brucella to susceptible cow [12], where α(B) is the pathogen density dependent component. Unlike [12, 20] using Holling-II response curve when the bacteria density is over the MID, a simple linear form is used in this paper. There are chiefly the following reasons. Firstly, the half-saturation pathogen density (H) is difficult to be determined in actual application. Secondly, the Brucella which is discharged to the nature by dairy cows and other livestock is hard to reach saturation. Lastly, that can get more practical conclusions with immunological threshold than [22]. The expression of α(B) in this paper is defined as(3)αB=0if0≤B≤M,aB-MifB>M.If the number of bacteria is more than MID, then the Brucella of environment can bring the indirect influence to the dairy cows. At this time,a is representative of maximum rate of indirect infection, and M is denoted MID of Brucella, B-M is representative of bacteria amount which can enter the body and take effect. If the number is less than or equal to MID, then the bacteria of environment have no effect on the susceptible and dairy cows were infected just by the direct transmission from infected to susceptible cow.So the model is a system of ordinary differential equation as follows:(4)S˙=A+qS+E+I-βSI-dS-αBS,E˙=βSI+αBS-dE-δE,I˙=δE-dI-cI,B˙=mE+I-w+ehB,where the description of other parameters which are not mentioned in Table 1 can be found in Table 2.Table 2
The description of new parameters in model (4).
Parameter Description Dimension M MID c e l l   l i t e r - 1 a Maximum rate of infection none β Indirect transmission rate noneRemark 1.
In this paper, we also use the conditiond>q like [22].


### 2.2. Forward Invariance

The first equation of (4) gives us that S(t)˙≥0, when S=0. In the same way, we can get E(t)˙≥0 and I(t)˙≥0. Hence, we know that S(t)≥0, E(t)≥0, and I(t)≥0 for t>0. As S˙+E˙+I˙≤A-(d-q)(S+E+I), we can get 0≤S+E+I≤A/d-q by using the character of S˙+E˙+I˙. Cause S˙+E˙+I˙=A-(d-q)(S+E+I) is a linear function [12]. In the same way, if B=0, then B(t)˙=0. Thus, the inequality B(t)≥0 for t>0 is true. If B(0)∈[0,m(E+I)/w+eh), then B(t)∈[0,m(E+I)/w+eh) for any t>0. To sum up, we get the following theorem.Theorem 2 (feasible region).
The set(5)Ω=S,E,I,B:0≤S+E+I≤Ad-q,0≤B<mAw+ehd-q,S≥0,E≥0,I≥0defines a forward invariant region of system (2).


### 2.3. Equilibria of System (4)

#### 2.3.1. Equilibria When Bacteria Are Less Than or Equal to MID

When0≤B≤M, then α(B)=0 and the model is(6)S˙=A+qS+E+I-βSI-dS,E˙=βSI-dE-δE,I˙=δE-dI-cI,B˙=mE+I-w+ehB.We solve the equation and can getE0=(S0,0,0,0), a disease-free and bacteria-free equilibrium where S0=A/d-q, and a endemic equilibrium E1=(S1,E1,I1,B1), where(7)S1=d+δd+cβδ,E1=d+cβδA-d-qd+δd+cβδd+δd+c-qd+c+δ,I1=βδA-d-qd+δd+cβd+δd+c-qd+c+δ,B1=md+c+δβδA-d-qd+δd+cw+ehβδd+δd+c-qd+c+δ.Remark 3.
S 1 > 0, E1>0, I1>0 and B1>0.According to the spectral radius theory [25, 26], we can directly obtain the basic reproduction number based on [22] when β2=0. That is (8)R0=ρFV-1=βδAd-qd+cd+δ.Then we can attain the existence of equilibria: if R0≤1, then there is only one equilibrium E0; if R0>1, then there are disease-free and bacteria-free equilibrium E0 and endemic equilibrium E1.


#### 2.3.2. Equilibria When Bacteria Are More Than MID

WhenB>M, the equation is(9)S˙=A+qS+E+I-βSI-dS-aB-MS,E˙=βSI+aB-MS-dE-δE,I˙=δE-dI-cI,B˙=mE+I-w+ehB.Obviously, there is no disease-free and bacteria-free equilibria because of M≥0. Thus, system (9) implies that(10)E=d+cw+ehBmd+c+δ,I=δw+ehBmd+c+δ,S=d+δEβI+αB=qE+I+AβI+αB+d-q.The equation ofS about B can be written as F(B)=0. Define F(B)=B2(w+eh)[d+cd+δ-qd+c+δ][βδ(w+eh)+am(d+c+δ)]+Bm(d+c+δ){-aM(w+eh)[(d+c)(d+δ)-q(d+c+δ)]-(w+eh)[βδA-(d-q)(d+δ)(d+c)]-am(d+c+δ)A}+aMAm2(d+c+δ)2. From the expression, we can know that F(B) is a quadratic function, in which the first and third coefficients are positive. Therefore, the existence of solution for F(B)=0 depends on the second coefficient and Δ. If the second coefficient is nonnegative, no matter Δ is, then F(B)=0 has no positive solution. If the second coefficient is negative, the positive solution depends on the sign of Δ: when Δ<0, there is no solution; when Δ=0, there is one positive solution; and when Δ>0, there are two positive solutions.Now we simplify them and can get the new second coefficientb1:(11)b1≔md+c+δ-βδw+eh-mad+c+δA+b11,where(12)b11≔-aMw+ehd+cd+δ-qd+c+δ+w+ehd-qd+δd+c, and the new Δ about A:(13)ΔA≔A2m2d+c+δ2-βδw+eh-amd+c+δ2+A-2m2d+c+δ2w+ehβδw+eh+amd+c+δaMd+cd+δ-qd+c+δ+d-qd+δd+c+m2d+c+δ2-aMw+ehd+cd+δ-qd+c+δ+w+ehd-qd+δd+c2.We can knowb1 is a liner function about A and Δ(A) is still a quadratic function which the first and third coefficients are positive and the second coefficient is less than zero. Then we calculate the Δ(A) and can know that the sign of Δ(Δ(A)) is always positive. So we can solve Δ(A)=0 and get A1, A2, where(14)A1=w+ehaMd+cd+δ-qd+c+δ+d-qd+cd+δ-2aMd+cd+δ-qd+c+δd-qd+cd+δβδw+eh+mad+c+δ,A2=w+ehaMd+cd+δ-qd+c+δ+d-qd+cd+δ+2aMd+cd+δ-qd+c+δd-qd+cd+δβδw+eh+mad+c+δ.By the analysis with the feature of function, we can get the conclusion as follows.Proposition 4.
The existence of positive solutions forF(B)=0:
(i) Whenb11≥0: if 0≤A<A2, then F(B)=0 does not get any positive solution; if A=A2, then F(B)=0 has one positive solution B1+; if A>A2, then F(B)=0 has two positive solutions B1+, B2+ (B1+<B2+).
(ii) Whenb11<0: if A1<A<A2, then F(B)=0 does not get any positive solution; if A=A1 or A=A2, then F(B)=0 has one positive solution B+ (=B1+=B2+); if 0≤A<A1 or A>A2, then F(B)=0 has two positive solutions B1+, B2+.Note theM~=(d-q)(d+c)(d+δ)/a((d+c)(d+δ)-q(d+c+δ)) (by the equation of b11=0) and we can put Proposition 4 to the chart as shown in Table 3.Table 3 F ( B ) = 0 has one positive solution B+ F ( B ) = 0 has two positive solutions B1+, B2+  M ≥ 0, A=A2 M ≥ 0, A>A2  M > M ~, A=A1 M > M ~, 0≤A<A1It just discusses the existence of positive solutions, but we need consider the value ofM, where only internal equilibria exist while positive solutions are more than M. Hence, we will get conclusions as follows.Theorem 5.
The existence of internal equilibria aboutB>M:(1)
IfM>0:
(i)
whenA>Af, there is an internal equilibrium E2+;(ii)
whenA=Af>Ag, there is an internal equilibrium E2+;(iii)
whenmax⁡Ag,A2<A<Af, there are two internal equilibria E1+, E2+;(iv)
whenA=A2>Ag, there is an internal equilibrium E+;(2)
IfM=0 and A>A2=Ag, then there is an internal equilibrium E2+.
where(15)Af=βδw+ehd+cd+δ-qd+c+δM+md+c+δd-qd+cd+δβδmd+c+δ,Ag=2βδw+eh2+amw+ehd+c+δd+cd+δ-qd+c+δM+w+ehmd+c+δd-qd+cd+δmd+c+δβδw+eh+mad+c+δ.Proof.
We have known the positive value ofB; now we need combine with the inequality B>M. The proof uses character of function F(B).
Firstly, we should know the value ofF(M). When B=M, we can get F(M)=βδ(w+eh)2((d+c)(d+δ)-q(d+c+δ))M2-m(d+c+δ)(w+eh)(βδA-(d-q)(d+c)(d+δ))M, and the relationship between F(M) and 0 is translated into the relationship between A and Af: (1)
WhenM>0: if A>Af, then F(M)<0, if A=Af, then F(M)=0, and if A<Af, then F(M)>0.(2)
WhenM=0, we can get F(M)=0.
Secondly, we need to know position ofM and symmetry axis of F(M). We put them into the formula g(M)≔[2βδ(w+eh)2+am(w+eh)(d+c+δ)]((d+c)(d+δ)-q(d+c+δ))M-m(d+c+δ)[(w+eh)(βδA-(d-q)(d+c)(d+δ))+ma(d+c+δ)A]. In the same way, we translate them into the relationship between A and Ag. It means that when M≥0: if A>Ag, then g(M)<0 (M is on the left of axis of symmetry); if A=Ag, then g(M)=0 (M is the axis of symmetry); and if A<Ag, then g(M)>0 (M is on the right of axis of symmetry).
Thirdly, we calculate the relationship ofAf, Ag, A1, A2. Then we can get the relationships among them as follows: (1)
The relationships ofAf>A1 and Ag>A1 are always true.(2)
When0<M≤Mg, there is Ag≤A2≤Af. When M>Mg, there is A2≤Af<Ag. When M=0, there is Ag=A2≤Af.
Finally, we use the above-mentioned relational expressions and combine with Proposition4; then, we can obtain the final chart as shown in Table 4.
Because positive solutions ofF(B) is equal to existences of internal equilibria of system. Thus, the existence of internal equilibria about B>M is proved.Table 4 F(B) F ( M ) g ( M ) Conclusions Solutions M ≥ 0, A=A2 – M ≥ 0, A>Ag M > 0, A=A2>Ag B +  M ≥ 0, A>A2 M > 0, A>Af – M > 0, A>Af B 2 + M = 0 M ≥ 0, A>Ag M = 0, A>A2=Ag B 2 + M > 0, A=Af M ≥ 0, A>Ag M > 0, A=Af>Ag B 2 + M > 0, A<Af M ≥ 0, A>Ag M > 0, max⁡A2,Ag<A<Af B 1 +, B2+Remark 6.
M g ≔ m 2 ( d + c + δ ) 2 a ( d - q ) ( d + c ) ( d + δ ) / β 2 δ 2 ( w + e h ) 2 ( ( d + c ) ( d + δ ) - q ( d + c + δ ) ) is the criterion when A2=Ag.After knowing internal equilibria whenB>M, we want to get the relationship between internal equilibria and R0, so we need to change the form above all. We had known the R0=βδA/(d-q)(d+c)(d+δ). Thus, we change A2 to P2≔βδA2/(d-q)(d+c)(d+δ), Ag to Pg≔βδAg/(d-q)(d+c)(d+δ), and Af to Pf≔βδAf/(d-q)(d+c)(d+δ). From the expression of these, we get R0=1<Pf when M>0 and R0=1=Pf when M=0.Theorem 7 (existence of equilibria in Figures2 and 3).
The equilibriumE0=(S0,0,0,0) always exists in Ω. Here are the existences of internal equilibria: (1)
If0<M<Mg:
(i)
IfP2>1: when 1<R0<P2, there only exists E1; when R0=P2, there exist E1, E+; when P2<R0<Pf, there exist E1, E1+, E2+; when R0=Pf, there exist E1, E2+; when R0>Pf, there exists E2+.(ii)
IfP2=1: when R0=P2, there exists E+; when 1<R0<Pf, there exist E1, E1+, E2+; when R0=Pf, there exist E1, E2+; when R0>Pf, there exists E2+.(iii)
IfP2<1: when R0=P2, there exists E+; when P2<R0≤1, there exist E1+, E2+; when 1<R0<Pf, there exist E1, E1+, E2+; when R0=Pf, there exist E1, E2+; when R0>Pf, there exists E2+.(2)
IfM≥Mg: when 1<R0≤Pf, there exists E1; when R0>Pf, there exists E2+.(3)
IfM=0 (in this case, Pf=1.): when R0>P2, there exists E2+.Figure 2
Density of the infected dairy cows as a function of the basic reproduction number with0<M<Mg. The common values of (a), (b), and (c) used for the parameters are as follows: w=3, d=2, c=0.8, δ=4, q=1, β=0.2, and a=2. The different values are (a) m=0.1, M=0.1; (b) m=0.0932, M=0.0468; and (c) m=0.05, M=0.0135. The forward bifurcation occurs at R0=1 in all the pictures and the backward bifurcation occurs at R0=1.2101 in picture (a), at R0=1.1056 in picture (b) and at R0=1.0567 in picture (c).
(a)
P 2 > 1 (b)
P 2 = 1 (c)
P 2 < 1Figure 3
The relationship between the basic reproduction number and infected dairy cows whileM≥Mg and M=0, where (2) and (3) of Theorem 9 holds. The common values of (a) and (b) used for the parameters are as follows: w=3, d=2, c=0.8, δ=4, q=1, β=0.2, a=2, and m=0.1. The different values of (a) and (b) are (a) M=0.5, (b) M=0. There is a forward bifurcation at R0=1 in (a). The bifurcation at R0=0.6383 in (b) is still forward.
(a)
M ≥ M g (b)
M = 0Theorem7 describes the existence of equilibria under the different parameters condition (Table 5), where disease-free and bacteria-free equilibrium E0 always exists. Besides, we give the stability of the equilibria, which will be proved in next part. We just give the stability of E0 and E1; others are too difficult to give the strict mathematical proofs because of lacking exact expression. Therefore, the local stability of E1+ and E2+ will be demonstrated numerically in Section 4.Table 5
Existence of endemic equilibria of model (4).
Condition 1 Condition 2 Condition 3 Endemic equilibria 0<M<Mg P 2 > 1 1 < R 0 < P 2 E 1 R 0 = P 2 E 1, E+ P 2 < R 0 < P f E 1, E1+, E2+ R0=Pf E 1, E2+ R0>Pf E2+ P 2 = 1 R 0 = P 2 E+ 1 < R 0 < P f E 1, E1+, E2+ R0=Pf E 1, E2+ R0>Pf E2+ P 2 < 1 R0=P2 E + P 2 < R 0 ≤ 1 E 1 +, E2+ 1 < R 0 < P f E 1, E1+, E2+ R0=Pf E 1, E2+ R0>Pf E2+ M≥Mg — 1 < R 0 ≤ P f E 1 R 0 > P f E 2 + M=0 — R 0 > P 2 E 2 +


## 2.1. The Dynamical Model

We classify the dairy cows into three compartments: the susceptible compartmentS, the exposed compartment E, and the infectious compartment I in Figure 1. And the Brucella in environment is denoted by B.Figure 1
A flow diagram of brucellosis transmission.One of the key differences of (2) is the incidence term α(B)S of indirect transmission of Brucella to susceptible cow [12], where α(B) is the pathogen density dependent component. Unlike [12, 20] using Holling-II response curve when the bacteria density is over the MID, a simple linear form is used in this paper. There are chiefly the following reasons. Firstly, the half-saturation pathogen density (H) is difficult to be determined in actual application. Secondly, the Brucella which is discharged to the nature by dairy cows and other livestock is hard to reach saturation. Lastly, that can get more practical conclusions with immunological threshold than [22]. The expression of α(B) in this paper is defined as(3)αB=0if0≤B≤M,aB-MifB>M.If the number of bacteria is more than MID, then the Brucella of environment can bring the indirect influence to the dairy cows. At this time,a is representative of maximum rate of indirect infection, and M is denoted MID of Brucella, B-M is representative of bacteria amount which can enter the body and take effect. If the number is less than or equal to MID, then the bacteria of environment have no effect on the susceptible and dairy cows were infected just by the direct transmission from infected to susceptible cow.So the model is a system of ordinary differential equation as follows:(4)S˙=A+qS+E+I-βSI-dS-αBS,E˙=βSI+αBS-dE-δE,I˙=δE-dI-cI,B˙=mE+I-w+ehB,where the description of other parameters which are not mentioned in Table 1 can be found in Table 2.Table 2
The description of new parameters in model (4).
Parameter Description Dimension M MID c e l l   l i t e r - 1 a Maximum rate of infection none β Indirect transmission rate noneRemark 1.
In this paper, we also use the conditiond>q like [22].


## 2.2. Forward Invariance

The first equation of (4) gives us that S(t)˙≥0, when S=0. In the same way, we can get E(t)˙≥0 and I(t)˙≥0. Hence, we know that S(t)≥0, E(t)≥0, and I(t)≥0 for t>0. As S˙+E˙+I˙≤A-(d-q)(S+E+I), we can get 0≤S+E+I≤A/d-q by using the character of S˙+E˙+I˙. Cause S˙+E˙+I˙=A-(d-q)(S+E+I) is a linear function [12]. In the same way, if B=0, then B(t)˙=0. Thus, the inequality B(t)≥0 for t>0 is true. If B(0)∈[0,m(E+I)/w+eh), then B(t)∈[0,m(E+I)/w+eh) for any t>0. To sum up, we get the following theorem.Theorem 2 (feasible region).
The set(5)Ω=S,E,I,B:0≤S+E+I≤Ad-q,0≤B<mAw+ehd-q,S≥0,E≥0,I≥0defines a forward invariant region of system (2).


## 2.3. Equilibria of System (4)

### 2.3.1. Equilibria When Bacteria Are Less Than or Equal to MID

When0≤B≤M, then α(B)=0 and the model is(6)S˙=A+qS+E+I-βSI-dS,E˙=βSI-dE-δE,I˙=δE-dI-cI,B˙=mE+I-w+ehB.We solve the equation and can getE0=(S0,0,0,0), a disease-free and bacteria-free equilibrium where S0=A/d-q, and a endemic equilibrium E1=(S1,E1,I1,B1), where(7)S1=d+δd+cβδ,E1=d+cβδA-d-qd+δd+cβδd+δd+c-qd+c+δ,I1=βδA-d-qd+δd+cβd+δd+c-qd+c+δ,B1=md+c+δβδA-d-qd+δd+cw+ehβδd+δd+c-qd+c+δ.Remark 3.
S 1 > 0, E1>0, I1>0 and B1>0.According to the spectral radius theory [25, 26], we can directly obtain the basic reproduction number based on [22] when β2=0. That is (8)R0=ρFV-1=βδAd-qd+cd+δ.Then we can attain the existence of equilibria: if R0≤1, then there is only one equilibrium E0; if R0>1, then there are disease-free and bacteria-free equilibrium E0 and endemic equilibrium E1.


### 2.3.2. Equilibria When Bacteria Are More Than MID

WhenB>M, the equation is(9)S˙=A+qS+E+I-βSI-dS-aB-MS,E˙=βSI+aB-MS-dE-δE,I˙=δE-dI-cI,B˙=mE+I-w+ehB.Obviously, there is no disease-free and bacteria-free equilibria because of M≥0. Thus, system (9) implies that(10)E=d+cw+ehBmd+c+δ,I=δw+ehBmd+c+δ,S=d+δEβI+αB=qE+I+AβI+αB+d-q.The equation ofS about B can be written as F(B)=0. Define F(B)=B2(w+eh)[d+cd+δ-qd+c+δ][βδ(w+eh)+am(d+c+δ)]+Bm(d+c+δ){-aM(w+eh)[(d+c)(d+δ)-q(d+c+δ)]-(w+eh)[βδA-(d-q)(d+δ)(d+c)]-am(d+c+δ)A}+aMAm2(d+c+δ)2. From the expression, we can know that F(B) is a quadratic function, in which the first and third coefficients are positive. Therefore, the existence of solution for F(B)=0 depends on the second coefficient and Δ. If the second coefficient is nonnegative, no matter Δ is, then F(B)=0 has no positive solution. If the second coefficient is negative, the positive solution depends on the sign of Δ: when Δ<0, there is no solution; when Δ=0, there is one positive solution; and when Δ>0, there are two positive solutions.Now we simplify them and can get the new second coefficientb1:(11)b1≔md+c+δ-βδw+eh-mad+c+δA+b11,where(12)b11≔-aMw+ehd+cd+δ-qd+c+δ+w+ehd-qd+δd+c, and the new Δ about A:(13)ΔA≔A2m2d+c+δ2-βδw+eh-amd+c+δ2+A-2m2d+c+δ2w+ehβδw+eh+amd+c+δaMd+cd+δ-qd+c+δ+d-qd+δd+c+m2d+c+δ2-aMw+ehd+cd+δ-qd+c+δ+w+ehd-qd+δd+c2.We can knowb1 is a liner function about A and Δ(A) is still a quadratic function which the first and third coefficients are positive and the second coefficient is less than zero. Then we calculate the Δ(A) and can know that the sign of Δ(Δ(A)) is always positive. So we can solve Δ(A)=0 and get A1, A2, where(14)A1=w+ehaMd+cd+δ-qd+c+δ+d-qd+cd+δ-2aMd+cd+δ-qd+c+δd-qd+cd+δβδw+eh+mad+c+δ,A2=w+ehaMd+cd+δ-qd+c+δ+d-qd+cd+δ+2aMd+cd+δ-qd+c+δd-qd+cd+δβδw+eh+mad+c+δ.By the analysis with the feature of function, we can get the conclusion as follows.Proposition 4.
The existence of positive solutions forF(B)=0:
(i) Whenb11≥0: if 0≤A<A2, then F(B)=0 does not get any positive solution; if A=A2, then F(B)=0 has one positive solution B1+; if A>A2, then F(B)=0 has two positive solutions B1+, B2+ (B1+<B2+).
(ii) Whenb11<0: if A1<A<A2, then F(B)=0 does not get any positive solution; if A=A1 or A=A2, then F(B)=0 has one positive solution B+ (=B1+=B2+); if 0≤A<A1 or A>A2, then F(B)=0 has two positive solutions B1+, B2+.Note theM~=(d-q)(d+c)(d+δ)/a((d+c)(d+δ)-q(d+c+δ)) (by the equation of b11=0) and we can put Proposition 4 to the chart as shown in Table 3.Table 3 F ( B ) = 0 has one positive solution B+ F ( B ) = 0 has two positive solutions B1+, B2+  M ≥ 0, A=A2 M ≥ 0, A>A2  M > M ~, A=A1 M > M ~, 0≤A<A1It just discusses the existence of positive solutions, but we need consider the value ofM, where only internal equilibria exist while positive solutions are more than M. Hence, we will get conclusions as follows.Theorem 5.
The existence of internal equilibria aboutB>M:(1)
IfM>0:
(i)
whenA>Af, there is an internal equilibrium E2+;(ii)
whenA=Af>Ag, there is an internal equilibrium E2+;(iii)
whenmax⁡Ag,A2<A<Af, there are two internal equilibria E1+, E2+;(iv)
whenA=A2>Ag, there is an internal equilibrium E+;(2)
IfM=0 and A>A2=Ag, then there is an internal equilibrium E2+.
where(15)Af=βδw+ehd+cd+δ-qd+c+δM+md+c+δd-qd+cd+δβδmd+c+δ,Ag=2βδw+eh2+amw+ehd+c+δd+cd+δ-qd+c+δM+w+ehmd+c+δd-qd+cd+δmd+c+δβδw+eh+mad+c+δ.Proof.
We have known the positive value ofB; now we need combine with the inequality B>M. The proof uses character of function F(B).
Firstly, we should know the value ofF(M). When B=M, we can get F(M)=βδ(w+eh)2((d+c)(d+δ)-q(d+c+δ))M2-m(d+c+δ)(w+eh)(βδA-(d-q)(d+c)(d+δ))M, and the relationship between F(M) and 0 is translated into the relationship between A and Af: (1)
WhenM>0: if A>Af, then F(M)<0, if A=Af, then F(M)=0, and if A<Af, then F(M)>0.(2)
WhenM=0, we can get F(M)=0.
Secondly, we need to know position ofM and symmetry axis of F(M). We put them into the formula g(M)≔[2βδ(w+eh)2+am(w+eh)(d+c+δ)]((d+c)(d+δ)-q(d+c+δ))M-m(d+c+δ)[(w+eh)(βδA-(d-q)(d+c)(d+δ))+ma(d+c+δ)A]. In the same way, we translate them into the relationship between A and Ag. It means that when M≥0: if A>Ag, then g(M)<0 (M is on the left of axis of symmetry); if A=Ag, then g(M)=0 (M is the axis of symmetry); and if A<Ag, then g(M)>0 (M is on the right of axis of symmetry).
Thirdly, we calculate the relationship ofAf, Ag, A1, A2. Then we can get the relationships among them as follows: (1)
The relationships ofAf>A1 and Ag>A1 are always true.(2)
When0<M≤Mg, there is Ag≤A2≤Af. When M>Mg, there is A2≤Af<Ag. When M=0, there is Ag=A2≤Af.
Finally, we use the above-mentioned relational expressions and combine with Proposition4; then, we can obtain the final chart as shown in Table 4.
Because positive solutions ofF(B) is equal to existences of internal equilibria of system. Thus, the existence of internal equilibria about B>M is proved.Table 4 F(B) F ( M ) g ( M ) Conclusions Solutions M ≥ 0, A=A2 – M ≥ 0, A>Ag M > 0, A=A2>Ag B +  M ≥ 0, A>A2 M > 0, A>Af – M > 0, A>Af B 2 + M = 0 M ≥ 0, A>Ag M = 0, A>A2=Ag B 2 + M > 0, A=Af M ≥ 0, A>Ag M > 0, A=Af>Ag B 2 + M > 0, A<Af M ≥ 0, A>Ag M > 0, max⁡A2,Ag<A<Af B 1 +, B2+Remark 6.
M g ≔ m 2 ( d + c + δ ) 2 a ( d - q ) ( d + c ) ( d + δ ) / β 2 δ 2 ( w + e h ) 2 ( ( d + c ) ( d + δ ) - q ( d + c + δ ) ) is the criterion when A2=Ag.After knowing internal equilibria whenB>M, we want to get the relationship between internal equilibria and R0, so we need to change the form above all. We had known the R0=βδA/(d-q)(d+c)(d+δ). Thus, we change A2 to P2≔βδA2/(d-q)(d+c)(d+δ), Ag to Pg≔βδAg/(d-q)(d+c)(d+δ), and Af to Pf≔βδAf/(d-q)(d+c)(d+δ). From the expression of these, we get R0=1<Pf when M>0 and R0=1=Pf when M=0.Theorem 7 (existence of equilibria in Figures2 and 3).
The equilibriumE0=(S0,0,0,0) always exists in Ω. Here are the existences of internal equilibria: (1)
If0<M<Mg:
(i)
IfP2>1: when 1<R0<P2, there only exists E1; when R0=P2, there exist E1, E+; when P2<R0<Pf, there exist E1, E1+, E2+; when R0=Pf, there exist E1, E2+; when R0>Pf, there exists E2+.(ii)
IfP2=1: when R0=P2, there exists E+; when 1<R0<Pf, there exist E1, E1+, E2+; when R0=Pf, there exist E1, E2+; when R0>Pf, there exists E2+.(iii)
IfP2<1: when R0=P2, there exists E+; when P2<R0≤1, there exist E1+, E2+; when 1<R0<Pf, there exist E1, E1+, E2+; when R0=Pf, there exist E1, E2+; when R0>Pf, there exists E2+.(2)
IfM≥Mg: when 1<R0≤Pf, there exists E1; when R0>Pf, there exists E2+.(3)
IfM=0 (in this case, Pf=1.): when R0>P2, there exists E2+.Figure 2
Density of the infected dairy cows as a function of the basic reproduction number with0<M<Mg. The common values of (a), (b), and (c) used for the parameters are as follows: w=3, d=2, c=0.8, δ=4, q=1, β=0.2, and a=2. The different values are (a) m=0.1, M=0.1; (b) m=0.0932, M=0.0468; and (c) m=0.05, M=0.0135. The forward bifurcation occurs at R0=1 in all the pictures and the backward bifurcation occurs at R0=1.2101 in picture (a), at R0=1.1056 in picture (b) and at R0=1.0567 in picture (c).
(a)
P 2 > 1 (b)
P 2 = 1 (c)
P 2 < 1Figure 3
The relationship between the basic reproduction number and infected dairy cows whileM≥Mg and M=0, where (2) and (3) of Theorem 9 holds. The common values of (a) and (b) used for the parameters are as follows: w=3, d=2, c=0.8, δ=4, q=1, β=0.2, a=2, and m=0.1. The different values of (a) and (b) are (a) M=0.5, (b) M=0. There is a forward bifurcation at R0=1 in (a). The bifurcation at R0=0.6383 in (b) is still forward.
(a)
M ≥ M g (b)
M = 0Theorem7 describes the existence of equilibria under the different parameters condition (Table 5), where disease-free and bacteria-free equilibrium E0 always exists. Besides, we give the stability of the equilibria, which will be proved in next part. We just give the stability of E0 and E1; others are too difficult to give the strict mathematical proofs because of lacking exact expression. Therefore, the local stability of E1+ and E2+ will be demonstrated numerically in Section 4.Table 5
Existence of endemic equilibria of model (4).
Condition 1 Condition 2 Condition 3 Endemic equilibria 0<M<Mg P 2 > 1 1 < R 0 < P 2 E 1 R 0 = P 2 E 1, E+ P 2 < R 0 < P f E 1, E1+, E2+ R0=Pf E 1, E2+ R0>Pf E2+ P 2 = 1 R 0 = P 2 E+ 1 < R 0 < P f E 1, E1+, E2+ R0=Pf E 1, E2+ R0>Pf E2+ P 2 < 1 R0=P2 E + P 2 < R 0 ≤ 1 E 1 +, E2+ 1 < R 0 < P f E 1, E1+, E2+ R0=Pf E 1, E2+ R0>Pf E2+ M≥Mg — 1 < R 0 ≤ P f E 1 R 0 > P f E 2 + M=0 — R 0 > P 2 E 2 +


## 2.3.1. Equilibria When Bacteria Are Less Than or Equal to MID

When0≤B≤M, then α(B)=0 and the model is(6)S˙=A+qS+E+I-βSI-dS,E˙=βSI-dE-δE,I˙=δE-dI-cI,B˙=mE+I-w+ehB.We solve the equation and can getE0=(S0,0,0,0), a disease-free and bacteria-free equilibrium where S0=A/d-q, and a endemic equilibrium E1=(S1,E1,I1,B1), where(7)S1=d+δd+cβδ,E1=d+cβδA-d-qd+δd+cβδd+δd+c-qd+c+δ,I1=βδA-d-qd+δd+cβd+δd+c-qd+c+δ,B1=md+c+δβδA-d-qd+δd+cw+ehβδd+δd+c-qd+c+δ.Remark 3.
S 1 > 0, E1>0, I1>0 and B1>0.According to the spectral radius theory [25, 26], we can directly obtain the basic reproduction number based on [22] when β2=0. That is (8)R0=ρFV-1=βδAd-qd+cd+δ.Then we can attain the existence of equilibria: if R0≤1, then there is only one equilibrium E0; if R0>1, then there are disease-free and bacteria-free equilibrium E0 and endemic equilibrium E1.


## 2.3.2. Equilibria When Bacteria Are More Than MID

WhenB>M, the equation is(9)S˙=A+qS+E+I-βSI-dS-aB-MS,E˙=βSI+aB-MS-dE-δE,I˙=δE-dI-cI,B˙=mE+I-w+ehB.Obviously, there is no disease-free and bacteria-free equilibria because of M≥0. Thus, system (9) implies that(10)E=d+cw+ehBmd+c+δ,I=δw+ehBmd+c+δ,S=d+δEβI+αB=qE+I+AβI+αB+d-q.The equation ofS about B can be written as F(B)=0. Define F(B)=B2(w+eh)[d+cd+δ-qd+c+δ][βδ(w+eh)+am(d+c+δ)]+Bm(d+c+δ){-aM(w+eh)[(d+c)(d+δ)-q(d+c+δ)]-(w+eh)[βδA-(d-q)(d+δ)(d+c)]-am(d+c+δ)A}+aMAm2(d+c+δ)2. From the expression, we can know that F(B) is a quadratic function, in which the first and third coefficients are positive. Therefore, the existence of solution for F(B)=0 depends on the second coefficient and Δ. If the second coefficient is nonnegative, no matter Δ is, then F(B)=0 has no positive solution. If the second coefficient is negative, the positive solution depends on the sign of Δ: when Δ<0, there is no solution; when Δ=0, there is one positive solution; and when Δ>0, there are two positive solutions.Now we simplify them and can get the new second coefficientb1:(11)b1≔md+c+δ-βδw+eh-mad+c+δA+b11,where(12)b11≔-aMw+ehd+cd+δ-qd+c+δ+w+ehd-qd+δd+c, and the new Δ about A:(13)ΔA≔A2m2d+c+δ2-βδw+eh-amd+c+δ2+A-2m2d+c+δ2w+ehβδw+eh+amd+c+δaMd+cd+δ-qd+c+δ+d-qd+δd+c+m2d+c+δ2-aMw+ehd+cd+δ-qd+c+δ+w+ehd-qd+δd+c2.We can knowb1 is a liner function about A and Δ(A) is still a quadratic function which the first and third coefficients are positive and the second coefficient is less than zero. Then we calculate the Δ(A) and can know that the sign of Δ(Δ(A)) is always positive. So we can solve Δ(A)=0 and get A1, A2, where(14)A1=w+ehaMd+cd+δ-qd+c+δ+d-qd+cd+δ-2aMd+cd+δ-qd+c+δd-qd+cd+δβδw+eh+mad+c+δ,A2=w+ehaMd+cd+δ-qd+c+δ+d-qd+cd+δ+2aMd+cd+δ-qd+c+δd-qd+cd+δβδw+eh+mad+c+δ.By the analysis with the feature of function, we can get the conclusion as follows.Proposition 4.
The existence of positive solutions forF(B)=0:
(i) Whenb11≥0: if 0≤A<A2, then F(B)=0 does not get any positive solution; if A=A2, then F(B)=0 has one positive solution B1+; if A>A2, then F(B)=0 has two positive solutions B1+, B2+ (B1+<B2+).
(ii) Whenb11<0: if A1<A<A2, then F(B)=0 does not get any positive solution; if A=A1 or A=A2, then F(B)=0 has one positive solution B+ (=B1+=B2+); if 0≤A<A1 or A>A2, then F(B)=0 has two positive solutions B1+, B2+.Note theM~=(d-q)(d+c)(d+δ)/a((d+c)(d+δ)-q(d+c+δ)) (by the equation of b11=0) and we can put Proposition 4 to the chart as shown in Table 3.Table 3 F ( B ) = 0 has one positive solution B+ F ( B ) = 0 has two positive solutions B1+, B2+  M ≥ 0, A=A2 M ≥ 0, A>A2  M > M ~, A=A1 M > M ~, 0≤A<A1It just discusses the existence of positive solutions, but we need consider the value ofM, where only internal equilibria exist while positive solutions are more than M. Hence, we will get conclusions as follows.Theorem 5.
The existence of internal equilibria aboutB>M:(1)
IfM>0:
(i)
whenA>Af, there is an internal equilibrium E2+;(ii)
whenA=Af>Ag, there is an internal equilibrium E2+;(iii)
whenmax⁡Ag,A2<A<Af, there are two internal equilibria E1+, E2+;(iv)
whenA=A2>Ag, there is an internal equilibrium E+;(2)
IfM=0 and A>A2=Ag, then there is an internal equilibrium E2+.
where(15)Af=βδw+ehd+cd+δ-qd+c+δM+md+c+δd-qd+cd+δβδmd+c+δ,Ag=2βδw+eh2+amw+ehd+c+δd+cd+δ-qd+c+δM+w+ehmd+c+δd-qd+cd+δmd+c+δβδw+eh+mad+c+δ.Proof.
We have known the positive value ofB; now we need combine with the inequality B>M. The proof uses character of function F(B).
Firstly, we should know the value ofF(M). When B=M, we can get F(M)=βδ(w+eh)2((d+c)(d+δ)-q(d+c+δ))M2-m(d+c+δ)(w+eh)(βδA-(d-q)(d+c)(d+δ))M, and the relationship between F(M) and 0 is translated into the relationship between A and Af: (1)
WhenM>0: if A>Af, then F(M)<0, if A=Af, then F(M)=0, and if A<Af, then F(M)>0.(2)
WhenM=0, we can get F(M)=0.
Secondly, we need to know position ofM and symmetry axis of F(M). We put them into the formula g(M)≔[2βδ(w+eh)2+am(w+eh)(d+c+δ)]((d+c)(d+δ)-q(d+c+δ))M-m(d+c+δ)[(w+eh)(βδA-(d-q)(d+c)(d+δ))+ma(d+c+δ)A]. In the same way, we translate them into the relationship between A and Ag. It means that when M≥0: if A>Ag, then g(M)<0 (M is on the left of axis of symmetry); if A=Ag, then g(M)=0 (M is the axis of symmetry); and if A<Ag, then g(M)>0 (M is on the right of axis of symmetry).
Thirdly, we calculate the relationship ofAf, Ag, A1, A2. Then we can get the relationships among them as follows: (1)
The relationships ofAf>A1 and Ag>A1 are always true.(2)
When0<M≤Mg, there is Ag≤A2≤Af. When M>Mg, there is A2≤Af<Ag. When M=0, there is Ag=A2≤Af.
Finally, we use the above-mentioned relational expressions and combine with Proposition4; then, we can obtain the final chart as shown in Table 4.
Because positive solutions ofF(B) is equal to existences of internal equilibria of system. Thus, the existence of internal equilibria about B>M is proved.Table 4 F(B) F ( M ) g ( M ) Conclusions Solutions M ≥ 0, A=A2 – M ≥ 0, A>Ag M > 0, A=A2>Ag B +  M ≥ 0, A>A2 M > 0, A>Af – M > 0, A>Af B 2 + M = 0 M ≥ 0, A>Ag M = 0, A>A2=Ag B 2 + M > 0, A=Af M ≥ 0, A>Ag M > 0, A=Af>Ag B 2 + M > 0, A<Af M ≥ 0, A>Ag M > 0, max⁡A2,Ag<A<Af B 1 +, B2+Remark 6.
M g ≔ m 2 ( d + c + δ ) 2 a ( d - q ) ( d + c ) ( d + δ ) / β 2 δ 2 ( w + e h ) 2 ( ( d + c ) ( d + δ ) - q ( d + c + δ ) ) is the criterion when A2=Ag.After knowing internal equilibria whenB>M, we want to get the relationship between internal equilibria and R0, so we need to change the form above all. We had known the R0=βδA/(d-q)(d+c)(d+δ). Thus, we change A2 to P2≔βδA2/(d-q)(d+c)(d+δ), Ag to Pg≔βδAg/(d-q)(d+c)(d+δ), and Af to Pf≔βδAf/(d-q)(d+c)(d+δ). From the expression of these, we get R0=1<Pf when M>0 and R0=1=Pf when M=0.Theorem 7 (existence of equilibria in Figures2 and 3).
The equilibriumE0=(S0,0,0,0) always exists in Ω. Here are the existences of internal equilibria: (1)
If0<M<Mg:
(i)
IfP2>1: when 1<R0<P2, there only exists E1; when R0=P2, there exist E1, E+; when P2<R0<Pf, there exist E1, E1+, E2+; when R0=Pf, there exist E1, E2+; when R0>Pf, there exists E2+.(ii)
IfP2=1: when R0=P2, there exists E+; when 1<R0<Pf, there exist E1, E1+, E2+; when R0=Pf, there exist E1, E2+; when R0>Pf, there exists E2+.(iii)
IfP2<1: when R0=P2, there exists E+; when P2<R0≤1, there exist E1+, E2+; when 1<R0<Pf, there exist E1, E1+, E2+; when R0=Pf, there exist E1, E2+; when R0>Pf, there exists E2+.(2)
IfM≥Mg: when 1<R0≤Pf, there exists E1; when R0>Pf, there exists E2+.(3)
IfM=0 (in this case, Pf=1.): when R0>P2, there exists E2+.Figure 2
Density of the infected dairy cows as a function of the basic reproduction number with0<M<Mg. The common values of (a), (b), and (c) used for the parameters are as follows: w=3, d=2, c=0.8, δ=4, q=1, β=0.2, and a=2. The different values are (a) m=0.1, M=0.1; (b) m=0.0932, M=0.0468; and (c) m=0.05, M=0.0135. The forward bifurcation occurs at R0=1 in all the pictures and the backward bifurcation occurs at R0=1.2101 in picture (a), at R0=1.1056 in picture (b) and at R0=1.0567 in picture (c).
(a)
P 2 > 1 (b)
P 2 = 1 (c)
P 2 < 1Figure 3
The relationship between the basic reproduction number and infected dairy cows whileM≥Mg and M=0, where (2) and (3) of Theorem 9 holds. The common values of (a) and (b) used for the parameters are as follows: w=3, d=2, c=0.8, δ=4, q=1, β=0.2, a=2, and m=0.1. The different values of (a) and (b) are (a) M=0.5, (b) M=0. There is a forward bifurcation at R0=1 in (a). The bifurcation at R0=0.6383 in (b) is still forward.
(a)
M ≥ M g (b)
M = 0Theorem7 describes the existence of equilibria under the different parameters condition (Table 5), where disease-free and bacteria-free equilibrium E0 always exists. Besides, we give the stability of the equilibria, which will be proved in next part. We just give the stability of E0 and E1; others are too difficult to give the strict mathematical proofs because of lacking exact expression. Therefore, the local stability of E1+ and E2+ will be demonstrated numerically in Section 4.Table 5
Existence of endemic equilibria of model (4).
Condition 1 Condition 2 Condition 3 Endemic equilibria 0<M<Mg P 2 > 1 1 < R 0 < P 2 E 1 R 0 = P 2 E 1, E+ P 2 < R 0 < P f E 1, E1+, E2+ R0=Pf E 1, E2+ R0>Pf E2+ P 2 = 1 R 0 = P 2 E+ 1 < R 0 < P f E 1, E1+, E2+ R0=Pf E 1, E2+ R0>Pf E2+ P 2 < 1 R0=P2 E + P 2 < R 0 ≤ 1 E 1 +, E2+ 1 < R 0 < P f E 1, E1+, E2+ R0=Pf E 1, E2+ R0>Pf E2+ M≥Mg — 1 < R 0 ≤ P f E 1 R 0 > P f E 2 + M=0 — R 0 > P 2 E 2 +


## 3. Stability of Equilibria

### 3.1. Local Stability ofE0 and E1

We calculate the Jacobian to analyze the local stability of each of the equilibria. When0≤B≤M, the Jacobian is(16)JS,E,I,B=q-d-βIqq-βS0βI-d+δβS00δ-d+c00mm-w+eh.Now considering E0, we use λ for eigenvalues. We can compute(17)det⁡λI-JE0=λ-q-dλ+w+ehλ+d+δλ+d+c-βδS0.Note h(λ)=(λ+d+δ)(λ+d+c)-βδS0. We can find that λ1=-(d-q)<0, λ2=-(w+eh)<0. Now we just consider the eigenvalues of h(λ). The Routh-Hurwitz coefficients of the h(λ) are(18)Δ1=d+c+d+δ>0,Δ2=d+c+d+δd+δd+c-βδS0.If R0<1, then Δ2>0; if R0>1, then Δ2<0. Thus the eigenvalues of h(λ) are negative when R0<1. So the disease-free and bacteria-free equilibria E0 are locally asymptotically stable while R0<1. If R0>1, then E0 is a saddle-point equilibrium.Now consideringE1. The characteristic polynomial of Jacobian of E1 is(19)det⁡λI-JE1=a0λ3+a1λ2+a2λ+a3=0, where(20)a0=1,a1=βδA-qd-qd+c+δd+δd+c-qd+c+δ+d+c+d+δ,a2=βδA-qd-qd+c+δd+δ+d+cd+δd+c-qd+c+δ,a3=βδA-d-qd+δd+c. The Routh-Hurwitz coefficients of the det⁡(λI-JE0) are(21)Δ1=a1,Δ2=a1a2-a3,Δ3=a3a1a2-a3. where a1>βδA-q(d-q)(d+c+δ)/(d+δ)(d+c)-q(d+c+δ)>βδA-(d-q)(d+δ)(d+c)/(d+δ)(d+c)-q(d+c+δ)>0 and a3>0, when R0>1.As for the conditiona1a2-a3, we have the following expression:(22)a1a2-a3=h1+h2+h3, where(23)h1=βδA-qd-qd+c+δβδA-qd-qd+c+δd+c+d+δ-qβδA-qd-qd+cd+δd+cd+δ-qd+c+δ2,h2=βδA-qd-qd+c+δd+c+d+δ2-d+c+d+δqβδA-qd-qd+cd+δd+cd+δ-qd+c+δ,h3=-βδA-d-qd+δd+c.Now, we magnify them:(24)h1>d-qd+cd+δ-qd-qd+c+δβδA-qd-qd+c+δd+c+d+δ-qβδA-qd-qd+cd+δd+δd+c-qd+c+δ2=d-qβδA-qd-qd+c+δd+c+d+δ-qβδA-d-qd+cd+δd+δd+c-qd+c+δ,h1+h2>βδA-qd-qd+c+δd+c+d+δd+δ+d+c+d-q-qβδA-qd-qd+cd+δd+δ+d+c+d-qd+δd+c-qd+c+δ,h1+h2+h3>βδA-qd-qd+c+δd+c+d+δd+δ+d+c+d-q-βδA-d-qd+cd+δ2dq-q2+d+cd+δd+δd+c-qd+c+δ>βδA-d-qd+cd+δd+c+d+δd+c+d+δ+d-q-βδA-d-qd+cd+δ2dq-q2+d+cd+δd+δd+c-qd+c+δ>βδA-d-qd+cd+δd+c2+d+δ+d-qd+δ+d+c+q2-2d2d+δd+c-qd+c+δ.We know thatd>q, so (d+c)(d+δ)>q(d+c+δ) and then a1a2-a3>0. Hence, the Routh-Hurwitz conditions are satisfied when R0>1 and E1 is locally asymptotically stable.Theorem 8 (local stability).
System (4) has equilibria as follows:
When0≤B≤M and R0<1, there is only one equilibrium E0 which is locally asymptotically stable. When 0≤B≤M and R0>1, E0 is unstable and a unique endemic equilibrium E1 exists which is locally asymptotically stable.We study the stability of equilibria in order to explain more epidemiological implications. Taking Figure2(a) as an example, we can get the following biological significance. If E0 is stable (R0<1), the brucellosis will fade as time goes by. If E1 is stable (1<R0<P2), disease will maintain the level of I1 all the time. With the scope of R0 (P2<R0<Pf), disease will continue at low-level state (I1) when the current status is in the attraction domain of E1 and will be going to stay at high-level state (I2+) while the current status exactly situate in the attraction domain of E2+, which will be demonstrated numerically in Section 4. Similarly, disease will be going to stay at the level of I2+ only if R0>Pf is true.


### 3.2. Global Stability ofE0 and E1

In this section, we will demonstrate the global stability of the endemic equilibrium of system (4) by using the Lyapunovfunction [27, 28].Theorem 9.
If0≤B≤M and R0<1, then E0 is globally asymptotically stable in Ω.Proof.
As we can see,B is not acting on the other subclasses. So we can define the Lyapunovfunction(25)VS,E,I,B=S-S0-S0ln⁡SS0+E+d+δδI.Then the derivative of V along system (6) is(26)
d
V
d
t
=
1
-
S
0
S
A
-
d
-
q
S
+
q
E
+
I
-
β
S
I
+
β
S
I
-
d
+
δ
E
+
d
+
δ
δ
δ
E
-
d
+
c
I
=
1
-
S
0
S
A
-
d
-
q
S
+
q
E
+
I
+
β
S
0
I
-
d
+
δ
E
+
d
+
δ
δ
δ
E
-
d
+
c
I
=
1
-
S
0
S
A
-
d
-
q
S
+
q
E
+
I
+
I
β
δ
A
-
d
-
q
d
+
c
d
+
δ
δ
d
-
q
.
According to the definition of the invariant region, we can know thatS≤S0. So dV/dt<0 is true and dV/dt=0 holds only if S=S0, when R0<1. dV/dt=0 holds only if S=S0 when R0=1. Thus the disease-free and bacteria-free equilibrium E0 is globally asymptotically stable by LaSalle′sInvariancePrinciple [29].Theorem 10.
If0≤B≤M, R0>1, and q=0, then E1 is a conditionally globally asymptotically stable in Ω.Proof.
In this proof, we use the same way to demonstrate the stability whenq=0. Similarly, B is not acting on the other subclasses, so we define the Lyapunovfunction as follows:(27)VS,E,I,B=S-S1-S1ln⁡SS1+E-E1-E1ln⁡EE1+d+δδI-I1-I1ln⁡II1. Then the derivative of V along system (6) is(28)dVdt=1-S1SA-dS-βSI+1-E1EβSI-d+δE+d+δδ1-I1IδE-d+cI. Let s=S/S1, e=E/E1, and i=I/I1. Then we can get(29)A-dS1-βS1I1=0,βS1I1-d+δE1=0,δE1-d+cI1=0. Combining with them, we can obtain that(30)dVdt=S12d-ds-ds+βS1I11-1s+i-sie+E1d+δ-d+δei+I1-d+δd+ciδ+d+cd+δδ=dS12-s-1s+βS1I13-1s-sie-ei.Combining with inequality of arithmetic and geometric means, we can know 2-s-1/s≤0 and (3-1/s-si/e-e/i)≤0. Thus, (31)dVdt≤0.
According to the definition of the invariant region, only ifs=e=i=1 means that S=S1, E=E1 and I=I1, then dV/dt=0. Similarly with Theorem 9, we can get the endemic equilibrium E1 which is globally asymptotically stable by LaSalle′sInvariancePrinciple in Ω when q=0.


## 3.1. Local Stability ofE0 and E1

We calculate the Jacobian to analyze the local stability of each of the equilibria. When0≤B≤M, the Jacobian is(16)JS,E,I,B=q-d-βIqq-βS0βI-d+δβS00δ-d+c00mm-w+eh.Now considering E0, we use λ for eigenvalues. We can compute(17)det⁡λI-JE0=λ-q-dλ+w+ehλ+d+δλ+d+c-βδS0.Note h(λ)=(λ+d+δ)(λ+d+c)-βδS0. We can find that λ1=-(d-q)<0, λ2=-(w+eh)<0. Now we just consider the eigenvalues of h(λ). The Routh-Hurwitz coefficients of the h(λ) are(18)Δ1=d+c+d+δ>0,Δ2=d+c+d+δd+δd+c-βδS0.If R0<1, then Δ2>0; if R0>1, then Δ2<0. Thus the eigenvalues of h(λ) are negative when R0<1. So the disease-free and bacteria-free equilibria E0 are locally asymptotically stable while R0<1. If R0>1, then E0 is a saddle-point equilibrium.Now consideringE1. The characteristic polynomial of Jacobian of E1 is(19)det⁡λI-JE1=a0λ3+a1λ2+a2λ+a3=0, where(20)a0=1,a1=βδA-qd-qd+c+δd+δd+c-qd+c+δ+d+c+d+δ,a2=βδA-qd-qd+c+δd+δ+d+cd+δd+c-qd+c+δ,a3=βδA-d-qd+δd+c. The Routh-Hurwitz coefficients of the det⁡(λI-JE0) are(21)Δ1=a1,Δ2=a1a2-a3,Δ3=a3a1a2-a3. where a1>βδA-q(d-q)(d+c+δ)/(d+δ)(d+c)-q(d+c+δ)>βδA-(d-q)(d+δ)(d+c)/(d+δ)(d+c)-q(d+c+δ)>0 and a3>0, when R0>1.As for the conditiona1a2-a3, we have the following expression:(22)a1a2-a3=h1+h2+h3, where(23)h1=βδA-qd-qd+c+δβδA-qd-qd+c+δd+c+d+δ-qβδA-qd-qd+cd+δd+cd+δ-qd+c+δ2,h2=βδA-qd-qd+c+δd+c+d+δ2-d+c+d+δqβδA-qd-qd+cd+δd+cd+δ-qd+c+δ,h3=-βδA-d-qd+δd+c.Now, we magnify them:(24)h1>d-qd+cd+δ-qd-qd+c+δβδA-qd-qd+c+δd+c+d+δ-qβδA-qd-qd+cd+δd+δd+c-qd+c+δ2=d-qβδA-qd-qd+c+δd+c+d+δ-qβδA-d-qd+cd+δd+δd+c-qd+c+δ,h1+h2>βδA-qd-qd+c+δd+c+d+δd+δ+d+c+d-q-qβδA-qd-qd+cd+δd+δ+d+c+d-qd+δd+c-qd+c+δ,h1+h2+h3>βδA-qd-qd+c+δd+c+d+δd+δ+d+c+d-q-βδA-d-qd+cd+δ2dq-q2+d+cd+δd+δd+c-qd+c+δ>βδA-d-qd+cd+δd+c+d+δd+c+d+δ+d-q-βδA-d-qd+cd+δ2dq-q2+d+cd+δd+δd+c-qd+c+δ>βδA-d-qd+cd+δd+c2+d+δ+d-qd+δ+d+c+q2-2d2d+δd+c-qd+c+δ.We know thatd>q, so (d+c)(d+δ)>q(d+c+δ) and then a1a2-a3>0. Hence, the Routh-Hurwitz conditions are satisfied when R0>1 and E1 is locally asymptotically stable.Theorem 8 (local stability).
System (4) has equilibria as follows:
When0≤B≤M and R0<1, there is only one equilibrium E0 which is locally asymptotically stable. When 0≤B≤M and R0>1, E0 is unstable and a unique endemic equilibrium E1 exists which is locally asymptotically stable.We study the stability of equilibria in order to explain more epidemiological implications. Taking Figure2(a) as an example, we can get the following biological significance. If E0 is stable (R0<1), the brucellosis will fade as time goes by. If E1 is stable (1<R0<P2), disease will maintain the level of I1 all the time. With the scope of R0 (P2<R0<Pf), disease will continue at low-level state (I1) when the current status is in the attraction domain of E1 and will be going to stay at high-level state (I2+) while the current status exactly situate in the attraction domain of E2+, which will be demonstrated numerically in Section 4. Similarly, disease will be going to stay at the level of I2+ only if R0>Pf is true.


## 3.2. Global Stability ofE0 and E1

In this section, we will demonstrate the global stability of the endemic equilibrium of system (4) by using the Lyapunovfunction [27, 28].Theorem 9.
If0≤B≤M and R0<1, then E0 is globally asymptotically stable in Ω.Proof.
As we can see,B is not acting on the other subclasses. So we can define the Lyapunovfunction(25)VS,E,I,B=S-S0-S0ln⁡SS0+E+d+δδI.Then the derivative of V along system (6) is(26)
d
V
d
t
=
1
-
S
0
S
A
-
d
-
q
S
+
q
E
+
I
-
β
S
I
+
β
S
I
-
d
+
δ
E
+
d
+
δ
δ
δ
E
-
d
+
c
I
=
1
-
S
0
S
A
-
d
-
q
S
+
q
E
+
I
+
β
S
0
I
-
d
+
δ
E
+
d
+
δ
δ
δ
E
-
d
+
c
I
=
1
-
S
0
S
A
-
d
-
q
S
+
q
E
+
I
+
I
β
δ
A
-
d
-
q
d
+
c
d
+
δ
δ
d
-
q
.
According to the definition of the invariant region, we can know thatS≤S0. So dV/dt<0 is true and dV/dt=0 holds only if S=S0, when R0<1. dV/dt=0 holds only if S=S0 when R0=1. Thus the disease-free and bacteria-free equilibrium E0 is globally asymptotically stable by LaSalle′sInvariancePrinciple [29].Theorem 10.
If0≤B≤M, R0>1, and q=0, then E1 is a conditionally globally asymptotically stable in Ω.Proof.
In this proof, we use the same way to demonstrate the stability whenq=0. Similarly, B is not acting on the other subclasses, so we define the Lyapunovfunction as follows:(27)VS,E,I,B=S-S1-S1ln⁡SS1+E-E1-E1ln⁡EE1+d+δδI-I1-I1ln⁡II1. Then the derivative of V along system (6) is(28)dVdt=1-S1SA-dS-βSI+1-E1EβSI-d+δE+d+δδ1-I1IδE-d+cI. Let s=S/S1, e=E/E1, and i=I/I1. Then we can get(29)A-dS1-βS1I1=0,βS1I1-d+δE1=0,δE1-d+cI1=0. Combining with them, we can obtain that(30)dVdt=S12d-ds-ds+βS1I11-1s+i-sie+E1d+δ-d+δei+I1-d+δd+ciδ+d+cd+δδ=dS12-s-1s+βS1I13-1s-sie-ei.Combining with inequality of arithmetic and geometric means, we can know 2-s-1/s≤0 and (3-1/s-si/e-e/i)≤0. Thus, (31)dVdt≤0.
According to the definition of the invariant region, only ifs=e=i=1 means that S=S1, E=E1 and I=I1, then dV/dt=0. Similarly with Theorem 9, we can get the endemic equilibrium E1 which is globally asymptotically stable by LaSalle′sInvariancePrinciple in Ω when q=0.


## 4. Numerical Results

In an epidemic model, the basic reproduction numberR0 is calculated and verified to be a threshold for the dynamics of the disease. The ultimate objective is to control the disease by making the basic reproduction number R0. In the previous analysis of this paper, we analyze the different equilibria in variant ranges of R0. Now, we should know the final infectious size of dairy cows depending on the R0.In Figure4(a), only if R0<min⁡{1,P2} (Figures 2(a), 2(b), 2(c), 3(a), and 3(b)), we can know that the number of infected dairy cows will obtain the same disease-free and bacteria-free equilibrium E0 no matter how much the initial values are and the relationship about M, Mg, and 0 is. In biology, this means only if the basic reproduction number demands the formula, the final infectious size of dairy cows is tending to zero. It is to say, the disease will vanish without any measures after a certain time.Figure 4
Time series of the infected dairy cows. The common values of (a) and (b) used for the parameters are as follows:w=3, d=2, c=0.8, δ=4, q=1, β=0.2, a=2, m=0.2. The different values of (a) and (b) are (a) A=10 and (b) A=50. The solid line or S1 means that I is entering system (6); the dotted line or S2 indicates that I is going into system (9).
(a)
R 0 < min ⁡ 1 , P 2 (b)
R 0 > max ⁡ P 2 , P fFigure4(b) demonstrates that when R0>max⁡{P2,Pf}, the infected dairy cows will obtain the endemic equilibrium E2+. In the meantime, different values of M have different equilibria: when M=0.5, I2+=15.3531; when M=0.1, I2+=15.9432; and when M=0, I2+=16.0625. This means that only if the basic reproduction number satisfies the above formula, brucellosis will never disappear and it would become epidemic disease. In addition, the more size of MID will lead to the less final infected dairy cows. It might be easy to understand that when the MID is heightened, the bacteria invading the body are diminished by immunology theory [10].Figures5 and 7 show that when min⁡1,P2<R0<Pf, the final size of infectious cows will be in bistable state. Figure 5 comes from Figure 2(a) when P2<R0<Pf and Figure 7 comes from Figure 2(c) when 1<R0<Pf, in which Figures 5(a) and 7(a) are time series plots of the number of Brucella and Figures 5(b) and 7(b) are time series plots of the number of infected dairy cows. Figures 6 and 8 show tendency of S, E, I. From them we can see that when R0 meets the relational expression, different initial values will eventually lead to two equilibria; one is E1 and the other is E2+ (E1<E2+). In biology, brucellosis will become endemic in dairy cows with two different stable states.Figure 5
The common values of (a) and (b) used for the parameters are as follows:w=3, d=2, c=0.8, δ=4, q=1, β=0.2, a=2, m=0.1, M=0.1, and A=25. Different colors are representative of different initial values, of which the dot line means that it enters system (6) and the full line means that it goes into system (9) as in Figure 4.
(a)
The tendency of bacteria (b)
The tendency of infected dairy cowsFigure 6
The tendency ofS, E, and I.Figure 7
The common values of (a) and (b) used for the parameters are as follows:w=3, d=2, c=0.8, δ=4, q=1, β=0.2, a=2, m=0.05, M=0.014, and A=22.
(a)
The tendency of bacteria (b)
The tendency of infected dairy cowsFigure 8
The tendency ofS, E, and I.


## 5. Conclusion and Discussion

Brucellosis is becoming more and more popular in recent years, since it has formed serious endemic disease and it cannot be thoroughly removed in northeast and northwest China. In this paper, we take account of immune influence of Brucella just like researching of cholera [12]. We build the new brucellosis model based on [22]. Compared with other studies, we are more concerned with the effects of immunology on Brucella. We just change the transmission rate from constant β2 to the piecewise continuous function between the research [12] and our paper. But the mathematical results we obtain are more complicated and more biological. The original work obtained the basic reproduction number and got two equilibria and then simply analyzed the stability state. In this paper, we take immunological threshold into consideration; then we get the whole mathematical analysis of equilibria in each case and do the locally and globally stable analysis about some equilibria. The MID of Brucella has not been tested by experiment, so we get the serious discussion with M>0 and M=0: when 0<M<Mg (the formation of Mg is given in this paper), we get the forward bifurcation occurring at R0=1 and backward bifurcation appearing at R0=P2; when M≥Mg, it just has forward bifurcation at R0=1; when M=0, the model is exactly same as in [22] and produces the forward bifurcation at R0=P2. By mathematical analysis of them, we can get the global stability E0 when R0<1 and a global stability E1 when R0>1. The character of E1,2+ is given by numerical simulations. From that, we know it will be in bistable state when the basic reproduction number demands max{1,P2}<R0<Pf. We make the basic reproduction number as the measuring standard in order to give suggestions of the brucellosis prevention [30, 31]. While the basic reproduction number is less than 1, brucellosis is becoming a disease-free and bacteria-free stable state. After a certain time, brucellosis will be removed. And while the basic reproduction number is more than 1 by recording to different relationship about M, Mg and 0, the brucellosis will be endemic disease which will hold a stable state with variant values of infected dairy cows. Using the analysis of this paper, we can reduce the number of infected dairy cows to prevent the spreading of epidemic disease, which can even obtain zero by adjusting the size of the basic reproduction number.Refs. [20, 21] using Holling-II functional response both show that MID is vital for the spread of cholera. Reducing the pathogen density to under the MID in an aquatic reservoir would be more effective measure to control cholera. Similarly, in order to control and prevent the spreading of brucellosis, we can place an emphasis on the incident rate because of the existence of MID. If we can adjust the initial value of Brucella to be below the MID and constrain growth and reproduction conditions of Brucella in environment, then we can obtain the incidence rate of indirect transmission of Brucella to susceptible cow to be zero all the time and the ultimate scale of brucellosis will decrease and even disappear.In further planning, the incidence termα(B) can be used in the formation by a Holling-II response which is more reasonable in biology than the linear form in this paper. In addition, we have not yet captured and processed the real data about brucellosis of dairy cows to get more effective conclusions. That is what we are going to do next. Furthermore, the real spreading of brucellosis is via multiple zones rather than one zone and the multizone patchy model can describe the transmission in multiple regions. Hence, with considering the network models [32, 33] or reaction-diffusion equations [34–38], conclusions can accord well with practice even more.



---

*Source: 6526589-2019-05-23.xml*