1. Introduction

Hepatitis B virus infection is a major public-health concern, both in developing countries and immigrant communities in developed countries. An estimated two billion people worldwide are infected with the virus, with about 378 million chronic carriers worldwide and approximately 620, 000 deaths each year [8]. Each year, about 4.5 million new infections occur, of which a quarter progress to liver disease [8]. Approximately 45% of the world's population live in areas with high prevalence of chronic HBV infection [16]. Prevalence ranges between 2% in low-endemic countries and 8% in high-endemic countries [8]. In high-endemic countries, the lifetime risk for HBV infection is greater than 60% [16].

HBV infection exhibits an acute infection stage and a chronic liver infection, characterized by persistent serum level of HBV surface antigen (HBsAg), IgG anti-core antigen (anti-HBc) and HBV DNA [4]. Acute disease usually occurs when the immune response is well preserved, while patients with an immunodeficiency are more likely to develop a chronic disease, in turn becoming a source for new infections [8]. Chronic infection may later develop into serious disease such as cirrhosis or liver cancer, causing major morbidity and mortality [14].

Safe and effective HBV vaccines have been commercially available since 1962 [16], being also the first vaccine to protect against cancer [14]. More than 150 countries have HBV vaccine immunization programs, with routine infant vaccination designated a high priority in all countries [16]. However, coverage in developing countries with high endemicity is limited due to high cost and social hurdles [8]. This is despite the cost-effectiveness of the vaccine in both high- and low-endemic countries [14]. In the developed world, the hepatitis B vaccine is one that was targeted by anti-vaccine campaigners, due to the presence of mercury-based thimerosal; however, by March 2001, thimerosal had been removed from all childhood vaccines [10]. An investigation of the literature established that no credible studies demonstrated a link between hepatitis B vaccination and autism [21]. Nevertheless, vaccination in adults in the United States decreased by 2.1% between 2012 and 2013, with 25% of adults $\geq 19$ years currently vaccinated [25].

The probability of becoming chronically infected is positively correlated to the age of the infected host [14,6,19,9], while progression of acute-stage and chronic-stage HBV is related to the time since infection. Furthermore, the infectivity of patients varies at different ages during the infectious period. Thus, mathematical models that can examine infection-age structure can be useful in investigating the consequences of infection age on HBV transmission dynamics and disease progression. Continuous age-structured models usually lead to partial differential equation (PDE) formulation. Population migration occurs at regional, national and global scales due to various factors, such as economic development. Labour workers transition from rural regions to large cities in both developing and industrialized countries, while immigrants and refugees migrate from developing countries to immigrant-receiving countries [20].

Although dynamical analysis of epidemic models with age structures is particularly challenging, there has been recent progress in global analysis [13,22,2,7,15,17,18]. Age-structured models have also been developed to study the epidemiology of HBV infection [6,19,27,29,26]. Medley et al. observed a feedback mechanism that determines the prevalence of HBV infection, using an ODE model to relate the rate of transmission, average age at infection and age-related probability of developing carriage following infection [19]. Based on sero-survey data in China, Zhao et al. [27] constructed an age-structured HBV model to evaluate the long-term effectiveness of vaccination programmes. Zou et al. in [29] proposed a full PDE model that incorporated multiple age structures to study the transmission dynamics of HBV and analyzed the existence and stability of the disease-free and endemic equilibria. In [26], age of infection and nontrivial vertical transmission were incorporated into the model to study the possible effects of variable infectivity on HBV dynamics.

Brauer and van den Driessche [3] studied an SIS model with a constant flow of infective immigrants into the infectious compartment. Such immigration models do not contain an infection-free equilibrium and consequently have no reproduction number. Guo and Li [11] generalized the immigration model to a high-dimensional SEIR model with constant immigration into each compartment, and the global stability of unique endemic equilibrium was shown by a global Lyapunov function. McCluskey [17,18] proved the global stability of the unique endemic equilibrium by a Lyapunov functional for an SEIR model with age structure consisting of a latent compartment, an infectious compartment and variable-age immigration terms. Guo and Li [12] demonstrated that the immigration inflow is vital to the dynamical behaviour: small variations in the immigration term into infectious compartment and/or latent compartment can profoundly change the resultant dynamics.

---

2. Model formulation

Based on the characteristics of HBV transmission and progression, we divide the population into five classes: susceptible, exposed, acute infection, chronic carriers and immunized. Let $S(t)$, $E(t)$ and $R(t)$ denote the population of susceptible, exposed and immunized individuals, respectively, at time $t$. Let $i(a, t)$ and $c(a, t)$ denote the densities of acute HBV infections and chronic HBV carriers with infection age $a$ at time $t$. The acute-infected population is thus $I(t) = \int_0^{a_1}i(a, t)da$, while $C(t) = \int_{a_1}^\infty c(a, t)da$ represents the total number of chronic HBV carriers at time $t$. The critical infection age $a_1$ is set at three months, representing the average duration of time in acute infection before the onset of chronic stage.

Vertical transmission from mother to child is also incorporated into our model. For simplicity, we assume all newborns are vaccinated with the same efficacy and that all the neonates who acquire HBV infection perinatally become chronic carriers, due to the high risk (up to 90%) of becoming chronic for this group. The model flow diagram is shown in Figure 1. Based on the natural history of HBV transmission and the assumptions of infection age, we have the following system of combined ordinary and partial differential equations:

$\begin{split} S'(t) = {}&\Lambda_S-b\omega\int_{a_1}^\infty v(a)c(a, t)da -(\mu_S+p)S(t)\\ & -\int_0^{a_1}\beta(a)i(a, t)S(t)da-\int_{a_1}^\infty \alpha\beta(a)c(a, t)S(t)da, \\ E'(t) = {}&\Lambda_E+ \int_0^{a_1}\beta(a)i(a, t)S(t)da \end{split}$

(1)

$\begin{split} &+\int_{a_1}^\infty \alpha\beta(a)c(a, t)S(t)da-(\mu_E+\sigma)E(t), \\ \dfrac{\partial i(a, t)}{\partial a}+ \dfrac{\partial i(a, t)}{\partial t} = {}&\Lambda_i(a)-(\mu_i(a)+\gamma_1(a))i(a, t), 0 <a\leq a_1\\ \dfrac{\partial c(a, t)}{\partial a}+ \dfrac{\partial c(a, t)}{\partial t} = {}&\Lambda_c(a)-(\mu_c(a)+\gamma_2(a)+\theta(a))c(a, t), a_1 <a <\infty \\ R'(t) = {}&\Lambda_R+ b(1-\omega)+\int_0^{a_1}(1-q(a))\gamma_1(a)i(a, t)da\\ &+\int_{a_1}^\infty \gamma_2(a)c(a, t)da -\mu_R R(t)+pS(t), \\ \end{split}$

with the boundary conditions

$\begin{align*} i(0, t)& = \sigma E(t), \\ c(a_1, t)& = \int_0^{a_1}q(a)\gamma_1(a)i(a, t)da+b\omega\int_{a_1}^\infty v(a)c(a, t)da \end{align*}$

and initial conditions

$S\left( 0 \right){\rm{ }} = {S_0}, \;\;\;E\left( 0 \right){\rm{ }} = {E_0}, \;\;\;i(a, 0){\rm{ }} = {i_0}(a), \;\;\;c(a, 0){\rm{ }} = {c_0}(a), \;\;\;R\left( 0 \right){\rm{ }} = {R_0}.$

The age-dependent function $\beta(a)$ describes the rate of infectiousness as disease progresses within an infected individual. The vertical transmission rate is

$\begin{align*} V& = \int_{a_1}^{\infty}\nu(a)c(a, t)da. \end{align*}$

The movement from acute individuals to carriers is given by

$\begin{align*} I_c& = \int_0^{a_1}q(a)\gamma_1(a)i(a, t)da. \end{align*}$

The recovery rate is

$\begin{align*} I_R& = \int_0^{a_1}(1-q(a))\gamma_1(a)i(a, t)da. \end{align*}$

The definitions of the other parameters in system (1) are listed in the Table 1. See also Zhang and Xu [26]. Note that unsuccessfully immunized birth is included in $\Lambda_S$ (i.e., $\Lambda_S>b\omega$) and the vertical transmission is given by $b\omega V$. From the definition of $v(a)$, we have $V <1$, so it follows that $\Lambda_S>b\omega V$.

In order to simplify the analysis, we made the following assumptions.

(H1) $v$, $\beta$, $\mu_i$, $\mu_c$, $\gamma_1$, $\gamma_2$, $q$, $\theta\in L^{\infty}(\mathbb{R}_{\geq 0}, \mathbb{R}_{\geq 0})$.

Let ${\underline \mu _i},{\underline \mu _c},{\rm{ }}\underline \beta$ be the essential infimums of $\mu_i, \mu_c, \beta$, respectively, and $\bar{\mu}_i, \bar{\mu}_c, \bar{q}, \bar{\gamma}_1, \bar{v}$ be the respective essential supremums of $\mu_i, \mu_c, q, \gamma_1$ and $v$.

(H2) ${\underline \mu _i},{\underline \mu _c}>0$ and $q, \gamma_1, v$ are Lipschitz continuous with Lipschitz coefficients $L_q, L_{\gamma_1}$ and $L_v$, respectively.

(H3) $\Lambda_i, \Lambda_c\in L^1(\mathbb{R}_{\geq 0}, \mathbb{R}_{\geq 0})$.

