Hopf bifurcation of the age-structured SIRS model with the varying population sizes

1.   Introduction

The recurrent epidemics refers to that the recovered individuals don't have the permanent immunity, and may be infected due to contact with the infectious individuals again^[1,2]. For epidemics without the permanent immunity and obvious latent infection period, we can structure the SIRS model on the basis of the SIR epidemic model established earlier by Kermack and McKendrick in 1927 to describe it^[3]. In fact, many papers have established and studied the dynamic behavior of SIRS epidemic model based on the above ideas^{[4,5,6,7,8,9,10,11,12]}. However, even if the acquired immunity is temporary, it also doesn't disappear immediately, and it still ensure that the recovered individuals have immunity for a certain period of time, and then the immunity will gradually disappear with time^{[13,14,15,16]}. Therefore, the model with immune age structure can well describe the disappearance of temporary acquired immunity over time^[17,18,19].

At present, most SIRS models with immune age structure assume that epidemics spread in populations with constant population complement or fixed size^[17,18,19]. However, it is unreasonable to assume that the population size is constant if the population growth or decrease is significant, or if the disease causes enough deaths to affect the population size^[5]. Size structure is considered to explore the effect of maturity size on the insect population ^[20].

In this paper, our purpose is to investigate the effect of the varying population sizes on the dynamical behavior of the age-structured SIRS model. For this reason, we structure a age-structured SIRS model with the varying population sizes as follows:

Here, we use $S(t)$ and $I(t)$ to represent the number of the suspectable and infectious individuals at time $t$, respectively, and $R(t, a)$ represents the density of the recovered individuals whose age of acquired immunity is $a$ at time $t$.

In particular, $k(a)$ represents the acquired immunity loss function of the recovered individuals with the acquired immunity age $a$, which interprets the recovered individuals leaving the recovered class and entering the susceptible class again because the recovered individuals lose their acquired immunity. In fact, the acquired immunity of the recovered individuals gradually disappear, Therefore, $k(a)$ follows the function:

and $k(a)\in L^\infty_+((0, +\infty), \mathbb{R})$, where $\tau > 0$ is the longest time for the recovered individuals to maintain the acquired immunity.

In addition, $b$ represents the birth rate of susceptible individuals, $d$ represents the natural mortality, $\beta$ represents the transmission rate, $\mu$ represents the disease-related mortality of infectious individuals, and $\delta$ represents the recovery rate.

Next, we will discuss the non-negativity of the solution and the boundedness of the system (1.1) in Section 2. The existence and stability of all the feasible equilibria in Sections 3 and 4, respectively, including the population extinction equilibrium $P_0$, the disease free equilibrium $P_1$, and the endemic equilibrium $P_*$. The conditions for the system (1.1) experiences Hopf bifurcation are given in Section 5. The numerical simulations and conclusions are presented in Section 6.

2.   Preliminaries and well-posedness

Theorem 2.1. Any solution $(S(t), I(t), R(t, a))$ of system (1.1) starting from any non-negative initial value $(S_0, I_0, R_0(\cdot))$ is non-negative when $t\geq0$, where $R_0(\cdot)\in L^1_+(0, +\infty)$.

Proof. Firstly, we prove the non-negativity of $R(t, a)$ for $t\geq0$. For this purpose, we use the characteristic line method to integrate the third equation of system (1.1), then

It implies that $R(t, a)$ is non-negative for $t\geq0$.

Secondly, we prove the non-negativity of $S(t)$ for $t\geq0$. Without loss of generality, for the positive initial value, we assume that there exists $T_1 > 0$, such that $S(T_1) = 0$, and $S(t) > 0$ when $t\in(0, T_1)$. Then, the first equation in system (1.1) means that $\frac{dS(T_1)}{dt} = \int_0^{+\infty}k(a)R(t, a)da > 0$. It contradicts $S(T_1) = 0$. Namely, $S(t)$ is non-negative for $t\geq0$.

Finally, we prove the non-negativity of $I(t)$ for $t\geq0$. By directly solving the second equation in system (1.1), we have

Obviously, $I(t)$ is non-negative for $t\geq0$ when $I_0\geq0$.

Summing up the above discussion, we know any solution of system (1.1) starting from any non-negative initial value is non-negative.

