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: wherewhile and can be computed by the formulas in [17, Equation (53)]. If the transversality conditionis satisfied, we can obtain local approximations of the bifurcation curves for the original system, as stated in the following theorem.
Theorem 16. Suppose that is an endemic equilibrium of (6) such that (57) and (68) hold, , and . Then the system admits the following behavior near . (1)There is a saddle-node bifurcation curve:(2)There is a Hopf bifurcation curve:(3)There is a homoclinic bifurcation curve:
Theorems 15 and 16 generalize the corresponding results by Li et al. [4, Theorems 4.1 & 4.2], which were proved for a model with bilinear incidence, by taking into account a saturated nonlinear incidence rate. It is worth noticing that our results require the hypotheses and , which were not present in [4]. When these conditions are not met, system (6) may undergo a bifurcation of codimension higher than 2, which exhibits the richer dynamics that our model has.
6.2. Codimension 3 BT Bifurcation
From Theorem 15, we know that the endemic equilibrium is a BT point of codimension 2 when (57) holds, and . If , then may be a cusp of codimension 3. In such case, we can compute the corresponding normal form for system (6) by following the same steps as in [18, 19], which we summarize as follows:(i) First, we take the Taylor expansion (58) of system (6) up to fourth-order terms and diagonalize the linear part using the transformation , , obtaining (ii) We apply a near-identity transformation to eliminate the non-resonant terms, obtaining the system (iii) Apply a near-identity transformation to eliminate the term and rescale the coordinates in the above system, which becomes (iv) Make and and eliminate the term by the transformation , , , obtaining where Alternatively, the coefficient can be computed by the formulas in [20].
Conditions and imply that system (6) is topologically equivalent to (75) near the endemic equilibrium , so we can conclude the following.
Theorem 17. Suppose that (57) holds. If , and , then is a BT point of codimension 3, and the system (6) localized at is topologically equivalent to (75).
In [21], it was shown that a generic unfolding with the parameters of the codimension 3 cusp singularity is equivalent toFor this system, the plane , excluding the origin in the parameter space, is a saddle-node bifurcation surface. When , system (77) has no equilibria, so the bifurcation surfaces are in the half-space . The bifurcation diagram of system (77) is a cone, which can be visualized by drawing its intersection with the 2-sphere: when is sufficiently small, as shown in Figure 3.
There are three bifurcation curves on : a Hopf bifurcation curve , a homoclinic bifurcation curve , and a double limit cycle bifurcation curve . The curve connects the two points and and is tangent to the curves and at these two points respectively. The curves and touch the boundary at the points and . In the neighborhood of and , one can find the unfolding of the cusp singularity of codimension 2. System (77) has a unique unstable limit cycle when the parameters lie between and in a small neighborhood of , and it has a unique stable limit cycle when the parameters lie between and in a small neighborhood of . Inside the curved triangle , the system has two limit cycles: the inner one is unstable and the outer one is stable. These two limit cycles coalesce at , where there exists a unique semi-stable limit cycle.
Remark 18. Together, Theorems 15–17 allow us to assert that the model will undergo a Bogdanov-Takens bifurcation at the endemic equilibrium when the conditions (57) and are met. The codimension of the bifurcation can be further determined by calculating the coefficient ( is a cusp of codimension 2 when and codimension 3 otherwise). All coefficients for the normal forms computed in this section can be expressed in terms of the original parameters of the model via a chain of substitutions (see formulas (67) and (59)); however, the resulting expressions would be too unwieldy to be useful. In practice, these coefficients can be calculated by numerical software packages like Matcont, which implement the computation of normal forms coefficients by more efficient methods based on combining the approximation of the center manifold with the normalization on it and using the Fredholm Alternative. Matcont also uses orthogonal collocation for the discretization of periodic and homoclinic orbits. For homoclinic orbits, it also uses the Lindstedt–Poincaré method to achieve a better rendering of these orbits compared to the regular perturbation method [22, 23]. Therefore, the analytical results developed here serve as the theoretical foundations, which is the mathematically solid part that gives support to what is calculated by numerical simulation software.
In this section and the previous one, we described the lengthy procedure required in order to compute the coefficients of the normal forms for the Hopf and codimension 2 and 3 Bogdanov-Takens bifurcations. In the next section, their numerical value will be actually obtained providing a set of numerical values of all the parameters involved in the model.
7. Numerical Continuation and Simulations
In this section, we will illustrate the analytical results obtained in the previous sections by performing numerical continuation of equilibria and numerical simulations. We will use the Matcont package [24] to carry out the numerical continuation of equilibria in order to obtain bifurcation curves of system (6) for different sets of parameters and locate points of interest like the codimension 2 BT points and the approximate location of a codimension 3 BT point. Matcont also allows us to perform the numerical continuation of limit cycles and homoclinic orbits. For some scenarios, the corresponding phase portraits will be plotted. The coefficients of the corresponding normal forms are also provided by Matcont, complementing the description of the long chain of substitutions carried out in the previous two sections in order to compute the numerical values of such coefficients.
Example 1. We consider first system (6) with the parameters from Table 1 and choose as a bifurcation parameter. Using numerical software, we can compute the endemic equilibria of the system as is varied and check that the function defined in Theorem 14 has a positive root (see Figure 4(a)). Similarly, it is readily seen that (43) and (46) hold for . Hence, all hypotheses of Theorem 14 are met, which implies that the model undergoes a Hopf bifurcation at . Since the coefficient is positive, the bifurcation is subcritical. The bifurcation diagram shown in Figure 4(b) reveals that, when , there is a unique limit cycle (represented by green curves in the diagram), which is stable. As becomes less than 0.5, an unstable limit cycle (blue curves in the diagram) appears around the endemic equilibrium . When is further reduced to 0.43, the two limit cycles collide and disappear via a double limit cycle bifurcation. The phase portrait of the system when can be seen in Figure 5: in this case, is unstable and all solutions approach the unique limit cycle. If we change to 0.48, the endemic equilibrium becomes stable, and we can see that there are two limit cycles (Figure 6).
With these parameters, we also get the bifurcation diagram in the plane shown in Figure 7. Note that , implying that there is a backward bifurcation at when , and for this value of the bifurcation parameter we have This is labeled as BP (branching point), the point where the unstable endemic equilibrium meets the disease-free equilibrium The red dotted curve corresponds to that unstable equilibrium point, which coalesces with the stable endemic equilibrium at LP (limit point, where a saddle-node bifurcation takes place). The continuous blue curve describes the movement of the stable endemic equilibrium as increases from , where The label on that curve indicates that after this point as we move along the blue curve. Label H indicates a subcritical Hopf bifurcation at that point on the blue curve, which occurs at . As Figure 4(b) shows, the outer stable limit cycle and the inner unstable limit cycle meet at and disappear for less than this critical value (saddle-node bifurcation of limit cycles). Associated with point labeled H, is the corresponding point in the bifurcation graph shown in Figure 8.
(a)
(b)
Example 2. If we fix the parameters , , , , and let both and vary, we can see by an application of Theorem 15 that one of the endemic equilibria undergoes a Bogdanov-Takens bifurcation at (denoted by ) and at (denoted by ). The corresponding normal form coefficients and , as defined in Section 6, can be computed numerically, and we obtain at and at . Since the coefficients are nonzero, this implies that both BT points are cusps of codimension 2.
If we plot the bifurcation curves on the plane, we obtain the diagram shown in Figure 8. Near the point, we can see a saddle-node bifurcation curve and a subcritical Hopf bifurcation curve (Figure 9(a)). The saddle-node bifurcation (magenta line) and neutral saddle curves (green line) connect the two Bogdanov-Takens points. Near , there exists also a supercritical Hopf bifurcation curve (Figure 9(b)).
(a)
(b)
Example 3. We will now plot the bifurcation curves on the plane for different values of , the other parameters being the same as in Example 2. Figure 10 shows the corresponding saddle-node bifurcation curves for taking several values between 0.5 and 0.0916. We can see that the two BT points get closer to each other as decreases, and they coalesce and become a BT point of codimension 3 at the critical value (which is approximately ), as stated in Theorem 17. At such point, the normal form coefficient vanishes, while the other coefficients become approximately and .
If we fix the parameters and , then a supercritical Hopf bifurcation occurs as is varied. Figure 11(a) shows the phase portrait of the system when and the endemic equilibrium is stable. When we increase the value of to 0.97, the endemic equilibrium loses its stability and a stable limit cycle appears around it, as we can see in Figure 11(b).
(a)
(b)
(a)
(b)
Example 4. In the following set of figures, the dynamics near the point shown in Figures 8 and 9(a) will be analyzed in close detail. Figure 12 shows a 3D bifurcation diagram that includes Figures 7 and 9(a). The axes are the control parameters and used to locate the Bogdanov-Takens points, as well as , the infected population. Note that the graph of Figure 7 is embedded in the vertical slice with the constant value of because this curve is generated varying only and keeping the mentioned value of constant. The region displayed is the one surrounding the point. BP−R1 is the bifurcation point where the unstable endemic equilibrium meets the disease-free equilibrium. R1 denotes that here , as it also occurs at the other points with this label along the bifurcation curves.
Example 5. Figure 13 shows the same region as Figure 12, but including the information of the homoclinic orbits arising from the BT bifurcation point. This information is shown as a blue vertical plane surface following the path of the homoclinic bifurcation curve in the - parameter plane. The vertical coordinate in this surface shows the varying values of the state variable in the homoclinic orbits. Figure 14 depicts the same region in the - parameter space, but the coordinate axes are , , and now. This graph allows us to notice the homoclinic orbits stacked as vertical slices of a 3D surface starting at the BT bifurcation point, partially obscured by this surface. The vertical slices are generated as we travel along the homoclinic bifurcation curve in the - parameter space.
The homoclinic bifurcation curve is displayed as the brown curve starting at the BT bifurcation point and to the left of the Hopf bifurcation curve of this BT point in Figure 15, a 2D figure displaying the region in - parameter space. This curve intersects the neutral saddle curve. At this parameter space point, the saddle point at the end of the corresponding homoclinic orbit is a neutral saddle. The curve ends at the other BT point that is not displayed in this picture. This curve is displayed until its intersection with the Hopf bifurcation curve. Note that at this point is where the backward bifurcation takes place (the endemic saddle point and the disease-free equilibrium collide here) and . Given that the homoclinic orbits corresponding to the part of this curve from a little bit before the intersection with the neutral saddle curve to its intersection with the Hopf bifurcation curve are packed very close together, they are not displayed in Figures 13, 14, and 15 in order to make them clearer. In the remaining part of the curve that is not displayed and that goes almost parallel to the Hopf bifurcation curve, the curve crosses a region of the parameter space where is equal to 1 or is very close to 1, and the corresponding homoclinic orbits are extremely close together; therefore they are also not displayed in Figures 13, 14, and 15. In Figure 15, the green curve is the saddle-node bifurcation curve. Note that in the left part of this green curve is the cusp bifurcation point also displayed in Figure 9(a). The magenta curve is the Hopf bifurcation and neutral saddle curve. The part of it starting from the BT point and to the right of the (brown) homoclinic bifurcation curve is the Hopf part, and the other part starting to the left of the BT point is the neutral saddle part, which was also displayed as a green curve in Figure 9(a). This part intersects the Hopf part as we travel from the BT point, but this intersection does not mean that the neutral saddle equilibrium point is colliding with a Hopf bifurcation point (this is evident in Figure 12). The correct interpretation of this intersection is that for this particular value of the pair there are a Hopf bifurcation point and a neutral saddle equilibrium point, but both located at different points in the - state space. Figure 15 also displays the graph of Figure 7 as the horizontal blue curve with constant .
Figure 16 is a 2D plot with coordinates that also depicts the same region. The curves shown in this picture are the paths of the corresponding equilibrium points (bifurcation curves) as and vary in parameter space. The blue region shows the graphs of the homoclinic curves as the control parameters and vary. In the case of the graph of Figure 7 included here (the curve in aquamarine and red colors), has a fixed value of 23.5 and varies. The aquamarine part is the path of the stable endemic equilibrium, and the red part is the path of the unstable endemic equilibrium. The points labeled as R1 denote when The point is at the same time a branching point and labeled as R1. A cusp point along the green curve is very close to that point; hence a CP label is also there. In all the previous figures mentioned in this analysis, all the remaining parameters are fixed with the values given in Table 1.
Figures 17, 18, and 19 display phase planes in three different scenarios. The first one (Figure 17) corresponds to a point near the end of the homoclinic bifurcation curve, with . Note that there is one homoclinic orbit separating the green solution that departs from very close to it towards a stable endemic equilibrium point and the blue solutions elsewhere. Figure 18 corresponds to , a point located the in - parameter space between the homoclinic bifurcation point and the Hopf bifurcation point. The green solution emanates from very close to an unstable limit cycle and wraps around towards a stable endemic equilibrium point. Solutions parting from outside the limit cycle are depicted in magenta. Solutions elsewhere are depicted in blue. Finally, Figure 19 has , a point located in the - parameter space to the right of the Hopf bifurcation point. Note the stability change of the previously stable endemic equilibrium point. All green solutions depicted diverge from this unstable endemic equilibrium point and there is no longer an unstable limit cycle. All remaining solutions are depicted in blue. Hence, the Hopf bifurcation is subcritical as we move from the scenario of Figure 18 to the scenario of Figure 19. In all these Figures 17–19, there are two other equilibria: the disease-free equilibrium, which is always a stable node, and an unstable endemic equilibrium point, which is a saddle point. Note that for the three scenarios, and case (3) of Theorem 2 holds here.
8. Conclusions
We generalized the SIR epidemic model with logistic growth and Holling type II treatment studied in [4] by studying a new model with a nonlinear saturated force of infection.
We proved that the local stability of the disease-free equilibrium depends only on the basic reproduction number : the DFE is stable when and unstable when . Our model may present a backward bifurcation under certain conditions, which implies that it may not be enough to reduce below unity in order to obtain the global asymptotic stability of the DFE. Theorem 11 shows that the disease is expected to become extinct when the basic reproduction number is sufficiently low, which can be achieved by decreasing the parameter or increasing . This implies that, just like in many other epidemic models, there is a threshold value for that guarantees the eradication of the disease, although this value becomes smaller and hence more difficult to achieve for epidemic control strategies due to the bifurcation dynamics.
We showed that, when , there is a unique endemic equilibrium, whose stability depends on inequalities (42) and (43). Our results indicate that the behavior of the model around endemic equilibria can become very complicated since variations in the parameters may induce different types of bifurcation.
As the disease transmission rate increases, our model undergoes a Hopf bifurcation under certain conditions, and both stable and unstable periodic solutions may emerge around an endemic equilibrium. When both and the treatment parameter are varied, a BT bifurcation of codimension 2 or 3 may occur. The use of numerical software can help us to corroborate our analytical results and to detect other types of complicated dynamics, such as cusp bifurcation and generalized Hopf bifurcation (see the points labeled CP and GH in Figure 8). The use of different 2D and 3D bifurcation diagrams allowed us to describe in some detail different aspects of the complicated dynamics due to the interplay of the nonlinear force of infection with a nonlinear treatment term. and play an important role as control parameters to detect the codimension 2 BT bifurcation points, and adding as an additional control parameter was required to locate the approximate position of the codimension 3 BT bifurcation point by running several times the required Matcont numerical routines used to locate the codimension 2 BT points, each run with a different fixed value of in a decreasing sequence until it disappears (Figure 10). Several phase portraits were also included to show some of the different scenarios that arise due to these complex dynamics.
As it happens with other epidemic models studied in the literature, such as [4], we can see that the nonlinear force of infection is an important factor that can lead to complex dynamics. It is noteworthy that our model has a similar mathematical structure to other models that arise in different application areas. For example, the predator–prey model studied in [6] is similar to (6) except that the saturated function would be replaced with and in the first and second equations of the system, respectively, and is interpreted as a Holling type II functional response. The authors in [6] studied the dynamics of that model and found the existence of Hopf and BT bifurcations, but only up to codimension 2. Another predator–prey model related to this is the one studied in [7], which uses the bilinear functional response and presents similarly varied dynamics. Our analysis for model (6) went further to explore the existence of codimension 3 BT bifurcation, which shows that our model may represent a more drastic variation of the population dynamics.
Given that the final purpose of our model is to eradicate the disease, it will be useful as long as it allows us to provide an effective control strategy. The numerical experiments performed in the previous section give us a hint of how to achieve this goal. From all the parameters of the model, we identified , , and as suitable control parameters. Varying has an effect to modify the separation of the two codimension 2 bifurcation points and it also has some role in modifying the saddle-node bifurcation curve, as Figure 10 shows. This bifurcation curve is very important towards an effective disease control strategy. This is evident by noting that the disease dies out in the region of the - plane to the left of the SN bifurcation curve because solutions converge to the disease-free equilibrium point as tends to . Hence, parameters and have a major role in controlling the disease. Parameter can be modified by implementing isolation strategies like quarantine or asking people to stay at home for as long as possible or avoiding public spaces in order to minimize exposure to infected individuals. The parameter measures the effectiveness of the treatment as the power of the medicines used to treat the infected individuals, and can be considered to be related to medicine availability. Hence, since the other parameters are related to the population, such as growth and mortality rates (natural and disease-induced), population carrying capacity, and disease recovery rate, they can be considered fixed and are determined by the disease and population susceptible to it, and thus the remaining parameters are the three previously considered as control parameters.
Therefore, the steps towards an effective control strategy are as follows: First, determine the population and disease-related parameters and set a value for according to treatment availability. Then, perform the numerical continuation of equilibria using Matcont to obtain the saddle-node bifurcation curve, and finally, choose a feasible pair of and such that the corresponding point in the - plane lies to the left of the bifurcation curve. Note that if for some reason it is not possible to choose suitable values of and that lay there, the numerical continuation of equilibria provided by Matcont can give us information about the long term behavior of the disease given a known initial population of susceptible and infected individuals. In some cases, such as the ones depicted in Figures 17 and 18, we can choose values of and such that and there are two endemic equilibria (one stable and one unstable), so the positive quadrant of the plane is divided into two regions by the stable manifold of the saddle endemic equilibrium. If our initial point lies between the unstable endemic equilibrium and the DFE, it could lay in the basin of attraction of the disease-free equilibrium and the disease will die out. In order to achieve this, the values of and should make the endemic equilibria to be as close to each other as possible.
For other values of the parameters, periodic solutions of the system may appear around endemic equilibria (see, for example, Figures 5 and 11(b)). The existence of stable limit cycles in certain regions of the - - parameter space can provide theoretical explanations for the appearance of periodic outbreaks of a disease (or oscillatory dynamics and coexistence of predator and prey populations, in the case where the system is interpreted as a predator–prey model). Numerical and theoretical analysis of the different types of bifurcations allows us to show that small variations in the model parameters can trigger a significant change in the long-term behavior of the subpopulations under study in both epidemic and ecological models.
Data Availability
The data used to support the findings of this study are included within the article.
Conflicts of Interest
The authors declare that there are no conflicts of interest regarding the publication of this paper.
Acknowledgments
This work was supported by Universidad Autónoma de Yucatán and Conacyt SNI under Grants nos. 15284 and 33365.