(H4) The constant parameters $\Lambda_S, b, \mu_S, \mu_E, \mu_R, \omega, \sigma, \alpha$ are positive, and $\Lambda_E, \Lambda_R, p$ are nonnegative.

(H5) The initial conditions satisfy $S_0\geq 0$, $E_0\geq 0$, $R_0\geq 0$, $i_0(a), c_0(a)\in L^1(\mathbb{R}_{\geq 0}, \mathbb{R}_{\geq 0})$.

(H6) The supports of $\beta, q\gamma_1$ and $\Lambda_i+\Lambda_c$ have positive measure. Let $\tilde{\Lambda}_i = \int_0^{a_1} \Lambda_i(a)da$ and $\tilde{\Lambda}_c = \int_{a_1}^{\infty}\Lambda_c(a)da$, then $\tilde{\Lambda}_i+\tilde{\Lambda}_c>0$.

(H7) essential infimum (support($\Lambda_i$))$<$ essential supremum (support($q\gamma_1$)),

essential infimum (support($\Lambda_j$))$<$ essential supremum (support($\beta$))

for $j = i$ or $c$.

To simplify expressions, we introduce the following notations

$\begin{align*} \pi_1(a)& = e^{-\int_0^a(\mu_i(s)+\gamma_1(s))ds}& a&\in[0, a_1], \\ \pi_2(a)& = e^{-\int_{a_1}^a(\mu_c(s)+\gamma_2(s)+\theta(s))ds}& a&\in[a_1, \infty), \end{align*}$

where $\pi_1(a)$ is the age-specific survival probability of being acutely infected and $\pi_2(a)$ is the age-specific survival probability of being a chronic carrier. Based on the boundary and initial conditions and methods in [24], integrating $i(a, t)$ and $c(a, t)$ equations in system (1) yields

$\label{a} i(a, t) = \begin{cases} i(0, t-a)\pi_1(a)+\int_0^a \Lambda_i(s)\frac{\pi_1(a)}{\pi_1(s)}ds, &t>a, ~ a\in[0, a_1], \\ i_0(a-t)\frac{\pi_1(a)}{\pi_1(a-t)}+\int_{a-t}^a \Lambda_i(s)\frac{\pi_1(a)}{\pi_1(s)}ds, &t\leq a, ~ a\in[0, a_1], \end{cases}$

(2)

and $c(a, t) =$

$\begin{cases} c(a_1, a_1+t-a)\pi_2(a)+\int_{a_1}^a \Lambda_c(s)\frac{\pi_2(a)}{\pi_2(s)}ds, &t+a_1>a, a\in[a_1, \infty), \\ c_0(a-t-a_1)\frac{\pi_2(a)}{\pi_2(a-t-a_1)}+\int_{a-t-a_1}^a \Lambda_c(s)\frac{\pi_2(a)}{\pi_2(s)}ds, &t+a_1\leq a, a\in[a_1, \infty). \end{cases}$

(3)

By classical existence and uniqueness results for functional differential equations, there exists a unique solution for the integro-differential system (1) in which $i(a, t)$ and $c(a, t)$ are substituted for the expressions (2) and (3), respectively.

For (2) and (3), it is easy to see that $i(a, t)$ and $c(a, t)$ remain nonnegative for any nonnegative initial value. Furthermore, if there exists a $t^*$ such that $S(t^*) = 0$ and $S(t)>0$ for $0 <t <t^*$, then, from the $S$ equation of (1), we have $S'(t^*) = \Lambda_S-b\omega\int_{a_1}^\infty v(a)c(a, t)da>0$, which implies that $S(t)\geq 0$ for all $t\geq0$, noting that unsuccessfully immunized birth $b\omega$ is included in $\Lambda_S$. Similarly, it can be shown that $E(t)\geq 0$ for all $t\geq 0$ and all nonnegative initial values. Let

$\mathcal{Y} = \mathbb{R}_{\geq 0}^2\times L^1(\mathbb{R}_{\geq 0}, \mathbb{R}_{\geq 0}) \times L^1(\mathbb{R}_{\geq 0}, \mathbb{R}_{\geq 0})\times\mathbb{R}_{\geq 0}$

be the state space of system (1). Then $\mathcal{Y}$ is positively invariant and there is a continuous semiflow defined by

$\Phi_t:\mathcal{Y}\rightarrow\mathcal{Y}.$

Denote $u(t) = (S(t), E(t), i(., t), c(., t), R(t))\in\mathcal{Y}$, which is endowed with the following norm:

$\| u\| = S+E+\int_0^{a_1}i(a, t)da+\int_{a_1}^\infty c(a, t) da+R.$

Notice that the variable $R$ does not appear in other equations in (1); thus the equation of $R$ can be ignored when studying the model dynamics, and the reduced system has the same dynamical behavior as the original system.

---

3. Boundedness

Let $N(t)$ denote the total population size at time $t$. Then we have

$N(t) = S(t)+E(t)+\int_0^{a_1}i(a, t)da+\int_{a_1}^\infty c(a, t)da+R(t) = \|u\|.$

In the following, we establish that $N(t)$ is bounded and the generated semiflow $\Phi_t$ is point dissipative.

From (2), we can obtain

$\begin{align*} \int_0^{a_1}i(a, t)da = {}&\int_0^t i(a, t)da+\int_t^{a_1}i(a, t)da \\ = {}&\int_0^t i(0, t-a)\pi_1(a)da+\int_0^t\int_0^a\Lambda_i(\tau)\frac{\pi_1(a)}{\pi_1(\tau)}d\tau da \\ &+\int_t^{a_1}i_0(a-t)\frac{\pi_1(a)}{\pi_1(a-t)}da+\int_t^{a_1}\int_{a-t}^a \Lambda_i(\tau)\frac{\pi_1(a)}{\pi_1(\tau)}d\tau da. \end{align*}$

Changing the order of integration for two double integrals in the above, we have

$\begin{align*} \int_0^t\int_0^a\Lambda_i(\tau)\frac{\pi_1(a)}{\pi_1(\tau)}d\tau da = {}&\int_0^t\int_{\tau}^t\Lambda_i(\tau)\frac{\pi_1(a)}{\pi_1(\tau)}da d\tau, \\ \int_t^{a_1}\int_{a-t}^a \Lambda_i(\tau)\frac{\pi_1(a)}{\pi_1(\tau)}d\tau da = {}&\int_0^t\int_t^{\tau+t} \Lambda_i(\tau)\frac{\pi_1(a)}{\pi_1(\tau)}da d\tau\\ & +\int_t^{a_1}\int_{\tau}^{\tau+t} \Lambda_i(\tau)\frac{\pi_1(a)}{\pi_1(\tau)}dad\tau\\ &-\int_{a_1-t}^{a_1}\int_{a_1}^{\tau+t} \Lambda_i(\tau)\frac{\pi_1(a)}{\pi_1(\tau)}da d\tau. \end{align*}$

Making substitutions $\tau = t-a$ and $\tau = a-t$ in the remaining two integrals gives

$\begin{align*} \int_0^{a_1}i(a, t)da = {}&\int_0^t i(0, \tau)\pi_1(t-\tau)d\tau+\int_0^{a_1-t}i_0(\tau)\frac{\pi_1(\tau+t)}{\pi_1(\tau)}d\tau \\ &+\int_0^{a_1}\int_{\tau}^{\tau+t} \Lambda_i(\tau)\frac{\pi_1(a)}{\pi_1(\tau)}dad\tau -\int_{a_1-t}^{a_1}\int_{a_1}^{\tau+t} \Lambda_i(\tau)\frac{\pi_1(a)}{\pi_1(\tau)}d\tau da. \end{align*}$

Thus differentiating $\int_0^{a_1}i(a, t)da$ with respect to $t$ leads to

$\begin{align*} \frac{d}{dt}\int_0^{a_1}i(a, t)da = {}&i(0, t)-i_0(a_1-t)\frac{\pi_1(a_1)}{\pi_1(a_1-t)}+\int_0^t i(0, \tau)\pi_1^{'}(t-\tau)d\tau\\ &+\int_0^{a_1-t}i_0(\tau)\frac{\pi_1^{'}(\tau+t)}{\pi_1(\tau)}d\tau+\int_0^{a_1}\Lambda_i(\tau)\frac{\pi_1(\tau+t)}{\pi_1(\tau)}d\tau\\ &-\frac{d}{dt}\int_{a_1-t}^{a_1}\int_{a_1}^{\tau+t} \Lambda_i(\tau)\frac{\pi_1(a)}{\pi_1(\tau)}da d\tau. \end{align*}$

We have

$\begin{align*} \int_0^{a_1-t}i_0(\tau)\frac{\pi_1^{'}(\tau+t)}{\pi_1(\tau)}d\tau = \int_t^{a_1}i_0(a-t)\frac{\pi_1^{'}(a)}{\pi_1(a-t)}da \end{align*}$

and

$\begin{align*} \int_0^t\int_0^a\Lambda_i(\tau)\frac{\pi_1^{'}(a)}{\pi_1(\tau)}d\tau da&+\int_t^{a_1}\int_{a-t}^a \Lambda_i(\tau)\frac{\pi_1^{'}(a)}{\pi_1(\tau)}d\tau da +\int_{a_1-t}^{a_1}\int_{a_1}^{\tau+t} \Lambda_i(\tau)\frac{\pi_1^{'}(a)}{\pi_1(\tau)}da d\tau \\ & = \int_0^{a_1}\int_{\tau}^{\tau+t} \Lambda_i(\tau)\frac{\pi_1^{'}(a)}{\pi_1(\tau)}da d\tau \\ & = \int_0^{a_1}\Lambda_i(\tau)\frac{\pi_1(\tau+t)-\pi_1(\tau)}{\pi_1(\tau)}d\tau\\ & = \int_0^{a_1}\Lambda_i(\tau)\frac{\pi_1(\tau+t)}{\pi_1(\tau)}d\tau-\int_0^{a_1}\Lambda_i(\tau)d\tau. \end{align*}$