At the end of this subsection, we discuss the boundedness of the system (1.1). To this end, we denote $\bar{R}(t) = \int_0^{+\infty}R(t, a)da$, and assume that $lim_{a\rightarrow +\infty}R(t, a) = 0$ since the age of an individual is limited in the real world. Then,

That is, $S(t)+I(t)+\bar{R}(t)\leq N_0e^{(b-d)t}$ with $N_0 = S(0)+I(0)+\bar{R}(0)$. It implies that the total population exponentially decays to $0$ when $b < d$, is constant $N_0$ when $b = d$, and gradually tends to positive infinity when $b > d$.

3.   The existence of the equilibria

It is obvious that system (1.1) always has a population extinction equilibrium $P_0(0, 0, 0)$. When $b = d$, system (1.1) only exists the disease free equilibrium $P_1(S^0, 0, 0)$, where $S^0$ is any positive constant. To discuss the existence of the endemic equilibrium $P_*(S_*, I_*, R_*(a))$, we need to solve the following equations:

The second equations of (3.1) implies that

By direct calculation, the third and the last equation of (3.1) implies that

Substituting the expression of $S_*$ and $R_*(a)$ into the first equation of (3.1), we have

Denote

then, the expression of $I_*$ can be rewritten as $I_* = \frac{b-d}{\beta\big(1-M\big)}$. It implies that $b > d$ and $M < 1$ together can guarantee $I_* > 0$. Therefore, we give the conclusion of the existence of the equilibria of system (1.1):

Theorem 3.1. (i) System (1.1) always exists a population extinction equilibrium $P_0 = (0, 0, 0)$.

(ii) If $b = d$, system (1.1) exists a disease free equilibrium $P_1 = (S^0, 0, 0)$, where $S^0$ is any positive constant;

(iii) If $b > d$ and $M < 1$, system (1.1) exists a endemic equilibrium $P_* = (S_*, I_*, R_*(a))$ with

4.   The stability of the equilibria

We will mainly discuss the stability of equilibria in this section, including the global stability of the population extinction equilibrium $P_0$ and disease free equilibrium $P_1$, and the local stability of the endemic equilibrium $P_*$ when $\tau = 0$.

4.1.   Stability of the population extinction equilibrium $P_0$

According to the discussion in Section 2, we know the total population will decrease exponentially and approach $0$ when $b < d$. Furthermore, the following result holds.

Theorem 4.1. If $b < d$, the population extinction equilibrium $P_0$ of system (1.1) is globally asymptotically stable.

Proof. In order to discuss the global stability of $P_0$, we construct the following Lyapunov function:

where $\varphi(a) = \int_a^{+\infty}k(\xi)e^{-\int_a^{\xi}(d+m(\theta))d\theta}d\xi]$ is a differential function of a on $\left[0, \infty\right)$, and $\varphi(0) < 1$. Then

Obviously, $b < d$ can ensure that $\frac{dV_0(t)}{dt} < 0$. Therefore, the Lyapunov stability theorem means that $P_0$ is globally asymptotically stable when $b < d$.

The conclusion of Theorem 4.1 is easy to understand. It means that once the birth rate $b$ of the population is less than the natural mortality $d$, the total population will inevitably decay to extinction as $t$ tends to $+\infty$.

4.2.   The stability of the disease free equilibrium $P_1$

When $b = d$, the total population always keep constant $N_0$. In this case, we study the stability of the disease free equilibrium $P_1$, and give the following stability conclusion.

Theorem 4.2. If $b = d$ and $S^0 < \frac{d+\mu+\delta}{\beta}$, then the disease free equilibrium $P_1$ of system (1.1) is globally asymptotically stable.

Proof. Let

Then, we have

Because of $0\leq S(t)\leq N_0 = S^0$, we have $V_1(t)\geq0$, and $\frac{dV_1(t)}{dt} < 0$ when $S^0 < \frac{d+\mu+\delta}{\beta}$. Therefore, the Lyapunov stability theorem means that $P_1$ is globally asymptotically stable.

Theorem 4.2 shows when the total population remains unchanged, if the number of susceptible people is small, the disease will not spread within the population and eventually die out.

4.3.   The stability of the endemic equilibria $P_*$

When $b > d$, Theorem 3.1(iii) shows that system (1.1) has unique endemic equilibrium $P_*$. In this case, we denote $k_*(a) = k$, $k > 0$. It implies that $M = \frac{k\delta e^{-d\tau}}{(d+k)(d+\mu+\delta)}\triangleq M(\tau)$, and satisfies $0 < M(\tau) < 1$. And then, we give the following conclusion.

