Dynamic analysis of the recurrent epidemic model

1.   Introduction

In the real world, many diseases, such as Influenza, hand-foot-mouth-disease, rotavirus, and so on, can occur the secondary infection after recovery ^[1,2,3]. We can classify such diseases as recurrent epidemics. The study of recurrent epidemics which has been attracting great attention for decades, is currently one of the hottest topics in the field of epidemiology. A great deal of significant progress has been achieved so far (for examples, see ^{[4,5,6,7,8,9,10,11,12]}), including derivation of the threshold for outbreak and persistence of recurrent infectious disease, and the dynamic analysis of the seasonal and non-seasonal behaviors.

At present, most of the existing deterministic models present periodic behaviors of recurrent infectious diseases by selecting seasonal parameters. Whether there are the deterministic models with non-seasonal parameters which can describe the possible periodic behaviors of recurrent infectious disease is worth exploring. Since the most important feature of recurrent infectious diseases is that the recovered individuals' immunity are temporary, the immunity age of the recovered individuals and the time required for loss of immunity must be characterized in the model. That is, the model should be built by combining immunity age and delayed immune loss.

We assume that recovered individuals lose their immunity after a period of recovery, and then, they become susceptible again. Let $m(a)$ be the immunity loss rate function of the recovered individual with immunity age $a$, and $\tau > 0$ be the minimal time that a recovered individual has immunity. If the immunity age $a$ is less than $\tau$, then the recovered individuals remain in the recovered class, which means $m(a) = 0$. If the immunity age $a$ is great than $\tau$, then the recovered individuals may lose their immunity and enter into the susceptible class at the rate $m(a) = m_* > 0$. Therefore, $m(a)$ follows the following function:

Apparently, $m(a)\in L^\infty_+((0, +\infty), \mathbb{R})$.

In this paper, we mainly consider the recurrent infectious diseases with temporary immunity. Therefore, we incorporate loss of immunity into epidemic model by establishing an SIRS model with age structure and delay. Denote $S(t)$ and $I(t)$ respectively as the number of susceptible and infectious individuals, and denote $R(t, a)$ as the number of recovered individuals with immunity age $a$ at time $t$. The model is given as follows:

with the initial condition

where $\Lambda$ is the recruitment rate of susceptible individuals, $\beta$ is the transmission rate, $\dfrac{1}{\alpha}$ is the half-saturation constant of infectious individuals, $d$ is the natural death rate, $\mu$ is the mortality rate due to disease, and $\gamma$ is the recovery rate. All parameters are positive and constant.

This paper is organized as follows. Some preliminary results and the well-posedness of system (1.2) are presented in Section 2. In Section 3, we prove the existence of equilibria, especially the existence and uniqueness of positive equilibrium, and linearize system (1.2) around the equilibrium $E_*$. In Section 4, we discuss the stability of both equilibrium, and analyze the existence of Hopf bifurcations when the stability of $E_*$ changes. We conclude and discuss our findings in Section 5 with some numerical examples given to illustrate our theoretical results.

2.   Preliminaries and Well-posedness

In this section, we will rewrite system (1.2) into an abstract equation on a suitable Banach lattice, establish the well-posedness result for the system and prove the nonnegativity and boundedness of solutions. At first, we collect some preliminaries on linear operators and $C_0$-semigroup theory and some notations to be used in this paper.

Let $L:D(L)\subset X\rightarrow X$ be a linear operator on a Banach space $X$. Denote the resolvent set of $L$ as $\rho(L)$. The spectrum of $L$ is $\sigma(L) = \mathbb{C}\; n\; \rho(L)$. The point spectrum of $L$ is the set

Definition 2.1. Let $L: D(L)\subseteq X \rightarrow X$ be a linear operator. If there exist real constants $M\geq1$ and $\omega\in\mathbb{R}$, such that $(\omega, +\infty)\subseteq\rho(L)$, and

Then the linear operator $(L, D(L))$ is called a Hille-Yosida operator.

For a Hille-Yosida operator one has the following perturbation result.

