Abstract
A new SEIRS epidemic model with nonlinear incidence rate and nonpermanent immunity is presented in the present paper. The fact that the incidence rate per infective individual is given by a nonlinear function and product of rational powers of two state variables, as well as the introduction of an epidemic-induced death rate, leads to a more realistic modeling of the physical problem itself. A stability analysis is performed and the features of Hopf bifurcation are investigated. Both the corresponding critical regions in the parameter space and their stability characteristics are presented. Furthermore, by using algorithms based on a new symbolic form as regards the restriction of an -dimensional nonlinear parametric system to the center manifold and the normal forms of the corresponding Hopf bifurcation, as well, the associated bifurcation diagram is derived, and finally various emerging limit cycles are numerically obtained by appropriate implemented methods.
1. Introduction
The realistic modeling of epidemic models constitutes an important issue of modern research as it can contribute to both a better understanding and more accurate modeling of the actual dynamics and the interrelation of the populations involved. Nonpermanent immunity leads to SEIRS or SIRS models which have been studied with respect to the effects of the epidemiological parameters, with bilinear (see, e.g., [1]) or nonlinear incidence rate (see [2, 3] and the references therein). In particular, in [3] the writers investigate the stability of both the disease-free equilibrium and the endemic one. Also, they determine conditions regarding the existence and stability of Hopf-bifurcated limit cycles with respect to the latter, concerning both SEIRS and SIRS models. The nonlinear incidence rate offers a deeper insight into the actual relation between the populations of susceptible and infective individuals. Furthermore, we introduce an additional death rate solely due to the disease, and hence the SEIRS model is enriched with new nonlinear terms.
Let us now present the specifics of the aforementioned model, described by the following 4D differential system:where , , , , denote the number of susceptible, exposed (incubating), infective, and recovered individuals and the total population, respectively, while represents the incidence rate per infective individual and , , , , , stand for the system parameters. As regards their physical meaning, denotes the birth rate, denotes the physical death rate, denotes the rate of loss of immunity, denotes the rate of incubation, is the additional death rate due to the epidemic, and denotes the recovery rate. Then by normalizing with respect to the total population which is considered constant and taking into account (2), the system becomeswithNow, by eliminating , system (3) is reduced to the following three-dimensional one:By settingwhere , are positive constants and , (5) takes the final formThe analysis is multiparametric in that the parameter space of the system is structured by three varying parameters. Thus in Section 2 a stability analysis of the system is performed, where the active parameters are determined and various graphical representations are obtained concerning the critical (with respect to Hopf bifurcations) values of the varying parameters, as well as the critical, noncritical, and stability regions in the parameter space considered. Then, based on a new proper symbolic form as regards the center manifold analysis (see [4]), the basic features and steps of which are generally presented in Section 3, effective algorithms are implemented, by using symbolic computational software, which result in the associated bifurcation portraits throughout the regions of the parameter space under consideration. These portraits are presented in Section 4, together with limit cycles corresponding to the cases resulting from the respective analysis and obtained by use of a custom orthogonal collocation method on finite elements. Finally, Appendices A and B include algebraic manipulations and formulae related to the analysis carried out throughout this work.
2. Stability Analysis-Hopf Bifurcation
Final reduced system (7) possesses two types of equilibria: a disease-free one, namely,and an endemic one of the form with obtained after some tedious algebraic manipulations (see Appendix A), aswhere , , , , and . We focus on the endemic equilibrium with given in (9)–(11), since it corresponds to persistence of the disease.
As regards the local stability of , taking into account (7), (10), and (11), the Jacobian matrix evaluated at becomeswhereThe associated characteristic equation isNow, by considering the well-known Routh-Hurwitz necessary and sufficient stability conditions, namely,related to the equilibrium , we conclude with the formulaewhere , , , are polynomials with respect to .
Moreover, by solving the equation (resulting from (14) for ) with respect to (after substitution of the right-hand side of (9) for , we finally evaluate the real roots of a 9th-degree polynomial numerically), we further evaluate and by substituting the obtained root of into (9) and (11), respectively. Then, taking into account the fact that , we arrive at by means of (11) in the case where . Thus ifthen represent the critical values .
Then, by solving (10) with respect to and considering varying parameters of the problem we obtain the critical valueprovided that , with being the aforementioned critical equilibrium of the system (7). Thus a critical surface is generated in the parameter space (we have and due to (9) we also have , with being fixed), defined over the area of the parameter plane , where the critical values of are obtained via (17) and (18); we call this area critical region. Moreover, we differentiate , given by (16), with respect to the active parameters , namely, where and , and also and denote the partial derivatives of and with respect to , provided in Appendix A (A.2). Then by introducing the critical equilibrium and taking into account the fact that is a numerically obtained root of a high degree polynomial, we conclude that by setting the right-hand side of (19) equal to zero, no explicit relation can be extracted involving the parameters of the system. Furthermore, for any fixed values of , we numerically computeeverywhere on the critical surface. Thus, considering (17) and (20), according to Liu criterion [5] has a pair of purely imaginary eigenvalues together with a negative real one on this surface; the transversality condition (see Appendix B, ) holds, as well. Hence, a Hopf bifurcation occurs at the critical equilibrium. Graphical representations of (evaluated by using (18)) versus (for different values of , with being fixed), versus (for different values of , with being fixed), and versus (for different values of , with being fixed) are presented in Figures 1(a), 1(b), and 1(c), respectively, while critical regions are obtained in the parameter plane for fixed values of in Figures 2(a) and 2(b).
(a)
(b)
(c)
(a)
(b)
We should note that variation of the values of fixed parameters does not affect the number of critical values as regards . Thus, in any case, the expression (17), under restrictions (18), yields zero or at most one critical value (). Additionally, we note that an increase in or gives rise to an expansion of the zero region, namely, the area of the parameter plane where no critical values of exist.
Furthermore, the status of the equilibrium points corresponding to the values of in the range is shown in Figures 3(a) and 3(b). More precisely, after the right-hand side of (9) has been substituted for , by using standard numerical computation routines, (10) is solved with respect to for any given value of and fixed values for the other parameters. Then and can be obtained by means of (9) and (11), respectively, and finally the coefficients , of the characteristic equation (14) are determined. Thus, concerning the critical pairs (the points of the critical region), we focus on the slope of the real part of the complex conjugate eigenvalues of the Jacobian as a function of , inside an interval (where a pair of complex eigenvalues exist), in order to determine the direction and stability of the occurring codimension 1 Hopf bifurcation (see Appendix B). As regards the zero regions (no ), one investigates the existence of equilibrium points, as well as their stability, as varies within the range .
(a)
(b)
As a result, we conclude that regardless of the parameter values, the real part of the complex eigenvalues transverses the -axis (at ) with negative slope, all over the critical surface (for lying into the aforementioned “complex” interval, the sign of the real eigenvalue is always negative). Moreover, we encounter an unstable endemic equilibrium depending on in the zero regions, for this active parameter lying inside an interval (we have no equilibria at ), with notably sensitive to -variation, in the sense that as increases, shifts rapidly towards unity, shrinking the (unstable) “equilibrium” - interval. Additionally, by increasing or , a “- equilibrium” region emerged inside the zero region (there exist no equilibria in the whole range ), getting larger as these two rational exponents (especially ) increase.
3. Analytical Formulae for a Parametric -Dimensional System
3.1. Reduction to an -Dimensional Coordinate Space
Now, we briefly discuss the basic features of a new formulation regarding a parametric -dimensional nonlinear system analysis. The procedure adopted is presented in detail in [4, Section ], resulting in the derivation of the two-dimensional restriction to the center manifold of the system, as well as in the fast numerical computation of the Lyapunov coefficients associated with simple or degenerate Hopf bifurcations of the system. Thus considering a smooth continuous-time three-parameter system with smooth dependence on the parameters, namely, (for SEIRS system (7) we have that ), then expansion around the equilibrium path yieldswith being the Jacobian matrix evaluated at this equilibrium path. The smooth vector function represents the nonlinear terms; that is, andwherewith . Additionally, if a critical triplet exists, where has a pair of pure imaginary eigenvalues , , while the real part of the rest of the eigenvalues is negative, then a Hopf bifurcation occurs at , leading to a family of limit cycles. This implies that has a pair of complex conjugate eigenvalues , in a region around , withWe also note that the classic Hopf theory additionally demands that, at , cross the imaginary axis with nonzero speed (transversality condition). Moreover we consider the normalized complex eigenvectors and of and , respectively, having the propertiesAlso, the two-dimensional parameter-dependent real eigenspace (corresponding to ) (spanned by ), is denoted by (which becomes critical at : ), and the -dimensional real eigenspace, corresponding to all eigenvalues of other than , is denoted by . Finally, the following Lemma ([6], Lemma ) holds.
Lemma 1. if and only if .
Then, based on the decomposition by means of the third and fourth relation of (26), the last equation yieldsThus taking into account (22) and the first and second equation of (26), as well as (27), differentiation of (28) results in the reduced form of system (22) in the -dimensional coordinate space :where the dot denotes differentiation with respect to and Systems (22) and (29)-(30) are in fact dimensionally equivalent, due to the two real constraints imposed on by means of Lemma 1.
3.2. A New Symbolic Representation of the System
By expressing the complex function in the form wherewith , , then (29) and (30) become wherewith as in (33), , , .
Now, by introducing the center manifold expressionwith , and , and combining (37), as well as the invariance relation , with (34) and (35), the analysis proceeds according to the steps figured out in [4, Subsection ]. Thus after extensive algebraic manipulations combined with computer assisted calculations carried out by use of the symbolic mathematical computational software Mathematica 7, the necessary (depending on the considered bifurcation codimension) center manifold coefficients, as well as the coefficients of the system restricted to the center manifold, are derived as explicit and implicit expressions of the parameters involved. Finally, the corresponding Lyapunov coefficients are numerically evaluated fast, not only at the critical parameter values, but throughout the whole parameter space, as well, and consequently the respective bifurcation portraits can be constructed.
4. Bifurcation Results-Discussion
After setting , , , , regarding the values of the fixed parameters, by following the procedure developed in [4] (briefly discussed in Section 3 of the present paper), we arrive at the bifurcation portrait of the system with respect to the equilibrium path, presented in Figure 4. Regarding the sign of the Lyapunov coefficient, it remains strictly negative on the whole critical surface . Thus stable limit cycles are generated through supercritical Hopf bifurcations, arising for taking values in the critical region and (the real part of the complex conjugate eigenvalues is a decreasing function of around ; see Section 2). Since the critical Lyapunov coefficient never becomes zero, the system undergoes solely a codimension 1 Hopf bifurcation, for any critical triplet of the active parameters (hence, in Appendix B, we only refer (after a general inspection) to the features of the third-order normal form of the planar equation).
The fact that the limit cycles bifurcated are stable means that the phenomenon is persistent; that is, the flows in a nearby neighbourhood of the limit cycle are attracted by the cycle itself, leading to the corresponding disequilibrium fluctuations defined by the periodic trajectory, as expected in a considerable number of epidemics. The bifurcation results are verified by the computation and presentation of one cycle for specific values of the parameters by use of a custom algorithm of orthogonal collocation on finite elements, shown in Figures 5(a), 5(b), and 5(c) and a family of limit cycles obtained for different values of the epidemic-induced parameter , shown in Figures 6(a), 6(b), and 6(c), where the corresponding and also the period of the periodic orbits increase with . The stability of the obtained cycles is additionally verified by numerical computation of the respective Floquet-multipliers and exponents. For the limit cycles presented in Figure 5, let the Floquet-multipliers be , with being the respective exponents and being the fundamental period of the cycle, computed as follows: and , respectively. Moreover we note that the same cycles are generated by the variable-step, variable-order Adams-Bashforth-Moulton predictor-corrector method of orders 1 to 12, which is the standard integrated Matlab routine used to solve nonstiff ODEs, noted as “ode113” in the graphs illustrated in Figures 5 and 6.
(a)
(b)
(c)
(a)
(b)
(c)
As regards the role and significance of the introduction of the additional epidemic-induced death rate, , apart from the fact that it contributes to a more accurate and realistic description of the occurrence of the epidemic, which leads to the de facto increase of the mortality rate of infected individuals, it also offers the opportunity for a more detailed and richer parameterisation as well as the effect of the epidemic on the dynamics and interrelation between the populations involved, especially in aggressive diseases. Evidently, the introduction of gives rise to new nonlinear terms involved in all four equations of the original system, as shown in (3), and has an important effect on the numeric value of the fundamental period of the bifurcating cycles, which would have been overseen, otherwise. Moreover, the introduction of the above-mentioned parameter could also contribute to making estimations of the additional resources needed, based on the changes in the duration of critical phases and stages during an epidemic cycle and the maximum number of infected individuals (i.e., estimating the additional cost and health resources needed), thus making it possible to manage and control the epidemics more effectively, efficiently, and with a better allocation of available resources.
Concerning the results presented in [3] with respect to SEIRS model, especially the ones related to Hopf bifurcation, firstly we would like to note that in the system investigated therein an equilibrium status between the birth and the death rate is considered, without taking into account the aforementioned pure epidemic-induced death rate introduced herein. Furthermore, the bifurcation analysis in [3] is performed in a two-dimensional parameter space. More precisely the respective diagrams are obtained with respect to the rational power and the parameter (called “contact number”; see [7]). We believe that the analysis and hence the results which are established in a higher dimension parameter space, like the three-dimensional one structured in the present work (with considered as the varying parameters of the system), enhance the parametric “resolution,” that is, the physical insight into the dynamics associated with the interaction of the parameters involved in the physical background. Therefore, taking into account these essential differences, we see that, in [3], depending on the parameter values and hence the sign of the first Lyapunov coefficient, either only a subcritical bifurcation (positive sign) or an alternation (due to multiple sign changes) of subcritical and supercritical (negative sign) bifurcations occurs, located by the corresponding critical curves in the parameter plane .
5. Conclusion
In the present paper a new SEIRS epidemic model with nonpermanent immunity, nonlinear incidence rate, and an additional disease-induced death rate is presented. Then as regards the characteristic equilibrium equation, by means of the Liu criterion we conclude with the determination of critical loci in the parameter space, where a Hopf bifurcation associated with the endemic equilibrium occurs. In addition, a three-dimensional parameter space structured by the recovery, incubation, and transmission rate forms the basis of the dynamic analysis, where we proceed by using a new symbolic form introduced by the writers in a previous work, applicable in any multidimensional and multiparameter system.
Therefore we believe that the adopted general form of the system combined with the dimension of the active parameter space satisfies the demand for a more realistic modeling and offers a deeper insight into the dynamic behaviour of the epidemic, with respect to the interaction of the populations involved. Additionally, due to the lengthy expressions derived in the treatment of the dynamics concerning the restriction of the system to the center manifold and the associated normal forms, as well, the algebraic platform mentioned above is proved to be appropriate in manipulating a large amount of analytical data. Thus constructing algorithms of an implicit structure, by use of chain numerical computations, all the quantities associated with the bifurcation analysis are evaluated fast and the bifurcation diagrams of the system can be obtained throughout the whole parameter space under consideration, for any values of the fixed parameters.
We should further note that variations of the fixed parameters change the shape and the size of the critical and zero regions, without affecting the qualitative profile of the results. In particular, increase in the rational exponents involved in the incidence rate, especially in that concerning the infective individual, results in a shrinkage of the critical region, while the zero region, as well as the nonequilibrium area inside that region (as regards the endemic one), expands towards lower values of the recovery rate. Finally, the algorithm of orthogonal collocation on finite elements with Legendre orthogonal polynomials is proved to be excellent in the fast and precise numerical computation of the bifurcated stable limit cycles.
Appendix
A. Algebraic Treatment of System (7)
A.1. Algebraic Manipulations Leading to the Final Form of System (7)
System (7), with , results inBy replacing the second equation of (A.1) with the sum of the last two ones, we getThen, by multiplying the second relation of (A.2) by and the last one by and adding them up we obtainBy further multiplying the last relation of (A.3) by and adding it to the second one, we obtainMoreover, by multiplying the first relation of (A.4) by and the third one by and adding them up we get Then by multiplying the second relation of (A.5) by and by adding the first two relations we finally obtainThe latter system yields (9)–(11).
A.2. Derivatives of with respect to the Varying Parameters
The derivatives of , with respect to the varying parameters of the problem , are involved in the transversality condition associated with the emerged Hopf bifurcation (see relation in Appendix B), as well as in the equivalent (according to Liu criterion) derivatives of (see (19)). Thus taking into account the fact that and are obtained from (10) and (9), respectively, asthen by differentiating (9) and (10), and by denoting the partial derivative of with respect to , and the partial derivatives of and with respect to the parameters , , and , respectively, the following formulae are obtained:where is given by (13) and
B. Normal Forms
B.1. Poincaré Normal Forms for the Planar Case
In this case the terms with in the first part of (33) and in (34) vanish. Moreover, (35) does not exist. Now, in order to obtain the desired normal form, we first introduce the transformationthe inversion of which is given by (see [6, Section ])Then by differentiating (B.2) with respect to and substituting the right-hand side of (29) for (with ), as well as the corresponding conjugate equation for , by keeping terms up to order with respect to , we obtainFinally substitution of the transformation (B.1) for (as well as the respective conjugate relation for ) in (B.3) yieldswhere terms up to order with respect to have been evaluated. Now, by setting equal to zero we define the coefficients of transformation (B.1), while the rest of the resonant terms of the odd order, , constitute the desirable normal form of (29). Furthermore, by substituting the obtained expressions of involved in () and then by means of computer calculations we arrive at the analytic formulae of the normal coefficientswhich are lengthy for .
B.2. Normal Form of the Third Order
In this case, the above procedure results in the third-order normal form:where represents the cubic resonant term, due to the term included in (we set ). Then by means of the linear time scalingand the nonlinear time reparameterisation (B.6) becomeswhere the parameter function represents the first Lyapunov coefficient. If there exist critical values of the parameters, wherethen by the linear local (around ) transformation , (B.9) yieldswith . According to [6] (Lemma , valid for both supercritical and subcritical case) the terms which are higher than the third order in (B.10) do not affect the bifurcation behaviour of the system near , and thus we consider the locally topologically equivalent system (B.10), where the terms have been dropped. Moreover, if the transversality condition holds at the critical values of the parameters, that is, then one limit cycle is bifurcated from the equilibrium at . By using polar coordinates , the equivalent system (B.10) (without the terms) gives , , corresponding to the well-known supercritical (, ) and subcritical (, ) cases of the Hopf bifurcation. We call a Hopf point of codimension 1 (H1 point).
In the case where , we proceed to higher order normal forms.
Conflicts of Interest
The authors declare that there are no conflicts of interest regarding the publication of this paper.
Acknowledgments
P. S. Douris is pleased to acknowledge his financial support from “Andreas Mentzelopoulos Scholarships for the University of Patras.”