Theorem 4.3. If $b > d$, then the endemic equilibrium $P_* = (S_*, I_*, R_*(a))$ of system (1.1) is locally asymptotically stable for $\tau = 0$.

Proof. Let $x(t) = S(t)-S_*$, $y(t) = I(t)-I_*$, and $z(t, a) = R(t, a)-R_*(a)$, then system (1.1) can be reorganized as

Suppose that system (4.3) has the solution in the form of $x(t) = x_0e^{\lambda t}$, $y(t) = y_0e^{\lambda t}$, and $z(t, a) = z_0(a)e^{\lambda t}$. Then, we obtain the following equations:

Direct calculation allows us to obtain the following characteristic equation

with

In the case where $\tau = 0$, $A_i(0) > 0$, $i = 2, 1, 0$, $B_0(0) > 0$, and

Then, according to the Routh-Hurwitz criterion, we can see that the real parts of all roots of $g(\lambda, 0) = 0$ are negative. It implies that $P_*$ is locally asymptotically stable when $b > d$ and $\tau = 0$.

If $b > d$, the total population gradually tends to $+\infty$ as $t$ tends to $+\infty$. In this case, Theorem 4.3 implies that system (1.1) does not exist the disease-free equilibrium, and the disease will certainly spread within the population. In addition, since the roots of $g(\lambda, \tau) = 0$ also depend on $\tau$ continuously, the roots of $g(\lambda, \tau) = 0$ may pass through the imaginary axis and enter the right side as $\tau$ increases, which may lead to the change in the stability of $P_*$. For this reason, we will study the stability of $P_*$ when $\tau > 0$ in the following.

5.   Hopf bifurcation

In the case where $b > d$, $k_*(a) = k$, and $\tau > 0$, since the coefficients $A_i(\tau)(i = 0, 1, 2)$ and $B_0(\tau)$ of the characteristic equation $g(\lambda, \tau) = 0$ are functions of $\tau$, we can use the method in ^[21] to study the stability of endemic equibrium $P_*$ and the conditions under which Hopf bifurcation occurs. To this end, the characteristic equation $g(\lambda, \tau) = 0$ can be rewritten as

with

It is clear that $A_i(\tau)(i = 0, 1, 2)$ and $B_0(\tau)$ are all continuous and differentiable respect to $\tau$, and satisfy

In addition, since $P_3(\lambda, \tau)$ and $Q_0(\lambda, \tau)$ are both differentiable in $\tau$ and analytic functions in $\lambda$, we can get

and

It implies that these assumptions (i)–(iv) in ^[21] hold.

Suppose $\lambda = i\omega(\tau)$, $\omega(\tau) > 0$ is the root of the characteristic of $g(\lambda, \tau) = 0$, then we know

Put $\Theta = \omega^2$, then (5.2) can be rewrite as

where $q_2(\tau) = A_2^2(\tau)-2A_1(\tau)$, $q_1(\tau) = A_1^2(\tau)-2A_2A_0(\tau)$, and $q_0(\tau) = A_0^2(\tau)-B_0^2(\tau) > 0$. Let $\Theta^* = \frac{-q_2(\tau)+\sqrt{q_2^2(\tau)-3q_1(\tau)}}{3}$ be the roots of $Q'(\Theta) = 0$ when $q_2^2(\tau)-3q_1(\tau)\geq0$, we give the results of the positive root of the equation $Q(\Theta) = 0$ in the following Lemma 5.1.

Lemma 5.1.

(i) If $q_1(\tau) < 0$, the equation $Q(\Theta) = 0$ has the positive real root when $Q(\Theta^*)\leq0$, andhas no positive root when $Q(\Theta^*) > 0$;

(ii) If $q_1(\tau)\geq0$, $q_2(\tau) < 0$ and $q_2^2(\tau)-3q_1(\tau)\geq0$, the equation $Q(\Theta) = 0$ has the positive real root when $Q(\Theta^*)\leq0$, and has no positive root when $Q(\Theta^*) > 0$;

(iii) If $q_1(\tau)\geq0$ and $q_2(\tau)\geq0$, the equation $Q(\Theta) = 0$ has no positive root.

