Bifurcation analysis of an age-structured SIRI epidemic model

Abstract: In this paper, an SIRI epidemic model with age of infection and the proliferation of susceptible individuals with logistic growth is investigated. By using the theory of integral semigroup and Hopf bifurcation theory for semilinear equations with non-dense domain, it is shown that if the threshold parameter is greater than unity, sufficient condition is derived for the occurrence of the Hopf bifurcation. Numerical simulations are carried out to illustrate the theoretical results.


Introduction
In the past decades, great attention has been paid by many researchers to SIR epidemic models [1][2][3][4][5][6][7], in which total host population is divided into three classes called susceptible (S ), infective (I) and removed (R), and the immunity that is obtained upon recovery is assumed to be permanent. In [8], for herpes viral infections, considering the fact of recovered individuals may relapse with reactivation of latent infection and reenter the infective class, Tudor proposed the following SIRI epidemic model: where S (t), I(t) and R(t) represent the number of susceptible individuals, infectious individuals and recovered individuals at time t, respectively. Assumptions made in the system (1.1) are homogeneous mixing and constant population size. The parameter A is the constant birth rate, µ is the natural death rate, β is the contact rate, and γ is the recovery rate from the infective class. It is assumed that an individual in the recovered class can revert to the infective class with a constant rate δ. Here δ > 0 implies that the recovered individuals would lose the immunity, and δ = 0 implies that the recovered individuals acquire permanent immunity.
We note that in system (1.1), the total population size is assumed to be constant. In reality, demographic features which allow the population size to vary should be incorporated into epidemiological models in some cases. For a relatively long-lasting disease or a disease with high death rate, the assumption of logistic growth input of the susceptible seems more reasonable [9]. In [10], Gao and Hethcote investigated a SIRS model with the standard incidence rate, and considered a demographic structure with density-dependent restricted population growth by the logistic equation. In [11], Li et al. studied a SIR epidemic model with logistic growth and saturated treatment, and the existence of the stable limit cycles was obtained. Rencently, there are growing interests in epidemiological models with demographic structures of logistic type [12][13][14][15][16].
It is worth noting that in the above models, the transmission coefficient is assumed to be constant, and the infected person has the same infectivity during their periodic infection. However, laboratory studies suggest that the infectivity of infectious individuals be different at the differential age of infection [17,18]. Further, it was reported in [19,20] that, age-structure of a population is an important factor which affects the dynamics of disease transmission. In [21], Magal et al. discussed an infection-age model of disease transmission, where both the infectiousness and the removal rate depend on the infection age. In [22], an age structured SIRS epidemic model with age of recovery is studied, and the existence of a local Hopf bifurcation is proved under certain conditions. In [23], Chen et al. considered an SIR epidemic model with infection age and saturated incidence, and established a threshold dynamics by applying the fluctuation lemma and Lyapunov functional.
Motivated by the works of Chen et al. [23], Gao and Hethcote [10], and Tudor [8], in the present paper, we are concerned with the effects of logistic growth and age of infection on the transmission dynamics of infectious diseases. To this end, we consider the following differential equation system: with the boundary condition and the initial condition is the set of all integrable functions from (0, ∞) into R + = [0, ∞). In system (1.2), S (t) represents the number of susceptible individuals at time t, i(a, t) represents the density of infected individuals with infection age a at time t, and R(t) is the number of individuals who have been infected and temporarily recovered at time t. All parameters in system (1.2) are positive constants, and their definitions are listed in Table 1. The rate of natural death γ The recovery rate of the infective individuals δ The rate at which recovered individuals lose immunity and return to the infective class a Age of infection, i.e., the time that has lapsed since the individual became infected β(a) Transmission coefficient of the infected individuals at age of infection a In the sequel, we further make the following assumption: For convenience, we assume that where e −(µ+γ)a is the probability of infected individual to survive to age a and τ > 0, β * > 0.
The organization of this paper is as follows. In Section 2, we formulate system (1.2) as an abstract non-densely defined Cauchy problem. In Section 3, we study the existence of feasible equilibria of system (1.2). In Section 4, the linearized equation and the characteristic equation of system (1.2) at the interior equilibrium are investigated, respectively. In Section 5, by analyzing corresponding characteristic equation, we discuss the existence of Hopf bifurcation. In Section 6, numerical examples are carried out to illustrate the theoretical results, and sensitivity analysis on several important parameters is carried out.
Accordingly, system (2.1) is equivalent to the following system: We now consider the following Banach space Define the linear operator L τ : D(L τ ) → X by .
Therefore, the linear operator L τ is non-densely defined due to Letting , system (2.3) is transformed into the following non-densely defined abstract Based on the results in [24] and [25], the global existence, uniqueness and positivity of solutions of system (2.4) are obtained.

Existence of feasible equilibria
In this section, we study the existence of feasible equilibria of system (2.4). Define the threshold parameter R 0 by ∈ X 0 is a equilibrium of system (2.4). Then we have By direct calculation, we obtain We derive from the first and second equations of Eq (3.1) that This, together with (3.4), yieldsR On substituting (3.6) into (3.3), we obtain that We therefore follows from (3.6) and (3.8) that Based on the discussions above, we have the following result.
In addition, if R 0 > 1, there exists a unique positive equilibrium Correspondingly, for system (1.2), we have the following result.

