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

Authors: M. P. Markakis; P. S. Douris
Journal: International Journal of Mathematics and Mathematical Sciences
 (2018)

---

## 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 ann-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.


---

## Body

## 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:(1)S˙=rN-dS+εR-hS,II,E˙=-d+σE+hS,II,I˙=-γ+d+αI+σE,R˙=γI-d+εR,(2)N˙=rN-dN-αI,where S, E, I, R, N denote the number of susceptible, exposed (incubating), infective, and recovered individuals and the total population, respectively, while h(S,I) represents the incidence rate per infective individual and r, d, ε, σ, α, γ stand for the system parameters. As regards their physical meaning, r denotes the birth rate, d 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 N=S+E+I+R which is considered constant and taking into account (2), the system becomes(3)x˙=r-rx+εz+αxy-h-x,yy,w˙=-r+σw+αwy+h-x,yy,y˙=σw-r+γ+αy+αy2,z˙=γy-r+εz+αyzwith(4)x=SN,w=EN,y=IN,z=RN,x+w+y+z=1.Now, by eliminating z=1-x-w-y, system (3) is reduced to the following three-dimensional one:(5)x˙=r+ε-r+εx-εw-εy+αxy-h-x,yy,w˙=-r+σw+αwy+h-x,yy,y˙=σw-r+γ+αy+αy2.By setting(6)h-x,y=βxmys-1,where β, m are positive constants and s>1, (5) takes the final form(7)x˙=r+ε-r+εx-εw-εy+αxy-βxmys,w˙=-r+σw+αwy+βxmys,y˙=σw-r+γ+αy+αy2.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.


## 2. Stability Analysis-Hopf Bifurcation

