Asymptotic behavior of an SEIR epidemic model with quadratic treatment

Abstract

Treatment is of great importance in fighting against infectious diseases. In this paper, we investigate a disease transmission model of SEIR type with quadratic treatment function. We consider r as the parameter which describes the societal effort to fight the infection. It is found that the system has four possible equilibria when r changes. Moreover, the existence and stability of equilibria of the model are discussed. Finally, the complex dynamics of the system are displayed by numerical simulations such as bi-stability, backward bifurcations.

1 Introduction

Kermack and Mckendrick [1] constructed a system of ordinary differential equations to study epidemiology in 1927. From then on, mathematical models have become important tools in analyzing the spread and control of epidemic diseases. Various epidemic models have been formulated and analyzed by many researchers (for example [2–12]).

In the course of disease control, the treatment is an important method to decrease the spread of diseases. Recently, many researchers have explored the role of treatment function T(I), which represents the number of infectious individuals of healing, within epidemic disease dynamic models. In classical epidemic models, the treatment rate is assumed to be proportional to the number of the infectives. That is, it is usually assumed that the capacity for the treatment of a disease in a community is a constant treatment rate 0≤r≤1, in this case the treatment function is rI, i.e., T(I)=rI.

For recent work, using an SIR model, Wang and Ruan [13] supposed that the treatment function has the form