Lemma 2.1 (see ^[13,14]). Let $(A, D(A))$ be a Hille-Yosida operator on a Banach space X and $B\in \mathscr{L}(X)$, $\mathscr{L}(X)$ denotes the set of all bounded linear operators on $X$, then the sum $C = A+B$ is a Hille-Yosida operator as well.

If $(L, D(L))$ is a Hille-Yosida operator on the Banach space $X$ and set

then the operator $(L_0, D(L_0))$ is called the part of $L$ in $X_0$ and we have that:

Lemma 2.2 (see ^[13,14]). If $(L, D(L))$ is a Hille-Yosida operator, then its part $(L_0, D(L_0))$ generates a $C_0$-semigroup $(T_0(t))_{t\geqslant 0}$ on $X_0$.

Now we set about to rewrite system (1.2) into an abstract evolution equation. Let

Define the linear operator $\mathcal{A}:D(\mathcal{A})\subseteq X\longrightarrow X$ by

with $D(\mathcal{A}) = \mathbb{R}^2\times \{0\}\times W^{1, 1}((0, +\infty), \mathbb{R})$. Then $\overline{D(\mathcal{A})} = \mathbb{R}^2\times\{0\}\times L^{1}((0, +\infty), \mathbb{R})$, which shows $D(\mathcal{A})$ is not dense in $X$. We also introduce a nonlinear map $\mathcal{F}: \overline{D(\mathcal{A})}\longrightarrow X$ given by

and let

Then we can reformulate system (1.2) as the following abstract Cauchy problem:

where $u_0 = (S_0, I_{0}, 0, R_0(a))^T$.

In general, it is difficult to find a strong solution for an abstract differential equation like (2.1). So, we solve (2.1) in integrated form

Set

We will show in Theorem 3.2 in the next section that $(\mathcal{A}, D(\mathcal{A}))$ is a Hille-Yosida operator and hence from Lemma 2.2 it generates a $C_0$-semigroup on the closure of its domain. As a result, we have the following well-posedness theorem for system (2.1).

Theorem 2.1. For any $u_0\in X_{0+}$, the system (1.2) represented by the integral equation (2.2) has a unique continuous solution with values in $X_{0+}$. Moreover, the map $\Phi: [0, +\infty) \times X_{0+}\mapsto X_{0+}$ defined by $\Phi(t, u_0) = u(t, u_0)$ is a continuous semi-flow, i.e. the map is continuous and satisfies that $\Phi(0, \cdot) = I$ and $\Phi(t, \Phi(s, \cdot)) = \Phi(t+s, \cdot)$.

Due to the biological interpretation of system (1.2), only non-negative solutions are meaningful to be considered. The following result reveals that the solutions $(S(t), I(t), R(t, a))$ of system (1.2) with non-negative initial value remain non-negative and bounded ultimately.

Theorem 2.2. All solutions of system (1.2) with non-negative initial value remain non-negative for all $t\geq 0$ and are ultimately bounded.

Proof. By using the second equation of system (1.2), we have

Integrating the third equation of system (1.2) along the characteristic line yields that

It is clear that $R(t, a)$ remains nonnegative for all nonnegative initial values.

In the following, we prove that $S(t)$ is nonnegative for $t\geq 0$. In fact, if there exists $t_1 > 0$ such that $S(t_1) = 0$, and $S(t) > 0$ for $\forall t\in(0, t_1)$, then by the first equation of system (1.2), we have $S'(t_1) = \Lambda+\int_{0}^{+\infty}m(a)R(t_1, a)da > 0$. It implies that $S(t)\geq0$ for all $t\geq0$.

Summarizing the above analysis, we know that any solution of system (1.2) with non-negative initial data remains nonnegative for all $t\geq 0$.

Next, we will show that the solutions of system (1.2) are ultimately bounded. Let $\bar{R}(t) = \int_0^{+\infty}R(t, a)da$, which represents the total number of recovered individuals at time t. Biologically, there exists a finite maximum age, so it is reasonable to assume that $\lim_{a\rightarrow+\infty}R(t, a) = 0$. Then from system (1.2), we have