If there is no positive real root for equation $Q(\Theta) = 0$, the stability of the endemic equilibrium $P_*$ will not change with the increase of $\tau$. It implies that the following conclusions are true.

Theorem 5.1. Suppose that $b > d$ and $M(\tau) < 1$.

(i) If $q_1(\tau) < 0$, the endemic equilibrium $P_*$ of system (1.1) is locally asymptotically stable for any $\tau > 0$ when $Q(\Theta^*) > 0$;

(ii) If $q_1(\tau)\geq0$ and $q_2(\tau) < 0$, the endemic equilibrium $P_*$ of system (1.1) is locally asymptotically stable for any $\tau > 0$ when $Q(\Theta^*) > 0$;

(iii) If $q_1(\tau)\geq0$ and $q_2(\tau)\geq0$, the endemic equilibrium $P_*$ of system (1.1) is locally asymptotically stable for any $\tau > 0$.

Theorem 5.1 implies that the stability of the endemic equilibrium $P_*$ does not change as $\tau$ increases when $b > d$, $M(\tau) < 1$, and the conditions in (i), (ii), and (iii) of Theorem 5.1 hold. Next, let's discuss the case where the equation $Q(\Theta) = 0$ has the positive real root. Without loss of generality, we assume that $\Theta_*$ is the root of $Q(\Theta) = 0$. That is, $\omega(\tau_*) = \sqrt{\Theta_*}$ is the positive real root of $G(\omega, \tau) = 0$. Thus, by Lemma 5.1, we can define a set

with

Namely, there exists $\omega = \omega(\tau) > 0$ such that $G(\omega, \tau) = 0$ when $\tau\in\Gamma$.

Define $\theta(\tau)\in[0, 2\pi]$, $\tau\in\Gamma$ as

The relation between $\omega(\tau)\tau$ and $\theta(\tau)$ implies that $\omega(\tau)\tau = \theta(\tau)+2l\pi$, $\tau\in\Gamma$, $l\in\mathbb{N}$. And then, we define $S_l(\tau)$ as

where $\omega_* > 0$ with $\omega_* = \omega(\tau_*)$ is the positive root of (5.2) if and only if $S_l(\tau_*)$ = 0 for some $\tau_*\in\Gamma$, $l\in\mathbb{N}$. Then, Theorem 2.2 in ^[21] guarantees that the characteristic equation (5.1) has a pair of simple conjugate pure imaginary roots when $\tau = \tau_*$, which are $\lambda_+(\tau_*) = i\omega(\tau_*)$, and $\lambda_-(\tau_*) = -i\omega(\tau_*)$, respectively. Furthermore, it crosses the imaginary axis from left to right if $F(\tau_*) > 0$ and crosses the imaginary axis from right to left if $F(\tau_*) < 0$, where

Based on the Hopf bifurcation theory of functional differential equations in ^[22], we know the following conclusion is true.

Theorem 5.2. Suppose that $b > d$, $Q(\Theta^*)\leq0$, and $Q'(\Theta_*)\neq0$.

(i) If $q_1(\tau) < 0$, the endemic equilibrium $P_*$ of system (1.1) is locally asymptotically stable for $\tau\in[0, \tau_*)$, and system (1.1) experiences the Hopf bifurcation at $\tau = \tau_*$;

(ii) If $q_1(\tau)\geq0$, and $q_2(\tau) < 0$, the endemic equilibrium $P_*$ of system (1.1) is locally asymptotically stable for $\tau\in[0, \tau_*)$, and system (1.1) experiences the Hopf bifurcation at $\tau = \tau_*$.

6.   Numerical simulations and conclusions

We have established and discussed a age-structured SIRS epidemic model with the varying population sizes in this paper. The purpose is to investigate the impact of the varying population sizes on the dynamic behavior of the age-structured epidemic model. In theoretical analysis, we study the boundedness and non-negativity of the solution of system (1.1), the existence and stability of the equilibria of system (1.1), including the population extinction equilibrium $P_0$, the disease free equilibrium $P_1$ and the endemic equilibrium $P_*$, and give the conditions that system (1.1) undergoes Hopf bifurcation.