$$T(I)=\left \{ \begin{array}{l@{\quad}l} r, & I>0,\\ 0, & I=0, \end{array} \right . $$

where r is a constant value. They discussed the stability of equilibria and proved the model exhibits Bogdanov-Takens bifurcation, Hopf bifurcation and Homo-clinic bifurcation. And they investigated the effect of the magnitude of parameter r on the behavior of the system. However, for the new outbreak disease such as SARS, etc., there is short of effective treatment techniques at the beginning of the outbreak. As the improving of the hospital’s treatment conditions, such as effective medicines, skillful techniques, the treatment rate will be increased. At last, the treatment rate of the disease won’t increase until a plateau is reached. Therefore, constant recovery is a poor description of real-world infections. Aiming at this kind of improvement of circumstance, Wang [14] incorporated the following piecewise linear treatment function into an SIR model:

$$T(I)=\left \{ \begin{array}{l@{\quad}l} rI, & 0\leq I\leq I_0,\\ rI_0, & I>I_0, \end{array} \right . $$

where rI ₀ is the capacity of treatment for infected population, and r is a positive constant. It is found that a backward bifurcation occurs if the capacity is small, and there exist bistable endemic equilibria if the capacity is low. In addition, Zhang and Liu [15] studied an SIR model with saturated treatment function:

$$T(I)=\frac{\beta I}{1+\alpha I}, $$

where α,β>0. The authors found sufficient conditions for global stability of endemic and disease-free equilibria, and they discussed the equilibrium effects of the magnitudes of the parameters in the treatment function.

However, in some extent, it’s possible that a society’s capacity to provide treatment for a severe epidemic will decrease. The reasons lead to the above mentioned condition are that equipments and supplies are run out or that healthcare assistants are subjected to victim to the disease. According to this, Eckalbar and Eckalbar [16] studied an SIR epidemic model with quadratic treatment. The treatment function is assumed to be of the form

$$T(I)=\max \bigl[rI-gI^2,0 \bigr], $$

where r is regarded as the treatment rate which is described the societal effort to fight the infection. It means that the treatment increases at a decreasing rate as the sub-population of infected rises. And at some finite number of infected individuals, society’s capacity for treating the infected reaches a peak and then begins to fall in that supplies are more and more diminishing or efficiency of health care resources is low. One of the important assumptions in this function is that g is quite small, perhaps on the order of r/N, where N is the population size. If g were precisely r/N, then T(I)=rI(1−I/N), and T(I)=0 for I=0 and for I=N. T(I) has the maximum value when I=N/2. The system is found to have as many as four equilibria, with possible bi-stability, backward bifurcations, and limit cycles.

We know that some diseases have a latent period, in which individuals are exposed to a disease but not yet infected. Motivated by the above work, in this paper we will consider an SEIR epidemic model with quadratic treatment. The model equations are given as follows:

$$ \left \{ \begin{array}{l} \displaystyle \frac{dS}{dt}=A-\beta SI-\mu S,\\\noalign{\vspace*{3pt}} \displaystyle \frac{dE}{dt}=\beta SI-\mu E-\delta E,\\\noalign{\vspace*{3pt}} \displaystyle \frac{dI}{dt}=\delta E-(\mu+\varepsilon+\eta)I-\max \bigl[rI-gI^2,0 \bigr],\\\noalign{\vspace*{3pt}} \displaystyle \frac{dR}{dt}=\eta I+\max \bigl[rI-gI^2,0 \bigr]-\mu R, \end{array} \right . $$

(1.1)

where A is the recruitment rate, μ is the natural death rate, β is the infection force, δ is the rate at which exposed individuals become infectious, η is the recovery rate for the infectious class, ε is the disease-related death rate, parameters r and g are as introduced as the above discussion. All parameters are assumed to be positive. If we define Ω={(S,E,I,R)|S≥0,E≥0,I≥0,R≥0}, it is obvious from inspection of the above system that if the initial value (S(0),E(0),I(0),R(0)) is in Ω, then (S(t),E(t),I(t),R(t)) will be in Ω at all t>0.

Note that the system has four dynamic variables, S, E, I, and R. We can observe that the variable R is described by differential equation dR/dt=ηI+max[rI−gI ²,0]−μR, but the variable R does not appear in the first three equations of the system (1.1). Hence we only need to consider the subsystem of (1.1) as follows:

$$ \left \{ \begin{array}{l} \displaystyle \frac{dS}{dt}=A-\beta SI-\mu S,\\ \noalign{\vspace*{3pt}} \displaystyle \frac{dE}{dt}=\beta SI-\mu E-\delta E,\\ \noalign{\vspace*{3pt}} \displaystyle \frac{dI}{dt}=\delta E-(\mu+\varepsilon+\eta)I-\max \bigl[rI-gI^2,0 \bigr]. \end{array} \right . $$

(1.2)

System (1.2) can be decomposed into two sub-systems according to whether max[rI−gI ²,0]=0 or rI−gI ². To this end, we define two domains: D ₁={(S,I,E)|rI−gI ²≥0,S,I,E≥0} and D ₂={(S,I,E)|rI−gI ²<0,S,I,E≥0}, and explore the two associated sub-systems

$$ \left \{ \begin{array}{l} \displaystyle \frac{dS}{dt}=A-\beta SI-\mu S,\\\noalign{\vspace*{3pt}} \displaystyle \frac{dE}{dt}=\beta SI-\mu E-\delta E,\\\noalign{\vspace*{3pt}} \vspace{1mm} \displaystyle \frac{dI}{dt}=\delta E-(\mu+\varepsilon+\eta)I- \bigl(rI-gI^2 \bigr), \end{array} \right . \quad (S, I, E)\in D_1, $$

(1.3)

and

$$ \left \{ \begin{array}{l} \displaystyle \frac{dS}{dt}=A-\beta SI-\mu S,\\\noalign{\vspace*{3pt}} \displaystyle \frac{dE}{dt}=\beta SI-\mu E-\delta E,\\\noalign{\vspace*{3pt}} \displaystyle \frac{dI}{dt}=\delta E-(\mu+\varepsilon+\eta)I, \end{array} \right . \quad (S, I, E)\in D_2. $$

(1.4)

The purpose of this paper is to show that quadratic treatment has significant impact on an SEIR epidemic model. And the remaining parts of this paper are organized as follows: in the next section, we investigate the equilibria of system (1.2). It is easy to check that the system has four possible equilibria, and the number of real and positive endemic equilibria depend upon the parameters chosen. We mainly discuss treatment parameter r which is a metric for the societal effort to fight the infection. Then, we discuss the stability of equilibria of the model in Sect. 3. We end the paper with a brief discussion in Sect. 4.

2 Equilibrium states

In the section, we want to investigate the equilibria of system (1.2), which can be obtained by setting the right side of each of the three differential equations equal to zero. Thus, any equilibrium P ^∗=(S ^∗,E ^∗,I ^∗) must satisfy the following algebraic system

$$ \left \{ \begin{array}{l} A-\beta SI-\mu S=0,\\ \beta SI-\mu E-\delta E=0,\\ \delta E-(\mu+\varepsilon+\eta)I-\max \bigl[rI-gI^2,0 \bigr]=0. \end{array} \right . $$

(2.1)

From this, it is easy to check that the system (1.2) have four possible equilibrium vectors. The possible equilibria are:

$$P_0^*=\bigl(S_0^*, E_0^*, I_0^* \bigr) \quad\mbox{iff} \quad\bigl(S_0^*, E_0^*, I_0^*\bigr)\in D_1, $$

$$P_1^*=\bigl(S_1^*, E_1^*, I_1^* \bigr) \quad\mbox{iff} \quad\bigl(S_1^*, E_1^*, I_1^*\bigr)\in D_1, $$

$$P_2^*=\bigl(S_2^*, E_2^*, I_2^* \bigr)\quad\mbox{iff} \quad\bigl(S_2^*, E_2^*, I_2^*\bigr)\in D_1, $$

$$P_3^*=\bigl(S_3^*, E_3^*, I_3^* \bigr)\quad\mbox{iff} \quad\bigl(S_2^*, E_2^*, I_2^*\bigr)\in D_2, $$

where

$$S_0^*=A/\mu, \qquad I_0^*=0, $$

$$S_1^*=\frac{(\mu+\delta)(\mu g+\beta\mu+\beta\varepsilon+\beta\eta+\beta r)+\sqrt{\varDelta}}{2\delta\beta^2}, $$

$$I_1^*=\frac{(\mu+\delta)(-\mu g+\beta\mu+\beta\varepsilon+\beta\eta+\beta r)-\sqrt{\varDelta}}{2g\beta(\mu+\delta)}, $$

$$S_2^*=\frac{(\mu+\delta)(\mu g+\beta\mu+\beta\varepsilon+\beta\eta+\beta r)-\sqrt{\varDelta}}{2\delta\beta^2}, $$

$$I_2^*=\frac{(\mu+\delta)(-\mu g+\beta\mu+\beta\varepsilon+\beta\eta+\beta r)+\sqrt{\varDelta}}{2g\beta(\mu+\delta)}, $$

$$S_3^*=\frac{(\mu+\delta)(\mu+\varepsilon+\eta)}{\delta\beta }, \qquad I_3^*= \frac{\delta\beta A-\mu(\mu+\delta)(\mu+\varepsilon+\eta)}{\beta(\mu+\delta)(\mu +\varepsilon+\eta)}, $$

and

$$E_i^*=\beta S_i^*I_i^*/(\mu+ \delta) \quad (i=1,2,3), $$

where

$$\varDelta=-4A g\delta\beta^2(\mu+\delta)+(\mu+\delta)^2 \bigl[\beta(\varepsilon+\eta+r)+\mu (g+\beta)\bigr]^2. $$

$P_{0}^{*}$ is the disease-free equilibrium, and $P_{1}^{*}$, $P_{2}^{*}$, and $P_{3}^{*}$ are the possible endemic equilibria. It follows from the second equation of (2.1) that E ^∗=βS ^∗ I ^∗/(μ+δ). It is clear that if S ^∗ and I ^∗ are fixed at an equilibrium, then E ^∗ is fixed. Thus, it’s not necessary to calculate $E_{i}^{*}$ to study the equilibria of system (1.2). Next, we only consider the following algebraic equations

$$ \left \{ \begin{array}{l} A-\beta SI-\mu S=0,\\ \delta\beta SI-(\mu+\delta)(\mu+\varepsilon+\eta)I-(\mu+\delta)\max\bigl[rI-gI^2,0\bigr]=0. \end{array} \right . $$

(2.2)

As described above we also define two domains: $D_{1}'=\{(S, I)|rI-gI^{2}\geq0, S, I\geq\nobreak 0\}$ and $D_{2}'=\{(S ,I)|rI-gI^{2}<0, S, I\geq 0\}$. The curves A−βSI−μS=0 and δβSI−(μ+δ)(μ+ε+η)I−(μ+δ)max[rI−gI ²,0]=0 meet or intersect, and denote that $G_{0}^{*}=(S_{0}^{*}, I_{0}^{*})\in D_{1}'$, $G_{1}^{*}=(S_{1}^{*}, I_{1}^{*})\in D_{1}'$, $G_{2}^{*}=(S_{2}^{*}, I_{2}^{*})\in D_{1}'$ and $G_{3}^{*}=(S_{3}^{*}, I_{3}^{*})\in D_{2}'$ are points where the two curves meet. As discussed in [16] that a real and weakly positive $G_{0}^{*}$ will always exist in $D_{1}'$, but the number of real and positive endemic equilibria depends upon the parameters chosen. Sometimes all three $G_{1}^{*}$, $G_{2}^{*}$, and $G_{3}^{*}$ will exist in real and positive forms, and sometimes none, one, or two of them will exist. Let iso-I={(S,I)|A−βSI−μS=0} and iso-S={(S,I)|δβSI−(μ+δ)(μ+ε+η)I−(μ+δ)max[rI−gI ²,0]=0}. Then, equilibria are found at points where iso-S meets iso-I. In particular, we give Fig. 1, which indicates four equilibria exist together.

Fig. 1

Iso-S and iso-I curves showing four equilibria with: {A=561.9677, β=0.1352458, μ=0.57253, δ=0.026936157, ε=0.14288351, η=0.5, r=2.9, g=0.24147661}. I=r/g divides the domain into $D_{1}'$ and $D_{2}'$

We remark that S=A/(βI+μ), E=βSI/(μ+δ) at any equilibrium. Therefore, if we verify that I is real and positive, so are S and E. In view of the existence of the fixed points, we have the following theorem which follows immediately from the results of Theorem 4.1 in [16].

Theorem 2.1

For system (1.2), the real and positive values for $P_{i}^{*}$ will exist in their appropriate domains under the following conditions:

• (1) $P_{0}^{*}$ is real and weakly positive, and $P_{0}^{*}$ will lie in D ₁ at all parameter values.

• (2) Equilibrium $P_{1}^{*}$ will be real and positive, if and only if

  • (2a) $\frac{\sqrt{A\delta}}{\sqrt{g(\mu+\delta)}}-\frac{\mu }{\beta}>0$ and

  • (2b) r ₀<r<r ₁, where

    $$r_0= {\small \begin{aligned} \frac{\sqrt{4Ag\beta^2\delta(\mu+\delta)}-\mu g(\mu+\delta)-\mu \beta(\mu+\delta)-\varepsilon\beta(\mu+\delta)-\eta\beta(\mu+\delta )}{\beta(\mu+\delta)} \end{aligned}} $$

    and

    $$r_1=\frac{-(\mu+\delta)(\mu^2+\mu\varepsilon+\mu\eta)+A\delta\beta}{\mu (\mu+\delta)}. $$

    (Note: $r_{1}-r_{0}=\frac{(\beta\sqrt{A\delta}-\mu\sqrt{g(\mu +\delta)})^{2}}{\mu\beta(\mu+\delta)}$, so r ₁>r ₀.)

  • (2c) Assuming (2a) and (2b) hold, real and positive $P_{1}^{*}$ will exist in D ₁ if any of the following hold:

    • (i) r ₀>0 and at r ₀, $P_{1}^{*}\in D_{1}$, or

    • (ii) r ₀>0 and at r ₀, $P_{1}^{*}\in D_{2}$, but r is restricted to r>r ₂, where $r_{2}=\frac{g(\delta A\beta-\mu(\mu+\delta)(\mu+\varepsilon+\eta))}{\beta(\mu+\delta)(\mu +\varepsilon+\eta)}$, or

    • (iii) r ₀<0 but $I_{1}^{*}$ is positive at r=0 and r is restricted to r>r ₂.

• (3) Equilibria $P_{2}^{*}$ will be real positive, and in D ₁ if and only if r>r ₀ and

  • (3a) r ₀≥0, and

    • (i) if $I_{2}^{*}$ at r ₀ is positive and less than r/g. And r∈[r ₀,r ₂], or

    • (ii) if $I_{2}^{*}$ at r ₀ is negative and r ₁<r<r ₂. Or

  • (3b) r ₀<0, $I_{2}^{*}$ at r ₀ is negative and r ₁<r<r ₂.

• (4) Real and positive values for $P_{3}^{*}$ will only obtain when r ₂>0 and r<r ₂.

From (2c)(iii) and (3b), we know that when r ₀<0, a real and positive $P_{1}^{*}$ or a real and positive $P_{2}^{*}$ can be obtained, but we cannot have both, and we will have neither if

$$(\mu+\delta)(\beta\mu+\beta\varepsilon+\beta\eta)+\phi>(\mu+\delta)\mu g>(\mu+\delta)(\beta\mu+\beta\varepsilon+\beta\eta)-\phi, $$

where $\phi=\sqrt{-4A g\delta\beta^{2}(\mu+\delta)+(\mu+\delta)^{2}[\beta(\varepsilon+\eta)+\mu (g+\beta)]^{2}}$.

Also notice that both $P_{2}^{*}$ and $P_{3}^{*}$ drop out of existence when r rises above r ₂. What is happening is that each solution is crossing out of its domain as r crosses from below to above r ₂.

From the above result, it is clear that changes of the parameter r not only lead to changes of the values of the possible endemic equilibria, but also affect the existence of them.

3 Equilibria stability

In this section, we discuss the locally asymptotical stability of the equilibria and global issues of disease-free equilibrium. Note that $P_{0}^{*}, P_{1}^{*}, P_{2}^{*}\in D_{1}$, and $P_{3}^{*}\in D_{2}$. Firstly, we shall consider the local asymptotic stability of equilibria which lie in D ₁.

The stability examination begins with the Jacobian matrix for sub-system (1.3):

$$ J(S,E,I) =\left ( \begin{array}{c@{\quad}c@{\quad}c} -\beta I-\mu& 0 & -\beta S \\ \beta I & -(\mu+\delta) & \beta S \\ 0 & \delta&-(\mu+\varepsilon+\eta+r)+2gI \end{array} \right ). $$

(3.1)

The characteristic equation of the Jacobian matrix of (1.2) corresponding to (S,E,I) is

$$\lambda^3+a_1\lambda^2+a_2 \lambda+a_3=0 $$

where

$$ \left \{ \begin{array}{l} a_1=(\mu+\delta)+(\mu+\varepsilon+\eta+r-2gI)+(\beta I+\mu), \\[2pt] a_2=(\mu+\delta)(\mu+\varepsilon+\eta+r-2gI) \\[2pt] \hphantom{a_2=}{} -\beta\delta S+(\beta I+\mu)\bigl[(\mu+\delta)+(\mu+\varepsilon+\eta+r-2gI)\bigr], \\[2pt] a_3=(\beta I+\mu)(\mu+\delta)(\mu+\varepsilon+\eta+r-2gI)-\beta\delta S\mu. \end{array} \right . $$

(3.2)

By the Routh-Hurwitz criterion, we know that (S,E,I) is locally asymptotically stable which requires a ₁>0, a ₁ a ₂−a ₃>0, a ₃>0.

Theorem 3.1

The disease-free equilibrium, $P_{0}^{*}=(S_{0}^{*}, E_{0}^{*}, I_{0}^{*})$ is locally asymptotically stable if and only if

$$r>r_1=\frac{-(\mu+\delta)(\mu^2+\mu\varepsilon+\mu\eta )+A\delta\beta}{\mu(\mu+\delta)}. $$

Proof

The characteristic equation of the Jacobian matrix of (1.2) corresponding to $P_{0}^{*}$ is

The root of the first factor is negative real −μ. It is easy to see the roots of λ ²+((μ+δ)+(μ+ε+η+r))λ+(μ+δ)(μ+ε+η+r)−βAδ/μ=0 are either negative real or complex conjugates having negative real part when (μ+δ)(μ+ε+η+r)−βAδ/μ>0, i.e., $r>\frac{-(\mu+\delta)(\mu^{2}+\mu\varepsilon+\mu\eta )+A\delta\beta}{\mu(\mu+\delta)}=r_{1}$. Therefore, $P_{0}^{*}$ is locally asymptotically stable if and only if r>r ₁. □

Theorem 3.2

Assuming a real and positive $P_{1}^{*}$ exists and β>g, then $P_{1}^{*}$ will be locally asymptotically stable.

Proof

For $P_{1}^{*}$, it follows from (3.2) that

and

and

where

where ψ=βμ+βε+βη+βr.

Since $(\mu+\delta)(\mu g+\beta\mu+\beta\varepsilon+\beta\eta+\beta r)\geq\sqrt{\varDelta}$, so

Thus, we have

therefore by Routh-Hurwitz criterion, we have proved that when β>g then $P_{1}^{*}$ is locally asymptotically stable. □

Similarly, by using the Routh-Hurwitz criterion, we can also obtain the conditions of the local stability of the equilibria $P_{2}^{*}$ and $P_{3}^{*}$.

Theorem 3.3

If a real and positive $P_{2}^{*}$ exists, then $P_{2}^{*}$ is unstable.

Theorem 3.4

If a real and positive $P_{3}^{*}$ equilibrium exists in D ₂, then $P_{3}^{*}$ is locally asymptotically stable at all possible positive parameter values.

In the following, we shall study global issues of disease-free equilibrium. Here we define Σ=S(t)+I(t)+E(t) and Φ={(S,E,I)|S≥0,E≥0,I≥0,S+E+I≤A/μ}, respectively.

Theorem 3.5

For system (1.2), if Σ=A1/μ, where A1≥A, then dΣ(t)≤0.

Proof

According to system (1.2), we have

$$\frac{dS}{dt}+\frac{dE}{dt}+\frac{dI}{dt}=A-\mu S-\mu E-\mu I-( \varepsilon+\eta)I-\max\bigl[rI-gI^2,0\bigr]. $$

We note that

$$\frac{d\varSigma}{dt}=\frac{dS}{dt}+\frac{dE}{dt}+ \frac{dI}{dt}, $$

and substituting S=A1/μ−E−I into above equation, we obtain

$$\frac{d\varSigma}{dt}=A-A1-(\varepsilon+\eta)I-\max\bigl[rI-gI^2,0\bigr] \leq0, $$

since A≤A1, I≥0, and max[rI−gI ²,0]≥0. □

Several comments on the implications of Theorem 3.5 are similar to Theorem 6.6 in [13]. That is, we obtain that every solution of system (1.2) with the initial condition eventually enters and remains in the region Φ. Also we have the following theorem.

Theorem 3.6

The disease-free equilibrium $P_{0}^{*}$ is globally asymptotically stable if r>max(r ₁,r ₂).

Proof

According to Theorem 3.5, we have that all points in Ω∖Φ will flow ultimately into Φ. Thus, without loss of generality, we assume that all of the positive semi-paths originating on the Φ. Obviously, Φ is positively invariant. And from Theorem 2.1, we know that $P_{0}^{*}$ is a unique equilibrium when r>max(r ₁,r ₂). Moreover, there are no closed paths around $P_{0}^{*}$ since the trajectories cannot cross E=0 and I=0. It follows that all positive semi-paths approach $P_{0}^{*}$ if r>max(r ₁,r ₂). □

4 Numerical simulation and discussion

Treatment is an important method to decrease the spread of diseases. Based on the ideas of Eckalbar and Eckalbar [16], in this paper, an SEIR epidemic model with quadratic treatment has been investigated. The key point of this paper is to analyze the impact of the changes of parameter r, which we have regarded as a indicator of the extent of the social mobilization for treatment. It has been shown that the system has four possible equilibria under some certain conditions. By analyzing the characteristic equation of the Jacobian matrix of (1.2) corresponding to equilibria, the conditions of stability are obtained. Finally, we have proved that the disease-free equilibrium of the system is globally asymptotically stable once r passes a threshold value r=max(r ₁,r ₂).

Our focus so far has been on the existence and asymptotic stability of the equilibria of the system (1.2). Now we will study the effects of the treatment rate r on the complexity of the system (1.2). There are several dynamic possibilities of system (1.2) as the value of r changes. Let A=9.5040669, β=0.0010958995, μ=0.027727999, δ=0.10034215, ε=0.1429976, η=0.001, and g=0.0060806619. Numerical analysis is done that can track the stability of equilibria of the equations, and the bifurcation diagram is plotted in Fig. 2. From Fig. 2, it can be observed that: (i) as r∈(0,r ₀) (r ₀=0.10000076), there are two equilibria $P_{0}^{*}$ and $P_{3}^{*}$. $P_{0}^{*}$ is unstable and $P_{3}^{*}$ is locally stable. (ii) As r increases and exceeds the critical value r ₀=0.10000076, the other two endemic equilibria $P_{1}^{*}$, $P_{2}^{*}$ become exist, and the endemic state $P_{1}^{*}$ is locally stable and $P_{2}^{*}$ is unstable. It follows that, bistable phenomenon can occur when r is between r ₀ and r ₂ (see Fig. 3). (iii) As r increases and exceeds the critical value r ₂=0.1098196, the two endemic equilibria $P_{2}^{*}$, $P_{3}^{*}$ vanish, and $P_{0}^{*}$ is unstable and $P_{1}^{*}$ is locally stable. (iv) As r increases and exceeds the other critical value r ₁=0.122579, the equilibrium $P_{1}^{*}$ vanishes and the stability of $P_{0}^{*}$ will be changed. $P_{0}^{*}$ becomes locally stable.

Fig. 2

The bifurcation diagram of population I with respect to r on [0,0.2]. I=r/g divides the domain into $D_{1}'$ and $D_{2}'$

Fig. 3

Bi-stability of the equilibria $P_{1}^{*}$ and $P_{3}^{*}$ with r=0.105, r ₀<r<r ₂

Note that as r=0.116 (r ₂<r<r ₁, i.e., the case (iii)), there are only two equilibria $P_{0}^{*}$ and $P_{1}^{*}$, and $P_{0}^{*}$ is unstable and $P_{1}^{*}$ is locally stable. By numerical simulation, Fig. 4 shows that all positive trajectories tend to $P_{1}^{*}$. So we guess the $P_{1}^{*}$ is globally asymptotically stable if r ₂<r<r ₁. Theoretical proofs will be considered in the future research.

Fig. 4

Positive trajectories from all initial points other than $P_{0}^{*}$ and its stable manifolds lead to $P_{1}^{*}$ with r=0.116 (r ₂<r<r ₁)