Therefore,

It follows that the omega limit set of system (1.2) is contained in the following bounded feasible region:

Obviously, this region is positively invariant with respect to system (1.2). It implies that the system (1.2) is ultimately bounded.

3.   Equilibria and linearized equations

In this section, we devote to linearizing the nonlinear system (1.2) around the equilibrium solutions. For this purpose, we firstly discuss the existence of equilibria of system (1.2).

It is clear that system (1.2) always has a disease free equilibrium $E_0 = (S^{(0)}, 0, 0)$ with $S^{(0)} = \dfrac{\Lambda}{d}$. In order to find the positive equilibrium $E_\ast = (S_*, I_*, R_*(a))$ of system (1.2), we have

Solving the third equation of (3.1), we get

By using the second equation of (3.1), we have

Substituting $R_\ast(a)$ and $S_\ast$ into the first equation of (3.1), we obtain

Let $R_0 = \dfrac{\beta S^{(0)}}{d+\mu+\gamma} = \dfrac{\Lambda\beta}{d(d+\mu+\gamma)}$. Then, $I_*$ can be rewritten as

which is positive if $R_0 > 1$. Hence, system (1.2) has a unique positive equilibrium when $R_0 > 1$.

Summarizing the above analysis, we have the following result.

Theorem 3.1. System (1.2) always has a disease free equilibrium $E_0 = (S^{(0)}, 0, 0)$. If $R_0 > 1$, there also exists a unique endemic equilibrium $E_* = (S_*, I_*, R_*(a))$ for fixed $a$.

In fact, each term in $R_0$ has clear epidemiological interpretation. $\dfrac{1}{d+\mu+\gamma}$ is the average infection period. $\beta$ denotes the transmission rate of an infectious individual. $S^{(0)}$ is the total number of susceptible individuals. Therefore, $R_0$ represents average new cases generated by a typical infectious member in the entire infection period. That is, $R_0$ is the basic reproduction number of system (1.2).

Let $S(t) = x(t)+\overline{S}$, $I(t) = y(t)+\overline{I}$, $R(t, a) = z(t, a)+\overline{R}(a)$, where $\overline{E} = (\overline{S}, \overline{I}, \overline{R}(a))$ is a steady state of system (1.2), and let $\tilde{u}(t) = (x(t), y(t), 0, z(t, a))$, $\bar{u} = (\overline{S}, \overline{I}, 0, \overline{R}(a))$. Then, system (2.1) is equivalent to the following Cauchy problem

By conducting direct computations one can obtain readily the linearized system of (2.1) around $\bar{u}$ as the following form

in which

Clearly, $D\mathcal{F}(\bar{u})$ is a compact bounded linear operator on $X$.

Denote $\Omega = \{\lambda\in\mathbb{C}:\; \mbox{Re}(\lambda) > -d\}$. We claim then prove the following statement.

Theorem 3.2. The operator $(\mathcal{A}, D(\mathcal{A}))$ is a Hille-Yosida operator.

Proof. For $(\phi, \varphi, \psi, \omega)\in X$, $(\tilde{\phi}, \tilde{\varphi}, 0, \tilde{\omega})\in D(\mathcal{A})$, $\lambda \in \Omega$, we have

It then follows that

Integrating the last equation of (3.4) with regard to the age variable $a$ and adding all the equations, we obtain that

Thus, we have

which shows that $(\mathcal {A}, D(\mathcal {A}))$ is a Hille-Yosida operator.

By Lemma 2.1 and Theorem 3.2, it follows that

Theorem 3.3. The operator $\mathcal{A}+D\mathcal{F}(\bar{u})$ is a Hille-Yosida operator.

Using Lemma 2.2, we further derive that

heorem 3.4. The part of $(\mathcal{A}, D(\mathcal{A}))$ and $\left(\mathcal{A}+D\mathcal{F}(\bar{u}), D\left(\mathcal{A}+D\mathcal{F}(\bar{u}\right))\right)$ generate $C_0$-semigroups $(\mathscr{S}(t))_{t\geq0}$ and $(\mathscr{T}(t))_{t\geq0}$, respectively, on space $X_0$.