In order to use the numerical simulations to illustrate the complex dynamical behaviors of system (1.1), let the maximum infection age be $500$. Firstly, we demonstrate the stability of $P_*$ in the case where $b > d$ and $\tau\geq0$. Taking $\tau = 0$, $b = 0.0852$, $d = 0.001$, $\mu = 0.01$, $\beta = 0.0011$, $\delta = 0.01$, and $k = 0.06$, Figure 1(a) shows that $P_*$ is stable focus. Choosing $\tau = 180,200,220$, respectively, and choosing $b = 0.0155$, $\beta = 0.001$, $\delta = 0.005$, $d = 0.0001$, $\mu = 0.001$, $k = 0.01$, then we have $q_1(\tau) > 0$, $q_2(\tau) > 0$, and Figure 1(b) displays $P_*$ is stable focus. When we fix $b = 0.0852$, $\beta = 0.001$, $\delta = 0.01$, $d = 0.001$, $\mu = 0.01$, $k = 0.01$, and $\tau = 180,200,220$, respectively, then $q_1(\tau) > 0$, $q_2(\tau) < 0$, and $Q(\Theta^*) > 0$, Theorem 5.1(ii) implies that $P_*$ is stable focus, see Figure 1(c). In addition, we take $b = 0.0155$, $\beta = 0.001$, $\delta = 0.005$, $d = 0.002$, $\mu = 0.005$, $k = 0.01$, and the value of $\tau$ is consistent with these in Figure 1(c), then $q_1(\tau) < 0$, and $Q(\Theta^*) > 0$. In this case, Figure 1(d) shows that $P_*$ is still stable focus, which is consistent with Theorem 5.1(i).

And then, we demonstrate that system (1.1) experiences the Hopf bifurcation as $\tau$ changes when $q_1(\tau)\geq0$, $q_2(\tau) < 0$, and $Q(\Theta^*) < 0$. Fixing $b = 0.00852$, $d = 0.001$, $\mu = 0.03$, $\beta = 0.002$, $\delta = 0.005$, $k = 0.009$, and taking $\tau = 180,200,220$, respectively, we have $q_1(\tau) > 0$, $q_2(\tau) < 0$, $Q(\Theta^*) < 0$, and $\tau_* = 163$. Figure 2 displays system (1.1) undergoes the Hopf bifurcation, where Figure 2(a) is the time series diagram, and Figure 2(b) is the phase diagrams.

Similarly, we also demonstrate that system (1.1) experiences the Hopf bifurcation when $q_1(\tau) < 0$ and $Q(\Theta^*) < 0$. For this purpose, we chose $b = 0.00852$, $d = 0.002$, $\mu = 0.01$, $\beta = 0.001$, $\delta = 0.005$, $k = 0.015$, and $\tau = 280,290,300$, respectively, then we have $q_1(\tau) < 0$, $Q(\Theta^*) < 0$, and $\tau_* = 273$. Figure 3 displays that system (1.1) experiences the Hopf bifurcation when $q_1(\tau) < 0$ and $Q(\Theta^*) < 0$, where Figure 3(a) is the time series diagram, and Figure 3(b) is the phase diagrams.

Finally, we take $k(a) = 0.1669e^{-1.319a}+0.0119e^{-0.000152a}$ to numerically study the complex dynamic behavior of system (1.1). Let $b = 0.00852$, $d = 0.001$, $\mu = 0.03$, $\beta = 0.002$, $\delta = 0.005$, and choose $\tau = 180,200,220$, respectively, Then, Figure 4 displays system (1.1) still undergoes the Hopf bifurcation.

In particular, fixing $b = 0.0155$, $d = 0.001$, $\mu = 0.1$, $\beta = 0.001$, $\delta = 0.05$, and using $\tau$ as the bifurcation parameter, we demonstrates system (1.1) may happen multiple stability switches (see Figure 5(a)). That is, the varying population sizes can lead to the age-structured SIRS epidemic model happens the complex dynamical behavior.

It is clear that $\tau$ cannot reach to the infinite since the coefficient of the characteristic equation $g(\lambda, \tau) = 0$ is related to $\tau$. The numerical simulation shows that the system (1.1) is ultimately stable when $\tau = 250$, see Figure 5(b), where the parameter values are the same as those in Figure 5(a).

Acknowledgments

This work is supported by National Natural Science Foundation of China 12071268 and 11971281, and by Natural Science Basic Research Plan in Shaanxi Province of China Grant 2022JM-038.

Conflict of interest

The authors declare there is no conflict of interest.