Mathematical Modeling and Analysis of Soft ComputingView this Special Issue
Stability and Bogdanov-Takens Bifurcation of an SIS Epidemic Model with Saturated Treatment Function
This paper introduces the global dynamics of an SIS model with bilinear incidence rate and saturated treatment function. The treatment function is a continuous and differential function which shows the effect of delayed treatment when the rate of treatment is lower and the number of infected individuals is getting larger. Sufficient conditions for the existence and global asymptotic stability of the disease-free and endemic equilibria are given in this paper. The first Lyapunov coefficient is computed to determine various types of Hopf bifurcation, such as subcritical or supercritical. By some complex algebra, the Bogdanov-Takens normal form and the three types of bifurcation curves are derived. Finally, mathematical analysis and numerical simulations are given to support our theoretical results.
In the mathematical modeling of epidemic transmission, there are several factors that substantially affect the dynamical behavior of the models. Incidence rate functions are seen as a major factor in producing the rich dynamics of epidemic models in many literatures (see [1–25]). In most classical models of epidemics, the incidence rate is taken to be mass action incidence with bilinear interaction, that is, , where is the probability of transmission per contact and and represent the number of susceptible and infected individuals, respectively. Epidemic models with such bilinear incidence rates usually have at most one endemic equilibrium, and then the diseases will be eradicated if the basic reproduction number is less than one and will persist otherwise. Besides, there are also many other types of incidence rate functions, such as nonlinear incidence rate, standard incidence rate, and saturated incidence rate. Recently, there are many studies that have demonstrated the nonlinear incidence rate which is one of the key factors that induce periodic oscillations in epidemic models (see [1, 2, 9, 10, 15]). Moreover, Liu et al.  introduced a nonlinear incidence rate of the form where means the infection force of the disease and measures the inhibition effect from the behavioral change of the susceptible individuals when the number of infective individuals increases. So we can see that the bilinear incidence rate is a special case of (1) with and or . Furthermore, Wang and Ruan in  studied the global dynamics of an SIRS model with the incidence function ; that is, ; they also showed that the SIRS epidemic model undergoes a Bogdanov-Takens bifurcation, that is, saddle-node, Hopf, and homoclinic bifurcations. To have a better understanding of the dynamics of the system, Tang et al.  calculated higher order Lyapunov values of the weak focus and reduced the system to a universal unfolding form for a cusp of codimension 3, and they gave the bifurcation surfaces and displayed all limit cycles and monoclinic loops of order up to 2. Wei and Cui explored an SIS epidemic model with standard incidence rate in . They took the incidence rate form and showed the dynamics and backward bifurcation of the SIS epidemic model.
Recently, in order to prevent and control the spread of the infectious diseases such as measles, tuberculosis, and flu, many mathematicians (see [7, 11, 16, 20, 23, 25–29]) have begun to investigate the role of treatment functions in epidemiological models. In some classical epidemic models, the treatment function is an important method to decrease the spread of the epidemiological diseases. Generally speaking, the treatment function of the infective individuals is assumed to be proportional to the number of the infective individuals. But every community should have a maximal capacity for the treatment of a disease and the resources for treatment should be very large. Therefore, it is very important to adopt a suitable treatment function. In , Wang and Ruan introduced a constant treatment function of diseases in an SIR model; that is, This means that they use the maximal treatment capacity to cure infective individuals so that the disease can be eradicated. They also found that the model undergoes saddle-node bifurcation, Hopf bifurcation, and Bogdanov-Takens bifurcation. Further, a piecewise linear treatment function was considered in ; that is, where and and are positive constants. This means that the treatment rate is proportional to the number of the infective individuals when the capacity of treatment has not been reached; otherwise it takes the maximal capacity of treatment . By considering the above treatment function, Wang  found that a backward bifurcation takes place in an SIR epidemic model. In , J. C. Eckalbar and W. L. Eckalbar constructed an SIR epidemic model with a quadratic treatment function; that is, , . They found that the system has as many as four equilibria, with possible bistability, backward bifurcations, and limit cycles.
Recently, saturated treatment function has been widely applied in many epidemic models. In particular, in paper , Zhang and Liu took a continuous and differentiable saturated treatment function , where , . stands for the cure rate and measures the extent of the effect of the infected being delayed for treatment. We can see that the treatment function when is small enough, whereas when is large enough. It is more realistic and it has the convenience of being continuous and differential compared to the previous ones. Furthermore, the authors in  found that is a critical threshold for disease eradication when this delayed effect for treatment is weak and a backward bifurcation will take place when this effect is strong. So, it is really important to adequately stress the interesting connection recently established between the choice of saturated treatment functions in epidemic models and the occurrence of backward bifurcation in the related system dynamics. In fact, recently, saturated-type treatment functions have been indicated as responsible for the occurrence of backward bifurcations for SIR [20, 23], for SIS [19, 21], and for SEIR [28, 29] models, supporting the general circumstance that saturated-type treatments can be one of the causes of backward bifurcations in epidemic models. In particular, in , such connection has also been shown and validated in a specific concrete disease-control setting.
In the real world, some infectious diseases do not confer immunity. Such infections do not have a recovered state and individuals become susceptible again after infection. This type of disease can be modeled by the SIS type. So the SIS epidemic model has been adopted by many mathematicians (see [8, 11, 18, 19, 21, 25]). And SIS models are appropriate for some bacterial agent diseases such as meningitis, plague, and venereal diseases and for protozoan agent diseases such as malaria and sleeping sickness.
Motivated by the above points, we will consider the following SIS model with bilinear incidence rate and saturated treatment function: where and denote the numbers of susceptible and infective individuals, respectively. Positive constant is the recruitment rate of the population. Positive constant is the nature death rate of population. The bilinear incidence rate is , where is positive. Positive constant is the natural recovery rate of infective individuals. Positive constant is the disease-related death rate. The saturated treatment function , where is positive and is nonnegative.
This paper focuses on the detailed dynamics analysis of the model (4). The local stability of these equilibria is investigated, which enables us to classify the types of model equilibria (e.g., attractor, saddle, or repeller). We show that the system has backward bifurcation and Bogdanov-Takens bifurcation (i.e., Hopf bifurcation, saddle-node bifurcation, and homoclinic bifurcation) under some certain conditions. Finally, the three bifurcation curves and the complicated global bifurcation phase portraits are derived by applying the Bogdanov-Takens normal form and the corresponding parameters which satisfy the conditions that ensure Bogdanov-Takens bifurcation exists.
The organization of this paper is as follows. In Section 2, we study the existence and local stability of equilibria and backward bifurcation. In Section 3, we investigate the global stability of the model. In Section 4, we give the supercritical and subcritical bifurcation under two different conditions in system (4). In Section 5, we show that the system (4) undergoes Bogdanov-Takens bifurcation under some certain conditions. In Section 6, some numerical simulations are displayed in detail. We close with a discussion in Section 7 on our mathematical results and epidemiological implications.
2. Equilibria and Backward Bifurcation
2.1. Disease-Free Equilibrium
Obviously, system (4) has a disease-free equilibrium . The Jacobian matrix of (4) at is By using the next generation matrix of (4), we get the basic reproduction number . has negative eigenvalues if . Then we have the following result.
2.2. Endemic Equilibria
An endemic equilibrium always satisfies In view of , we get and substitute it into the second equation of (6). When , we obtain Then we have an equation of the form with This equation may admit positive solution
Obviously, if , then , if , then , and if , then . From (8), it is obvious that we have the following results.
Theorem 2. The following results hold. Let . Equation (8) is a linear equation with a unique solution . Then the system (4) has a unique endemic equilibrium when and has no endemic equilibrium when .Let . If , system (4) has a unique endemic equilibrium when and no endemic equilibrium when .Let . If , system (4) has a unique endemic equilibrium when , no endemic equilibrium when , and two endemic equilibria and when . When and , one has which is equivalent to
We define the right-hand side of (11) as
Therefore, according to the qualitative approach recently proposed by  which is based on the analysis of the equilibria curve in the neighboring of the critical threshold , we have the following theorem.
In order to verify the bifurcation curve (the graph of as a function of ) in Figure 3, we think of as a variable with the other parameters as constant. Through implicit differentiation of (8) with respect to , we get From (13), we know that the sign of is opposite to that of . And from the definition of we know that decreases when increases. It implies that the bifurcation curve has positive slope at equilibrium values with and negative slope at equilibrium values with . If there is no backward bifurcation at , then the unique endemic equilibrium for satisfies and the bifurcation curve has positive slope at all points where . If there is a backward bifurcation at , then there is an interval on which there are two endemic equilibria given by The bifurcation curve has negative slope at the smaller one and positive slope at the larger one. Thus the bifurcation curve is as shown in Figure 3.
Under the conditions of Theorem 3, if a backward bifurcation takes place, we can see from Figure 3 that there is a critical value at the turning point. In this case, the disease will not die out when . However, the disease will die out when . Therefore, the critical value can be taken as a new threshold for the control of the disease.
In the following, we give an explicit criterion of a backward bifurcation at .
For convenience, we define
Corollary 4. When , system (4) has a backward bifurcation at .
Proof. When , The condition is equivalent to From (17) and (18), we get , which reduces to It means that is big enough to lead a backward bifurcation with two endemic equilibria when . Therefore, the proof is complete.
Next, we consider the local stability of the unique endemic equilibrium when .
Theorem 5. When and , the unique endemic equilibrium is locally asymptotically stable.
Proof. Firstly, from Theorem 2, we can know that system (4) has a unique endemic equilibrium when . Moreover, the Jacobian matrix of system (4) is From the second equation of (6), we have From (21), the Jacobian matrix reduces to We obtain In fact, there is . Since and , we have So we get . The trace of is given by In the same way as the above calculation of (24), we have So we get . The proof is complete.
Now we consider the case that there are two endemic equilibria and ; let be the Jacobian matrix at , .
Theorem 6. The endemic equilibrium is a saddle whenever it exists.
Proof. Since and , we have . Thus From of Theorem 2, we can know that if , , and , then exists. Hence, we can get and . It follows that there exists a unique such that where On the other hand, where In the following we will show that . Since , we have that is, Therefore, Obviously, we have . From (30), one has . So we get . Hence the endemic equilibrium is a saddle. The proof is complete.
In order to explore the stability of the endemic equilibrium , define
Theorem 7. If , then endemic equilibrium is locally asymptotically stable; if , then endemic equilibrium is unstable, where .
Proof. Consider By carrying out arguments similar to that of Theorem 6, we have . Therefore, . In addition, we have and then , where Using the expression of and , one has where is a first degree polynomial of . Since , . From the expression of , we have Thus, is locally asymptotically stable if and is unstable if . The proof is complete.
From the above discussion, we have the following conclusion. If two endemic equilibria and exist, the stable manifolds of the saddle split into two regions. The disease is persistent in the upper region and dies out in the lower region (see Figure 4).
3. Global Analysis
Firstly, we consider the global stability of the disease-free equilibrium . Let be the total population size. Now we note that the equation for total population is given by . It follows that . Let which is positive invariant with respect to system (4).
Theorem 8. If , the disease-free equilibrium is globally asymptotically stable; that is, the disease dies out.
Proof. Suppose . From the of Theorem 2, we know that the model has no endemic equilibrium. From the corollary of Poincaré-Bendixson theorem , we know that there is no periodic orbit in as there is a disease-free equilibrium in . Since is a bounded positively invariant region and is the only equilibrium in , the local stability of implies that every solution initiating in approaches . Thus, the disease-free equilibrium is globally asymptotically stable. The proof is complete.
Now we analyze the global dynamics of the endemic equilibrium when .
Theorem 9. If and , the system (4) has no limit cycles.
Proof. We use Dulac theorem to exclude the limit cycle. Let and take the Dulac function According to , we can get Hence, the system (4) has no limit cycles. The proof is complete.
Therefore, we obtain the global result of the unique endemic equilibrium.
Theorem 10. If and , the unique endemic equilibrium is globally asymptotically stable (see Figure 5).
4. Hopf Bifurcation
In this section, we study the Hopf bifurcation of system (4). From the above discussion, we know that there is no closed orbit surrounding or because the -axis is invariant with respect to system (4) and is always a saddle. Therefore, Hopf bifurcation can only occur at . Set where , , and are defined by (50) and (52).
Theorem 11. System (4) undergoes a Hopf bifurcation if . Moreover, if , there is a family of stable period orbits of system (4) as decreases from 0; that is, a supercritical Hopf bifurcation occurs; if , there is a family of unstable period orbits of system (4) as increases from 0; that is, a subcritical Hopf bifurcation occurs.
Proof. The proof of Theorem 7 shows that if and only if , and when exists. Therefore, the eigenvalues of are a pair of pure imaginary roots if and only if . The direct calculations show that
By Theorem 3.4.2 in , is the Hopf bifurcation point for system (4).
To be concise in notations, rescale (4) by . For simplicity, we still use instead of . Then we obtain Let and ; then (48) becomes where Let denote the origin of plane. Since satisfies (6), we obtain From the proof of Theorem 7, it follows that is always positive. It is easy to verify that if and only if . Set and let and ; then the normal form of (48) for Hopf bifurcation reads where Now, we evaluate the first Lyapunov coefficient of system (4) as follows: where denotes and so forth. Then Obviously, the sign of is determined by . By Theorem and () in , the rest of the claims in Theorem 11 are valid. The proof is complete.
5. Bogdanov-Takens Bifurcations
The purpose of this section is to study the Bogdanov-Takens bifurcation of (4). Now, we assume the following:,.Then (6) admits a unique positive equilibrium , where The Jacobian matrix of system (4) at is By (58), we have because of Furthermore, implies that Thus, and imply that the Jacobian matrix has a zero eigenvalue with multiplicity 2. This suggests that (4) may admit a Bogdanov-Takens singularity.
Proof. Using the transformation of and , system (4) becomes where is a smooth function of at least of the third order and Set , . Then (62) is transformed into where are smooth functions of at least of the third order and In order to obtain the canonical normal form, we perform the transformation of variables by Then, we obtain where are smooth functions of at least of the third order. Note that and . It follows from [33–35] that (4) admits a Bogdanov-Takens bifurcation.
In the following, we will find the versal unfolding in terms of the original parameters in (4). In this way, we will know the approximate saddle-node, Hopf, and homoclinic bifurcation curves. We choose and as bifurcation parameters. Fix , , , , and . Let and , where and are parameters which vary in a small neighborhood of the origin.
Suppose that , , , , , , and satisfy and . Consider the following system: By the transformations of and , system (68) becomes where is a smooth function of at least of the third order and Making the change of variables , and rewriting , as and , respectively, we have where , is a smooth function of , , , and at least of the third order, and Next, introduce a new time variable by . Rewriting as , we obtain Let , and rename and as and ; we have where , is a smooth function of , , , and at least of the third order, and Now, we assume that and when are small. Set and rewrite as ; we can get where , is a smooth function of , , , and at least of the third order, and Making the final change of variables by , , and and then denoting them again by , , and , respectively, we obtain where and is a smooth function of , , , and at least of the third order.
We substitute values of , , , , , , and for the above system and these values satisfy conditions and . And we obtain the following equations: where and in a small neighborhood of . Let And after simple calculation we obtain that Thus, and are regular maps in a small neighborhood of . By the Bogdanov and Takens bifurcation theorems , we obtain the following conclusion.
Theorem 13. Suppose that , , , , , , and satisfy , , , and when are small. Then (4) admits the following bifurcation behavior. (1)There is a saddle-node bifurcation curve , −, .(2)There is a Hopf bifurcation curve , .(3)There is a homoclinic bifurcation curve , −−−.
6. Numerical Simulations
When , , , , , , and , . By applying PPLANE8 and Photoshop software, the (, )-plane near the origin is divided into 4 regions by these bifurcation curves as shown in Figure 6. Fix and decrease from 0; our conclusions are summarized as follows.(a)When (, ) lies in region I, there is no positive equilibrium, which implies that the positive orbits of (4) meet the positive -axis in finite time, and therefore the disease disappears.(b)There is a saddle-node point when (, ) lies on curve .(c)When (, ) crosses into region II, two positive equilibria which are an unstable focus and a saddle appear.(d)When the parameters lie on the curve H, there is also a saddle and an unstable focus. An unstable limit cycle appears with the parameters crossing H into III.(e)A homoclinic cycle appears as the parameters passing III into HL through the homoclinic bifurcation.(f)The homoclinic loop breaks with the parameter crossing HL into IV, and a saddle and a stable focus appear.
Furthermore, using PPLANE8, some numerical simulations of system (68) are depicted in Figures 7–10. When , that is , there is a unique positive equilibrium , which is a cusp of codimension 2 (Figure 7). When