In order to establish the stability results for system (1.2), we will analyze the compactness of the generated $C_0$-semigroups. Firstly, we introduce the definition of quasi-compactness for a semigroup below.

Definition 3.1 (cf. ^[15]). A $C_0$-semigroup $\left(\mathcal{T}(t)\right)_{t\ge 0}$ is called quasi-compact if $\mathcal{T}(t) = \mathcal{T}_1(t)+\mathcal{T}_2(t)$ with the operator families $\mathcal{T}_1(t)$ and $\mathcal{T}_2(t)$ satisfying that

(i) $\mathcal{T}_1(t)\rightarrow0$, as $t\rightarrow+\infty$,

(ii) $\mathcal{T}_2(t)$ is eventually compact, that is, there is $t_0 > 0$, such that $\mathcal{T}_2(t)$ is compact for all $t > t_0$.

For a quasi-compact $C_0$-semigroup, one has that

Lemma 3.1 (cf. ^[15]). Let $\left(\mathcal{T}(t)\right)_{t\ge 0}$ be a quasi-compact $C_0$-semigroup and $(B, D(B))$ its infinitesimal generator. Then $e^{\delta t}\|\mathcal{T}(t)\| \rightarrow 0$, as $t\rightarrow +\infty$ for $\delta > 0$ if and only if all eigenvalues of $B$ have strictly negative real part.

By the Hille-Yosida estimate in the proof of Theorem 3.2, we have $\|\mathscr{S}(t))\|\leq e^{-\xi t}$. Furthermore, $D\mathcal{F}(\bar{u})\mathscr{S}(t) : X_0 \rightarrow X$ is compact for every $t > 0$. Since

it is seen that $(\mathscr{T}(t))_{t\geq0}$ is quasi-compact. Then by Lemma 3.1 we deduce that, for some $\eta > 0$, $e^{\eta t} \|\mathscr{T}(t)\|\rightarrow 0$ as $t\rightarrow+\infty$ whenever all the eigenvalues of $(\mathcal{A}+D\mathcal{F}(\bar{u}))$ have negative real part.

From the above arguments we can now make the following conclusion.

Theorem 3.5. The solution semi-flow $\Phi(t, u_0)$ of system (1.2), defined as in Theorem 2.1, satisfies the following properties.

(i) If all the eigenvalues of $(\mathcal{A}+D\mathcal{F}(\bar{u}))$ have strictly negative real part, then the steady state $\bar{u}$ is locally asymptotically stable.

(ii) If, however, at least one eigenvalue of $(\mathcal{A}+D\mathcal{F}(\bar{u}))$ has strictly positive part, then the steady state $\bar{u}$ is unstable.

4.   The stability and Hopf bifurcation

Based on the preceding analysis, in this section, we will firstly discuss the global stability of the disease free equilibrium $E_0 = (S^{(0)}, 0, 0)$. Then, we study the stability of the endemic equilibrium $E_* = (S_*, I_*, R_*(a))$. At last, we analyze the existence of the Hopf bifurcation when $E_*$ is unstable.

4.1.   Stability of disease-free equilibrium

Theorem 4.1. If $R_0 < 1$, then disease free equilibrium $E_0 = (S^{(0)}, 0, 0)$ of system (1.2) is globally asymptotically stable. While if $R_0 > 1$, $E_0$ is unstable.

Proof. Let $x(t) = S(t)-S^{(0)}$, $y(t) = I(t)$, $z(t, a) = R(t, a)$. Linearizing system (3.3) at $E_0$ turns out to be the following system:

Solving the second and the third equations of (4.1), we obtain

The second equation of (4.2) implies that $\lim_{t\rightarrow+\infty}y(t) = 0$ when $R_0 < 1$. Furthermore, we have $\lim_{t\rightarrow+\infty}z(t, a) = 0$. Substituting $\lim_{t\rightarrow+\infty}y(t) = 0$ and $\lim_{t\rightarrow+\infty}z(t, a) = 0$ into the first equation of (4.1) and solving for $x(t)$, we obtain that $\lim_{t\rightarrow +\infty}x(t) = 0$ as well. It implies that $E_0$ is locally asymptotically stable when $R_0 < 1$. In addition, when $R_0 > 1$, we have $y(t)\rightarrow\infty$ as $t\rightarrow+\infty$, which implies that $E_0$ is unstable when $R_0 > 1$.

