An age-structured epidemic model with boosting and waning of immune status

1.   Introduction

The dynamics of immune status among host individuals plays a crucial role in the spread of infectious diseases. Individuals which have just recovered from disease obtain any immunity, however, its level of effectiveness is not necessarily time-constant but can vary as time evolves. An individual's immunity may decay as time goes by, while it could be regained by boosting, which means the immunity enhancement by continued or intermittent exposure to infectious agent.

The boosting and waning dynamics of immune status and its epidemiological consequences have been studied by many authors. Historically speaking, it should be firstly noted that Kermack and McKendrick ^[2] had already developed a class-age structured epidemic model allowing reinfection of recovered individuals from the 1930s. Reinfection of recovered individuals can be seen as a consequence of waning immunity, and it can cause a backward bifurcation of endemic steady state (^[3,4]). If the loss of immunity of recovered individuals implies the direct reversion from the recovered class to the susceptible class, the basic system is formulated as the well-known SIRS model, which has been long studied in mathematical epidemiology. For example, Bhattacharya and Adler ^[5] developed a recovery-age dependent SIRS model to show that loss of immunity could produce oscillations. On the other hand, Nakata, et al. ^[6] considered an infection-age structured SIRS model to examine global stability conditions for an endemic equilibrium. Okuwa, et al. ^[7] studied a chronological age-structured SIRS model to establish its endemic threshold results.

