Abstract
There is a wide range of works that have proposed mathematical models to describe the spread of infectious diseases within human populations. Based on such models, researchers can evaluate the effect of applying different strategies for the treatment of diseases. In this article, we generalize previous models by studying an SIR epidemic model with a nonlinear incidence rate, saturated Holling type II treatment rate, and logistic growth. We compute the basic reproduction number and determine conditions for the local stability of equilibria and the existence of backward bifurcation and Hopf bifurcation. We also show that, when the disease transmission rate and treatment parameter are varied, our model undergoes a Bogdanov-Takens bifurcation of codimension 2 or 3. Simulations of the solutions and numerical continuation of equilibria are carried out to generate 2D and 3D bifurcation diagrams, as well as several related phase portraits that illustrate our results. Our work shows that incorporating these factors into epidemic models can lead to very complex dynamics.
1. Introduction
Treatment is an important and effective method that can be used for controlling the spread of infectious diseases within a population. In epidemic models, the treatment function represents the probability that each infected individual has of receiving treatment against the disease at a given time. Many classical models incorporated treatment rates that were assumed to be proportional to the number of infectives. Ideally, if we face a curable infectious disease, all infected people should receive treatment. In practice, however, it is known that communities have limited medical resources for the treatment of diseases, so providing this service puts great pressure on public health infrastructure. Hence, it is more realistic to consider saturated type treatment functions, which tend to a finite limit as the number of infectives increases.
The authors in [1] considered an SIR epidemic model with bilinear incidence rate and piecewise treatment function: This means that the treatment rate is proportional to the number of infectives when the capacity of treatment is not reached; otherwise, it takes the maximal capacity.
In [2], Zhang and Liu studied a model with nonlinear incidence rate and saturated treatment function The authors studied there the local stability of equilibria and showed that backward bifurcation may occur when the number of infectives is high. This model was further explored in [3], where Zhou and Fan modified the treatment function to the following form and determined sufficient conditions for the existence of backward bifurcation, as well as the existence, stability, and direction of Hopf bifurcation.
Li et al. studied in [4] an SIR model with logistic growth rate, bilinear incidence rate and a saturated treatment function of the form . They studied the local stability of the disease-free and endemic equilibria and showed that the system exhibits backward bifurcation, Hopf bifurcation, and Bogdanov-Takens bifurcation of codimension 2.
The study of nonlinear equations in the modeling of biological processes has gained attention in recent years due to the fact that many systems in nature present inherently nonlinear dynamics. Recent research has shown that the use of nonlinear, saturated functions in both epidemic models [2–5] and ecological (predator–prey) models [6, 7] can lead to the appearance of very complex dynamics from the mathematical viewpoint, such as several types of bifurcation, bistability, heteroclinic and homoclinic orbits, and periodic oscillations. Although these functions make the models more difficult to analyze, they have been proved to represent more accurately the processes that we want to describe, so nonlinear models are worthy of being studied in greater detail.
In this paper, we will modify the SIR model introduced in [4] by using the nonlinear saturated incidence rate studied in [2, 3]. The rest of the paper is organized as follows. We begin by presenting our model and establishing the basic properties of solutions in Section 2. Then, we analyze the different cases for the existence of equilibrium points in Section 3. In Section 4, we determine necessary and sufficient conditions for the stability of the disease-free and endemic equilibria, and we show that a backward bifurcation may occur. In Section 5, we prove that our model may also present a Hopf bifurcation. In Section 6, we compute the normal form of Bogdanov-Takens bifurcation of codimension 2 and 3 for our model. In Section 7, we carry out numerical simulations of the solutions of the model and perform numerical continuation of equilibria to generate 2D and 3D bifurcation diagrams, which illustrate the complex dynamics that can arise under certain conditions. Finally, we summarize our results and discuss their implications for epidemic control strategies in Section 8.
2. Description of the Model
Based on the above, we propose here the following model:where the variables , , denote the the number of susceptible, infected, and recovered individuals at time , respectively. The parameter is the intrinsic growth rate of susceptible population, while represents the carrying capacity of population; represents the natural death rate, is the disease-induced death rate, and is the natural recovery rate of infected population. The incidence rate is of saturated type: is the maximal disease transmission rate, while the factor models the “psychological” or inhibition effect when the number of infectives is large. The parameter is the maximal treatment rate for each individual per time unit, and is a constant that measures the saturation effect caused by infected population being delayed for treatment ( is called the half-saturation constant). We will assume that is nonnegative, while all other parameters are positive constants.
We can notice that the variable does not appear in the first two equations of system (5). Therefore, we can omit the third equation and focus on the following subsystem of (5):Notice that we have written for simplicity. From now on, we will use the parameter to account for the total removal rate of infectives.
It is easily verified that all solutions of (6) with nonnegative initial conditions remain nonnegative for .
Lemma 1. Let . Then the region is a positively invariant and attracting set for system (6).
Proof. From the first equation of (6), we have Hence, a standard comparison argument yields .
Let . Then Thus, we have for sufficiently large, so all solutions of (6) are ultimately bounded and enter the region .
3. Existence of Equilibria
It can be easily verified that system (6) always has a disease-free equilibrium (DFE), which we will denote by , and it is given by
We can compute the basic reproduction number of the model by means of the next-generation matrix method. Following the procedure described in [8], we obtain
We will now determine some conditions for the existence of positive endemic equilibria. For that, we begin by making the right-hand side of (6) equal to zero. Assuming that , we obtainSolving for in (13) and substituting the resulting expression in (14), we obtain, after simplification, the following cubic equation in :whereFor each positive root of the polynomial , we can make so is a positive equilibrium of (6). Thus, the number of endemic equilibria of the model is the same as the number of positive roots of .
It is clear that the coefficient is always positive. On the other hand, we have when , and when . The discriminant of the cubic polynomial is given byand it has the property that has one real root and two complex roots when , three real roots (at least two of which are equal) when , and three distinct real roots when . Note that and the discriminant introduced by Tartaglia and Cardano have opposite sign [9]. Let and be the critical points of , i.e., the points where its derivative equals zero. By calculation, we obtain Assume first that . Then . If is complex or , then is increasing on , which implies that it has no positive roots. This allows us to make a graphical analysis of said roots, similar to Theorem 0.1 in [5], from which we get the following theorem.
Theorem 2. Let be given by (18) and assume . (1)If , then system (6) has no endemic equilibria.(2)If , then there is a unique endemic equilibrium when and no endemic equilibria otherwise.(3)If , then there are two endemic equilibria when and no endemic equilibria otherwise.
Next, for the case when , we have . We can see that is increasing on whenever is complex or , in which case there are no negative roots. In this case, we also have the following result.
Lemma 3. Assume . If , then .
Proof. Assume and let . Then Recall that and
Based on this, we obtain the following result.
Theorem 4. Assume . Then system (6) has a unique endemic equilibrium.
Proof. Let be given by (18). If , them has one real root and two complex conjugate roots. Using the fact that and , it follows that the real root is positive.
If , then has three real roots, being at least two of them equal and located either at or at . Note that if the three real roots of were positive, this would imply due to the fact that is minus the sum of the roots of . Lemma 3 would imply , and hence and , the roots of , would be real and of opposite sign. The fact that and would emphasize that would be the negative one and a double root of , contradicting that the three roots of are positive. Therefore, at least one root of is negative in this case, and it must be its double root. Hence, the remaining root of must be positive, and is negative.
If , then the three roots of are real and different from each other. The previous discussion shows that , implying that two roots are negative and the remaining root is positive.
Note that in all the previous three cases there is only one positive root. Therefore, system (6) has a unique endemic equilibrium and the theorem is proved.
We will now consider the case . In this case, we obtain , so (15) becomes . Hence, the number of endemic equilibria depends on the number of positive roots of the quadratic equation:We can make a graphical analysis of the roots of by looking at the signs of and . For that, we will first prove the following result.
Lemma 5. Assume . If , then .
Proof. We will proceed by contradiction. Assume , , and . Then , and must satisfy in order to have .
Note that and .
If , then Therfore, being , this implies if we assume , which contradicts our assumption.
Theorem 6. Assume . (1)If , then (6) has a unique endemic equilibrium.(2)If and , then (6) has no endemic equilibria.
Proof. If , then the roots of are real of opposite sign. This proves case (1) of the theorem. If and , then the roots of are real non-positive or complex conjugate, proving case (2). Lemma 5 shows that the case with and never happens.
4. Stability of Equilibria
In this section, we will study the local and global stability of equilibria for system (6).
4.1. Stability of the DFE and Backward Bifurcation
The Jacobian matrix of system (6) evaluated at the disease-free equilibrium is given by so the characteristic equation is Thus, there are two eigenvalues: and . It is clear that is always negative, while when and when . If , we know by Hartman-Grobman’s theorem that the solutions of (6) and its linearization are qualitatively equivalent near , so we can conclude the following.
Theorem 7. The disease-free equilibrium is a stable node when , and it is a saddle when .
According to the above theorem and Theorems 2 and 4, we can see that is locally stable whenever there are no endemic equilibria. Furthermore, we can establish the global stability of as follows.
Theorem 8. If (6) has no endemic equilibria, then the disease-free equilibrium is globally asymptotically stable.
Proof. If there is no endemic equilibrium, then is the only equilibrium in the region , which was proved to be positively invariant for system (6). By the Poincaré–Bendixson theory, every solution trajectory starting in will approach either an equilibrium or a closed orbit contained in . Since is located in the boundary of said region, there cannot exist a closed orbit totally contained in enclosing . This implies that all solutions of the system with initial conditions in must approach as tends to infinity.
In the case when , one of the eigenvalues of vanishes, so we cannot apply Hartman-Grobman’s theorem. In such case, we can resort to the center manifold theory.
Theorem 9. LetIf and , then is a saddle-node; i.e., the neighborhood of consists of a parabolic sector and two hyperbolic sectors. If , the hyperbolic sectors are on the half-plane above the -axis; if , the parabolic sector is on the upper half-plane.
Proof. Assume . Making the change of variables , , we obtain that system (6) is equivalent to which has an equilibrium point at the origin. The Jacobian matrix of this system at is with eigenvalues 0 and , and respective eigenvectors and . Using the change of variables , , we obtain the system in diagonal form:where and . We now use the local center manifold theorem (see [10, 11]) to compute the expansion of the center manifold and substitute in the first equation of (30). This yields the flow in the center manifold, which is given by the following equation: Then the dynamics of (30) near the origin are determined by the quadratic term of the above equation, provided that is nonzero. Since , then , so Thus, the equilibrium at the origin corresponds to a saddle-node whenever . Using and to go back to the original coordinates, we can see that has stable directions on the positive and negative -axes and, when is close to zero, we have for and for , hence the theorem.
If we consider as a bifurcation parameter, we can calculate a critical value such that when and when . The Jacobian matrix of the system evaluated at when has a simple zero eigenvalue, whose corresponding right eigenvector and left eigenvector are given by We can now apply Theorem 4.1 from [12] to calculate the bifurcation constants and . Let , where is the function defined by the right-hand side of (6). Taking into account that and , we have The nonzero second partial derivatives of evaluated at are so we obtain and . Since is always positive, the type of bifurcation that occurs at when the basic reproduction number crosses unity depends only on the sign of . Thus, from [12, Theorem 4.1], we can conclude the following result.
Theorem 10. If , then system (6) undergoes a backward bifurcation at when crosses unity, while if , the bifurcation is forward.
An example of the bifurcation diagram for the model when can be seen in Figure 1, which depicts the number of infected individuals at equilibria as varies. We can see that there is a critical value such that for there exist two endemic equilibria: one stable and one unstable. As crosses the value , the two endemic equilibria coalesce at a limit point LP and disappear via a saddle-node bifurcation. When , an endemic equilibrium switches stability with the DFE at a point BP and becomes negative for .
The phase portrait when can be seen in Figure 2(a) using the parameters , , , , , , and . In this case, the DFE is globally asymptotically stable. We can also see that the origin is a saddle point: this can be verified via the linearization, but we omit this analysis since the trivial equilibrium is not interesting from an epidemiological point of view.
(a)
(b)
With the parameters , , , , , , and , we obtain the phase portrait shown in Figure 2(b) with two endemic equilibria. The first one, a saddle point, has an approximate value of , and the second one, a stable node, has an approximate value of . For this set of parameters, and . Theorem 10 indicates that if we increase in order to obtain , we have a backward bifurcation at the disease-free equilibrium .
We can calculate an expression for the critical value at which the two endemic equilibria collide with each other as follows. Assume that and . Since there are two positive equilibria for , this corresponds to case (3) of Theorem 2. These two equilibria coalesce when the discriminant given by (18) equals zero. The condition is equivalent to which can be viewed as a quadratic equation in . Since , we must take the positive root of this equation, which is given by where , and . Then, when , we can solve for in the equation and substitute, which yields
Since there is no endemic equilibrium when , it follows from Theorem 8 that the DFE will be globally asymptotically stable in such case. Hence, we can assert the following theorem, which could also be proved alternatively by the method of Lyapunov functionals.
Theorem 11. Assume , where is given by (39). Then the disease-free equilibrium of (6) is globally asymptotically stable.
4.2. Stability of Endemic Equilibria
The Jacobian matrix evaluated at an endemic equilibrium is given by and the characteristic equation is , where By the Routh–Hurwitz criteria, we know that all roots of the characteristic equation have negative real parts if and only if and . This allows us to determine necessary and sufficient conditions for local stability of the endemic equilibrium, as stated in the following theorem.
Theorem 12. Let be an endemic equilibrium of (6). Then is locally asymptotically stable if and only if the following inequalities hold:
Moreover, we can notice that the characteristic equation has one positive and one negative root whenever , so we can conclude the following result.
Theorem 13. If is an endemic equilibrium of (6) andthen is a saddle point.
5. Hopf Bifurcation
Since the disease-free equilibrium lies on the boundary of the invariant region , it is clear that there cannot exist a closed orbit surrounding . When inequality (44) holds, the endemic equilibrium is a saddle and there is no possibility of a Hopf bifurcation. Thus, we will only consider the existence of Hopf bifurcation around an endemic equilibrium when (43) holds. In the following, we will write to emphasize the dependence of the endemic equilibrium on the parameter , while we regard all other parameters as fixed.
Theorem 14. Let be an endemic equilibrium of (6) and define Assume that is a positive root of such that (43) holds when andThen system (6) undergoes a Hopf bifurcation around at .
Moreover, let be given by (56). If , there exists a family of stable periodic orbits as decreases from (supercritical bifurcation); if , there exists a family of unstable periodic orbits as increases from (subcritical bifurcation).
Proof. From the discussion in Section 4.2, we know that the characteristic polynomial at has a pair of purely imaginary roots if and only if (43) holds and . Then we know by Theorem 12 that switches stability as crosses the critical value . By calculation, we have By Theorem 3.4.2 from [13], we see that is a Hopf bifurcation point for (6) when (46) holds.
Next, set . By rewriting as , we obtain the following system equivalent to (6):By the change of variables , , we can transform (48) into where Let denote the origin of the - plane. Then since (43) holds by hypothesis. Let and make , . Then the normal form of Hopf bifurcation for (6) is where Letwhere denotes , etc. An explicit expression for can be computed using the above expressions for and and substituting . Using the software Maple [14], we obtainBy Theorem 3.4.2 and in [13], the proof is complete.
Theorem 14 can be viewed as a generalization of Theorem 3.1 of [4] that takes into account the nonlinearity of the incidence rate. Notice that the critical value at which the Hopf bifurcation takes place is given implicitly as a root of the function . Due to the complex conditions that are required to ensure the existence of Hopf bifurcation, we do not give an explicit expression for the value of in terms of the system parameters. However, its approximate value can be computed by numerical methods, as we will see in Section 7.
6. Bogdanov-Takens Bifurcation
From the analysis in Section 4.2, we know that the Jacobian matrix of an endemic equilibrium has a double zero eigenvalue if and only ifWhen this occurs, the endemic equilibrium may present a Bogdanov-Takens (BT) bifurcation. We will now investigate the conditions under which BT bifurcation takes place.
6.1. Codimension 2 BT Bifurcation
Let be an endemic equilibrium of (6) such that (57) holds. Making the change of variables , , we can transform system (6) intowhere
Condition (57) implies that and . It is clear that , so and . Using this and applying the change of variables , , we obtain the system where Finally, to obtain the normal form of BT bifurcation, we perform the near-identity transformation , . This yields the system For the BT bifurcation to be nondegenerate, we introduce the assumptions and , obtaining thus the following theorem [15, 16].
Theorem 15. Suppose that is an endemic equilibrium of (6) such that (57) holds, , and . Then is a cusp of codimension 2, that is, a Bogdanov-Takens singularity.
In the following, we will calculate the approximate saddle-node, Hopf, and homoclinic bifurcation curves by determining the versal unfolding of system (6). For that, we will use and as the bifurcation parameters. Let and , where and are close to zero. Suppose that satisfies (57) for system (6) with and . Then we obtain the perturbed system:By the transformation , , we can rewrite (63) as whereand , are smooth functions in at least of third order. We can now obtain the generic normal form: where