In the following, we construct the Lyapunov function to investigate the global stability of the disease free equilibrium $E_0$ of system (1.2). Define a Lyapunov function

where $\theta > 0$ is constant, which will be determined later, and $\alpha(a) = \int_a^{+\infty}m(\xi)e^{-\int_a^{\xi}(d+m(\eta))d\eta}d\xi$ is a differential function of $a$ on $[0, \infty)$, then the derivative of $V(t)$ along the solution of system is given by

By using of $\Lambda = dS^{(0)}$, we obtain

Because of

we have

It is clear that $S-\frac{\Lambda}{d}\le0$. If $R_0 < 1$, assuming that there exists $\varepsilon > 0$ such that $\theta\alpha(0) = \frac{(d+\mu+\gamma)(1-R_0)}{\gamma}-\varepsilon > 0$, we obtain

Therefore, $R_0 < 1$ ensures that the positive-definite function $V(t)$ has negative derivative $\dfrac{dV(t)}{dt}$. Moreover, the strict equality $\frac{dV(t)}{dt} = 0$ holds if and only if $S(t) = \frac{\Lambda}{d}$, $I(t) = 0$ and $R(t, a) = 0$. Thus, the singleton ${E_0}$ is the largest invariant subset of $\frac{dV(t)}{dt} = 0$. By the LaSalle's invariant principle, the disease free equilibrium $E_0$ is globally attractive. Therefore, $E_0$ is globally asymptotically stable when $R_0 < 1$.

4.2.   Hopf bifurcation around endemic equilibrium $E_*$

Next, we are interested in the stability of endemic equilibrium $E_* = (S_*, I_*, R_*(a))$ and the existence of Hopf bifurcations around $E_*$ when its stability changes. In order to show the local stability of $E_*$, we linearize system (1.2) around $E_*$ and put it as follows:

To analyze the asymptotic behavior of $E_\ast$, we look for solutions of the form $x(t) = x_0e^{\lambda t}$, $y(t) = y_0e^{\lambda t}$, $z(t, a) = z_0(a)e^{\lambda t}$. Thus, we have the following eigenvalue problem:

Solving the second and the third equation in (4.3), we obtain

Substituting $z_0(a)$ and $x_0$ into the first equation of (4.3), we get the following characteristic equation:

where

It is easy to see that

In addition, if $\tau = 0$, then

where $\tilde{b}_i = b_i|_{\tau = 0}, \; i = 0, \; 1, \; 2, \; 3, \; \; \tilde{a}_0 = a_0|_{\tau = 0}$. It is easy to check that $\tilde{b}_i > 0$ for $i = 0, \; 1, \; 2, \; 3$, and $\tilde{a}_0+\tilde{b}_0 > 0$. Therefore, by using of the Routh-Hurwitz criterion, we know if

then all roots of (4.5) have negative real parts, which implies that the endemic equilibrium $E_*$ is locally asymptotically stable. That is, the following result holds:

Theorem 4.2. If $R_0 > 1$, $\tau = 0$, and $\tilde{b}_1\tilde{b}_2 > \tilde{b}_3(\tilde{a}_0+\tilde{b}_0)$, then the endemic equilibrium $E_*$ of system (1.2) is locally asymptotically stable.

In fact, the roots of $f(\lambda, \tau) = 0$ depend on $\tau$ continuously, and the roots may pass through the imaginary axis and enter the right side as $\tau$ increasing. In the following, we discuss the case where $\tau > 0$.

Let $\lambda = i\omega\; (\omega > 0)$ be purely imaginary roots of $f(\lambda, \tau) = 0$. By submitting $\lambda = i\omega$ into $f(\lambda, \tau) = 0$ and separating the real and imaginary parts, we have

