Stability Analysis of an Age-Structured SEIRS Model with Time Delay

This paper is concerned with the stability of an age-structured susceptible–exposed– infective–recovered–susceptible (SEIRS) model with time delay. Firstly, the traveling wave solution of system can be obtained by using the method of characteristic. The existence and uniqueness of the continuous traveling wave solution is investigated under some hypotheses. Moreover, the age-structured SEIRS system is reduced to the nonlinear autonomous system of delay ODE using some insignificant simplifications. It is studied that the dimensionless indexes for the existence of one disease-free equilibrium point and one endemic equilibrium point of the model. Furthermore, the local stability for the disease-free equilibrium point and the endemic equilibrium point of the infection-induced disease model is established. Finally, some numerical simulations were carried out to illustrate our theoretical results.


Introduction
In recent years, the study of epidemiology has been a vital problem in ecology. The research of population dynamics has developed rapidly, and many mathematical models have been used to analyze various infectious diseases. Many results have been established in the stability analysis of different epidemic models. The first susceptible-infective-recovered (SIR) epidemic model about disease transmission was established by Kermack and McKendrick in 1927 [1]. Since then, the population dynamics of infectious diseases have attracted the attention of scientists. In 2012, the field of mathematical biology was expanded, particularly in the context of the spread of infectious diseases by Fred Brauer et al. [2]. Nowadays, there are some research work devoted to study the stability of steady states of the SIR, SIRS, SEIR, etc. models [3][4][5][6][7]. It is well known, in the spread of infectious diseases, some infective individuals of population are immune after being recovered (e.g., measles, smallpox, mumps, and others). Meanwhile, some recovered individuals have no immunity (e.g., AIDS, hyperthyroidism, lupus erythematosus, and others), who will return to the susceptible population and continue to be infected. In fact, the probability of becoming infected is different among different individuals, which may depend on the type of infectious diseases and the status of individuals. Therefore, it is necessary to discuss the SEIRS model, which can more clearly describe the spread of infectious diseases in real life.
Time delay is ubiquitous and can be applied in many epidemiology related studies [8,9]. For example, measles has an incubation period of 8-13 days and the incubation period of canine madness is a few months or several years after infection. Sharma et al. [10] developed a five compartmental infection model to describe the spread of avian influenza A (H7N9) virus with two discrete time delays. In addition, Xu et al. [11] analyzed the stability of a SIRS model with time delay. Similarly, Shu [12] and

An Age-Structured SEIRS Model with Time Delay
Motivated by the referred works [28], we discuss two age stages of each subpopulation of an age-structured SEIRS model with time delay in this paper, which include immature stage and mature stage. In the immature stage, individuals can be born, grow up, die, or survive until the maximum age a 1 , but, in this case, the individuals are not proliferate until the maximum age a 1 . The age a 1 is considered to be the maximum age in the immature stage and it is also considered to be the initial age in the mature stage. In the mature stage, individuals have reached maturity (a > a 1 ) and can grow up, proliferate, die, or survive to the maximum age A. Here, we consider susceptible, exposed, infective, and recovered individuals of two age stages in an age-structured SEIRS model. Let S(a, t), E(a, t), I(a, t), and R(a, t) be the distribution densities of susceptible, exposed, infective ,and recovered individuals of age x at time t, and they take place in the domain Ω = {(a, t) : 0 ≤ a ≤ A, 0 ≤ t ≤ T}. The integrals N s (t) = A 0 S(a, t)da, N e (t) = A 0 E(a, t)da, N i (t) = A 0 I(a, t)da, and N r (t) = A 0 R(a, t)da are considered as the number of susceptible, exposed, infective and recovered individuals, respectively. The total population N(t) is: N(t) = N s (t) + N e (t) + N i (t) + N r (t). What needs to be added is that S(a, t), E(a, t), I(a, t), and R(a, t) should belong to L 1 (Ω) because we assume that the initial total population is limited. Then, the following differential equations with time delay are