Noticing that $\pi_1'(a) = -(\mu_i(a)+\gamma_1(a))\pi_1(a)$, it follows that

$\begin{split} \frac{d}{dt}\int_0^{a_1}i(a, t)da = {}&i(0, t)-i_0(a_1-t)\frac{\pi_1(a_1)}{\pi_1(a_1-t)}+\int_0^t i(0, \tau)\pi_1^{'}(t-\tau)d\tau \\ &+\int_t^{a_1}i_0(a-t)\frac{\pi_1^{'}(a)}{\pi_1(a-t)}da+\int_0^t\int_0^a\Lambda_i(\tau)\frac{\pi_1^{'}(a)}{\pi_1(\tau)}d\tau da \\ &+\int_t^{a_1}\int_{a-t}^a \Lambda_i(\tau)\frac{\pi_1^{'}(a)}{\pi_1(\tau)}d\tau da+\int_0^{a_1}\Lambda_i(\tau)d\tau\\ &-\int_{a_1-t}^{a_1} \Lambda_i(\tau)\frac{\pi_1(\tau+t)}{\pi_1(\tau)}d\tau \\ = {}&i(0, t)-i_0(a_1-t)\frac{\pi_1(a_1)}{\pi_1(a_1-t)}+\tilde{\Lambda}_i -\int_{a_1-t}^{a_1} \Lambda_i(\tau)\frac{\pi_1(\tau+t)}{\pi_1(\tau)}d\tau \\ &-\int_0^{a_1}(\mu_i(a)+\gamma_1(a))i(a, t)da. \end{split}$

(4)

A similar calculation for (3) yields

$\begin{align*} \int_{a_1}^{\infty}c(a, t)da& = \int_{a_1}^{a_1+t}c(a, t)da+\int_{a_1+t}^{\infty}c(a, t)da\\ = {}&\int_0^tc(a_1, \tau)\pi_2(a_1+t-\tau)d\tau+\int_0^{\infty}c_0(\tau)\frac{\pi_2(t+a_1+\tau)}{\pi_2(\tau)}d\tau \\ &+\int_{a_1}^{\infty}\int_{\tau}^{\tau+t}\Lambda_c(\tau)\frac{\pi_2(a)}{\pi_2(\tau)}dad\tau, \end{align*}$

and hence

$\begin{split} \frac{d}{dt}\int_{a_1}^{\infty}c(a, t)da = {}&c(a_1, t)+\int_{a_1}^{a_1+t}c(a_1, a_1+t-a)\pi_2^{'}(a)da \\ &+\int_{a_1+t}^{\infty}c_0(a-t-a_1)\frac{\pi_2^{'}(a)}{\pi_2(a-t-a_1)}da\\ &+\int_{a_1}^{a_1+t}\int_{a_1}^a\Lambda_c(\tau)\frac{\pi_2^{'}(a)}{\pi_2(\tau)}d\tau da \\ &+\int_{a_1+t}^{\infty}\int_{a-t-a_1}^a\Lambda_c(\tau)\frac{\pi_2^{'}(a)}{\pi_2(\tau)}d\tau da+\int_{a_1}^{\infty}\Lambda_c(\tau)d\tau \\ = {}&c(a_1, t)-\int_{a_1}^{\infty}(\mu_c(a)+\gamma_2(a)+\theta(a))c(a, t)da+\tilde{\Lambda}_c. \end{split}$

(5)

Combining the equations gives

$\begin{align*} N'(t) = {}&\Lambda_S+\Lambda_E+\tilde{\Lambda}_i+\tilde{\Lambda}_c+\Lambda_R +b(1-\omega)-i_0(a_1-t)\frac{\pi_1(a_1)}{\pi_1(a_1-t)}\\ &-\int_{a_1-t}^{a_1} \Lambda_i(\tau)\frac{\pi_1(\tau+t)}{\pi_1(\tau)}d\tau -\mu_SS(t)-\mu_EE(t)-\mu_RR(t)\\ &-\int_0^{a_1}\mu_i(a)i(a, t)da-\int_{a_1}^{\infty}(\mu_c(a)+\theta(a))c(a, t)da\\ \leq{}&\Lambda^*-\mu^*N(t), \end{align*}$

where $\Lambda^* = \Lambda_S+\Lambda_E+\tilde{\Lambda}_i+\tilde{\Lambda}_c+\Lambda_R +b(1-\omega)$ and $\mu^* = \min\{\mu_S, \mu_E, \mu_R, {\underline \mu _i},{\underline \mu _c}\}$. Letting $N^* = \frac{\Lambda^*}{\mu^*}$, we have

$\displaystyle{\limsup\limits_{t\to\infty}} N(t)\leq N^*$

which implies that all solutions of system (1) are ultimately bounded. Moreover, when $N(t)>\frac{N^*}{\mu^*}$, we have $\frac{dN(t)}{dt} <0$, which implies that all solutions are uniformly bounded. Therefore, the solution semiflow $\Phi_t:\mathcal{Y}\rightarrow\mathcal{Y}$ is point dissipative. It follows that the set

$\Omega = \left\{u(t) = (S(t), E(t), i(., t), c(., t), R(t))\in\mathcal{Y}:\| u\| \leq \max\left\{N^*, N(0)\right\}\right\}$

is positively invariant and absorbing under the semiflow $\Phi_t$ on $\mathcal{Y}$.

Lemma 3.1. The unique solution semiflow $\Phi_t$ of system (1) is uniformly bounded and point dissipative in $\mathcal{Y}$.

Lemma 3.2. There exists $T, \delta>0$ such that $i(0, t), c(a_1, t)>\delta$ for all $t>T$.

Proof. First, since $\lim\limits_{t\to\infty}\sup N(t)\leq N^*$, then, for any $\epsilon >0$, there exists $t_1\geq 0$ such that $\int_0^{a_1}i(a, t)da <N^*+\epsilon$ and $\int_{a_1}^{\infty}c(a, t)da <N^*+\epsilon$ for all $t\geq t_1$. Thus, for sufficiently large $t$, from (1), we have

$\begin{align*} S'(t)\geq{}&\Lambda_S-b\omega\int_{a_1}^{\infty}v(a)c(a, t)da-\bar{\beta}S(t)\left(\int_0^{a_1}i(a, t)da+\alpha\int_{a_1}^{\infty}c(a, t)da\right)\\ &-(\mu_S+p)S(t)\\ \geq{}&\Lambda_S-b\omega\int_{a_1}^{\infty}v(a)c(a, t)da-S(t)\left[\bar{\beta}(N^*+\epsilon)(1+\alpha)+(\mu_S+p)\right], \end{align*}$

from which it follows that

$\liminf\limits_{t\to\infty} S(t)\geq\dfrac{\Lambda_S-b\omega\int_{a_1}^{\infty}v(a)c(a, t)da}{\bar{\beta}(N^*+\epsilon)(1+\alpha)+\mu_S+p}\equiv M_S.$

Next we have

$\begin{align*} E'(t)&\geq\Lambda_E+\underline \beta S(t)(\int_0^{a_1}i(a, t)da+\alpha\int_{a_1}^{\infty}c(a, t)da)-(\mu_E+\sigma)E(t)\\ &\geq\Lambda_E+\underline \beta M_S(\int_0^{a_1}i(a, t)da+\alpha\int_{a_1}^{\infty}c(a, t)da)-(\mu_E+\sigma)E(t), \end{align*}$

which implies

$\liminf\limits_{t\to\infty} E(t)\geq\frac{\Lambda_E+\underline \beta M_S(\int_0^{a_1}i(a, t)da+\alpha\int_{a_1}^{\infty}c(a, t)da)}{\mu_E+\sigma}\equiv M_E.$

Therefore there exists $t_2>t_1$ such that

$\label{e} i(0, t) = \sigma E(t)\geq\frac{\sigma M_E}{2}$

(6)

for all $t\geq t_2$. Note that (H7) implies that there exist $T\geq t_2$ and $\delta>0$ such that

$\begin{align*} \int_0^t\int_0^aq(a)\gamma_1(a)\Lambda_i(\tau)\frac{\pi_1(a)}{\pi_1(\tau)}d\tau da\geq\delta \end{align*}$

for all $t>T$, so it follows that

$\label{f} \begin{split} c(a_1, t)&\geq\int_0^{a_1}q(a)\gamma_1(a)i(a, t)da\\ &\geq\int_0^t\int_0^aq(a)\gamma_1(a)\Lambda_i(\tau)\frac{\pi_1(a)}{\pi_1(\tau)}d\tau da\\ &\geq\delta. \end{split}$

(7)

Then, when $t\geq T$, the result follows from (6) and (7).

---

4. Smoothness and attractiveness

Lemma 4.1. The semi-flow $\Phi:\mathbb{R}_{\geq 0}\times\mathcal{Y}\rightarrow\mathcal{Y}$ is asymptotically smooth if there are maps $\Theta, \Psi:\mathbb{R}_{\geq 0}\times\mathcal{Y}\rightarrow\mathcal{Y}$ such that