which yields that

where $p_3 = b_3^2, \; \; p_2 = b_2^2-2b_1b_3, \; \; p_1 = b_1^2-2b_0b_2, \; \; p_0 = b_0^2-a_0^2$. It is clear that $p_3 > 0$, and $p_0 > 0$.

Put $\Theta = \omega^2$, (4.8) turns out to be

Let $F(\Theta) = 3p_3\Theta^2+2p_2\Theta+p_1$. When $p_2^2-3p_1p_3 < 0$, we know $F(\Theta) = 0$ has no real roots. When $p_2^2-3p_1p_3\geq0$, we know $F(\Theta) = 0$ has two real roots, which are $\Theta_1 = \dfrac{-p_2-\sqrt{p_2^2-3p_1p_3}}{3p_3}$, and $\Theta_2 = \dfrac{-p_2+\sqrt{p_2^2-3p_1p_3}}{3p_3}$, respectively. The following lemma gives the results on the positive root of the equation $Q(\Theta) = 0$.

Lemma 4.1.

(i) If $p_2^2-3p_1p_3 < 0$, then $Q(\Theta) = 0$ has no positive root;

(ii) If $p_2^2-3p_1p_3\geq0$ and $\Theta_2\leq0$, then $Q(\Theta) = 0$ has no positive root;

(iii) If $p_2^2-3p_1p_3\geq0$, $\Theta_2 > 0$, and $Q(\Theta_2) > 0$, then $Q(\Theta) = 0$ has no positive root;

(iv) If $p_2^2-3p_1p_3\geq0$, then $Q(\Theta) = 0$ has positive roots if and only if $\Theta_2 > 0$ and $Q(\Theta_2)\leq0$.

If $Q(\Theta) = 0$ does not have a positive root, then the stability of $E_\ast$ will not change as $\tau$ increasing. Therefore, the following result holds:

Theorem 4.3. Assume that $R_0 > 1$, $\tau > 0$, and $\tilde{b}_1\tilde{b}_2 > \tilde{b}_3(\tilde{a}_0+\tilde{b}_0)$.

(i) If $p_2^2-3p_1p_3 < 0$, then the endemic equilibrium $E_*$ of system (1.2) is locally asymptotically stable;

(ii) If $p_2^2-3p_1p_3\geq0$ and $\Theta_2\leq0$, then the endemic equilibrium $E_*$ of system (1.2) is locally asymptotically stable;

(iii) If $p_2^2-3p_1p_3\geq0$, $\Theta_2 > 0$, and $Q(\Theta_2) > 0$, then the endemic equilibrium $E_*$ of system (1.2) is locally asymptotically stable.

While if $Q(\Theta) = 0$ has positive root, then the stability of $E_\ast$ may change when $\tau$ passes through some specific values. Let $\Theta_\ast$ be the positive real root of (4.9). Then $\omega_\ast = \sqrt{\Theta_\ast}$ is the only one positive real root of (4.8), so $f(\lambda, \tau) = 0$ with $\tau = \tau_k$, $k = 0, 1, 2, \cdots$, has a pair of purely imaginary roots $\pm i\omega_*$, where

for $k = 0, \; 1, \; 2, \; \cdots$, and $c = \dfrac{-b_3\omega^3_*+b_1\omega_*}{a_0}$.

Differentiating both sides of $f(\lambda, \tau) = 0$ with respect to $\tau$ yields

From (4.7) and the fact that $sign\left\{Re\left[\left(\frac {d\lambda}{d\tau}\right)^{-1}\big|_{\lambda = i\omega_\ast}\right]\right\} = sign\left\{\frac {d Re(\lambda)}{d\tau}\big|_{\tau = \tau_k}\right\}$, we have

The transversality condition holds and a Hopf bifurcation occurs at $\tau = \tau_k$, $k = 0, 1, 2, \cdots$ when $Q'(\Theta_*)\neq0$. According to the Hopf bifurcation theorem for functional differential equations ^[17], we have the following result.

Theorem 4.4. Assume $R_0 > 1$.