∂I ∂t
∂R ∂t with initial conditions and boundary conditions q(a, t)θ 3 (a, t)I(a, t)da, t ∈ (0, T), whereθ i (a, t)(i = 1, 2, 3, 4) denote fertility rates of females of each subpopulation of age a;θ i (a, t) = θ i (a, t), if a ∈ [a 1 , A];θ i (a, t) = 0, if a / ∈ [a 1 , A];α i (a, t) are natural death rates of each subpopulation of age a; σ is the provided parameter: σ = 1 if the individuals die when they produce and σ = 0 if the individuals continue to survive when they produce; γ 1 (a, t) is a transmission coefficient which describes the varying probability of infectiousness and it is related to a great many social, environmental, and epidemiological factors; β(a, a , t) is the contact rate between infected population (age a ) and susceptible population (age a) per unit time; τ (τ > 0) is a fixed incubation period of infection; γ 2 (a, t) is the conversion rate from exposed population to infected population of age a; γ 3 (a, t) is the recovery rate of age a; ρ(a, t) is the conversion rate from recovered population losing immunity to the susceptible population of age a; µ i (i = 1, 2, 3, 4) are reproductive rates of each subpopulation in the proliferating stage; p(a, t) is part of exposed individuals which procreate new exposed individuals of age a; and, similarly, q(a, t) is part of infective individuals which procreate new infective individuals of age a.