$\Phi(t, u) = \Theta(t, u)+\Psi(t, u)$

and, for any bounded closed set $C\subset\mathcal{Y}$ that is forward invariant under $\Phi$, the following holds:

$(a)$ $\lim\limits_{t\rightarrow\infty} {\rm{diam}}\, \Theta(t, C) = 0$,

$(b)$ there exists $t_C\geq 0$ such that $\Psi(t, C)$ has compact closure for each $t\geq t_C$.

This is a special case of Theorem 2.46 in [23]. Note that $\Phi_t = \Phi(t, u_0) = (S(t), E(t), i(., t), c(., t), R(t))$, where $u(0) = (S_0, E_0, i_0(a), c_0(a), R_0)$. Then, for $t\geq 0$, we define two flows $\Psi$ and $\Theta$ on $\mathcal{Y}$ so that $\Phi = \Psi+\Theta$. Let $\Psi(t, u_0) = (S(t), E(t), \tilde{i}(., t), \tilde{c}(., t), R(t))$ and $\Theta(t, u_0) = (0, 0, \tilde{i}_0(., t), \tilde{c}_0(., t), 0)$, where

$\begin{align*} \tilde{i}(a, t)& = \begin{cases} i(0, t-a)\pi_1(a), &t>a, \;\;a\in[0, a_1], \\ 0, &t\leq a, \;\;a\in[0, a_1], \end{cases}\\ \tilde{c}(a, t)& = \begin{cases} c(a_1, a_1+t-a)\pi_2(a), &t+a_1>a, a\in[a_1, \infty), \\ 0, &t+a_1\leq a, a\in[a_1, \infty), \end{cases}\\ \tilde{i}_0(a, t)& = \begin{cases} \int_0^a \Lambda_i(s)\frac{\pi_1(a)}{\pi_1(s)}ds, &t>a, \;\;a\in[0, a_1], \\ i_0(a-t)\frac{\pi_1(a)}{\pi_1(a-t)}+\int_{a-t}^a \Lambda_i(s)\frac{\pi_1(a)}{\pi_1(s)}ds, &t\leq a, \;\;a\in[0, a_1], \end{cases}\\ \tilde{c}_0(a, t)& = \begin{cases} \int_{a_1}^a \Lambda_c(s)\frac{\pi_2(a)}{\pi_2(s)}ds, &t+a_1>a, \\ c_0(a-t-a_1)\frac{\pi_2(a)}{\pi_2(a-t-a_1)}+\int_{a-t-a_1}^a \Lambda_c(s)\frac{\pi_2(a)}{\pi_2(s)}ds, &t+a_1\leq a, \end{cases} \end{align*}$

for $a\in[a_1, \infty)$. Let $C\subset\mathcal{Y}$ be bounded, with bound $K>N^*$. For $j = 1, 2$, denote $u_0^j = (S_0^j, E_0^j, i_0^j, c_0^j, R_0^j)\in C$, with corresponding solutions

$\Phi(t, u_0^j) = (S^j(t), E^j(t), i^j(., t), c^j(., t), R^j(t)).$

In the following, we evaluate the distance between $\Theta(t, u_0^1)$ and $\Theta(t, u_0^2)$, which corresponds to condition (a) in Lemma 4.1.

Noting that $i_0^j(a), c_0^j(a)\in L^1(\mathbb{R}_{\geq 0}, \mathbb{R}_{\geq 0})$ for $j = 1, 2$, let $\|.\|_1$ denote the standard norm on $L^1$. Since

$\tilde{i}_0^1(a, t)-\tilde{i}_0^2(a, t) = \begin{cases} 0, &t>a, a\in[0, a_1], \\ (i_0^1(a-t)-i_0^2(a-t))\frac{\pi_1(a)}{\pi_1(a-t)}&t\leq a,\;\; a\in[0, a_1], \end{cases}$

then

$\begin{align*} \|\tilde{i}_0^1(a, t)-\tilde{i}_0^2(a, t)\|_1& = \int_t^{a_1}\left|i_0^1(a-t)-i_0^2(a-t)\right|\frac{\pi_1(a)}{\pi_1(a-t)}da\\ & = \int_0^{a_1-t}\left|i_0^1(\tau)-i_0^2(\tau)\right|\frac{\pi_1(\tau+t)}{\pi_1(\tau)}d\tau\\ & = \int_0^{a_1-t}\left|i_0^1(\tau)-i_0^2(\tau)\right|e^{-\int_{\tau}^{\tau+t}(\mu_i(s)+\gamma_1(s))ds}d\tau\\ &\leq e^{-{\underline \mu _i}t}\int_0^{a_1}\left|i_0^1(\tau)-i_0^2(\tau)\right|d\tau\\ &\leq e^{-{\underline \mu _i}t}(\| i_0^1\|_1+\|i_0^2\|_1)\\ &\leq 2Ke^{-{\underline \mu _i}t}, \end{align*}$

which approaches zero as $t\rightarrow \infty$.

Similarly, $\|\tilde{c}_0^1(a, t)-\tilde{c}_0^2(a, t)\|_1\leq 2Ke^{-{\underline \mu _c}t}$, which leads to

$\|\Theta(t, u_0^1)-\Theta(t, u_0^2)\|\leq2K(e^{-{\underline \mu _i}t}+e^{-{\underline \mu _c}t})$

for all $t\geq 0$. Since $u_0^j$ is chosen arbitrarily in $C$, it follows that

$\text{diam}\, \Theta(t, C)\leq2K(e^{-{\underline \mu _i}t}+e^{-{\underline \mu _c}t})$

so condition (a) in Lemma 4.1 is satisfied.

In the following, we verify condition (b) in Lemma 4.1, which can be alternatively proved by verifying four conditions in Lemma 4.2 below. This shows that $\tilde{i}(a, t)$ and $\tilde{c}(a, t)$ remain in a subset of $L^1_{\geq 0}$, which has compact closure and is independent of $u_0$.

Lemma 4.2. $($Theorem B.2 in [23] for the case $\mathcal{S} = L^1(\mathbb{R}_{\geq 0}, \mathbb{R}_{\geq 0}))$ A set $\mathcal{S}\subseteq L^1_+(0, \infty)$ has compact closure if and only if the following conditions hold:

1. $\displaystyle\sup\limits_{f\in\mathcal{S}}\int_0^{\infty}f(a)da <\infty$,

2. $\displaystyle\lim\limits_{r\rightarrow\infty}\int_r^{\infty}f(a)da\rightarrow 0$ uniformly for $f\in\mathcal{S}$,

3. $\displaystyle\lim\limits_{h\rightarrow 0^+}\int_0^{\infty}\left|f(a+h)-f(a)\right|da\rightarrow 0$ uniformly for $f\in\mathcal{S}$,

4. $\displaystyle\lim\limits_{h\rightarrow 0^+}\int_0^hf(a)da\rightarrow 0$ uniformly for $f\in\mathcal{S}$.

Conditions 1, 2 and 4 are easy to show, since

$0\leq\tilde{i}(a, t) = \left\{ \begin{aligned} &i(0, t-a)\pi_1(a), &t>a, \;\;a\in[0, a_1]\\ &0, &t\leq a,\;\; a\in[0, a_1] \end{aligned} \right\} \leq\sigma Ke^{-{\underline \mu _i}a}$

and $0\leq\tilde{c}(a, t) =$

$\left\{ \begin{aligned} &c(a_1, a_1+t-a)\pi_2(a)&t+a_1>a, a\in[a_1, \infty)\\& 0, &t+a_1\leq a, a\in[a_1, \infty) \end{aligned} \right\}\leq (\bar{q}\bar{\gamma}_1+b\omega \bar{v}) Ke^{-{\underline \mu _c}a}.$

Next we show that Condition 3 is also satisfied. As $h$ tends to $0^+$, without loss of generality, we can assume $h\in(0, t)$. Thus

$\label{g} \begin{split} \int_0^{a_1}\left|\tilde{i}(a+h, t)-\tilde{i}(a, t)\right|da = {}&\int_0^{t-h}\left|i(0, t-a-h)\pi_1(a+h)-i(0, t-a)\pi_1(a)\right|da\\ ~&+\int_{t-h}^{t}\left|0-i(0, t-a)\pi_1(a)\right|da\\ ~\leq{}&\int_0^{t-h}\left|i(0, t-a-h)\pi_1(a+h)-i(0, t-a)\pi_1(a)\right|da\\&+\sigma Kh\\ ~\leq{}&\sigma Kh+\int_0^{t-h}i(0, t-a-h)\left|\pi_1(a+h)-\pi_1(a)\right|da\\ ~&+\int_0^{t-h}\left|i(0, t-a-h)-i(0, t-a)\right|\pi_1(a)da\\ ~\leq{}&\sigma Kh+\sigma K\int_0^{t-h}\left|\pi_1(a+h)-\pi_1(a)\right|da\\ ~&+\int_0^{t-h}\left|i(0, t-a-h)-i(0, t-a)\right|\pi_1(a)da. \end{split}$

(8)

Noting that $|E'(t)|$ is bounded by $L_E = \Lambda_E+\bar{\beta}K^2(1+\alpha)+(\mu_E+\sigma)K$, it follows that $E(\cdot)$ is Lipschitz for $t\geq 0$ with coefficient $L_E = L_E(K)$ and

$\left|i(0, t-a-h)-i(0, t-a)\right|\pi_1(a) = \sigma\left|E(t-a-h)-E(t-a)\right|\pi_1(a),$

which leads to