It is more challenging to consider synergistic consequences of boosting and waning of immunity. Heffernan and Keeling ^[8] discussed how vaccination can have a range of unexpected consequences as it reduces the natural boosting of immunity as well as reducing the number of naive susceptibles. Kababi{\'r}e, et al. ^[9] examined a model for a vector-borne disease with immunity decay and immunity boosting of human hosts. Leung, et al. ^[10] studied how immune boosting and cross immunity can influence the timing and severity of epidemics. It was also shown that the epidemiological patterns of an infectious disease may change considerably when the duration of vaccine-acquired immunity differs from that of infection-acquired immunity ^[11]. Dafilis, et al. ^[12] showed that allows for boosting of immunity may naturally give rise to undamped oscillatory behavior for biologically realistic parameter choices. Diekmann, et al. ^[13] developed a mathematical model to describe the distribution of immune status shaped by continuous waning and occasional boosting. Epidemiological consequences of waning of vaccine-derived immunity are also discussed in ^[14], which topic should be paid more attention in the context of COVID-19.

In a series of papers from the 1980s, Aron (^[15,16,17]) developed mathematical models for malaria to consider the effect of immunity boosting by reinfection, and discussed that the boosting mechanism may explain the age-specific prevalence of acute malaria and the functional relationship between the rate of reversion from the immune class to the susceptible class and the force of infection at the endemic steady state, although those pioneering works have been so far almost neglected. In those models, the host population are divided into three classes: susceptibles (uninfected individuals); infecteds (the individuals with severe and symptomatic infection); immunes (the individuals with mild, subclinical, or asymptomatic infection). The recovered individuals are assumed to be partially susceptible and infective, and their immune status can be boosted by reinfection. Since subclinical individuals are assumed to become completely susceptible again with a recovery-age dependent reversion rate, so the Aron model can be seen as an extension of the well-known SIRS epidemic model.

In ^[4], Inaba showed that the Aron model can be more generally formulated as an age-structured SIRS model, where asymptomatic infecteds are distributed by the recovery-age (or immunity clock), that is, time elapsed since the moment of recovery from the symptomatic infectious class. In Inaba's formulation, the boosting effect is expressed by a boundary condition such that the immunity clock is reset to zero by boosting (reinfection). Under this formulation, the functional relationship between the reversion rate and the force of infection is then naturally induced, and Aron's result is obtained as a special case in which the probability density function for the loss of immunity is given by a delta function. If the immunity level is waning monotonically as time evolves, it is inversely proportional to recovery-age, that is, the reset of the immunity clock means the boosting of immunity.

By using the immunity clock mechanism, we could consider the dynamics of continuous change of individual immune status in the spread of infectious diseases, which is more advantageous and natural than the traditional compartment model formulation. For another kind of structured boosting models, the reader may refer to ^[18], ^[19] and ^[20], where immune individuals are structured by a general immunity level parameter. Although the immunity clock mechanism may be more restrictive than the general immunity level parameter dynamics, the clock mechanism is theoretically much more tractable.

For the Aron–Inaba formulation, a crucial assumption is that newly boosted individuals have the same immunity level as individuals who have just recovered from symptomatic infection, but it would be a restrictive assumption in order to consider the incomplete boosting. Therefore, in this paper, we extend the Aron–Inaba model so that the effect of boosting is that the immunity clock can be reset to any time less than the recovery-age at which reinfection occurs. If the immunity level is monotone decreasing with respect to the recovery-age, our assumption implies that newly boosted (reinfected) individuals could get, with a given probability, any level of immunity less than the maximum level that is gained by recovery from symptomatic infection status.

In the following, we first give the well-posedness result of our basic system. Next, we investigate the initial invasion condition which is formulated by the local stability of disease-free steady state. Thirdly we consider the existence of endemic steady states. Finally, we provide a bifurcation analysis of endemic steady states. Based on Lyapunov–Schmidt type arguments, we determine the direction of bifurcation that endemic steady states bifurcate from the disease-free steady state when the basic reproduction number passes through the unity. We give a necessary and sufficient condition for backward bifurcation to occur.

2.   Basic model and preliminary observations

We consider a closed population with the constant total size $N$ divided into three parts: susceptibles, clinical (symptomatic) infectives and subclinical (asymptomatic) infectives. The subclinical class is structured by variable $\tau$, called the recovery-age, which is the time elapsed since the recovery from clinical infectious class. Let $S(t)$ and $I(t)$ be the size of susceptible and clinical infectives at time $t$, respectively, and let $r(t, \tau)$ be the age density function of subclinical infectious class. By extending the Aron–Inaba model ^[4], we formulate the following age-structured epidemic model:

where $\mu > 0$ is the natural death rate, $q > 0$ is the recovery rate from clinical infection, $\gamma(\tau)$ is the reversion rate at recovery-age $\tau$. The force of infection $\lambda$ is given by

where $\beta_1$ and $\beta_2(\tau)$ are transmission coefficients reflecting the infectivity for clinical individual and subclinical individual at recovery-age $\tau$, respectively. For simplicity, we assume that the same force of infection $\lambda$ is applied to both susceptibles and subclinical individuals.

As a technical but biologically reasonable assumption, we suppose that $\gamma$, $\beta_2 \in L^\infty (0, \infty)$ and denote their upper bound as $\gamma^\infty$ and $\beta_2^\infty$, respectively. Moreover, we also define the survival function $\Gamma$ induced from the reversion rate $\gamma$ as follows:

The parameter $p(\tau, \sigma)$ is a key to formulate our novel idea, which means the proportion that subclinical individuals with recovery-age $\sigma$ reset their recovery-age to $\tau < \sigma$ when their immunity is boosted by reinfection. Then we assume the following condition:

so that the number of subclinical population conserves through the boosting, that is,

holds. Furthermore, we assume $p^\infty: = \sup_{\tau < \sigma}p(\tau, \sigma) < \infty$. For instance, following types of function satisfy the condition (2.4):

1. $p(\tau, \sigma) = e^{-(\sigma-\tau)}$.

2.

where $s > 0$ is a given constant.

3. $p(\tau, \sigma) = \frac{1}{1+\sigma}$.

4. $p(0, \sigma) = 0$ and $p(\cdot, \sigma)$ is an arbitrary probability density function whose support is $[0, \sigma]$.

Note that the original Aron model corresponds to the special case that $p(0, \sigma) = 1$ and $p(\tau, \sigma) = 0$ for all $\tau > 0$ (^[4], chapter 8). As is mentioned above, if the immunity level is monotone decreasing with respect to recovery-age, our assumption implies that, by through the reset of recovery-age, the immunity level of reinfected individuals could be boosted to any level of immunity less than the highest immunity level gained by recovery from the clinical status. Finally it is easily checked that

is constant for all $t\geq 0$, since we assume there is no disease-induced death rate.

As a preliminary study for our basic system (2.1), we here consider the linearized system at the diesease-free steady state $(S, I, r) = (N, 0, 0)$. Let $(Y, z(\cdot))$ be the perturbation of the infected population densities $(I, r)$ from the equilibrium state $(0, 0)$. Then the linearized system at the disease-free steady state is described as follows:

which is already introduced in ^[4], because the boosting effect is the second order effect, so the parameter $p$ does not appear in (2.7), and the effect of immunity boosting can be neglected as long as we consider the initial invasion phase.

Let us introduce the renewal equation for the density of newly infecteds denoted by $v(t): = N\lambda(t)$. The linearized equation for clinically infected population density in (2.7) implies

Integrating the second equation in (2.7) along the characteristic line, we obtain

Then we obtain the following integral equation:

where the initial data is given by

Changing the order of integrals, we have the following renewal equation:

where the net reproduction kernel $\Phi$ is defined as

Then the basic reproduction number $R_0$ is given by

which is the expected number of secondary cases produced by an infected individual during its entire period of infectiousness in a completely susceptible population.

Note that the basic reproduction number is traditionally defined as the spectral radius of the next generation operator (^[4,21,22]). Since the density of newly infecteds in the initial invasion phase is described by the renewal integral equation as (2.11), its next generation operator is given by

which is acting on the state space of infected individuals. In our case, $K$ is a scalar, so $K$ itself gives $R_0$.

From the epidemiological threshold principle, we can expect that if $R_0 > 1$, the disease can invade into the host population, while it cannot if $R_0 < 1$. We will give a mathematical justification for this invasion principle for our system in section 4.

3.   Well-posedness of the basic model

Here we prove the existence, uniqueness and positivity of solution semiflow of the system (2.1). If $I$ and $r(\cdot)$ are once determined, $S$ is given by the relation (2.6). So it is sufficient to consider the $(I, r)$-system as follows:

where the state space of $(I, r)$ is given by $\left\{(I, r) \in \mathbb R_+ \times L^1_+(\mathbb R_+): I+\|r\|_{L^1} \le N\right\}$.

In order to formulate the basic system (3.1) as a semilinear abstract Cauchy problem (^[23,24,25]). Let us define the extended state space $Z$ as $Z = \mathbb R \times \mathbb R \times L^1(0, \infty)$ and its closed subspace $Z_0 = \{ 0 \} \times \mathbb R \times L^1(0, \infty)$, where the first element corresponds to the boundary value of $r$, and the second and third parts give the states of $I$ and $r$.

Let $A$ be a differential operator on $X: = \mathbb R \times L^1(\mathbb R_+)$ such as

and define the linear operator $\mathcal A$ acting on $Z$ such that

Next define the nonlinear operator $\mathcal F: Z_0 \to Z$ as

where $\lambda[\cdot]: \mathbb R_+ \times L^1(\mathbb R_+) \to \mathbb R_+$ is a functional defined as

Then it is easy to see that $\mathcal F$ is a Lipschitz continuous nonlinear perturbation, that is, there is a number $L > 0$ such that $\|\mathcal Fu-\mathcal Fv\|_Z \le L\|u-v\|_Z$ for all $u, v \in Z_0$.

Now formally the system (3.1) can be written as a semilinear abstract Cauchy problem in $Z$:

From biological meaning, we are interested in the solution semiflow in $C_0 \subset Z_0$, where

Then $\mathcal A$ is not densely defined in $Z$, but its domain is dense in $Z_0$, that is, $\overline{D(\mathcal A)} = Z_0$. Moreover, we can observe that $\mathcal A$ is a resolvent positive, Hille-Yosida type operator as follows:

Lemma 1. 1. $\{ \zeta\in\mathbb C \mid \Re \zeta > -\mu \}\subset\rho(A)$.

2. For $\zeta > -\mu$, the Hille-Yosida estimation holds as

3. For $\zeta > 0$, it holds that $\zeta(\zeta-\mathcal A)^{-1}(C)\subset C_0$.

Proof. Let $\Re \zeta > -\mu$. Consider the following resolvent equation

Then $v_2(0) = u_3$, $\zeta v_1+\mu v_1 = u_1$ and $\zeta v_2(\tau)+(\mu+\gamma(\tau))v_2(\tau)+v_2^\prime (\tau) = u_2 (\tau)$ hold. These imply

Hence for $\zeta > -\mu$ and $u\in Z$

From the above, it is obvious that $\zeta(\zeta-\mathcal A)^{-1}(C)\subset C_0$ for given $\zeta > 0$.

Subsequently we can observe the quasi-positivity of the nonlinear perturbation term $\mathcal F$:

Lemma 2. For sufficiently small $\epsilon > 0$, $(I+\epsilon \mathcal F)(C_0)\subset C$.

Proof. Let $\psi = \begin{pmatrix} 0\\ \begin{pmatrix} \psi_1\\ \psi_2(\cdot) \end{pmatrix} \end{pmatrix} \in C_0,$ then we have

Since the estimation

holds for $\psi\in C_0$, if we take $\epsilon > 0$ such as

then $(I+\epsilon\mathcal F)(C_0)\subset Z_+$. Suppose $0 < \epsilon < \epsilon_0$. Due to the positivity of each element, we obtain

where the condition (2.4) is used. Moreover, since $1-\epsilon \lambda [\psi] > 1-\frac{\lambda [\psi]}{ \max\{\beta_1, \beta_2^\infty\} }\geq 0$, the following estimation holds:

which implies $(I+\epsilon \mathcal F)(C_0)\subset C$.

In order to construct a semiflow in $C_0$, we rewrite (3.6) as follows:

where $\epsilon$ is a small positive number such that $I+\epsilon \mathcal F$ maps $C_0$ into $C$. Hence it is expected that the solution semiflow of (3.9) could be obtained by a positive iteration for the variation of constants formula, but the operator $\mathcal A$ has a non-dense domain in $Z$, $\mathcal A-1/\epsilon$ is not an infinitesimal generator of a strongly continuous semigroup. Thus we seek a solution in a weak sense.

For a given $T > 0$, $u: [0, T) \to C_0$ is an integral solution of (3.9) if $u$ is continuous on $[0, T)$ and it satisfies

and

Note that (3.10) implies that $u(t) \in Z_0 = \overline{D(\mathcal A)}$ and the integral solution becomes the classical solution if and only if $u$ is differentiable on $(0, T)$ ^[26].

Let $\mathcal A_0$ be the part of $\mathcal A$ in $Z_0$, that is,

Then the following holds:

Lemma 3. $\overline{D(\mathcal A_0)} = Z_0$ holds and $\mathcal A_0$ generates a strongly continuous semigroup $\{ \mathcal T_0 (t) \}_{t\geq 0}$ on $Z_0$ and $\mathcal T_0(t) C_0 \subset C_0$ for all $t > 0$.

Proof. It follows from definition that $\begin{pmatrix} 0\\ \psi \end{pmatrix}\in D(\mathcal A_0)$ if and only if $\psi \in D(A)$ with $\psi_2(0) = 0$. Then we have $\overline{D(\mathcal A_0)} = Z_0$. Define a differential operator $A_0$ such that

where

Then $A_0$ is a generator of $C_0$-semigroup $T_0(t)$, $t \ge 0$ on $\mathbb R \times L^1(\mathbb R_+)$ such that

where $U(t)$, $t \ge 0$ denotes the translation semigroup on $L^1$ defined by

and we have $D(\mathcal A_0) = \{0\} \times D(A_0)$. Thus $\mathcal A_0$ is densely defined on $Z_0$ and it generates a $C_0$-semigroup on $Z_0$ such that

Since the semiflow $T_0(t)$, $t > 0$ makes a closed subset $\{\psi = (\psi_1, \psi_2) \in \mathbb R_+ \times L^1_+(\mathbb R_+): |\psi_1|+\|\psi_2\|_{L^1} \le N\}$ positively invariant, we have for $t > 0$, $(\mathcal T_0 (t)) (C_0)\subset C_0$.

Proposition 1. For a given $T > 0$ and $u(0) \in C_0$, a positive function $u \in C(0, T; C_0)$ is an integral solution for $(3.9)$ if and only if it is the positive continuous solution of the extended variation of constants formula on $C_0$ as

Proof. First we assume that $u \in C(0, T; C_0)$ is an integral solution of (3.9). Applying the resolvent operator $\zeta(\zeta-(\mathcal A-\frac{1}{\epsilon}))^{-1}$ with $\zeta > 0$ to the both sides of (3.9), we obtain

where $u_\zeta(t): = \zeta(\zeta-(\mathcal A-\frac{1}{\epsilon}))^{-1}u(t)$. Then it is easily seen that $u_\zeta \in D(\mathcal A_0) \cap C_0$. Since $\zeta(\zeta-(\mathcal A-\frac{1}{\epsilon}))^{-1}(\mathcal A-\frac{1}{\epsilon})$ is bounded, it follows that

On the other hand, due to the continuity of $(I+\epsilon \mathcal F)u(t)$, integral and the resolvent operator can be interchanged in the third part of the right hand side in (3.18). Therefore we arrive at the following equation on $D(\mathcal A_0) \cap C_0$:

Then it follows that $u_\zeta$ is a classical solution of the Cauchy problem on $Z_0$ as

Applying the standard variation of constants formula to (3.19), we have

Therefore, it follows from $\lim_{\zeta \to \infty} u_\zeta = u$ that we have (3.17). Conversely, we assume that $u(t)$ is the positive continuous solution of (3.17). Consider a linear inhomogeneous Cauchy problem on $Z$ as:

Since the generator $\mathcal A-\frac{1}{\epsilon}$ satisfies the Hille–Yosida estimate, it follows from Da Prato–Sinestrari theorem ^[27] that (3.22) has a unique integral solution $v(t) \in Z$ such that

Applying the same argument as above, we obtain the extended variation of constants formula again as:

which shows that $u(t) = v(t)$ for all $t \in [0, T)$. Then the positive continuous solution of (3.17) is the integral solution of the basic system (3.9).

From the above arguments, we obtain the well-posedness result for the basic system (3.6) as follows:

Proposition 2. The basic model $(3.6)$ has a unique, global, positive integral solution.

Proof. From the above observation, it is sufficient to show that the extended variation of constants formula (3.17) has a unique continuous solution for some $T > 0$. Define an operator $\mathcal H$ on $C(0, T;C_0)$ by

Since the right hand side of (3.25) is a linear convex combination of elements in $C_0$, so $\mathcal H$ maps $C_0$ into itself. Therefore, using the local Lipschitz continuity of the bounded perturbation $\mathcal F$ and the contraction mapping principle to $\mathcal H$, we can find a small $T > 0$ such that a positive local solution for the extended variation of constants formula (3.17) uniquely exists. Moreover, it follows from the fact $\|\mathcal T_0(t)\| \le 1$ and $\|\zeta (\zeta-(\mathcal A-\frac{1}{\epsilon}))^{-1}\| \le 1$ that

Hence we know that the fixed point of $\mathcal H$ grows at most exponentially as time evolves, so the local solution can be extended to a global solution.

From the above argument, we know that for each $u(0) \in C_0$, (3.9) has a unique integral solution $u: \mathbb R_+ \to C_0$, then $\mathcal T(t)u(0) = u(t)$ defines the nonlinear semigroup $\mathcal T(t)$, $t \ge 0$ on $C_0$. Here we consider the stability of steady states of $\mathcal T(t)$, $t \ge 0$, which are time-independent solutions of (3.11)

4.   Disease-free steady state and its stability

Here we examine the initial invasion phase of epidemic, that is, we consider the local stability of the disease-free equilibrium $u = 0$ of $\mathcal T(t)$, $t \ge 0$. Since $\mathcal F$ is continuously F\'rechet differentiable at a steady state $u = u^*$, the nonlinear semigroup $\mathcal T(t)$, $t \ge 0$ is also Fréchet differentiable at $u = u^*$ and its derivative at $u = u^*$ is the strongly continuous linear semigroup $\mathcal S(t)$, $t \ge 0$ generated by the part of $\mathcal B: = \mathcal A+\mathcal F'[u^*]$ in $Z_0$. Let $\mathcal B_0$ be the part of $\mathcal B$ in $Z_0$. Then the growth bound $\omega_0(\mathcal B)$ and the essential growth bound $\omega_1(\mathcal B)$ are defined as

if these exist, where $\alpha[A]$ is the measure of noncompactness of an operator $A$. Let $s(\mathcal B_0): = \sup_{\zeta\in\sigma(\mathcal B_0)} \Re \zeta$ be the spectral bound of $\mathcal B_0$. Then the following holds (^[28,29]):

If $\omega_0(\mathcal B_0) < 0$, we can conclude that the steady state $u = u^*$ is locally asymptotically stable, in the sense that for $\omega \in (\omega_0(\mathcal B_0), 0)$ there exists $\delta > 0$ such that $\|\mathcal T(t)u-u^*\|_Z \le e^{\omega t}\|u-u^*\|_Z$ for all $u \in C_0$ with $\|u-u^* \|_Z \le \delta$. On the other hand, $u^*$ is unstable if at least one eigenvalue of $\mathcal B_0$ has strictly positive real part (Theorem 4.2, Corollary 4.3 in ^[25]).

Here let us consider the stability of the disease-free steady state for (3.6). The linearized system of (3.6) at $u = 0$ is given as follows:

where the Fréchet derivative of $\mathcal F$ at $0\in Z$ is denoted by $D\mathcal F[0]$ and is calculated as

Since the range of $D\mathcal F[0]$ is finite, it is a compact operator. In the following we investigate the spectrum of the part $\mathcal B_0$ of the linearized operator $\mathcal B = \mathcal A+D\mathcal F[0]$.

Lemma 4. Let $P_\sigma(\mathcal B_0)$ be the point spectrum of $\mathcal B_0$ and let $\Sigma: = \{z \in \mathbb C: \Re z > -\mu\}$. Then it follows that $\sigma(\mathcal B_0) \cap \Sigma = P_\sigma(\mathcal B_0) \cap \Sigma$.

Proof. Consider the resolvent equation

Then we have

Solving the above ODE, we have

Moreover, formally we can calculate $v_1$ as

where

Inserting the expression of (4.6) into (4.5), $v_2$ can be expressed by $u_1$, $u_2$, so we can formally calculate $(\zeta-\mathcal B_0)^{-1}u$. It is clear that $(\zeta-\mathcal B_0)^{-1}u \in D(\mathcal A)$ if $u \in Z_0$ and $\Re \zeta > -\mu$ and $\Delta(\zeta)\neq 1$, so $\{ \zeta\in\mathbb C: \Delta(\zeta) \neq 1, \Re\zeta > -\mu \} \subset \rho(\mathcal B_0)$. Then we have $\Omega: = \{ \zeta\in\mathbb C : \Delta(\zeta) = 1, \Re \zeta > -\mu \} \supset \sigma(\mathcal B_0) \cap \Sigma$. On the other hand, if $\zeta \in \Omega$, we can calculate an eigenfunction for $\mathcal B_0$ associated with $\zeta$ as $(0, v_1, qv_1 e^{-(\mu+\zeta)\tau}\Gamma(\tau))$ where $v_1$ is any nonzero constant, so we have $\zeta \in P_\sigma(\mathcal B_0)$. Thus we have $P_\sigma(\mathcal B_0) \cap \Sigma \supset \Omega \supset \sigma(\mathcal B_0) \cap \Sigma$, but trivially it follows that $P_\sigma(\mathcal B_0) \cap \Sigma \subset \sigma(\mathcal B_0) \cap \Sigma$, hence we obtain $\sigma(\mathcal B_0) \cap \Sigma = \Omega = P_\sigma(\mathcal B_0) \cap \Sigma$.

Lemma 5. If $\Delta(-\mu) \ge 1$, $\Omega = \{ \zeta\in\mathbb C : \Delta(\zeta) = 1, \Re \zeta > -\mu \}$ has a dominant real normal eigenvalue $\zeta_0$ and it follows that ${\rm sign}(\zeta_0) = {\rm sign}(\Delta(0)-1)$.

Proof. Since $\zeta \mapsto \Delta(\zeta)$ is an analytic function in $\Sigma$, each $\zeta \in \Omega$ is a normal eigenvalue. $\Delta(\zeta)$ is monotone decreasing from $\Delta(-\mu)$ to zero as $\zeta$ moves from $-\mu$ to $\infty$. Then there is a unique real root $\zeta_0 \in [-\mu, \infty)$ of $\Delta(\zeta) = 1$ if $\Delta(-\mu) \ge 1$, and it follows that ${\rm sign}(\zeta_0) = {\rm sign}(\Delta(0)-1)$. In this case, $\zeta_0$ is the dominant eigenvalue, that is, $\Re \zeta < \zeta_0$ for any $\zeta \in \Omega \setminus \{\zeta_0\}$. In fact, if $\zeta \in \Omega$ with $\Re \zeta \le \zeta_0$ and $\Im \zeta \neq 0$, we have $\Delta(\zeta_0) = 1 = |\Delta(\zeta)| < \Delta(\Re \zeta)$ holds, so $\Re \zeta < \zeta_0$ is obtained because of the monotonicity of $\Delta$ on $\mathbb R$.

Proposition 3. If $R_0 > 1$, then the disease-free steady state is unstable, while if $R_0 < 1$, then it is locally asymptotically stable.

Proof. It follows from (3.16) that $\omega_1(\mathcal A_0) \le -\mu$ and $D\mathcal F[0]$ is a compact perturbation, we have $\omega_1((\mathcal A + D\mathcal F [0])_0) = \omega_1(\mathcal A_0) \leq -\mu$ (^[28,30]). Therefore, ${\rm sign}(\omega_0(\mathcal B_0)) = {\rm sign}(s(\mathcal B_0))$. From the above Lemma 5, we know that $s(\mathcal B_0) < 0$ if $\Delta(0) < 1 \le \Delta(-\mu)$. If $\Delta(-\mu) < 1$, there is no eigenvalue in $\Sigma$ because $|\Delta(\zeta)| \le \Delta(\Re \zeta) \le \Delta(-\mu) < 1$, so we have $s(\mathcal B_0) < 0$ if $\Delta(0) < 1$. On the other hand, if $\Delta(0) > 1$, it follows from Lemma 5 that there exists a dominant positive eigenvalue. Thus we can conclude that the disease-free steady state is unstable if $\Delta(0) > 1$, while if $\Delta(0) < 1$, it is locally asymptotically stable. Finally we notice that $\Delta(0) = R_0$.

5.   Existence and non-existence of endemic steady states

Let $(S^*, I^*, r^*(\cdot))$ be a steady state. Then we have

It is easy to see that $S^*$ and $I^*$ are expressed by $\lambda^*$ and $r^*$ as follows:

Solving the equation of $r^*$ in (5.1), we obtain

Substituting the above expression into the boundary condition of $r^*$, it follows that

where for $\sigma < \eta$, we define

Then $p(\zeta, \eta)L(\sigma, \eta; \lambda^*)$ gives the probability density that asymptomatic individuals with recovery-age $\sigma$ survive to the age of $\eta$ and then returns to recovery-age $\zeta$. Since $\int_\sigma^\infty p(0, \eta)L(\sigma, \eta; \lambda^*)d\eta < 1$, we have

Then we obtain an integral equation with respect to $r^*$ as

which can be expressed as follows:

where

The integral equation (5.8) can be seen as the abstract equation with respect to $u \in X_+: = L^1_+(\mathbb R_+)$:

where $f_{\lambda^*}(\tau) = f(\tau; \lambda^*)$, $V_{\lambda^*}$ and $F_{\lambda^*}$ are bounded linear operators from $X_+$ into itself given by

Since $V_{\lambda^*}$ is a Volterra operator, its spectral radius is zero, so the positive operator $(I-V_{\lambda^*})^{-1}$ exists. Thus the equation (5.9) is rewritten as

Since the operator $((I-V_{\lambda^*})^{-1}F_{\lambda^*})$ is a Volterra operator again, the solution of (5.10) is obtained by the iteration:

Substituting the expression (5.11) into (5.2), we obtain the following one-dimensional equation with respect to unknown $\lambda^*$:

where

From the above characteristic equation, the following statement is obtained.

Proposition 4. If $R_0 > 1$, then there exists at least one endemic steady state.

Proof. Let $H(\lambda^*)$ be the right hand side of (5.12). First, we prove $H(0) = R_0$. For this purpose the value of $\frac{r^*(\tau)}{\lambda^*}$ at $\lambda^* = 0$ needs to be calculated. Since the operator $F_0$ is a zero operator, all of the terms for $k\geq 1$ in (5.11) are zero. Then we obtain

from which it follows that $H(0) = R_0$. Next we prove $\lim_{\lambda^*\rightarrow \infty}H(\lambda^*) = 0$. First, observe that

By using the equality obtained by dividing (5.9) by $\lambda^*$ the following representation follows:

where each integral of the right hand side of (5.13) can be estimated as follows:

Then $\lim_{\lambda^*\rightarrow \infty} H(\lambda^*) = 0$ holds. From the above, because of the continuity of $H$, it turns out that the equation $1 = H(\lambda^*)$ has at least one positive solution if $R_0 > 1$.

As is seen in the next section, in our basic model (2.1), $R_0 < 1$ does not necessarily imply non-existence of endemic steady states. However, we can formulate a sufficient condition for the global stability of the disease-free steady state implying the non-existence of endemic steady states. Let us define a time-dependent reproduction number $R_e(t)$ at time $t$ as

which is a kind of the effective reproduction number, because it is defined in a partially immunized host population. For the use and the definition of effective reproduction number, the reader may consult ^[4].

From the above definition, we have

Then we can expect that the disease will be eradicated if the effective reproduction number $R_e(t)$ is suppressed below the unity and it holds that

which is a criteria for disease prevention.

Proposition 5. Suppose that $(5.16)$ holds. Then the disease-free steady state is globally asymptotically stable, that is, there is no endemic steady state.

Proof. It follows from (5.15) that

from which we know that (5.16) implies that $\lim_{t \to \infty}I(t) = 0$. Next let $U(t): = \int_{0}^{\infty}r(t, \tau)d\tau$. Observe that

Then we have

which implies that $\lim_{t \to \infty}U(t) = 0$. Then $(S(t), I(t), r(t, \tau))$ converges to $(N, 0, 0)$ as $t \to \infty$, so the disease-free steady state is globally asymptotically stable.

Corollary 1. If

then the disease-free steady state is globally asymptotically stable.

Proof. Observe that

where we know that $S(N-S) \le N^2/4$. If (5.17) holds, we obtain

which implies (5.16).

6.   Bifurcation of endemic steady states from the disease-free steady state

We here consider the bifurcation problem of endemic steady states from the disease-free steady state when the basic reproduction number $R_0$ crosses the unity. Especially, the direction of bifurcation is biologically important because in the case of backward bifurcation, there exist multiple endemic steady states even when $R_0 < 1$, so $R_0 < 1$ is not a sufficient condition to eradicate the disease, whereas if forward bifurcation occurs at $R_0 = 1$ and there exists no endemic steady state for $R_0 < 1$, lowering $R_0$ less than unity becomes the almighty policy for disease eradication. Moreover, we can conclude from the Factorization Theorem ^[31] that the forwardly bifurcated solution is locally asymptotically stable, while the backwardly bifurcated solution is unstable as long as $|R_0-1|$ is small enough.

First, according to ^[1], we apply the Lyapunov-Schmidt type arguments to the one-parameter bifurcation model. Let us consider the nonlinear operator $\mathcal G \colon Z_0 \times \mathbb R \rightarrow Z$ with a parameter $\bar{\beta}$ defined as

for $u = {}^T(0, u_1, u_2(\cdot))\in Z_0$, where

where $\beta_{10}$ and $\beta_{20}$ are normalized transmission coefficients such that

Then it follows that $\bar{\beta} = R_0$ under the above assumption.

Then the system (2.1) is rewritten as the following abstract differential equation:

This equation has a trivial equilibrium $0\in Z_0$, which corresponds to the disease-free steady state. Let $D_u \mathcal G [u, \bar{\beta}]$ and $D^2_{uu} \mathcal G [u, \bar{\beta}]$ be the first and second order derivative of $\mathcal G$ with respect to $u$, respectively. Similarly, let $D^2_{u \bar{\beta}} \mathcal G [u, \bar{\beta}]$ be the second order derivative of $\mathcal G$ with respect to $u$ and $\bar{\beta}$. Note that $D_u \mathcal G [u, \bar{\beta}]$ is a linear map, while $D^2_{uu} \mathcal G [u, \bar{\beta}]$ and $D^2_{u \bar{\beta}} \mathcal G [u, \bar{\beta}]$ are bilinear maps. In order to study the direction of the bifurcation, we apply the following theorem:

Theorem 1 (^[1]). Suppose that the following conditions hold:

$A1$ Let $\mathcal L = D_u \mathcal G [0, 1]$ be the linearization around the zero equilibrium, evaluated at a critical value of parameter $\bar{\beta} = 1$, such that $\mathcal L$ is a closed operator with a simple isolated eigenvalue zero and remaining eigenvalues having negative real part. Let $\hat v_0$ be the unique (up to a constant) positive solution of $\mathcal L \xi = 0$.

$A2 \; \mathcal G(u, \bar{\beta})\in C^2(U_0\times I_0, Z)$ for some neighborhood $U_0 \times I_0$ of $(0, 1)$.

$A3$ Let $Z^*$ be the dual of $Z$ and $\langle \cdot, \cdot \rangle$ be the duality pairing between $Z$ and $Z^*$. Assume $\hat v_0^* \in Z^*$ is the unique (up to a constant) positive vector satisfying $\langle \mathcal L u, \hat v_0^* \rangle = 0$, $\forall u \in Z_0$, that is, $\dim(\ker \mathcal L^*) = 1$, where $\mathcal L^*$ is the adjoint of $\mathcal L$, and $\ker \mathcal L^* = {\rm span}\{\hat v_0^*\}$.

$A4 \; \langle D^2_{u \bar{\beta}} \mathcal G[0, 1](\hat v_0, 1), \hat v_0^* \rangle \neq 0$.

Then, the direction of the bifurcation is determined by the numbers

If $b > 0$, the bifurcation is backward if and only if $a > 0$ and forward if only if $a < 0$.

For our system, the derivatives of the nonlinear mapping $\mathcal G$ are calculated as the following:

for $\xi = {}^T(0, \xi_1, \xi_2(\cdot))\in Z_0$ and $\eta = {}^T(0, \eta_1, \eta_2(\cdot))\in Z_0$. We define the linear operator $\mathcal L : = D_u\mathcal G[0, 1]$. If $\zeta \in \mathbb C$ is an eigenvalue of $\mathcal L$, there exists $\xi = {}^T(0, \xi_1, \xi_2(\cdot))\in Z_0\setminus \{ 0 \}$ such that

Then the eigenvalue $\zeta$ is given as a root of the characteristic equation:

If $\bar{\beta} = 1$, we can use the similar discussion as in the Section 4 to show that the equation (6.7) has a single solution zero and each of the remaining solutions has negative real part. Solving $\mathcal L \xi = 0$, the positive eigenvector $\hat{v}_0 \in Z$ corresponding to the eigenvalue $0$ is obtained as

Next we calculate the adjoint operator of $\mathcal L$ and its eigenfunctional corresponding to the dominant eigenvalue zero. For $\xi = (0, \xi_1, \xi_2(\cdot)) \in Z_0$ and $\psi = {}^T (\psi_2(0), \psi_1, \psi_2(\cdot)) \in Z^* = \mathbb R \times \mathbb R \times L^\infty$,

hence the following representation is obtained:

where $Z^*$ denotes the adjoint space of $Z$. In order to find the eigenvector of $\mathcal L^*$ corresponding to the eigenvalue zero, we consider $\mathcal L^* \psi = 0$, that is,

By our assumption (6.3), we have

Then it follows that

Solving the ordinary differential equation in (6.10), we have

Therefore the unique (up to constant) eigenfunctional $\hat{v}_0^* \in Z^*$ of $\mathcal L^*$ corresponding to eigenvalue zero is obtained as

From the above, it is obvious that $b : = \langle D^2_{u \bar{\beta}} \mathcal G[0, 1] (\hat v_0, 1), \hat v_0^* \rangle > 0$, hence the bifurcation is backward if $a > 0$, while if $a < 0$, the bifurcation is forward.

Now calculate the value of $a: = \langle D^2_{uu} \mathcal G[0, 1](\hat v_0, \hat v_0), \hat v^*_0 \rangle$. Since

it follows that the value of $a$ satisfies

In order to see the meaning of the bifurcation condition by the parameter $a$, we adopt a more intuitive argument to find the direction of the bifurcation based on the characteristic equation (5.12). The direction of the bifurcation is determined by the sign of $H'(0)$ when $R_0 = 1$. In fact, if $R_0 < 1$ but $|R_0-1|$ is small enough, $H'(0) > 0$ implies the existence of positive solution for the equation (5.12) (Figure 1).

Then we can conclude that if $H'(0) > 0$, the bifurcation is backward, while it is forward if $H'(0) < 0$ at $R_0 = 1$. Now we prove the equivalence of both conditions.

Proposition 6. Suppose that $R_0 = 1$. Then it holds that ${\rm sign}(a) = {\rm sign}(H'(0))$.

Proof. First we calculate the derivative

using the equality (5.11). For the purpose of calculating the derivative at zero, we may neglect the second or higher terms of $\lambda^*$, hence

Since

holds, it follows that

Next,

and

are obtained. Finally we arrive at the following expression:

Now let us prove that

Let

where

Note that the following holds:

Then each $H_i$ can be calculated as follows:

Therefore we have

which implies the equality (6.18).

In order to see under what kind of biological conditions the backward bifurcation could occur, let us define a function:

Then $F(\eta)$ is the expected number of reproduction (reproductive value) for sub-clinically infected individuals at recovery-age $\eta$, especially $F(0)$ gives the reproduction number of newly recovered individuals. Note that our condition (6.3) is expressed as

Moreover we define

Then $G(\eta)$ denotes the average reproduction number of sub-clinical individuals boosted at recovery-age $\eta$.

Proposition 7. Under the assumption $(6.3)$, backward bifurcation of endemic steady states occurs at $R_0 = 1$ from the disease-free steady state if and only if

Proof. From the above arguments so far, it is sufficient to show that

It follows from (6.20) that

where

Therefore it is easy to see that (6.25) holds.

Corollary 2. Suppose that

Then the forward bifurcation of endemic steady states occurs at $R_0 = 1$ from the disease-free steady state. Especially one of the following typical conditions is satisfied, then (6.26) holds:

1. $F$ is constant,

2. $\sup_{\sigma \ge 0}F(\sigma) \le 1$,

3. $p(0, \eta) = 1$ and $F(0) \le 1+F(\eta)$ for all $\eta$.

Proof. From (6.25), it is clear that (6.26) is a sufficient condition in order to satisfy $H'(0) < 0$. If $F$ is constant, it follows from (2.4) that $G(\eta)-F(\eta) = 0$, hence we have (6.26). Next, from (2.4), we have $G(\eta) \le \sup_{\sigma \ge 0}F(\sigma)$, so (6.26) holds if $\sup_{\sigma \ge 0}F(\sigma) \le 1$. Finally, if $p(0, \eta) = 1$ for all $\eta$, that is, all newly reinfected individual has recovery age zero, it follows that $G(\eta) = F(0)$, so again (6.26) holds, if $F(0) \le 1+F(\eta)$ for all $\eta$.

Corollary 3. Suppose that $\int_\sigma^\infty p(\sigma, \eta) d\eta \leq 1$ for all $\sigma > 0$. Then, the endemic steady state bifurcates forwardly from the disease-free steady state at $R_0 = 1$.

Proof. Since

we have $\beta_{10} = \mu+q -q F(0)$. Since $\beta_{10}$ is nonnegative, $F(0)$ should satisfy

Moreover, since $p(0, \eta) = 1 - \int_0^\eta p (\sigma, \eta) d\sigma$ for all $\eta > 0$, we have

for all $\eta > 0$. From (6.27) and (6.28), we have

Moreover, we have

Since

we have from (6.29) and (6.30) that

Thus, the inequality (6.24) does not hold.

Corollary 4. Suppose that $F$ is monotone on $[0, \infty)$. Then, the endemic steady state bifurcates forwardly from the disease-free steady state at $R_0 = 1$.

Proof. First, suppose that $F$ is monotone nonincreasing on $[0, \infty)$. Since $p(0, \eta) = 1-\int_0^\eta p(\sigma, \eta) d\sigma$ for all $\eta > 0$, we have

Hence, we have

Thus, the inequality (6.24) in Proposition 7 does not hold. Next, suppose that $F$ is monotone nondecreasing on $[0, \infty)$. We then have

Hence, we have

Thus, (6.24) in Proposition 7 does not hold.

Proposition 8. $F(0) \le1$ if and only if

And the condition $(6.31)$ is satisfied if $\sup_{\tau \ge 0}\beta_{20}(\tau) \le \beta_{10}$. Then if $p(0, \eta) = 1$ for all $\eta \ge 0$, the condition $(6.31)$ is a sufficient condition to show that the forward bifurcation of endemic steady states occurs at $R_0 = 1$ from the disease-free steady state.

Proof. It follows from (6.22) that $\beta_{10}-\mu = q(1-F(0))$. Then $F(0) \le 1$ implies $\beta_{10} \ge \mu$. Hence we have

so $F(0) \le 1$ implies (6.31). Conversely suppose that (6.31) holds. If we assume that $F(0) > 1$, it follows that $\beta_{10}-\mu < 0$, which implies

which contradicts $F(0) > 1$, then it must hold that $F(0) \le 1$. Therefore, the condition (6.31) implies $F(0) \le 1$, and it also implies the condition 3 of Corollary 2, so the bifurcation is forward. Finally note that if $\sup_{\tau \ge 0}\beta_{20}(\tau) \le \beta_{10}$, then we can observe that

Then the condition (6.31) is satisfied.

From the above propositions, we know that age-dependent variability of the transmission coefficient for subclinical status is needed to produce possibility of backward bifurcation of endemic steady states at $R_0 = 1$ from the disease-free steady state. Yang and Nakata ^[32] proved that if $p(0, \eta) = 1$ for all $\eta \ge 0$ and the condition (6.31) is satisfied, there exists a unique endemic steady state exists if and only if $R_0 > 1$. Therefore, if all newly boosted infecteds have the recovery-age zero, subcritical endemic steady state does not exist if there is no enhancement of the infectivity such that $\sup_{\tau \ge 0}\beta_{20}(\tau) > \beta_{10}$ among subclinical infecteds.

7.   Discussion

In this paper, we have extended the Aron epidemic model so that the immunity clock can be reset to any time less than the recovery-age at which reinfection occurs. If we assume that the immunity level is determined by the immunity clock, our model assumption allows that newly boosted (reinfected) individuals could get any level of immunity, which is a most flexible extension for the epidemic with boosting and waning of immune status.

We have established the mathematical well-posedness of our formulation and have shown that the initial invasion condition and the endemicity can be characterized by the basic reproduction number $R_0$. Our focus is the bifurcation analysis of endemic steady states. Based on Lyapunov–Schmidt type arguments, we have determined the direction of bifurcation that endemic steady states bifurcate from the disease-free steady state when the basic reproduction number passes through the unity. We have given a necessary and sufficient condition for backward bifurcation to occur. Although we have not yet discovered biologically reasonable, concrete numerical examples of subcritical endemic steady states for our model^*, we have cleared what kind of mathematical conditions are needed to create the backward bifurcation from the disease-free steady state under the reset mechanism of the immunity clock.

^*In numerical simulation, we can use a standard method such as Gauss-Laguerre quadrature to calculate numerical integrals on infinite intervals.

Mathematically speaking, even in case that the forward bifurcation occurs, it is still a difficult problem to show conditions for the uniqueness of endemic steady state. It is a future challenge to investigate the global dynamics of the epidemic model with boosting and waning of immune status. From application point of view, it would be also important to extend our recovery-age dependent model to take into account the chronological-age structure of host population. As is suggested in ^[4], the chronological-age dependent model is essential to understand the age-specific prevalence dynamics for the epidemic with boosting and waning of immune status.

Finally, we can again stress importance to consider the waning and boosting dynamics in the real-world application. As is strongly expected in the pandemic of COVID-19, we believe that vaccination is the most effective tool to control the pandemic. However, if the effect of vaccination is not necessarily lifelong and it will wane as time evolves, so we may well have to boost the immune status by repeating vaccination to keep the herd immunity. To make an appropriate immunization policy, we may need to know the waning and boosting dynamics of immune status under any vaccination policy. It is a future challenge to extend our model by taking into account possible vaccination policy.

Acknowledgement

Kento Okuwa and Hisashi Inaba are supported by the Japan Society for the Promotion of Science (JSPS), Grant-in-Aid for Scientific Research (C) (No.16K05266). Toshikazu Kuniya is supported by JSPS Grant-in-Aid for Early-Career Scientists (No.19K14594).

Conflict of interest

There is no conflict of interest.