Traveling Wave Solution
In this section, we mainly use the method of characteristics [2,18,[25][26][27][28] for first-order hyperbolic partial differential equations (Equation (1)). Then, the system in Equation (1) is reduced to nonlinear delayed integro-differential equations along the characteristics curve a − t = constant [25,26]. Let u = a − t, the following system with time delay is obtained: In general, we divide time Then, Ω can be grouped into two sets: We also need to define an auxiliary set for k = 1, ..., K: denotes a whole part of real number x. The following hypotheses are given: (H 1 ): S 0 (a), E 0 (a) and I 0 (a, 0) are non-negative and continuous when a ∈ [ 0, A]; when a → A − 0, S 0 (A) = 0, E 0 (A) = 0 and I 0 (A, 0) = 0.
For convenience, we assume Hypotheses (H 1 )-(H 4 ) are satisfied. Then, we solve the system in Equation (4) by using the well-known steps method [8,9]. For the first step h = 1, i.e., for the first time interval t ∈ [0, τ], the solution to the initial value problem in Equation (4) can be obtained according to the known values of function I(x, t − τ). Repeating the pervious step for h = 2, 3, 4, ..., H, for the time interval t ∈ [(h − 1)τ, hτ], we can obtain the solution to the initial value problem in Equation (4) for the whole time interval t ∈ [0, T]: where The exact solution to Equations (11)-(14) of the system in Equation (4) can be expressed in terms of the original variables. Let k = 1 and u = a − t; we have the solution of the system in Equation (1) in the sets Ω where W 1 (u, t 0 , t, τ), W 2 (u, t 0 , t), W 3 (u, t 0 , t), and W 4 (u, t 0 , t) are shown in Equations (11)-(14), respectively. F (k) (u), G (k) (u), and Q (k) (u) are given by defining functions F n (u) can be defined according to the recurrent algorithm as follows: Proof. By analogy with [28], if S 0 (a), E 0 (a), I 0 (a, 0), and R 0 (a) satisfy compatibility conditions (i.e., Hypothesis H 4 ), then two parts of the traveling wave in Equation (15) can be combined continuously. If the parameters of the system in Equation (1)

Stability Analysis of System
Consider the nonlinear autonomous system in Equation (1) where the following parameters are constants: θ i (a, t) = θ i0 (i = 1, 2, 3, 4), α i (a, t) = α i0 = α i + σθ i0 (i = 1, 2, 3, 4), γ 1 (a, t) = γ 10 , γ 2 (a, t) = γ 20 , γ 3 (a, t) = γ 30 , β(a, a , t) = β 0 , p(a, t) = p 0 , q(a, t) = q 0 , and ρ(a, t) = ρ 0 , where α i0 , θ i0 , γ i0 , β 0 , p 0 , q 0 , and ρ 0 are positive constants. In this paper, a partial differential equation (Equation (1)) with the initial-boundary values in Equations (2) and (3) reduced to a nonlinear ordinary differential equation. Take the maturity age a 1 → 0, which is not an essential simplification. Integrating Equation (1) in regard to age a from 0 to A, and using the real conditions S(A, t) = 0, E(A, t) = 0, I(A, t) = 0 and R(A, t) = 0, the initial value problem of the nonlinear ordinary differential equation autonomous system that describes the population dynamics of the number of population N s (t), N e (t), N i (t), and N r (t) is: The basic reproduction number [31] is defined as According to the analysis of [28], if R 0 > 1, it can lead to the outbreak of infectious diseases. Hence, the following theorem is obtained. Theorem 2. If γ 10 = 0, the system in Equation (20) has only a disease-free equilibrium point H 0 =(0, 0, 0, 0); if R 0 > 1 and the parameters of the system in Equation (20)  Proof. When γ 10 = 0, that is no infectious diseases, there exists one disease-free equilibrium point H 0 =(0, 0, 0, 0). When the fertility rate and the death rate of susceptible population are satisfied by where R 1 is a dimensionless index for the existence of the disease-free equilibrium point. It can be seen that R 1 > 1 presents that death rate of susceptible population is less than their reproductive rate. We can take it further into consideration with the dynamic behavior of susceptible population. The endemic equilibrium point On the basis of the parameters of the system in Equation (1) where R 2 , R 3 , and R 4 are dimensionless indexes for the existence of the endemic equilibrium point H * . It can be known that R 2 < 1 presents that the death rate and conversion rate of exposed population outweigh their reproductive rate. The density of exposes cannot increase indefinitely because of the balance of death rate, conversion rate, and reproductive rate. From the biological point of view, higher death rate and conversion rate are the consequence of infectious disease. R 3 < 1 denotes that reproductive rate of infected population is less than their death rate and conversion rate.
First, the local stability of the disease-free equilibrium point H 0 is analyzed.
Next, we study the stability of the endemic equilibrium point H * by linearizing the nonlinear autonomous system and calculating characteristic equations. Accordingly, we can obtain polynomial equations which can analyze the stability of system. After linearization of the system in Equation (20), we get Consider the form of the exponential solution of the linearized system in Equation (32) as where λ is a parameter. Substituting these solutions into Equation (32), we have Denote M 1 = µ 1 θ 10 − α 10 where Theorem 4. For the system in Equation (20) is locally asymptotically stable when τ = 0.
Proof. When τ = 0, the characteristic equation (Equation (34)) becomes It is evident that D 2 > 0, D 3 > 0 and D 4 > 0 when N * s , N * e , N * i , N * r satisfy the nonnegativity conditions in Equation (23)- (30). The following equations are also satisfied: Then, according to the Routh-Hurwitz criterion, the roots of the characteristic equations of system are all negative real. Thus, the endemic equilibrium point of the autonomous system in Equation (20) is asymptotically stable with τ = 0.

Theorem 5. For the system in Equation
the results can be obtained: (i) When 0 < τ < τ 0 , the endemic equilibrium point H * =(N * s , N * e , N * i , N * r ) is locally asymptotically stable.
(ii) When τ > τ 0 , the endemic equilibrium point Proof. When τ > 0, suppose that Equation (34) has a purely imaginary root λ = iω(ω > 0). The characteristic equation (Equation (34)) converses to the following form: Separating the real and imaginary parts yields two corresponding equations: Thus, Let x = ω 2 . Squaring Equations (37) and (38) and adding them together, we can get According to the previous restrictions in Equation (23)-(30), we can also get Suppose that In accordance with Descartes' rule of signs, Equation (41) has only one changed sign. Thus, it only has one positive root. Let x * be the small unique positive root and it always exists. The unknown parameter ω of Equations (37) and (38) is defined as ±iω 0 = ±i √ x 0 . Combining with Equations (37) and (38), the form of time delay τ k is gained: Therefore, when τ ∈(0, τ 0 ), all roots of Equation (34) have strictly negative real parts. The endemic equilibrium point H * =(N * s , N * e , N * i , N * r ) of the system in Equation (20) is locally asymptotically stable. When τ = τ 0 , the roots of Equation (34) have strictly negative real parts except for ±iω 0 . When τ > τ 0 , the endemic equilibrium point H * = (N * s , N * e , N * i , N * r ) of the system in Equation (20) is unstable. Then, differentiating both sides of Equation (34) with respect to τ, we obtain Further, one leads to where Since conditions in Equation (42) where Hence, based on the properties of the Hopf bifurcation discussed in [32], the transversal condition holds and a Hopf bifurcation occurs at τ = τ 0 . The proof is complete.

Simulation
To affirm the stability analysis above, we numerically simulated the disease-free equilibrium point H 0 and the endemic equilibrium point H * of the system in Equation (20). We are more concerned with the indicators and conditions for outbreaks of "local" population or local asymptotic stability of equilibrium points. All programs were developed using the software Matlab R2016a (see Program S1).

Concluding Remarks
In this article, the theory of an age-structured SEIRS model with time delay is analyzed. The model is based on the delayed nonlinear partial differential equation of initial-boundary value problems. The traveling wave solution of the system in Equation (1) is obtained using the method of characteristic and the recurrent algorithm. Then, we can obtain the existence and uniqueness of the continuous traveling wave solution of system according to hypotheses. The age-structured SEIRS model with time delay is reduced to a nonlinear ordinary differential equation under some insufficient simplifications. This allows us to obtain some sufficient conditions of existence of two equilibrium points of an age-structured SEIRS system: R 1 is a dimensionless index for the existence of the disease-free equilibrium point H 0 . R 2 , R 3 , and R 4 are dimensionless indexes for the existence of the endemic equilibrium point H * . From the biological point of view, the endemic equilibrium point H * only exists in the case of high values of death and conversion rate of exposed and infected population. The disease-free equilibrium point H 0 and the endemic equilibrium point H * are given. The disease-free equilibrium point H 0 = (0, 0, 0, 0) is locally asymptotically stable if R 1 < 1, R 2 < 1, and R 3 < 1. The stability of the endemic equilibrium point H * = (N * s , N * e , N * i , N * r ) with τ = 0 and τ > 0 are analyzed: for R 0 > 1, if the condition in Equation (36) holds, the endemic equilibrium point H * is locally asymptotically stable when τ = 0; if the conditions in Equations (43) and (44) hold, the endemic equilibrium point H * is locally asymptotically stable when time delay 0 < τ < τ 0 ; if the conditions in Equations (43)-(44) hold, the endemic equilibrium point H * is unstable when τ satisfies τ > τ 0 ; Hopf bifurcation occurs at τ = τ k (k = 0, 1, 2, ...). When time delay exceeds the critical value τ 0 , the system in Equation (20) loses its stability and Hopf bifurcation occurs. In this case, the susceptible, exposed, infected, and recovered population in the model will coexist in an oscillating mode and infectious diseases will get out of control. It can be seen that time delay has an important effect on the spread of infectious diseases. Therefore, we should shorten the time delay as much as possible in order to predict and eliminate infectious diseases. Our further research is using some bifurcation control strategies to control the occurrence of the Hopf bifurcation so as to control the occurrence of infectious diseases.
In general, this study provides the practical understanding of the different dynamic behaviors of an age-structured susceptible-exposed-infected-recovered-susceptible model with time delay, which is helpful for us to understand the application of epidemiology better in real life.