$\label{h} \int_0^{t-h}\left|i(0, t-a-h)-i(0, t-a)\right|\pi_1(a)da\leq \sigma L_E\frac{h}{{\underline \mu _i}}.$

(9)

Note that $\pi_1(a)$ is a positive decreasing function with $1$ as its supremum. Thus

$\begin{split} \int_0^{t-h}\left|\pi_1(a+h)-\pi_1(a)\right| & = \int_0^{t-h}\pi_1(a)da-\int_0^{t-h}\pi_1(a+h)da\\ & = \int_0^{t-h}\pi_1(a)da-\int_h^{t}\pi_1(a)da\\ & = \int_0^{t-h}\pi_1(a)da+\int_0^{h}\pi_1(a)da-\int_0^t\pi_1(a)da\\ & = \int_0^{h}\pi_1(a)da-\int_{t-h}^{t}\pi_1(a)da\\ &<\int_0^{h}\pi_1(a)da\\ &\leq h. \end{split}$

(10)

Substituting (9) and (10) into (8), we have

$\label{j} \int_0^{a_1}\left|\tilde{i}(a+h, t)-\tilde{i}(a, t)\right|da\leq 2\sigma Kh+\sigma L_Eh/{\underline \mu _i}.$

(11)

The constant in (11) is dependent on $K$ but independent of $u_0$. Thus (11) holds for all $u_0\in C$, which implies that $\tilde{i}(a, t)$ satisfies Condition 3 in Lemma 4.2. It remains in a pre-compact subset $C_K^i$ of $L_{\geq 0}^1$. A similar result can be obtained for $\tilde{c}(a, t)$.

Lemma 4.3. Both $\tilde{i}(a, t)$ and $\tilde{c}(a, t)$ remain in a pre-compact subset of $L_{\geq 0}^1$.

Proof. The result for $\tilde{i}(a, t)$ follows from (11). Note that

$\begin{split} &\int_{a_1}^{\infty}\left|\tilde{c}(a+h, t)-\tilde{c}(a, t)\right|da\\ = &\int_{a_1}^{a_1+t-h}\left|c(a_1, a_1+t-a-h)\pi_2(a+h)-c(a_1, a_1+t-a)\pi_2(a)\right|da\\ &+\int_{a_1+t-h}^{a_1+t}\left|0-c(a_1, a_1+t-a)\pi_2(a)\right|da\\ = &\int_{a_1}^{a_1+t-h}c(a_1, a_1+t-a-h)\left|\pi_2(a+h)-\pi_2(a)\right|da\\ &+\int_{a_1}^{a_1+t-h}\left|c(a_1, a_1+t-a-h)-c(a_1, a_1+t-a)\right|\pi_2(a)da\\ &+\int_{a_1+t-h}^{a_1+t}\left|0-c(a_1, a_1+t-a)\pi_2(a)\right|da. \end{split}$

(12)

Recalling the boundary conditions on $c(a_1, t)$, it follows that

$\label{l} \left|0-c(a_1, a_1+t-a)\pi_2(a)\right|\leq(\bar{q}\bar{\gamma}_1+b\omega \bar{v})K.$

(13)

Similar to (10), we have

$\label{m} \int_{a_1}^{a_1+t-h}c(a_1, a_1+t-a-h)\left|\pi_2(a+h)-\pi_2(a)\right|da\leq(\bar{q}\bar{\gamma}_1+b\omega \bar{v})Kh.$

(14)

For the second integral in (12), we have

$\begin{align*} \big|c(a_1, a_1+t-a&-h)-c(a_1, a_1+t-a)\big|\\ \leq{}&\int_0^{a_1}\left|q(a)\gamma_1(a)i(a, a_1+t-a-h)-q(a)\gamma_1(a)i(a, a_1+t-a)\right|da\\ &+b\omega\int_{a_1}^{\infty}\left|v(a)c(a, a_1+t-a-h)-v(a)c(a, a_1+t-a)\right|da, \end{align*}$

since

$\begin{split} \int_0^{a_1}q(a)&\gamma_1(a)i(a, t+h)da-\int_0^{a_1}q(a)\gamma_1(a)i(a, t)da\\ = {}&\int_0^hq(a)\gamma_1(a)i(a, t+h)da+\int_h^{a_1}q(a)\gamma_1(a)i(a, t+h)da\\ &-\int_0^{a_1}q(a)\gamma_1(a)i(a, t)da\\ = {}&\int_0^hq(a)\gamma_1(a)i(0, t+h-a)\pi_1(a)da+\int_0^hq(a)\gamma_1(a)\int_0^a\Lambda_i(\tau)\frac{\pi_1(a)}{\pi_1(\tau)}d\tau da\\ &+\int_h^{a_1}q(a)\gamma_1(a)i(a, t+h)da-\int_0^{a_1}q(a)\gamma_1(a)i(a, t)da\\ \leq{}& \bar{q}\bar{\gamma}_1\sigma Kh+\bar{q}\bar{\gamma}_1\tilde{\Lambda}_i h+\int_0^{a_1-h}q(\tau+h)\gamma_1(\tau+h)i(\tau+h, t+h)d\tau\\ &-\int_0^{a_1}q(a)\gamma_1(a)i(a, t)da\\ = {}&\bar{q}\bar{\gamma}_1(\sigma K+\tilde{\Lambda}_i)h+\int_0^{a_1-h}q(\tau+h)\gamma_1(\tau+h)i(\tau, t)\frac{\pi_1(\tau+h)}{\pi_1(\tau)}d\tau\\ & +\int_0^{a_1-h}q(\tau+h)\gamma_1(\tau+h)\int_{\tau}^{\tau+h}\Lambda_i(s)\frac{\pi_1(\tau+h)}{\pi_1(s)}dsd\tau\\ & -\int_0^{a_1}q(a)\gamma_1(a)i(a, t)da. \end{split}$

(15)

Observing the fact that

$\begin{align*} i(\tau+h, t+h) = i(\tau, t)\frac{\pi_1(\tau+h)}{\pi_1(\tau)}+\int_{\tau}^{\tau+h}\Lambda_i(s)\frac{\pi_1(\tau+h)}{\pi_1(s)}ds \end{align*}$

and

$\begin{split} \int_0^{a_1-h}q(\tau+h)\gamma_1(\tau+h)&\int_{\tau}^{\tau+h}\Lambda_i(s)\frac{\pi_1(\tau+h)}{\pi_1(s)}dsd\tau\\ & \leq \bar{q}\bar{\gamma}_1\int_0^{a_1-h}\int_{\tau}^{\tau+h}\Lambda_i(s)dsd\tau\\ & = \bar{q}\bar{\gamma}_1\left(\int_0^h\int_0^s\Lambda_i(s)d\tau ds+\int_0^h\int_h^{a_1}\int_{s-h}^{s}\Lambda_i(s)d\tau ds\right)\\ & = \bar{q}\bar{\gamma}_1\left(\int_0^hs\Lambda_i(s)ds+\int_h^{a_1}h\Lambda_i(s)ds\right) \end{split}$

$\begin{split}& \leq\bar{q}\bar{\gamma}_1h\int_0^{a_1}\Lambda_i(s)ds\\ & = \bar{q}\bar{\gamma}_1\tilde{\Lambda}_ih, \end{split}$

(16)

and

$1-\bar{\mu}_ih\leq e^{-\bar{\mu}_ih}\leq e^{-\int_a^{a+h}\mu_i(s)ds}\leq 1$

for $h$ sufficiently small, substituting (16) into (15), we have

$\begin{split} \int_0^{a_1}&\left|q(a)\gamma_1(a)i(a, t+h)-q(a)\gamma_1(a)i(a, t)\right|da\\ \leq{}&\bar{q}\bar{\gamma}_1(\sigma K+2\tilde{\Lambda}_i)h+\int_0^{a_1-h}\left|q(a+h)\gamma_1(a+h)\frac{\pi_1(a+h)}{\pi_1(a)}-q(a)\gamma_1(a)\right|i(a, t)da\\ \leq{}&\bar{q}\bar{\gamma}_1(\sigma K+2\tilde{\Lambda}_i)h+\int_0^{a_1-h}q(a+h)\gamma_1(a+h)\left|e^{-\int_a^{a+h}\mu_i(s)ds}-1\right|i(a, t)da\\ ~~ &+\int_0^{a_1-h}\left|q(a+h)\gamma_1(a+h)-q(a)\gamma_1(a)\right|i(a, t)da\\ \leq{}&\bar{q}\bar{\gamma}_1(\sigma K+2\tilde{\Lambda}_i+\bar{\mu}_iK)h+L_{\gamma_1}^{q}Kh. \end{split}$

(17)

Here $L_{\gamma_1}^{q}$ is the Lipschitz coefficient of $q(a)\gamma_1(a)$, which is dependent on $L_q, L_{\gamma_1}, \bar{q}$ and $\bar{\gamma}_1$, since $q(a)$ and $\gamma_1(a)$ are both Lipschitz continuous on $[0, a_1]$. By similar calculations, we can obtain

$\label{o} \int_{a_1}^{\infty}\left|v(a)c(a, t+h)-v(a)c(a, t)\right|da \leq\bar{v}(K+2\tilde{\Lambda}_c+\bar{\mu}_cK)h+L_vKh.$

(18)

Substituting (13)-(18) into (12) gives

$\int_{a_1}^{\infty}|\tilde{c}(a+h, t)-\tilde{c}(a, t)|da\leq Lh,$