Final reduced system (7) possesses two types of equilibria: a disease-free one, namely,(8)Σ0=1,0,0,and an endemic one of the form Σ1(x0,y0,z0) with y0≠0 obtained after some tedious algebraic manipulations (see Appendix A), as(9)x0=1σαy0-κ0α2y03-ακ4y02+κ3+σεy0-σκ0,(10)βσx0my0s-1=α2y02-ακ4-εy0+κ3-εκ2,(11)w0=1σ-αy02+κ2y0,where κ0=r+ε, κ1=r+ε+σ, κ2=r+γ+α, κ3=κ1κ2, and κ4=κ1+κ2. We focus on the endemic equilibrium Σ1=x0,y0,w0 with x0,y0,w0 given in (9)–(11), since it corresponds to persistence of the disease.As regards the local stability ofΣ1, taking into account (7), (10), and (11), the Jacobian matrix evaluated at Σ1 becomes(12)J0=αy0-κ0-mσAy0x0-ε-ε+αx0-sσAmσAy0x0αy0-κ1+εασ-αy02+κ2y0+sσA0σ2αy0-κ2,where(13)y0=y0γ,r,α,ε,σ,β,m,s,x0=x0y0;γ,r,α,ε,σ,A=Ay0;γ,r,α,ε,σ=α2y02-ακ4-εy0+κ3-εκ2.The associated characteristic equation is(14)λ3+B2λ2+B1λ+B0=0.Now, by considering the well-known Routh-Hurwitz necessary and sufficient stability conditions, namely,(15)B0>0,B1>0,B1B2-B0>0,related to the equilibrium Σ1, we conclude with the formulae(16)B0=1σx0P00+P01x0,B1=1σx0P10+P11x0,B3=B1B2-B0=1σ2x02P30+P31x0+P32x02,where Pij=fij(y0;γ,r,α,ε,σ,m,s), i=0,1,3, j=0,1,2, are polynomials with respect to y0.Moreover, by solving the equationB3=0 (resulting from (14) for λ=iω) with respect to y0 (after substitution of the right-hand side of (9) for x0, we finally evaluate the real roots of a 9th-degree polynomial numerically), we further evaluate x0 and w0 by substituting the obtained root of B3 into (9) and (11), respectively. Then, taking into account the fact that κ2>α, we arrive at w0>0 by means of (11) in the case where 0<y0<1. Thus if(17)0<x0<1,0<y0<1,w0<1,B0>0,B1>0then (x0,y0,w0) represent the critical values (xcr0,ycr0,wcr0).Then, by solving (10) with respect to β and considering (γ,σ,β) varying parameters of the problem we obtain the critical value(18)βcr=α2ycr02-ακ4-εycr0+κ3-εκ2σxcr0mycr0s-1provided that 0<βcr<1, with xcr0,ycr0 being the aforementioned critical equilibrium of the system (7). Thus acritical surfaceβcr=βcr(γ,σ) is generated in the parameter space (γ,σ,β) (we have ycr0=ycr0(γ,σ) and due to (9) we also have xcr0=xcr0(γ,σ,ycr0(γ,σ)), with r,α,ε,m,s being fixed), defined over the area of the parameter plane (γ,σ), where the critical values of β are obtained via (17) and (18); we call this areacritical region. Moreover, we differentiate B3, given by (16), with respect to the active parameters p(p=γ,σ,β), namely, (19)B3,p=1σ2x02P30,p+P31,px0+P32,px02+P31+2P32x0x0p-2B3σx0δx0+σx0p,δ=1,p=σ0,p=γ,βP3j,p=f3j,p+f3j,y0y0p,j=0,1,2,where □,p=∂□/∂p and □,y0=∂□/∂y0, and also x0p and y0p denote the partial derivatives of x0 and y0 with respect to p, provided in Appendix A (A.2). Then by introducing the critical equilibrium (xcr0(ycr0),ycr0) and taking into account the fact that ycr0 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 r,α,ε,m,s, we numerically compute(20)B3,pγ,σ,βcr,xcr0,ycr0≠0,everywhere on thecritical surface. Thus, considering (17) and (20), according to Liu criterion [5] J0 has a pair of purely imaginary eigenvalues together with a negative real one on this surface; thetransversality condition (see Appendix B, C1-2) holds, as well. Hence, a Hopf bifurcation occurs at the critical equilibrium. Graphical representations of βcr (evaluated by using (18)) versus γ (for different values of σ, with r,α,ε,m,s being fixed), versus σ (for different values of γ, with r,α,ε,m,s being fixed), and versus s (for different values of γ, with r,α,ε,σ,m being fixed) are presented in Figures 1(a), 1(b), and 1(c), respectively, whilecritical regions are obtained in the parameter plane (γ,σ) for fixed values of r,α,ε,m,s in Figures 2(a) and 2(b).Figure 1
The critical value,βcr, as a function of the parameters: (a) γ for σ=0.25 (magenta) and σ=0.75 (green), (b) σ for γ=0.01 (magenta), γ=0.05 (blue), and γ=0.075 (green), and (c) s for σ=0.5 and γ=0.01 (magenta), γ=0.06 (blue), γ=0.12 (green), and γ=0.17 (red).
(a)
(b)
(c)Figure 2
Critical regions (magenta – one critical value for β) andzero regions (black – no critical value for β), for r=10-4, α=10-4, ε=10-3, and m=0.5 and (a) s=1.2 and (b) s=1.6.
(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 (0<βcr<1). Additionally, we note that an increase in s or m gives rise to an expansion of thezero 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) and 3(b). More precisely, after the right-hand side of (9) has been substituted for x0, by using standard numerical computation routines, (10) is solved with respect to y0 for any given value of β and fixed values for the other parameters. Then x0 and w0 can be obtained by means of (9) and (11), respectively, and finally the coefficients Bi, i=0,1,2 of the characteristic equation (14) are determined. Thus, concerning thecritical pairs (γ,σ) (the points of thecritical 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-d1,βcr+d2),d1,d2>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 thezero regions (no βcr), one investigates the existence of equilibrium points, as well as their stability, as β varies within the range (0,1).Figure 3
Stability features regarding thecritical (green) and thezero (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 r=10-4, α=10-4, ε=10-3, and m=0.5 and (a) s=1.2 and (b) s=1.6.
(a)
(b)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 thecritical 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 thezero regions, for this active parameter lying inside an interval (δ,1),δ=δ(γ,r,α,ε,σ,m,s)>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 s or m, a “0- equilibrium” region emerged inside thezero region (there exist no equilibria in the whole range (0,1)), getting larger as these two rational exponents (especially s) increase.


## 3. Analytical Formulae for a Parametricn-Dimensional System

### 3.1. Reduction to an(n+2)-Dimensional Coordinate Space

Now, we briefly discuss the basic features of a new formulation regarding a parametricn-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, (21)dxdt=fx;a,x∈Rn,a=a1,a2,a3T∈R3,f:Rn+3⟶Rn,f∈C∞(for SEIRS system (7) we have that (a1,a2,a3)=(γ,σ,β)), then expansion around the equilibrium path yields(22)dξdt=J0aξ+Fξ;a,ξ=x-x0a,with J0(a)=Dxf(x0(a);a) being the Jacobian matrix evaluated at this equilibrium path. The smooth vector function F:Rn+3→Rn represents the nonlinear terms; that is, F=Ox2 and(23)Fξ;a=∑m=2∞1m!Fmξ,…,ξ︸m;a,where(24)Fmu1,…,um;a=∑j1,…,jm=1n∂mfx;a∂xj1⋯∂xjmx=x0u1j1⋯umjm,m=2,3,…with ui=(ui1,…,uin)T∈Cn,i=1,…,m. Additionally, if a critical triplet a0=a10,a20,a30 exists, where J0 has a pair of pure imaginary eigenvalues λ1,2=±iω0, ω0>0, while the real part of the rest of the eigenvalues is negative, then a Hopf bifurcation occurs at a0, leading to a family of limit cycles. This implies that J0 has a pair of complex conjugate eigenvalues λa, λ-a in a region around a0, with(25)λa=μa+iωa,μa0=0,ωa0=ω0>0.We also note that the classic Hopf theory additionally demands that, at a0, λ(a) cross the imaginary axis with nonzero speed (transversality condition). Moreover we consider the normalized complex eigenvectors qa∈Cn and pa∈Cn of J0 and J0T, respectively, having the properties(26)J0q=λq,J0Tp=λ-p,p,q=∑i=1np-iqi=1,p,q-=0.Also, the two-dimensional parameter-dependent real eigenspace (corresponding to λ,λ-) (spanned by {Req(a),Imq(a)}), is denoted by T(a) (which becomes critical at a0: Tc=T(a0)), and the (n-2)-dimensional (n>2) real eigenspace, corresponding to all eigenvalues of J0 other than λ,λ-, is denoted by Tsu(a). Finally, the following Lemma ([6], Lemma 5.4) holds.Lemma 1.
y
∈
T
s
u if and only if 〈p,y〉=0.Then, based on the decomposition(27)x-x0=ξ=zq+z-q-+y,x∈Rn,z∈C1,y∈Rn∩Tsu,zq+z-q-∈T,by means of the third and fourth relation of (26), the last equation yields(28)z=p,ξ,y=ξ-p,ξq-p-,ξq-.Thus 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 (n+2)-dimensional coordinate space (z,y):(29)z˙=λaz+gz,z-,y;a,(30)y˙=J0ay+Fzq+z-q-+y;a-gz,z-,y;aq-g-z,z-,y;aq-,where the dot (·) denotes differentiation with respect to t and (31)gz,z-,y;a=p,Fzq+z-q-+y;a.Systems (22) and (29)-(30) are in fact dimensionally equivalent, due to the two real constraints imposed on y by means of Lemma 1.


### 3.2. A New Symbolic Representation of the System

By expressing the complex functiong(z,z-,y) in the form (32)g=12g002+16g003+⋯+g101+12g102+16g103+⋯z+g011+12g012+16g013+⋯z-+g2+g3+⋯,where(33)gmz,z-,y=∑k+l=m1k!l!gkl+gkl1+12gkl2+16gkl3+⋯zkz-l,gijm+y=p,Fi+j+m+q,…,q︸i,q-,…,q-︸j,y,…,y︸m+,gkl=gkl0with gij(m+)∈C, m=2,3,…,i+j=0,1,2,…,m+=0,1,2,…, then (29) and (30) become (34)z˙=12g002+16g003+⋯+λ+g101+12g102+16g103+⋯z+g011+12g012+16g013+⋯z-+g2+g3+⋯,(35)y˙=J0ay+12r002+16r003+⋯-12g002+16g003+⋯q-12g-002+16g-003+⋯q-+r101+12r102+16r103+⋯-g101+12g102+16g103+⋯q-g-011+12g-012+16g-013+⋯q-z+r011+12r012+16r013+⋯-g011+12g012+16g013+⋯q-g-101+12g-102+16g-103+⋯q-z-+r2+r3+⋯,where(36)rmz,z-,y=∑k+l=m1k!l!rkl+rkl1+12rkl2+16rkl3+⋯-gkl+gkl1+12gkl2+16gkl3+⋯q-g-lk+g-lk1+12g-lk2+16g-lk3+⋯q-zkz-l,rijm+y=Fi+j+m+q,…,q︸i,q-,…,q-︸j,y,…,y︸m+,rkl=rkl0with gij(m+) as in (33), rij(m+)∈Cn, m=2,3,…, i+j=0,1,2,…,m+=0,1,2,….Now, by introducing the center manifold expression(37)y=Vz,z-=∑m=2∞∑k+l=m1k!l!wklzkz-l,with wkl∈Cn, k+l=2,3,… and 〈p,wkl〉=0, and combining (37), as well as the invariance relation y˙=Vzz˙+Vz-z-˙, with (34) and (35), the analysis proceeds according to the steps figured out in [4, Subsection 2.3]. 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.


## 3.1. Reduction to an(n+2)-Dimensional Coordinate Space

Now, we briefly discuss the basic features of a new formulation regarding a parametricn-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, (21)dxdt=fx;a,x∈Rn,a=a1,a2,a3T∈R3,f:Rn+3⟶Rn,f∈C∞(for SEIRS system (7) we have that (a1,a2,a3)=(γ,σ,β)), then expansion around the equilibrium path yields(22)dξdt=J0aξ+Fξ;a,ξ=x-x0a,with J0(a)=Dxf(x0(a);a) being the Jacobian matrix evaluated at this equilibrium path. The smooth vector function F:Rn+3→Rn represents the nonlinear terms; that is, F=Ox2 and(23)Fξ;a=∑m=2∞1m!Fmξ,…,ξ︸m;a,where(24)Fmu1,…,um;a=∑j1,…,jm=1n∂mfx;a∂xj1⋯∂xjmx=x0u1j1⋯umjm,m=2,3,…with ui=(ui1,…,uin)T∈Cn,i=1,…,m. Additionally, if a critical triplet a0=a10,a20,a30 exists, where J0 has a pair of pure imaginary eigenvalues λ1,2=±iω0, ω0>0, while the real part of the rest of the eigenvalues is negative, then a Hopf bifurcation occurs at a0, leading to a family of limit cycles. This implies that J0 has a pair of complex conjugate eigenvalues λa, λ-a in a region around a0, with(25)λa=μa+iωa,μa0=0,ωa0=ω0>0.We also note that the classic Hopf theory additionally demands that, at a0, λ(a) cross the imaginary axis with nonzero speed (transversality condition). Moreover we consider the normalized complex eigenvectors qa∈Cn and pa∈Cn of J0 and J0T, respectively, having the properties(26)J0q=λq,J0Tp=λ-p,p,q=∑i=1np-iqi=1,p,q-=0.Also, the two-dimensional parameter-dependent real eigenspace (corresponding to λ,λ-) (spanned by {Req(a),Imq(a)}), is denoted by T(a) (which becomes critical at a0: Tc=T(a0)), and the (n-2)-dimensional (n>2) real eigenspace, corresponding to all eigenvalues of J0 other than λ,λ-, is denoted by Tsu(a). Finally, the following Lemma ([6], Lemma 5.4) holds.Lemma 1.
y
∈
T
s
u if and only if 〈p,y〉=0.Then, based on the decomposition(27)x-x0=ξ=zq+z-q-+y,x∈Rn,z∈C1,y∈Rn∩Tsu,zq+z-q-∈T,by means of the third and fourth relation of (26), the last equation yields(28)z=p,ξ,y=ξ-p,ξq-p-,ξq-.Thus 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 (n+2)-dimensional coordinate space (z,y):(29)z˙=λaz+gz,z-,y;a,(30)y˙=J0ay+Fzq+z-q-+y;a-gz,z-,y;aq-g-z,z-,y;aq-,where the dot (·) denotes differentiation with respect to t and (31)gz,z-,y;a=p,Fzq+z-q-+y;a.Systems (22) and (29)-(30) are in fact dimensionally equivalent, due to the two real constraints imposed on y by means of Lemma 1.


## 3.2. A New Symbolic Representation of the System

By expressing the complex functiong(z,z-,y) in the form (32)g=12g002+16g003+⋯+g101+12g102+16g103+⋯z+g011+12g012+16g013+⋯z-+g2+g3+⋯,where(33)gmz,z-,y=∑k+l=m1k!l!gkl+gkl1+12gkl2+16gkl3+⋯zkz-l,gijm+y=p,Fi+j+m+q,…,q︸i,q-,…,q-︸j,y,…,y︸m+,gkl=gkl0with gij(m+)∈C, m=2,3,…,i+j=0,1,2,…,m+=0,1,2,…, then (29) and (30) become (34)z˙=12g002+16g003+⋯+λ+g101+12g102+16g103+⋯z+g011+12g012+16g013+⋯z-+g2+g3+⋯,(35)y˙=J0ay+12r002+16r003+⋯-12g002+16g003+⋯q-12g-002+16g-003+⋯q-+r101+12r102+16r103+⋯-g101+12g102+16g103+⋯q-g-011+12g-012+16g-013+⋯q-z+r011+12r012+16r013+⋯-g011+12g012+16g013+⋯q-g-101+12g-102+16g-103+⋯q-z-+r2+r3+⋯,where(36)rmz,z-,y=∑k+l=m1k!l!rkl+rkl1+12rkl2+16rkl3+⋯-gkl+gkl1+12gkl2+16gkl3+⋯q-g-lk+g-lk1+12g-lk2+16g-lk3+⋯q-zkz-l,rijm+y=Fi+j+m+q,…,q︸i,q-,…,q-︸j,y,…,y︸m+,rkl=rkl0with gij(m+) as in (33), rij(m+)∈Cn, m=2,3,…, i+j=0,1,2,…,m+=0,1,2,….Now, by introducing the center manifold expression(37)y=Vz,z-=∑m=2∞∑k+l=m1k!l!wklzkz-l,with wkl∈Cn, k+l=2,3,… and 〈p,wkl〉=0, and combining (37), as well as the invariance relation y˙=Vzz˙+Vz-z-˙, with (34) and (35), the analysis proceeds according to the steps figured out in [4, Subsection 2.3]. 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 settingr=10-4, α=10-4, ε=10-3, m=0.5, s=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 wholecritical surfaceβcr=βcr(γ,σ). Thus stable limit cycles are generated through supercritical Hopf bifurcations, arising for (γ,σ) taking values in thecritical 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 l1 never becomes zero, the system undergoes solely acodimension 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).Figure 4
Bifurcation portrait as regards the supercritical Hopf bifurcation forr=10-4, α=10-4, ε=10-3, m=0.5, and s=1.2. Negative sign (blue region) of the first Lyapunov coefficient at thecritical surfaceβcr=βcr(γ,σ), in the parameter plane γ,σ (no βcr in the gray region).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 Figures5(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 βcr and also the period T 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 μi=eλi·T,i=1,2,3, with λi=(1/T)·ln(μi) being the respective exponents and T being the fundamental period of the cycle, computed as follows: μi={1,0.82,1.1×10-17} and λi={0,-4.7×10-4,-0.095}, 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.Figure 5
Stable cycle generated by supercritical Hopf bifurcation inxy(SI), wy(EI) planes and in 3D(SEI), respectively, for (r,α,ε,m,s,γ,σ,β)=(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): (x0,w0,y0)=(0.073948,0.001986,0.023535). Period: T=410.572533 days. ode113: Adams-Bashforth-Moulton PECE solver of orders 1 to 12 (Matlab).
(a)
(b)
(c)Figure 6
Stable cycles generated by supercritical Hopf bifurcation inxy(SI), wy(EI) planes and in 3D(SEI), respectively, for (r,ε,m,s,γ,σ,β)=(10-4,10-3,0.5,1.2,0.0548,0.24,0.46825), α∈[5⋅10-5,10-3], determined by a custom orthogonal collocation on finite elements algorithm. minβcr=0.51825 (α=5⋅10-5). Periods: Tmin=410.414990 days (α=5⋅10-5); Tmax=441.550487 days (α=10-3). ode113: Adams-Bashforth-Moulton PECE solver of orders 1 to 12 (Matlab).
(a)
(b)
(c)As regards the role and significance of the introduction of the additional epidemic-induced death rate,a, 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 a 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 s and the parameter c=βσ/[(σ+r)(γ+r)] (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 thefirst 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 correspondingcritical curves in the parameter plane (s,c).


## 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 thecritical andzero 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 thecritical region, while thezero 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.



---

*Source: 1467235-2018-01-17.xml*