Characteristic equation
In this section, we investivate the linearized equation of (2.4) around the positive equilibriumx * (a), and the characteristic equation of (2.4) atx * (a), respectively.
Making the change of variable y(t) Accordingly, the linearized equation of system (4.1) around the origin is After that, system (4.1) can be rewritten as where the linear operator L := L τ + τDF(x * (a)) and F (y(t)) = τF(y(t) +x * (a)) − τF(x * (a)) − τDF(x * (a))y(t) satisfying In the following, we give the characteristic equation of (2.4) at the positive equilibrium. By means of the method used in [26], we obtain the following lemma.
Hence, we have the following result.
Using a similar argument as in [27], it is easy to verify the following result.  with ∆(λ) and K(λ, ψ) given by (4.6) and (4.7).

Existence of Hopf bifurcation
In this section, by applying Hopf bifurcation theory [26], we are concerned with the existence of Hopf bifurcation for the Cauchy problem (2.4) by regarding τ as the bifurcation parameter.
without loss of generality, we may assume that Eq (5.6) has two positive roots denoted respectively as z 1 and z 2 . Then Eq (5.4) has two positive roots ω k = √ z k (k = 1, 2). Further, from (5.3), we have where k = 1, 2; j = 0, 1, · · · , and Based on the above discussion, we have the following result.
Therefore ζ = iω k is a simple root of Eq (5.1).
Noting that when τ = 0, the equilibrium x * (a) of (2.4) is locally asymptotically stable if R 0 > 1, from what has been discussed above, we have the following results.

Numerical simulations
In this section, numerical simulations will be given to illustrate the theoretical results in Section 5. Further, sensitivity analysis is used to quantify the range of variability in threshold parameter and to identify the key factors giving rise to threshold parameter, which is helpful to design treatment strategies.

Dynamics of system (1.2)
In this section, we give a numerical example to illustrate the main results in Section 5. The intrinsic growth rate of susceptible populations r (per year) 0.1 Assumed The rate of natural death µ (per year) 1/70 [29,30] The rate at which recovered individuals lose immunity and return to the infective class δ (per year) 0.01 Assumed The recovery rate of the infective individuals γ (per year) 0.3 [31]  Based on the research works of tuberculosis [29,30,32], parameter values of system (1.2) are summarized in Table 2 In this case, system (1.2) has an endemic steady state E * (0.6070, 0.0879τe −0.3143τa , 3.4534). By calculation, we obtain ω 0 = 0.0565 and τ 0 = 14.3556. By Theorem 5.1, we see that the endemic steady state E * is locally asymptotically stable if τ ∈ [0, τ 0 ) and is unstable if τ > τ 0 . Further, system (1.2) undergoes a Hopf bifurcation at E * when τ = τ 0 . Numerical simulation illustrates the result above (see, Figures 1 and 2). Remark. From Figure 1, we see that when the bifurcation parameter τ is less than the critical value τ 0 , the endemic steady state E * of system (1.2) is locally asymptotically stable. From Figure 2, we observe that E * losses its stability and Hopf bifurcation occurs when τ crosses τ 0 to the right (τ > τ 0 ). This implies that the age, i.e., infection period τ is the key factor that causes the endemic steady state E * to become unstable and the appearance of Hopf bifurcation.

Sensitivity analysis
Sensitivity analysis is used to quantify the range of variables in threshold parameter and identify the key factors giving rise to threshold parameter. In [33,34], Latin hypercube sampling (LHS) is found to be a more efficient statistical sampling technique which has been introduced to the field of disease modelling. LHS allows an un-biased estimate of the threshold parameter, with the advantage that it requires fewer samples than simple random sampling to achieve the same accuracy.
By analysis of the sample derived from Latin hypercube sampling, we can obtain large efficient data in respect to different parameters of R 0 . Figure 3 shows the scatter plots of R 0 in respect to K, µ, δ and γ, which implies that δ is a positively correlative variable with R 0 , while µ is a negatively correlative variable. But the correlation between K, γ and R 0 is not clear. In [34], Marino et al. mentioned that Partial Rank Correlation Coefficients (PRCCs) provide a measure of the strength of a linear association between the parameters and the threshold parameter. Furthermore, PRCCs are useful for identifying the most important parameters. The positive or negative of PRCCs respectively denote the positive or negative correlation with the threshold parameter, and the sizes of PRCCs measure the strength of the correlation. As can been seen in Figure 4, K, δ and γ are positively correlative variables with R 0 while µ is negatively correlative variables. By selecting different parameter values, we can explore the influence of the parameters µ and δ on the numbers of infected individuals at time t, which is denoted as I(t). As shown in Figure 5, increasing the natural death rate µ and decreasing the rate at which recovered individuals return to the infective class will have a positive impact on I(t) to some extent, which means that the influence of µ and δ on I(t) is consistent with that on R 0 .