where $L = 2(\bar{q}\bar{\gamma}_1+b\omega \bar{v})K+\bar{q}\bar{\gamma}_1(\sigma K+2\tilde{\Lambda}_i+\bar{\mu}_iK) +\bar{v}(K+2\tilde{\Lambda}_c+\bar{\mu}_cK)+(L_{\gamma_1}^{q}+L_v)K$ is independent of $u_0$. Hence condition (3) of Lemma 4.2 is satisfied for $\tilde{c}(a, t)$. This implies that $\tilde{c}(a, t)$ remains in a pre-compact subset $C_K^c$ of $L_{\geq 0}^1$.

Consequently, $\Phi_t(C)\subseteq [0, K]^2\times C_K^i\times C_K^c\times [0, K]$, which has compact closure in $\mathcal{Y}$. It follows that $\Phi_t$ has compact closure. By Lemma 4.1, Lemma 4.2 and Lemma 4.3, we have the following theorem.

Theorem 4.4. The solution semiflow $\Phi_t$ is asymptotically smooth, and there exists a compact global attractor $\mathcal{A}$ for $\Phi_t$.

This result is implied by Lemma 3.1, Theorem 4.4 and Theorem 1.1.3 in [28].

---

5. Equilibria and global stability

In this section, we will investigate the existence of the endemic equilibrium and its global stability, in the special case of $v(a) = 0$, using a Lyapunov functional as in [26]. Consider the reduced model

$\begin{split} S'(t)& = \Lambda_S-\int_0^{a_1}\beta(a)i(a, t)S(t)da-\int_{a_1}^\infty \alpha\beta(a)c(a, t)S(t)da-(\mu_S+p) S, \nonumber \end{split}$

$\label{reducedmodel} \begin{split} E'(t)& = \Lambda_E+\int_0^{a_1}\beta(a)i(a, t)S(t)da+\int_{a_1}^\infty \alpha\beta(a)c(a, t)S(t)da-(\mu_E+\sigma)E, \\ \frac{\partial i(a, t)}{\partial a}&+ \frac{\partial i(a, t)}{\partial t} = \Lambda_i(a)-(\mu_i(a)+\gamma_1(a))i(a, t), 0 <a\leq a_1 \\ \frac{\partial c(a, t)}{\partial a}&+ \frac{\partial c(a, t)}{\partial t} = \Lambda_c(a)-(\mu_c(a)+\gamma_2(a)+\theta(a))c(a, t), a_1 <a <\infty \end{split}$

(19)

with boundary conditions

$\begin{align*} i(0, t)& = \sigma E(t) \\ c(a_1, t)& = \int_0^{a_1}q(a)\gamma_1(a)i(a, t)da \end{align*}$

and initial conditions

$S\left( 0 \right){\rm{ }} = {S_0}, \;\;\;E\left( 0 \right){\rm{ }} = {E_0}, \;\;\;i(a, 0){\rm{ }} = {i_0}(a), \;\;\;c(a, 0){\rm{ }} = {c_0}(a).$

For convenience, denote

$\begin{align*} W_1& = \int_0^{a_1}\beta(a)\pi_1(a)da, \\ W_2& = \int_0^{a_1}q(a)\gamma_1(a)\pi_1(a)da, \\ W_3& = \int_{a_1}^\infty \beta(a)\pi_2(a)da, \\ W_4& = \int_{a_1}^\infty v(a)\pi_2(a)da, \\ W_5& = \int_0^{a_1}q(a)\gamma_1(a)\int_0^{a}\Lambda_i(\tau)\frac{\pi_1(a)}{\pi_1(\tau)}d\tau da, \\ W_6& = \int_{a_1}^\infty v(a)\int_{a_1}^{a}\Lambda_c(\tau)\frac{\pi_2(a)}{\pi_2(\tau)}d\tau da, \\ W_7& = \int_0^{a_1}\beta(a)\int_0^{a}\Lambda_i(\tau)\frac{\pi_1(a)}{\pi_1(\tau)}d\tau da, \\ W_8& = \int_{a_1}^{\infty}\beta(a)\int_{a_1}^{a}\Lambda_c(\tau)\frac{\pi_2(a)}{\pi_2(\tau)}d\tau da, \\ W_9(a)& = \int_0^{a}\frac{\Lambda_i(\tau)}{\pi_1(\tau)}d\tau, \\ W_{10}(a)& = \int_{a_1}^a\frac{\Lambda_c(\tau)}{\pi_2(\tau)}d\tau \end{align*}$

and

$V_1 = \frac{\Lambda_S+\Lambda_E}{\mu_S+p}, \;\;\;\;V_2 = \frac{\mu_E+\sigma}{\mu_S+p}, \;\;\;\;V_3 = \frac{b\omega(c^*(a_1)W_4+W_6)}{\mu_S+p}.$

Let $P^* = (S^*, E^*, i^*(a), c^*(a), R^*)$ denote the endemic equilibrium, with $S^*>0, E^*>0, R^*>0,$ where $i^*(a)$ and $c^*(a)$ satisfy

$\begin{align*} i^*(a)& = \sigma E^*\pi_1(a)+W_9\pi_1(a), \\ c^*(a)& = c^*(a_1)\pi_2(a)+W_{10}\pi_2(a). \end{align*}$

Substituting $i^*(a)$ and $c^*(a)$ into the boundary conditions yields

$\begin{align*} c^*(a_1) = {}&\sigma E^*W_2+\int_0^{a_1}q(a)\gamma_1(a)W_9(a)\pi_1(a)da\\ &+b\omega\int_{a_1}^\infty v(a)(c^*(a_1)\pi_2(a)+W_{10}(a)\pi_2(a))da. \end{align*}$

It follows that

$c^*(a_1) = \dfrac{\sigma E^*W_2+W_5+b\omega W_6}{1-b\omega W_4}$

Solving the first two equations of system (1) in terms of $S$ and $E$ gives

$\begin{align*} S^* = \frac{1}{\mu_S+p}\left[(\Lambda_S+\Lambda_E)-(\mu_E+\sigma)E^*-b\omega(c^*(a_1)W_4+W_6)\right] = V_1-V_2E^*-V_3. \end{align*}$

Noting that $b\omega W_4$ is the number of infants infected via vertical transmission, we have $1- b\omega W_4 >0$ and hence $c^*(a_1)>0$. Thus $i^*(a), c^*(a)$ and $S^*$ can be expressed in terms of $E^*$, and are all positive for $E^*\in \left(0, \frac{\Lambda_S+\Lambda_E-b\omega(c^*(a_1)W_4+W_6)}{\mu_E+\sigma}\right)$ when they exist.

Since

$\begin{align*} ~~\, \int_0^{a_1}\beta(a)i^*(a)da+&\alpha\int_{a_1}^\infty\beta(a)c^*(a)da\\ = {}&\int_0^{a_1}\beta(a)\sigma E^*\pi_1(a)da+\int_0^{a_1}\beta(a)W_9(a)\pi_1(a)da\\ ~~ &+\alpha\int_{a_1}^\infty\beta(a)c^*(a_1)\pi_2(a)da+\alpha\int_{a_1}^\infty\beta(a)W_{10}(a)\pi_2(a)da\\ = {}&\sigma E^*W_1+W_7+\alpha c^*(a_1)W_3+\alpha W_8, \end{align*}$

we have

$(\sigma E^*W_1+W_7+\alpha c^*(a_1)W_3+\alpha W_8)(V_1-V_2E^*-V_3) = (\mu_E+\sigma)E^*-\Lambda_E.$

Let

$f(E^*) = A_0(E^*)^2+A_1E^*+A_2.$

Then $f(E^*) = 0$ is a quadratic equation in $E^*$ with coefficients

$\begin{align*} A_0& = -V_2(\sigma W_1+\frac{\alpha\sigma W_2W_3}{1-b\omega W_4})\\ A_1& = (V_1-V_3)\left(\sigma W_1+\frac{\alpha\sigma W_2W_3}{1-b\omega W_4}\right)\\ &-V_2\left(\frac{\alpha W_3(W_5+W_6)}{1-b\omega W_4}+W_7+\alpha W_8\right)-(\mu_E+\sigma)\\ A_2& = \Lambda_E+(V_1-V_3)(W_7+\alpha W_8). \end{align*}$

Noting that $A_0 <0$ and $A_2>0$ (by assumptions H4 and H7), there exists a unique positive equilibrium $P^*$ when age-dependant immigration terms have positive measure.

In the following, based on the techniques used in [26,18], a Lyapunov functional is constructed to prove the global stability of the unique positive equilibrium $P^*$ when $v(a) = 0$.

Theorem 5.1. When $v(a) = 0$, $\Lambda_S>0$, and $\Lambda_i+\Lambda_c$ have positive measure, system (1) has a unique, globally asympototically stable endemic equilibrium $P^*$.

Proof. For $t\in \mathbb{R}$, let $X(t) = (S(t), E(t), i(\cdot, t), c(\cdot, t))$ be a solution trajectory to the reduced system $(19)$. For $x>0, x^*>0$, define

$G(x, x^*) = x-x^*-x^*\ln\frac{x}{x^*} ~~~~\hbox{and}~~~~ g(x) = x-1-\ln x.$

It is obvious that $G(x, x^*) = x^*g(\frac{x}{x^*})$, and they share the properties $g(x)\geq 0, -g(x) = 1-x+\ln x\leq 0$ for any $x>0$. Let