(i) If $p_2^2-3p_1p_3\geq0$, $\Theta_2 > 0$, and $Q(\Theta_2)\leq0$, then the endemic equilibrium $E_*$ of system (1.2) is asymptotically stable for all $\tau\in[0, \tau_0)$ under condition (4.6);

(ii) If $p_2^2-3p_1p_3\geq0$, $\Theta_2 > 0$, $Q(\Theta_2)\leq0$, and $Q'(\Theta_*)\neq0$, then system (1.2) undergoes Hopf bifurcation at $E_*$ when $\tau = \tau_k$, $k = 0, 1, 2, \cdots$.

5.   Numerical simulations and conclusion

In this paper, we have proposed and analyzed an SIRS model with age structure for recurrent infectious disease by incorporating temporary immunity. A delayed differential equation system can be derived from this model and the delay corresponds to the time that recovery individuals lose their immunity. The purpose of this article is to explore the conditions switching between periodic and non-periodic behavior of recurrent infectious diseases which have the temporary immunity.

On the dynamic behavior analysis of system (1.2), we showed the well-posedness, gave the basic reproduction number $R_0$, and proved that $R_0 = 1$ is the threshold that determines whether the epidemic persists or not by studying the stability of both disease free equilibrium $E_0$ and the endemic equilibrium $E_*$. The disease free equilibrium $E_0$ is globally asymptotically stable if $R_0 < 1$ and is unstable if $R_0 > 1$. Moreover, the disease persists in the later case, in the sense that infected individuals survive above a certain number for any initial infection numbers. We also proved the existence of Hopf bifurcation around the endemic equilibrium $E_*$ when $E_*$ is unstable.

In order to display our conclusions more intuitively, we will use Matlab to demonstrate the nonlinear dynamics behavior of system (1.2). We denote the numbers of recovery individuals at time $t$ as $R(t) = \int_0^{+\infty}R(t, a)da$. Numerically, we set the maximum immunity age as $100$.

Firstly, we illustrate that the disease free equilibrium $E_0$ is globally asymptotically stable when $R_0 < 1$. The parameter values are chosen as $\Lambda = 1$, $d = 0.007$, $\mu = 0.0025$, $\alpha = 0.1$, $\gamma = 0.9$, $m_* = 0.1$, and $\beta = 0.005$. Accordingly, we obtain $R_0 < 1$. Since $R_0$ is independent of $\tau$, $R_0$ does not change with time delay $\tau$. When $\tau = 0$, the solutions of system (1.2) with three different initial values all approach to $E_0$ as $t$ trends to infinity(see Figure 1(a)). When $\tau = 10$, something similar happens(see Figure 1(b)). It means if we control the basic regeneration number $R_0$ to be less than unity, the disease will die out in population. In this case, we don't need to worry about the recurrence of the epidemic disease.

Secondly, we illustrate that the stability of the endemic equilibrium $E_*$ when $R_0 > 1$. Here, we only change the values of $\gamma$ and $\beta$, and keep all other parameter values same as theses in Figure 1. In the case where $\tau = 0$, we take $\gamma = 0.9$ and $\beta = 0.01$, then, $R_0 > 1$, and $\tilde{b}_1\tilde{b}_2 > \tilde{b}_3(\tilde{a}_0+\tilde{b}_0)$. Figure 2(a) shows that the solutions of system (1.2) with three different initial values all approach to $E_*$ as $t$ trends to infinity. In the case where $\tau = 10 > 0$, we take three different pairs of value for $\beta$ and $\gamma$, which are $(\beta, \gamma) = (0.05, 0.9)$, $(\beta, \gamma) = (0.03, 0.1)$, and $(\beta, \gamma) = (0.03, 1)$, respectively. Accordingly, we have (ⅰ) $p_2^2-3p_1p_3 < 0$, (ⅱ) $p_2^2-3p_1p_3 > 0$ and $\Theta_2 < 0$, (ⅲ) $p_2^2-3p_1p_3 > 0$, $\Theta_2 > 0$, and $Q(\Theta_2) > 0$. Theorem 4.3 indicates that the endemic equilibrium $E_*$ of system (1.2) is asymptotically stable under these conditions, which is consistent with the results shown in Figure 2(b). That is, under these conditions of Theorem 4.2 and Theorem 4.3, no matter how the delay $\tau$ changes, the system (1.2) doesn't have the periodic behavior even if the disease persists in the population. Accordingly, the distribution of recovery individuals with respect to immunity age at the endemic equilibrium $E_*$, $R_*(a)$ is shown in Figure 3(a), which corresponding to the second solution line in Figure 2(a), and the distributions with both immunity age and time, $R(t, a)$ is shown in Figure 3(b).

