Hopf Bifurcation Analysis of a New SEIRS Epidemic Model with Nonlinear Incidence Rate and Nonpermanent Immunity

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 n-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.


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 ( By setting where ,  are positive constants and  > 1, (5)  ( The 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.
As regards the local stability of Σ 1 , taking into account ( 7), (10), and (11), the Jacobian matrix evaluated at Σ 1 becomes International Journal of Mathematics and Mathematical Sciences 3 where  0 =  0 (, , , , , , , ) ,  0 =  0 ( 0 ; , , , , ) ,  =  ( 0 ; , , , , ) The associated characteristic equation is Now, by considering the well-known Routh-Hurwitz necessary and sufficient stability conditions, namely, related to the equilibrium Σ 1 , we conclude with the formulae where   =   ( 0 ; , , , , , , ),  = 0, 1, 3,  = 0, 1, 2, are polynomials with respect to  0 .Moreover, by solving the equation  3 = 0 (resulting from (14) for  = ) with respect to  0 (after substitution of the right-hand side of (9) for  0 , we finally evaluate the real roots of a 9th-degree polynomial numerically), we further evaluate  0 and  0 by substituting the obtained root of  3 into (9) and (11), respectively.Then, taking into account the fact that  2 > , we arrive at  0 > 0 by means of (11) in the case where 0 <  0 < 1.Thus if then ( 0 ,  0 ,  0 ) represent the critical values ( 0 cr ,  0 cr ,  0 cr ).Then, by solving (10) with respect to  and considering (, , ) varying parameters of the problem we obtain the critical value provided that 0 <  cr < 1, with  0 cr ,  0 cr being the aforementioned critical equilibrium of the system (7).Thus a critical surface  cr =  cr (, ) is generated in the parameter space (, , ) (we have  0 cr =  0 cr (, ) and due to (9) we also have  0 cr =  0 cr (, ,  0 cr (, )), 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  3 , given by (16), with respect to the active parameters  ( = , , ), namely, where ◻ , = ◻/ and ◻ , 0 = ◻/ 0 , and also  0  and  0  denote the partial derivatives of  0 and  0 with respect to , provided in Appendix A (A.2).Then by introducing the critical equilibrium ( 0 cr ( 0 cr ),  0 cr ) and taking into account the fact that  0 cr 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 compute everywhere on the critical surface.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 (0 <  cr < 1).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 (0, 1) is shown in Figures 3(a numerical computation routines, (10) is solved with respect to  0 for any given value of  and fixed values for the other parameters.Then  0 and  0 can be obtained by means of ( 9) and (11), respectively, and finally the coefficients   ,  = 0, 1, 2 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 ( cr − 1 ,  cr + 2 ),  1 ,  2 > 0 (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  cr ), one investigates the existence of equilibrium points, as well as their stability, as  varies within the range (0, 1).
As a result, we conclude that regardless of the parameter values, the real part of the complex eigenvalues transverses the -axis (at  cr ) 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 (, 1),  = (, , , , , , ) > 0 (we have no equilibria at (0, )), 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 "0-equilibrium" region emerged inside the zero region (there exist no equilibria in the whole range (0, 1)), getting larger as these two rational exponents (especially ) increase.

Analytical Formulae for a Parametric
-Dimensional System 3.1.Reduction to an ( + 2)-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 2], 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 ( 1 ,  2 ,  3 ) = (, , )), then expansion around the equilibrium path yields with  0 () =   ( 0 (); ) being the Jacobian matrix evaluated at this equilibrium path.The smooth vector function  : R +3 → R  represents the nonlinear terms; that is,  = (‖‖ 2 ) and where with   = (   1 , . . .,    )  ∈ C  ,  = 1, . . ., .Additionally, if a critical triplet  0 = ( 10 ,  20 ,  30 ) exists, where  0 has a pair of pure imaginary eigenvalues  1,2 = ± 0 ,  0 > 0, while the real part of the rest of the eigenvalues is negative, then a Hopf bifurcation occurs at  0 , leading to a family of limit cycles.This implies that  0 has a pair of complex conjugate eigenvalues (), () in a region around  0 , with We also note that the classic Hopf theory additionally demands that, at  0 , () cross the imaginary axis with Then, based on the decomposition 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)   +  (3)

Bifurcation Results-Discussion
After setting  = 10 −4 ,  = 10 −4 ,  = 10 −3 ,  = 0.5,  = 1.2 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 Σ 1 equilibrium path, presented in Figure 4. Regarding the sign of the Lyapunov coefficient, it remains strictly negative on the whole critical surface  cr =  cr (, ).Thus stable limit cycles are generated through supercritical Hopf bifurcations, arising for (, ) taking values in the critical region and  <  cr (, ) (the real part of the complex conjugate eigenvalues is a decreasing function of  around  cr ; see Section 2).Since the critical Lyapunov coefficient  1 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) 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 where the corresponding  cr 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   =    ⋅ ,  = 1, 2, 3, with   = (1/) ⋅ ln(  ) being the respective exponents and  being the fundamental period of the cycle, computed as follows:   = {1, 0.82, 1.1 × 10 −17 } and   = {0, −4.7 × 10 −4 , −0.095}, respectively.Moreover we note that the same cycles are generated by the variablestep, variable-order Adams-Bashforth-Moulton predictorcorrector 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.
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 Figure 5: Stable cycle generated by supercritical Hopf bifurcation in  (),  () planes and in 3D (), respectively, for (, , , , , , , ) = (10 −4 , 10 −4 , 10 −3 , 0.5, 1.2, 0.042, 0.5, 0.328525) ( cr = 0.338525), determined by a custom orthogonal collocation on finite elements algorithm.Unstable endemic equilibrium (red marker): ( 0 ,  0 ,  0 ) = (0.073948, 0.001986, 0.023535).Period:  = 410.572533days.ode113: Adams-Bashforth-Moulton PECE solver of orders 1 to 12 (Matlab).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 (, ).

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.
In the case where  1 ( 0 ) = 0, we proceed to higher order normal forms.

Figure 3 :
Figure 3: Stability features regarding the critical (green) and the zero (black and red) regions.Green denotes the negative slope of the real part of the complex eigenvalues versus , around  cr , red indicates the existence of unstable endemic equilibria for  ∈ (, 1), and black outlines the pairs (, ) for which no endemic equilibria arise in the whole range (0, 1) for .The parameters are defined as  = 10 −4 ,  = 10 −4 ,  = 10 −3 , and  = 0.5 and (a)  = 1.2 and (b)  = 1.6.
), 5(b), and 5(c) and a family of limit cycles obtained for different values of the epidemicinduced parameter , shown in Figures 6(a), 6(b), and 6(c),
Symbolic Representation of the System.By expressing the complex function (, , ) in the form