$\begin{align*} L_1(t)& = G(S(t), S^*), \\ L_2(t)& = G(E(t), E^*), \\ L_3(t)& = \int_0^{a_1}\varphi(a)G(i(a, t), i^*(a))da, \\ L_4(t)& = \int_{a_1}^\infty\psi(a)G(c(a, t), c^*(a))da, \end{align*}$

where

$\begin{align*} \varphi(a)& = \int_a^{a_1}\left[\beta(s)S^*+\psi(a_1)q(s)\gamma_1(s)\right]\frac{\pi_1(s)}{\pi_1(a)}ds, \\ ~~~ \psi(a)& = \alpha\int_a^\infty\beta(s)S^*\frac{\pi_2(s)}{\pi_2(a)}ds. \end{align*}$

Noting that

$\begin{align*} ~~~ \varphi(0)& = \int_0^{a_1}[\beta(s)S^*+\psi(a_1)q(s)\gamma_1(s)]\pi_1(s)ds, \\ ~~~ \psi(a_1)& = \alpha\int_{a_1}^\infty \beta(s)S^*\pi_2(s)ds = \alpha S^*W_3, \\ ~~~ \varphi'(a)& = -(\beta(a)S^*+\psi(a_1)q(a)\gamma_1(a))+(\mu_i(a)+\gamma_1(a))\varphi(a), \\ ~~~ \psi'(a)& = -\alpha\beta(a)S^*+(\mu_c(a)+\gamma_2(a)+\theta(a))\psi(a), \end{align*}$

and

$\begin{align*} \dfrac{di^*(a)}{da}& = \Lambda_i(a)-(\mu_i(a)+\gamma_1(a))i^*(a), \\ \dfrac{i(a, t)}{i^*(a)}& = \dfrac{i(0, t-a)+W_9(a)}{i^*(0)+W_9(a)}, \end{align*}$

we have

$\begin{split} &\int_{t-a_1}^{t}\varphi(t-\tau)i^*(t-\tau)\frac{d}{dt}g\left(\frac{i(t-\tau, t)}{i^*(t-\tau)}\right)d\tau \\ & = \int_{t-a_1}^{t}\varphi(t-\tau)i^*(t-\tau)g'\left(\frac{i(t-\tau, t)}{i^*(t-\tau)}\right)\frac{d}{dt}\left(\frac{i(t-\tau, t)}{i^*(t-\tau)}\right)d\tau\\ & = \int_{t-a_1}^{t}\varphi(t-\tau)i^*(t-\tau)g'\left(\frac{i(t-\tau, t)}{i^*(t-\tau)}\right)\\ & \times\frac{W'_9(t-\tau)(i^*(0)+W_9(t-\tau))-(i(0, \tau)+W_9(t-\tau))W'_9(t-\tau)}{(i^*(0)+W_9(t-\tau))^2}d\tau\\ & = \int_0^{a_1}\varphi(a)i^*(a)g'\left(\frac{i(a, t)}{i^*(a)}\right)\frac{W'_9(a)(i^*(0)+W_9(a))-(i(0, t-a)+W_9(a))W'_9(a)}{(i^*(0)+W_9(a))^2}da \end{split}$

$\label{1} \begin{split}& = \int_0^{a_1}\varphi(a)i^*(a)g'\left(\frac{i(a, t)}{i^*(a)}\right)\frac{W'_9(a)}{i^*(0)+W_9(a)}\left(1-\frac{i(a, t)}{i^*(a)}\right)da\\ & = \int_0^{a_1}\varphi(a)\Lambda_i(a)\left(1-\frac{i^*(a)}{i(a, t)}\right)\left(1-\frac{i(a, t)}{i^*(a)}\right)da\\ & = \int_0^{a_1}\varphi(a)\Lambda_i(a)\left(2-\frac{i^*(a)}{i(a, t)}-\frac{i(a, t)}{i^*(a)}\right)da, \end{split}$

(20)

and

$\begin{split} ~~~ &\int_{t-a_1}^{t}\varphi(t-\tau)\frac{di^*(t-\tau)}{dt}g\left(\frac{i(t-\tau, t)}{i^*(t-\tau)}\right)d\tau\\ & = \int_{t-a_1}^{t}\varphi(t-\tau)\left[\Lambda_i(t-\tau)-(\mu_i(t-\tau)+\gamma_1(t-\tau))i^*(t-\tau)\right]g\left(\frac{i(t-\tau, t)}{i^*(t-\tau)}\right)d\tau\\ & = \int_0^{a_1}\varphi(a)\Lambda_i(a)g\left(\frac{i(a, t)}{i^*(a)}\right)da -\int_0^{a_1}\varphi(a)(\mu_i(a)+\gamma_1(a))i^*(a)g\left(\frac{i(a, t)}{i^*(a)}\right)da. \end{split}$

(21)

Combining (20) with (21), we get

$\label{3} \begin{split} ~~\int_{t-a_1}^{t}&\varphi(t-\tau)\left[i^*(t-\tau)\frac{d}{dt}g\left(\frac{i(t-\tau, t)}{i^*(t-\tau)}\right) +\frac{d}{dt}(i^*(t-\tau))g\left(\frac{i(t-\tau, t)}{i^*(t-\tau)}\right)\right]d\tau\\ = {}&\int_0^{a_1}\varphi(a)\Lambda_i(a)\left[2-\frac{i^*(a)}{i(a, t)}-\frac{i(a, t)}{i^*(a)}+g\left(\frac{i(a, t)}{i^*(a)}\right)\right]da\\ &-\int_0^{a_1}\varphi(a)(\mu_i(a)+\gamma_1(a))i^*(a)g\left(\frac{i(a, t)}{i^*(a)}\right)da\\ = {}&\int_0^{a_1}\varphi(a)\Lambda_i(a)\left(1-\frac{i^*(a)}{i(a, t)}+\ln\frac{i^*(a)}{i(a, t)}\right)da\\ &-\int_0^{a_1}\varphi(a)(\mu_i(a)+\gamma_1(a))i^*(a)g\left(\frac{i(a, t)}{i^*(a)}\right)da\\ = {}&-\int_0^{a_1}\varphi(a)\Lambda_i(a)g\left(\frac{i^*(a)}{i(a, t)}\right)da -\int_0^{a_1}\varphi(a)(\mu_i(a)+\gamma_1(a))i^*(a)g\left(\frac{i(a, t)}{i^*(a)}\right)da\\ \equiv{}& I_i\leq 0, \end{split}$

(22)

since $g(x)\geq 0$ for any $x>0$. Similarly,

$\label{4} \begin{split} ~~ \int_{-\infty}^{t-a_1}&\psi(t-\tau)\left [c^*(t-\tau)\frac{d}{dt}g\left(\frac{c(t-\tau, t)}{c^*(t-\tau)}\right)+\frac{d}{dt}(c^*(t-\tau))g\left(\frac{c(t-\tau, t)}{c^*(t-\tau)}\right)\right] d\tau\\ = {}&\int_{a_1}^{\infty}\psi(a)\Lambda_c(a)\left(1-\frac{c^*(a)}{c(a, t)}+\ln\frac{c^*(a)}{c(a, t)}\right)da\\ ~~&-\int_{a_1}^{\infty}\psi(a)(\mu_c(a)+\gamma_2(a)+\theta(a))c^*(a)g\left(\frac{c(a, t)}{c^*(a)}\right)da\\ = {}&-\int_{a_1}^{\infty}\psi(a)\Lambda_c(a)g\left(\frac{c^*(a)}{c(a, t)}\right)da\\ &-\int_{a_1}^{\infty}\psi(a)(\mu_c(a)+\gamma_2(a)+\theta(a))c^*(a)g\left(\frac{c(a, t)}{c^*(a)}\right)da\\ \equiv{}& I_c\leq 0. \end{split}$

(23)

The terms (22) and (23) are a result of incorporating immigration; compare with the model in [26].

The time derivatives of functions $L_1, L_2, L_3$ and $L_4$ along (1) are

$\begin{align*} L_1'(t) = {}&-\frac{\mu_S+p}{S}(S-S^*)^2+S^*\left(\int_0^{a_1}\beta(a)i^*(a)da+\alpha\int_{a_1}^{\infty}\beta(a)c^*(a)da\right)\\ ~~~~~~~~~~ &+S^*\left(\int_0^{a_1}\beta(a)i(a, t)da+\alpha\int_{a_1}^{\infty}\beta(a)c(a, t)da\right)\\ ~~~~~~~~~~ &-S\left(\int_0^{a_1}\beta(a)i(a, t)da+\alpha\int_{a_1}^{\infty}\beta(a)c(a, t)da\right)\\ ~~~~~~~~~~ &-\frac{(S^*)^2}{S}\left(\int_0^{a_1}\beta(a)i^*(a)da+\alpha\int_{a_1}^{\infty}\beta(a)c^*(a)da\right), \\ L_2'(t) = {}&S\left(\int_0^{a_1}\beta(a)i(a, t)da+\alpha\int_{a_1}^{\infty}\beta(a)c(a, t)da\right)\\ ~~~~~~~~~~ &-\frac{E^*}{E}S\left(\int_0^{a_1}\beta(a)i(a, t)da+\alpha\int_{a_1}^{\infty}\beta(a)c(a, t)da\right)\\ ~~~~~~~~~~ &-(\sigma+\mu_E)(E-E^*)+\Lambda_E\left(1-\frac{E^*}{E}\right), \\ L_3'(t) = {}&\frac{d}{dt}\int_0^{a_1}\varphi(a)i^*(a)g\left(\frac{i(a, t)}{i^*(a)}\right)da\\ ~~~~~~~ = {}&\frac{d}{dt}\int_{t-a_1}^{t}\varphi(t-\tau)i^*(t-\tau)g\left(\frac{i(t-\tau, t)}{i^*(t-\tau)}\right)d\tau\\ ~~~~~~~ = {}&\varphi(0)G(i(0, t), i^*(0))+\int_0^{a_1}\varphi'(a)G(i(a, t), i^*(a))da+I_i, \end{align*}$