Thirdly, we demonstrate the case when $R_0 > 1$ and $Q(\Theta) = 0$ has the positive root. According to Theorem 4.4, the endemic equilibrium $E_*$ is asymptotically stable for $\tau\in[0, \tau_0)$, and periodic solutions occur as the stability of $E_*$ changes, which means that Hopf bifurcation appears. Setting $\Lambda = 1$, $d = 0.0006$, $\mu = 0.00003$, $\gamma = 3$, $\alpha = 0.000001$, $\beta = 0.02$, and $m_\ast = 0.0085$, the conditions $R_0 > 1$, $p_2^2-3p_1p_3 > 0$, $\Theta_2 > 0$ and $Q(\Theta_2) < 0$ are satisfied. The critical time value is $\tau_0 = 10$.

When time delay is small, say $\tau = 0, 2, 4, 6, 8, 10$, the solutions shown in Figure 4 all approaches to the endemic equilibrium $E_*$ through damped oscillations. While if we choose $\tau = 12, 14, 16, 18, 20$, Figure 5 displays periodic solutions with different periods and amplitudes. That is, with increasing of time delay, the stable endemic equilibrium is replaced by stable periodic orbits. Accordingly, when $\tau = 12$, we displays the distributions $R(t, a)$ with both immunity age and time in Figure 6.

In fact, when Hopf bifurcation exists, it is global continuation. That is, Hopf bifurcation always exits for any $\tau\in[\tau_k, \tau_{k+1}]$. We show this numerically by the following Figure 7, where $\tau = 12.2, 12.4, 12.6, 12.8, 13$, and other parameter values are the same as those of Figure 5. The system presents periodic behaviors for all these chosen time delays.

Based on our analysis, we know that both non-periodic and periodic behaviors are possible when the disease persists in population. These findings are consistent with the results in ^[7,11]. It means that the dynamical behaviors of recurrent epidemics is dependent on the parameters of system (1.2). In particular, immune age and time delay is an important effects on the transmission of recurrent epidemics. Once $R_0 > 1$, we can control how disease spreads through controlling the conditions theorem 4.3 and theorem 4.4. If the parameters of system (1.2) satisfy these conditions of theorem 4.3, or the conditions of theorem 4.4(ⅰ), the recurrent infectious disease is going to be a steady state as $t$ goes to infinity. While if the parameters of system (1.2) satisfy these conditions of theorem 4.4(ⅱ), the recurrent infectious disease will persist in the population in the form of periodic oscillations.

Although the simple delay differential equations model may show the switching between periodic and non-periodic behavior, it cannot express the immunity age of the recovery individual. In fact, the immunity age of the individual is important since the immunity age of the individual can be used to describe whether the recovered individual has the immunity or not. The model (1.2) may more accurately describe how the recurrent infectious disease spreads. In addition, our results also implies that model (1.2) may be closer to the transmission mechanism of the recurrent infectious disease. We suspect that the seasonality of recurrent epidemics is also related to age structure and the delay.

Acknowledgments

We would like to thank the referees very much for the careful review and the valuable comments to this manuscript which improve it greatly.

This work is supported by National Natural Science Foundation of China(grant 11301314, 11671142, and 11371087), by Natural Science Basic Research Plan in Shaanxi Province of China grant 2019JM-081, and by Natural Science Foundation of Shaanxi Provincial Department of Education in China grant 18JK0092.

Conflict of interest

The authors have declared that no competing interests exist.