and

$\begin{align*} L_4'(t)& = \dfrac{d}{dt}\int_{a_1}^{\infty}\psi(a)c^*(a)g\left(\frac{c(a, t)}{c^*(a)}\right)da\\ ~~~~~~~ & = \dfrac{d}{dt}\int_{-\infty}^{t-a_1}\psi(t-\tau)c^*(t-\tau)g\left(\frac{c(t-\tau, t)}{c^*(t-\tau)}\right)d\tau\\ ~~~~~~~ & = \psi(a_1)G(c(0, t), c^*(0))+\int_{a_1}^{\infty}\psi'(a)G(c(a, t), c^*(a))da+I_c. \end{align*}$

We construct the following Lyapunov functional:

$L(t) = L_1(t)+L_2(t)+L_3(t)+L_4(t).$

We need to show that the time derivative along the solutions of (1) satisfies

$\frac{d}{dt}L(t) = L'_1(t)+L'_2(t)+L'_3(t)+L'_4(t)\leq 0.$

Note that

$\begin{align*} i(0, t)-i^*(0)& = \sigma(E-E^*), \\ c(a_1, t)-c^*(a_1)& = \int_0^{a_1}q(a)\gamma_1(a)[i(a, t)-i^*(a)]da, \\ (\sigma+\mu_E)E^*-\Lambda_E& = \int_0^{a_1}\beta(a)S^*i^*(a)da+\alpha\int_{a_1}^\infty\beta(a)S^*c^*(a)da. \end{align*}$

We have

$\begin{align*} (\sigma+\mu_E)E^*-\Lambda_E = {}&\int_0^{a_1}\beta(a)S^*i^*(0)\pi_1(a)da+\int_0^{a_1}\beta(a)S^*W_9(a)\pi_1(a)da\\ ~~&+\alpha\int_{a_1}^\infty\beta(a)S^*\pi_2(a)(c^*(a_1)+W_{10}(a))da\\ = {}&\int_0^{a_1}\beta(a)S^*i^*(0)\pi_1(a)da+\int_0^{a_1}\beta(a)S^*W_9(a)\pi_1(a)da\\ ~~&+\alpha\int_{a_1}^\infty\beta(s)S^*\pi_2(s)\int_0^{a_1}q(a)\gamma_1(a)(i^*(0)+W_9(a))\pi_1(a)dads\\ ~~&+\alpha\int_{a_1}^\infty\beta(s)S^*\pi_2(s)W_{10}(s)ds. \end{align*}$

Then

$\varphi(0)i^*(0)\ln\frac{i(0, t)}{i^*(0)} = ((\sigma+\mu_E)E^*-\Lambda_E)\ln\frac{i(0, t)}{i^*(0)}-W\ln\frac{i(0, t)}{i^*(0)}$

and

$\varphi(0)(i(0, t)-i^*(0)) = (\sigma+\mu_E)(E-E^*)-\Lambda_E\left(\frac{E}{E^*}-1\right)-W\left(\frac{E}{E^*}-1\right),$

where

$\begin{align*} W = {}&\int_0^{a_1}\beta(a)S^*W_9(a)\pi_1(a)da\\ &+\alpha\int_{a_1}^\infty\beta(s)S^*\pi_2(s)\int_0^{a_1}q(a)\gamma_1(a)W_9(a)\pi_1(a)dads\\ &+\alpha\int_{a_1}^\infty\beta(s)S^*\pi_2(s)W_{10}(s)ds. \end{align*}$

We also have

$\begin{align*} ~~\int_0^{a_1}q(a)&\gamma_1(a)i^*(a)\left[1-\frac{i(a, t)c^*(a_1)}{i^*(a)c(a_1, t)}\right]da \\ & = \int_0^{a_1}q(a)\gamma_1(a)i^*(a)da-\frac{c^*(a_1)}{c(a_1, t)}\int_0^{a_1}q(a)\gamma_1(a)i(a, t)da \\ & = c^*(a_1)-\frac{c^*(a_1)}{c(a_1, t)}c(a_1, t)\\ & = 0 \end{align*}$

and

$\begin{align*} ~~\frac{E^*}{E}S\bigg(\int_0^{a_1}\beta(a)&i(a, t)da+\alpha\int_{a_1}^\infty\beta(a)c(a, t)da\bigg)\\ & = \frac{i^*(0)}{i(0, t)}S\left(\int_0^{a_1}\beta(a)i(a, t)da+\alpha\int_{a_1}^\infty\beta(a)c(a, t)da\right)\\ & = \int_0^{a_1}\beta(a)S^*i^*(a)\frac{Si^*(0)i(a, t)}{S^*i(0, t)i^*(a)}da+\\ &\alpha\int_{a_1}^\infty\beta(a)S^*c^*(a)\frac{Si^*(0)c(a, t)}{S^*i(0, t)c^*(a)}da. \end{align*}$

Combining all these time derivatives of $L_i(t)$ yields

$\begin{align*} &L'(t) = -\frac{\mu_S+p}{S}(S-S^*)^2\\ ~~~~~~~~~ &+\int_0^{a_1}\beta(a)S^*i^*(a)\left[1-\frac{S^*}{S}+\ln\frac{S^*}{S}+1-\frac{i(a, t)Si^*(0)}{i^*(a)S^*i(0, t)}+\ln\frac{Si^*(0)i(a, t)}{S^*i^*(a)i(0, t)}\right]da \\ ~~~~~~~~~ &+\alpha\int_{a_1}^\infty\beta(a)S^*c^*(a)\left[1-\frac{S^*}{S}+\ln\frac{S^*}{S}+1-\frac{Sc(a, t)i^*(0)}{S^*c^*(a)i(0, t)}+\ln\frac{Sc(a, t)i^*(0)}{S^*c^*(a)i(0, t)}\right]da \\ ~~~~~~~~~ &+\alpha S^* W_3\int_0^{a_1}q(a)\gamma_1(a)i^*(a)\left[1-\frac{i(a, t)c^*(a_1)}{i^*(a)c(a_1, t)}+\ln\frac{i(a, t)c^*(a_1)}{i^*(a)c(a_1, t)}\right]da\\ ~~~~~~~~~ &+\Lambda_E\left(2-\frac{E^*}{E}-\frac{E}{E^*}\right)+W\left(1-\frac{E}{E^*}+\ln\frac{i(0, t)}{i^*(0)}\right)+I_i+I_c. \end{align*}$

Noting that $\frac{i(0, t)}{i^*(0)} = \frac{E}{E^*}$, we see that $L'(t)\leq 0$ and $L'(t) = 0$ holds, implying that $S(t) = S^*$ and

$\frac{i(a, t)}{i^*(a)} = \frac{i(0, t)}{i^*(0)} = \frac{c(a, t)}{c^*(a)} = \frac{c(a_1, t)}{c^*(a_1)}.$

It can be verified that the largest invariant set where $L'(t) = 0$ is the singleton $\{P^*\}$. Therefore, by Lyapunov-LaSalle asymptotic stability theorem, the unique endemic equilibrium $P^*$ is globally asymptotically stable when it exists.

---

6. Discussion

Many epidemiological models assume that all infected individuals are equally infectious during their infectivity period, which is reasonable for some diseases, such as influenza. However, the infection of HBV is a dynamic process characterized by replicative and non-replicative phases based on virus-host interaction so that the infectivity of HBV individuals varies at different age of infection. It follows that age structure of the host population is an important factor for the dynamics of HBV transmission. In general, there are two different age structures in disease models: biological age and infection age. In this paper, according to the characteristics of HBV, we formulated a PDE model incorporating infection-age structure, as well as immigration age, into all compartments, in order to describe the possible effects of variable infectivity and immigration on the transmission dynamics.

The model studied in this paper is a refined version of the one investigated in [26], with age-dependent immigration into both acute and chronic infection stages. Immigration models do not have infection-free equilibria and hence have no reproduction number [17,18,11,3]. When immigration is introduced into the infectious compartment, the model has a unique globally asymptotically stable equilibrium. Age structure is only applicable to two infectious classes in this paper: the acute and chronic stages. In a more general setting, it can be added to any compartment [29].

In the case that the perinatal-related parameter $v(a)$ is absent and age-related immigration is measurable, the unique endemic equilibrium is globally stable, and the disease will always persist in the endemic level. Epidemiologically, $v(a)>0$ indicates vertical transmission of hepatitis B from mother to baby, a complex process where the detailed mechanism remains unclear.

Relevant simulations in [12] suggest that immigration inflow is vital to the amplitude of the endemic equilibrium: small variations of the rate of immigration to the infectious and/or latent compartments can cause abrupt changes of quantity of endemicity. The feasible epidemiological solution is to screen immigrants with high risk of disease as a first step to control infection from the source. Finally, the scale of the problem and the availability of a safe and effective vaccine means that many more people should be protected from HBV than currently are. Outmoded and ignorant attitudes towards vaccines are costing a great deal of lives that could be saved with the application of evidence-based research.

---

Acknowledgments

For citation purposes, please note that the question mark in "Smith?" is part of his name.

---