1. Introduction

Determining the threshold dynamics of infection-free and infection equilibrium in viral infection model has made great progress in the last decades [1,2,4,8,9,10,16,17,24,31,36,32,33,34,35,40,44]. A key insight in this progress is that if threshold value (named, the basic reproduction number) is less than one then the infection-free equilibrium is globally asymptotically stable otherwise the endemic equilibrium attracts all solutions (is globally asymptotically stable) whenever it exists. One method adopted here is due to the classical Volterra type Lyapunov function, which was discovered by Volterra [30]. These confirmed global stability properties of steady states for within host virus model establish our understanding the virus dynamical behaviors, that is, whether the viruses die out or not.

Even large discrete and continuous delay differential equations of viral infection models have been successfully treated by Volterra type Lyapunov function, (ⅰ) nonlinear incidence rate functions [8,24,31]; (ⅱ) discrete delays [10,16,35] and finite distributed delays [17,31,36], and infinite distributed delays [8,24,34]; (ⅲ) immune responses [24,35,44,39,41]; and (ⅳ) additional infection processes [31]. It is still a hot topic in in-host model to determine how these factors affect the virus dynamical behaviors. We also refer the reader to see these citations for more references.

Recently, age-structured viral infection model has attracted much attentions of researchers. HIV latency remains a major obstacle to viral elimination. Although HIV-1 replication can be controlled by antiretroviral therapy in suppress the plasma viral load to below the detection limit, the time spent in this progress may last half life of months or years [21]. Virus persisting in reservoirs, such as latency infected CD4+ T cells, may the reason that long-term low viral load persistence in patients on antiretroviral therapy and keeping the virus from being eliminated. These latency infected CD4+ T cells are not affected by immune responses but can produce virus once activated by relevant antigens.

Some recent studies reveals the decay dynamics of the latent reservoir. For example, a model has been developed by Muller et al. [11] to describe the heterogeneity of latent cell activation. An ordinary differential equations (ODEs) model has been studied by Kim and Perelson [7] to include decreasing activation of latently infected cells. Activation of these latently infected cells needs specificity antigen. A recent study by Strain et al. [22] reveals that the dynamics of latently infected CD4+ T cells are often heterogeneous. They argued that cells specific to frequently encountered antigens are activated soon while cells specific to rare antigens need more time to be activated. Thus, the activation rate depends on the time spent since the cell is latently infected (that is the time elapsed since the establishment of latency), which we refer as latency age for short. A recent paper by Alshorman et al. [1] introduced a latency age model to mathematically analyze the dynamics of the latent reservoir under combination therapy. They give an affirmed answers that the long-term activation rate of latently infected cells plays an important role in determining the dynamics.

Taking into account the picture that the mortality rate and viral production rate of infected cells may depend on the infection age of cells, Nelson et al. [15], Huang et al. [4] and Wang et al. [37] have studied age-structured model of HIV infection by considering age to be a continuous variable rather than be constant in ODEs models. These assumptions lead to a hybrid system of ODEs and partial differential equations (PDEs) formulation and allow us to have a good understanding on productively infected cells. Together with the infinite-dimensional nature of system, this formulation creates some mathematical difficulties in establishing the existence of a global compact attractor, even in other epidemic models (see some relevant references for our discussion on age-structured models, [6,43,42,25,27,26,14]).

Denote by $T(t), e(a, t), i(a, t), V(t)$ the concentration of uninfected CD4+ T cells at time $t$, the concentration of latently infected T cells with latency age $a$ at time $t$, the concentration of productively infected cells, and the concentration of virions in plasma at $t$, respectively. The parameter $h$ is the production rate of uninfected CD4+ T cells, $d$ is the per capita death rate of uninfected cells, and $\beta$ is the infection rate of the target cell by virus. $c$ is the viral clearance rate.

The following assumptions are a compromise between generality and simplicity.

Assumption 1.1 (ⅰ) There is a small fraction ($f$) of infected cells lead to latency and that the remaining become productively infected cells [1,40].

(ⅱ) When latently infected cells are activated to become productively infected cells, an age-dependent remove rate $-\theta_1(a)$ is used to illustrate the decreasing effect of the pool size of latent reservoir. The integral term $\int_0^\infty \xi(a)e(a,t) da$ describes he total number of productively infected cells gained per unit time from the activation of latently infected cells, where $\xi(a)$ denotes the activation rate of latently infected T cells with latency age $a$. Biologically, we omit the proliferation rate [25] and the death rate of latently infected cells, which is assumed to be included in the removal rate.

(ⅲ) We assume that production rate of viral particles $p(b)$ and the death rate of productively infected cells $\theta_2(b)$ wiht infection age $b$ are two continuous functions of age (the time passed since infection), see, e.g., [15,4,32].

Biologically, (ⅰ) of Assumption 1.1 comes from the evidences that a very small fraction of CD4+ T cell infection leads to HIV latency. They don't produce new virus unless activated by antigens, please see [18,19,20]. (ⅱ) of Assumption 1.1 based on the fact the latently infected cell population is very likely to be heterogeneous [22,23]. Cells specific to frequently encountered antigens may be preferentially activated and quickly removed from the reservoir. It may depends on the time elapsed since latent infection and affect the activation rate, that is the reason why we are interested in the latent infection age. (ⅲ) of Assumption 1.1 it is known that viral proteins and unspliced viral RNA accumulate within the cytoplasm of an infected cell, and thus, they actually ramps up [3,12,29]. Therefore, infection age should be incorporated into the model.

In this paper, we introduce the following HIV infection model with latency and infection age,

$\left\{ \begin{array}{lll} \displaystyle \frac{ dT(t)}{ dt} =h -d T(t)-\beta T(t)V(t), \\ \displaystyle \left(\frac{\partial}{\partial t}+\frac{\partial}{\partial a}\right)e(a,t) =-\theta_1(a)e(a,t),\\ \displaystyle \left(\frac{\partial}{\partial t}+\frac{\partial}{\partial b}\right)i(b,t) =-\theta_2(b)i(b,t),\\ \displaystyle \frac{ dV(t)}{ dt} =\int_0^\infty p(b)i(b,t) db-cV(t), \end{array}\right.$

(1)

with boundary and initial conditions

$\begin{align*} \left\{ \begin{array}{lll} \displaystyle e(0,t) = f\beta T(t)V(t),\\ \displaystyle i(0,t) = (1-f)\beta T(t)V(t)+\int_0^\infty \xi(a)e(a,t) da,\\ \displaystyle T(0) =T_0 \ge 0 ,\ e(a,0) =e_0(a)\in L^1_{+}(0,\infty),\\ \displaystyle i(b,0) =i_0(b)\in L^1_{+}(0,\infty),\ V(0) =V_0 \ge 0, \end{array}\right. \end{align*}$

where $L_+^1(0, \infty)$ is the set of all integrable nonnegative functions on $\mathbb{R}_+: = [0,\infty)$.

Mathematically, for the ease of simplicity, we make the following assumptions.

Assumption 1.2 (ⅰ) $h$, $d$, $\beta$, $c>0$;

(ⅱ) For $1 = 1,2$, $\theta_i(\cdot), p(\cdot), \xi(\cdot)\in L^{\infty}_+ (0,\infty)$ satisfy the conditions:

$\bar{\theta_i}: = ess\sup\limits_{a\in[0,\infty)}\theta_i(a)<\infty,\ \ \bar{p}: = ess\sup\limits_{a\in[0,\infty)}p(a)<\infty,\ \ \bar{\xi}: = ess\sup\limits_{a\in[0,\infty)}\xi(a)<\infty,$

(ⅲ) $p(\cdot)$, $\xi(\cdot)$ are Lipschitz continuous on $\mathbb{R}_{+}$ with Lipschitz constants $M_p$, $M_\xi$ respectively;

(ⅳ) There exists $\mu_0 \in \left( c,d \right]$ such that $\theta_1(a), \theta_2(b) \geq \mu_0$ for all $a, b\geq 0$;

(ⅴ) There exists a maximum age $b^+>0$ for the viral production such that $p(b) > 0$ for $b\in (0,b^+)$ and $p(b) =0$ for $b>b^+$.

Our goal of the present paper is to adopt previous model in [32,33] by incorporating the latency age for infected cells as discussed in [1], and to study the threshold dynamics of infection-free and infection equilibrium in viral infection model subject to latently age and infection age. We will show the existence of a compact attractor of all compact sets of nonnegative initial data and use the Lyapunov function to show that this attractor is the singleton set containing the endemic equilibrium. Roughly speaking, if the basic reproduction number is less than one then the infection-free equilibrium is globally asymptotically stable otherwise the endemic equilibrium attracts all solutions with active infection at some time.

The remaining part of this paper is organized as follows. In Section 2, we present some preliminary results including model formulation (equivalent integrated semigroup formulation and Volterra formulation), properties of solutions and existence of equilibria. Then we show the asymptotic smoothness of $\Phi(t,X_0)$ of orbits in Section 3, where we arrive at a key result on the existence of global compact attractor. In section 4, we prove that system (1) is the uniformly persistent. The Section 5 is devoted to local stability analysis of the infection-free equilibrium and the infection equilibrium. Then we establish their global attractivity in Section 6 by constructing Lyapunov functions.

---

2. Preliminary

For ease of notations, we introduce the following notations:

$\Omega(a) = e^{-\int_0^a\theta_1(\tau)d\tau}\ \ \mbox{and}\ \ \Gamma(b) = e^{-\int_0^b \theta_2(\tau)d\tau}\ \ \mbox{for}\ \ a, b\ge 0.$

Biologically, $\Omega(a)$ is the probability of an infected cell staying in the latent state until age $a$. $\Gamma(b)$ is typically interpreted as the probability that an infected cell can survive to age $b$.

Then, by (ⅱ) and (ⅴ) of Assumption 1.2, we have that for all $a, b\ge 0$,

$\begin{align}\label{E3} &\ 0\leq\Omega(a)\leq \textrm{e}^{-\mu_0a}\ \ \mbox{and}\ \ 0\leq\Gamma(b)\leq \textrm{e}^{-\mu_0b}, \notag\\ &\ \frac{d\Omega(a)}{da} = -\theta_1(a)\Omega(a)\ \ \mbox{and}\ \ \frac{d\Gamma(b)}{db} = -\theta_2(b)\Gamma(b). \end{align}$

(2)

---

2.1. Integrated semigroup formulation

Following the line of [27], we reformulate the model (1) as a semilinear Cauchy problem. Taking into account the boundary conditions, we consider the following state space,

$\mathcal{X} = \mathbb{R}\times\mathbb{R}\times L^1((0,\infty), \mathbb{R})\times\mathbb{R}\times L^1((0,\infty), \mathbb{R})\times\mathbb{R},$

$\mathcal{X}_+ = \mathbb{R}_+\times\mathbb{R}_+\times L_+^1((0,\infty), \mathbb{R})\times\mathbb{R}_+\times L_+^1((0,\infty),\mathbb{R})\times\mathbb{R}_+,$

endowed with the usual product norm, and set

$\mathcal{X}_0 = \mathbb{R} \times \{0\} \times L^1((0,\infty), \mathbb{R})\times \{0\} \times L^1((0,\infty), \mathbb{R})\times\mathbb{R},$

$\mathcal{X}_{0+} = \mathcal{X}_0\cap\mathcal{X}_+.$

We consider the linear operator $\mathcal{A} : Dom(\mathcal{A})\subset \mathcal{X}\rightarrow\mathcal{X}$ defined by

$\begin{equation*} \mathcal{A}\begin{pmatrix} T\\ \begin{pmatrix} 0\\ e\\ \end{pmatrix}\\ \begin{pmatrix} 0\\ i\\ \end{pmatrix}\\ V \end{pmatrix} = \begin{pmatrix} -dT\\ \begin{pmatrix} -e(0)\\ -e'-\theta_1(a)e\\ \end{pmatrix}\\ \begin{pmatrix} -i(0)\\ -i'-\theta_2(a)i\\ \end{pmatrix}\\ -cV \end{pmatrix} \end{equation*}$

with

$\textrm{Dom}(\mathcal{A}) = \mathbb{R}\times \{0\} \times W^{1,1}((0,\infty), \mathbb{R})\times \{0\} \times W^{1,1}((0,\infty), \mathbb{R})\times\mathbb{R},$

where $W^{1,1}$ is a Sobolev space. Note that $\overline{\textrm{Dom}(\mathcal{A})} = \mathcal{X}_0$ is not dense in $\mathcal{X}$.

Define nonlinear operator $\mathcal{F} : \overline{\textrm{Dom}(\mathcal{A})} \rightarrow\mathcal{X}$ by

$\mathcal{F}\begin{pmatrix} T\\ \begin{pmatrix} 0\\ e\\ \end{pmatrix}\\ \begin{pmatrix} 0\\ i\\ \end{pmatrix}\\ V \end{pmatrix} =\begin{pmatrix} h-\beta TV\\ \begin{pmatrix} f\beta TV\\ 0_{L_1}\\ \end{pmatrix}\\ \begin{pmatrix} (1-f)\beta TV+M\\ 0_{L_1}\\ \end{pmatrix}\\ N \end{pmatrix}.$

where

$M(t) =\int_0^{\infty} \xi(a)e(a, t)da,\ \ N(t) =\int_0^{\infty} p(b)i(b, t)db.$

(3)

Then by setting $u(t) = \left( T(t), \begin{pmatrix} 0\\ e(\cdot,t)\\ \end{pmatrix}, \begin{pmatrix} 0\\ i(\cdot,t)\\ \end{pmatrix}, V(t)\right)^T$, we can reformulate system (1) with the boundary and initial conditions as the following abstract Cauchy problem

$\frac{du(t)}{dt} = \mathcal{A}u(t)+\mathcal{F}(u(t))\ \ \mbox{for}\ \ t\geq 0\ \ \mbox{and}\ \ u(0)\in \mathcal{X}_{0+}.$

If any initial value $X_0 = (T_0, e_0(\cdot), i_0(\cdot), V_0)\in \mathcal{Y}$ satisfies the coupling equations

$e(0,0) = f\beta T_0V_0$

and

$i(0,0) = (1-f)\beta T_0V_0+\int_0^{\infty}\xi(a)e_0(a)da,$

then (1) is well-posed under Assumption 1.2 due to Iannelli [6] and Magal [14]. Denote

$\mathcal{Y} = \mathbb{R}_+\times L_+^1(0,\infty)\times L_+^1(0,\infty)\times \mathbb{R}_+$

with the norm

$\|(x, \varphi, \psi, y)\|_{\mathcal{Y}} = |x|+\|\varphi\|_{L^1}+\|\psi\|_{L^1}+|y|\ \ \mbox{for}\ (x, \varphi, \psi, y)\in \mathcal{Y}.$

In fact, for such solutions, it is not difficult to show that $(T(t), e(\cdot, t), i(\cdot, t), V(t))\in \mathcal{Y}$ for each $t\ge 0$. In the sequel, we always assume that the initial values satisfy the coupling equations.

Using the results presented in [14,27], thus we can get a continuous solution semi-flow $\Phi:\mathbb{R}_+\times \mathcal{Y}\to\mathcal {Y}$ defined by

$\Phi \left( t,X_0\right) = \Phi_t(X_0): =\left( T(t), e(\cdot, t), i(\cdot, t), V(t) \right), t\geq 0,\ \ X_0 \in \mathcal{Y}.$

The precise result is the following proposition.

Proposition 1. For system (1), there exists a unique strongly continuous semiflow ${\Phi:\ \mathcal{X}_{0+}\rightarrow \mathcal{X}_{0+}}_{t\geq0}$ such that for each $x_0 \in \mathcal{X}_{0+}$, the operator $x\in C([0, \infty), \mathcal{X}_{0+})$ defined by $x =\Phi(t)x_0$ is a mild solution of (1), that is, it satisfies

$\int_0^{t}x(s)ds\in \textrm{Dom}(\mathcal{A}),\ \mbox{and}\ x(t) =x_0 +\mathcal{A}\int_0^tx(s)ds+ \int_0^t \mathcal{F}(x(s))ds, \forall t\geq0.$

---

2.2. Volterra formulation

According to the Volterra formulation (see Webb [42] and Iannelli [6]), integrating the second and third equations of (1) along the characteristic lines $t-a = \mathrm{const.}$ and $t-b = \mathrm{const.}$ respectively yields

$e(a,t) = \left\{ {\begin{array}{*{20}{l}} {f\beta T(t - a)V(t - a)\Omega (a) = e(0,t - a)\Omega (a),}&{{\rm{if}}\;\;t > a,}\\ {{e_0}(a - t)\frac{{\Omega (a)}}{{\Omega (a - t)}},}&{{\rm{if}}\;\;t \le a;} \end{array}} \right.$

(4)

and

$\begin{array}{*{20}{l}} {i(b,t) = }\\ {\left\{ {\begin{array}{*{20}{l}} {\left[ {(1 - f)\beta T(t - b)V(t - b) + M(t - b)} \right]\Gamma (b) = i(0,t - b)\Gamma (b),\;}&{{\rm{if}}\;t > b,}\\ {{i_0}(b - t)\frac{{\Gamma (b)}}{{\Gamma (b - t)}},}&{{\rm{if}}\;t \le b.} \end{array}} \right.} \end{array}$

(5)

Thus system (1) can be rewritten as the following Volterra-type equations,

$\left\{ \begin{array}{ll} \displaystyle \frac{ dT(t)}{ dt} =h -d T(t)-\beta T(t)V(t), \\ \displaystyle \frac{ dV(t)}{ dt} =\int_0^t p(b)\Gamma(b)((1-f)\beta T(t-b)V(t-b)+M(t-b)) db\\ \displaystyle +\int_t^{\infty}p(b)i_0(b-t)\frac{\Gamma(b)}{\Gamma(b-t)}db-cV(t), \end{array}\right.$

where $M(t-b) = \int_0^{t-b} \xi(a)e(0,t-b-a)da+\int_{t-b}^{\infty}\xi(a)e_0 (a-t-b)\frac{\Omega(a)}{\Omega(a-t-b)}da$.

---

2.3. Boundedness of solutions

Proposition 2. Define

$\begin{array}{l} \Xi : = \;\left\{ {{X_0} = ({T_0},{e_0},{i_0},{V_0}) \in {\cal Y}\;|\;{T_0} + \left\| {{e_0}(a)\left\| {_{{L^1}}} \right.} \right. \le \frac{h}{{{\mu _0}}},} \right.\\ \quad \quad \quad \quad \quad \quad \quad \quad \quad \quad \quad \quad \quad \quad {T_0} + \left\| {{e_0}(a)} \right.\left\| {_{{L^1}}} \right. + \left\| {{i_0}} \right.(b)\left\| {_{{L^1}}} \right. \le \frac{h}{{{\mu _1}}},\;\\ \quad \quad \quad \quad \quad \quad \quad \quad \quad \quad \quad \quad \quad \quad \left. {{V_0} \le \frac{{\bar ph}}{{c{\mu _0}}} + \frac{{h\bar p\bar \xi }}{{c\mu _0^2}},\;{{\left\| {{X_0}} \right\|}_{\cal Y}} \le \frac{h}{{{{\tilde \mu }_0}}}} \right\}, \end{array}$

where $\tilde{\mu}_0 \ : =\ \frac{\mu_0}{ 1+\frac{\bar{\xi}}{\mu_0}+\frac{\bar{p}}{c}+\frac{\bar{p}\bar{\xi}}{c\mu_0}}$, ${\mu}_1 \ : =\ \frac{\mu_0}{ 1+\frac{\bar{\xi}}{\mu_0}}$. Then $\Xi$ is a positively invariant subset for $\Phi$, that is,

$\Phi (t,{X_0}) \in \Xi \;\;for\;all\;\;t \ge 0\;\;and\;\;{X_0} \in \Xi .$

Moreover, $\Phi$ is point dissipative and $\Xi$ attracts all points in $\mathcal{Y}$.

Proof. By (4) and changes of variables, we have

$\begin{align*} \|e(\cdot,t)\|_{L^1} = &\ \int_0^t e(0,t-a)\Omega(a) da+\int_t^\infty e_0(a-t)\frac{\Omega(a)}{\Omega(a-t)} da\\ =&\ \int_0^t e(0,\sigma)\Omega(t-\sigma) d\sigma+ \int_0^\infty e_0 (\tau)\frac{\Omega(t+\tau)}{\Omega(\tau)} d\tau. \end{align*}$

We derivative this equality,

$\frac{{d{{\left\| {e( \cdot ,t)} \right\|}_{{L^1}}}}}{{dt}} = e(0,t)\Omega (0) + \int_0^t e (0,\sigma )\frac{{d\Omega (t - \sigma )}}{{dt}}d\sigma + \int_0^\infty {\frac{{{e_0}(\tau )}}{{\Omega (\tau )}}} \frac{{d\Omega (t + \tau )}}{{dt}}d\tau .$

By (2) and changing of variables, we have

$\begin{align*} \frac{d\|e(\cdot,t)\|_{L^1}}{dt} =&\ e(0,t)\Omega(0)-\int_0^t e(0,\sigma)\theta_1(t-\sigma)\Omega(t-\sigma)d\sigma\\ &\ -\int_0^{\infty} \frac {e_0(\tau)}{\Omega(\tau)}\theta_1(t+\tau)\Omega(t+\tau)d\tau \nonumber\\ = &\ e(0,t)\Omega(0)-\int_0^\infty \theta_1(a)e(a,t) da. \end{align*}$

By the first equation in (1), (ⅳ) of Assumption 1.2, and use $f<1$,

$\begin{array}{*{20}{l}} {\frac{{d(T(t) + {{\left\| {e( \cdot ,t)} \right\|}_{{L^1}}})}}{{dt}} = h - dT(t) - \beta T(t)V(t) + f\beta T(t)V(t) - \int_0^\infty {{\theta _1}} (a)e(a,t)da}\\ {\quad \quad \quad \quad \quad \quad \quad \quad \quad \le h - {\mu _0}(T(t) + {{\left\| {e( \cdot ,t)} \right\|}_{{L^1}}}).} \end{array}$

We integrate this differential inequality and obtain the a priori estimate,

$T(t)+ \|e(\cdot,t)\|_{L^1} \leq \frac{h}{\mu_0} -e^{-\mu_0 t}\left\{ \frac{h}{\mu_0} - ( T_0+ \|e(\cdot,t)\|_{L^1}) \right\}, t\geq 0.$

(6)

This implies that

$T(t)+\|e(\cdot,t)\|_{L^1}\leq \frac{h}{\mu_0}.$

Similarly, we have

$\frac{ d\|i(\cdot, t)\|_{L^1}}{ dt} = i(0, t)\Gamma(0)-\int_0^\infty \theta_2(b)i(b,t) db.$

We add the two equations,

$\begin{array}{l} \frac{{d(T(t) + {{\left\| {e( \cdot ,t)} \right\|}_{{L^1}}}) + {{\left\| {i( \cdot ,t)} \right\|}_{{L^1}}})}}{{dt}} = h - dT + \int_0^\infty \xi (a)e(a,t)da\\ \quad \;\quad \quad \quad \quad \quad \quad \quad \quad \quad \quad \quad - \int_0^\infty {{\theta _1}} (a)e(a,t)da - \int_0^\infty {{\theta _2}} (b)i(b,t)db. \end{array}$

and since (ⅱ) and (ⅳ) of Assumption 1.2, obtain the estimates

$\frac{{d(T(t) + {{\left\| {e( \cdot ,t)} \right\|}_{{L^1}}}) + {{\left\| {i( \cdot ,t)} \right\|}_{{L^1}}})}}{{dt}} \le h + \bar \xi \frac{h}{{{\mu _0}}} - {\mu _0}(T(t) + {\left\| {e( \cdot ,t)} \right\|_{{L^1}}} + {\left. {i( \cdot ,t)} \right\|_{{L^1}}}).$

We integrate this differential inequality and obtain the a priori estimate,

$\begin{align}\label{Ii} &\ \|T(t)+\|e(\cdot, t)\|_{L^1})+ \|i(\cdot, t)\|_{L^1}\nonumber\\ \leq&\ \frac{h+\bar{\xi}\frac{h}{\mu_0}}{\mu_0}-e^{-\mu_0 t}\left\{\frac{h+\bar{\xi}\frac{h}{\mu_0}}{\mu_0}-(T(t)+ \| e(\cdot, t)\|_{L^1}+i(\cdot, t)\|_{L^1})\right\}, \ \ t\geq 0. \end{align}$

(7)

This implies $\|i(\cdot,t)\|_{L^1}\leq \frac{h}{\mu_0}+\frac{\bar{\xi} h}{\mu_0^2}$. Further, since (ⅱ) of Assumption 1.2, we have

$\frac{ dV(t)}{ dt} \leq \bar{p} \| i(\cdot,t)\|_{L^1}-c V(t)\leq \bar{p}\left(\frac{h}{\mu_0}+\frac{\bar{\xi}h}{\mu_0^2}\right) -c V(t).$

It follows (ⅳ) of Assumption 1.2, we have that

$\begin{align}\label{p} V(t) \leq&\ \frac{\bar{p}(\frac{h}{\mu_0}+\frac{\bar{\xi}h}{\mu_0^2})}{c}-e^{-c t}\left\{\frac{\bar{p}h}{c\mu_0}+ \frac{h\bar{p}\bar{\xi}}{c\mu_0^2} -V_0 \right\}\nonumber\\ \leq&\ \frac{\bar{p}h}{c\mu_0}+\frac{h\bar{p}\bar{\xi}}{c\mu_0^2}-e^{-\mu_0 t}\left\{ \frac{\bar{ph}}{c\mu_0}+\frac{h\bar{p}\bar{\xi}}{c\mu_0^2} -V_0 \right\} \end{align}$

(8)

Consequently, from (6), (7) and (8), we conclude that if $X_0\in \Xi$, then for $t\ge 0$,

$\begin{array}{*{20}{c}} {{{\left\| {{\Phi _t}({X_0})} \right\|}_{\cal Y}}}&{ \le \;\left( {1 + \frac{{\bar \xi }}{{{\mu _0}}} + \frac{{\bar p}}{c} + \frac{{\bar p\bar \xi }}{{c{\mu _0}}}} \right)\frac{h}{{{\mu _0}}}} \end{array}$

(9)

$\begin{array}{*{20}{l}} {\quad \quad \quad \quad \quad \quad - {e^{ - {\mu _0}t}}\left\{ {\left( {1 + \frac{{\bar \xi }}{{{\mu _0}}} + \frac{{\bar p}}{c} + \frac{{\bar p\bar \xi }}{{c{\mu _0}}}} \right)\frac{h}{{{\mu _0}}} - {X_0}{_{\cal Y}}} \right\}}\\ {\;\quad \quad \quad \quad = \frac{h}{{{{\tilde \mu }_0}}} - {e^{ - {\mu _0}t}}\left\{ {\frac{h}{{{{\tilde \mu }_0}}} - {{\left\| {{X_0}} \right\|}_{\cal Y}}} \right\} \le \frac{h}{{{{\tilde \mu }_0}}}.} \end{array}$

(10)

In summary, we have shown that $\Xi$ is positively invariant with respect to $\Phi$. Lastly, it follows from (9) that $\limsup\limits_{t\rightarrow\infty}\|\Phi_{t}(X_0)\|_{\mathcal{Y}} \leq \frac{h}{\tilde{\mu}_0}$ for any $X_0 \in \mathcal{Y}$, that is, $\Phi$ is point dissipative and $\Xi$ attracts all points in $\mathcal{Y}$. This completes the proof.

As a consequence of Proposition 2, we have the following result.

Proposition 3. Let $A \ge \frac{h}{\tilde{\mu}_0}$ be given. If $X_0\in \mathcal{Y}$ satisfying $\|X_0\|_{\mathcal{Y}} \leq A$, then the following statements hold for all $t\geq 0$.

(ⅰ) $T(t)$, $\|e(\cdot,t)\|_{L^1}$, $\|i(\cdot,t) \|_{L^1}$, $V(t)\leq A$;

(ⅱ) $M(t)\leq \bar{\xi}A$ and $N(t)\leq \bar{p}A$;

(ⅲ) $e(0,t)\leq f\beta A^2$ and $i(0,t)\leq (1-f)\beta A^2+\bar{\xi}A$.

---

2.4. Existence of equilibria

System (1) always has an infection-free equilibrium

$P^0 =(T^0, e^0(a), i^0(b), V^0): = (\frac{h}{d}, 0, 0, 0).$

The equations for an equilibrium are obtained from (1) by setting the time derivatives equal to $0$ with boundary conditions, that is, infection equilibrium $P^* = (T^*, e^*(\cdot), i^*(\cdot), V^*)\in \mathcal{Y}$ of (1) satisfies

$\left\{ \begin{array}{ll} \displaystyle h -d T^*-\beta T^*V^* = 0,\\ \displaystyle \frac{ d}{ da}e^*(a) = -\theta_1(a)e^*(a), \\ \displaystyle \frac{d}{da}i^*(b) = -\theta_2(b)i^*(b), \\ \displaystyle \int_0^\infty p(b)i^*(b) db = cV^*,\\ \displaystyle e^*(0) = f\beta T^*V^*,\\ \displaystyle i^*(0) = (1-f)\beta T^*V^*+\int_0^\infty\xi(a) e^*(a) da. \end{array} \right.$

(11)

Denote

$K = \int_0^\infty \xi(a)\Omega(a) da,\ \ J = \int_0^\infty p(b)\Gamma(b) db.$

Biologically, $K$ is the total number of infected cells activated by latency infected cells. $J$ accounts for the total number of virus particles produced by an infected cell during its life-span, i.e., the burst size.

We define basic reproduction number, $\Re_0$ of (1) as

$\Re_0 =\frac{f\beta T^0 K J}{c} +\frac{(1-f)\beta T^0J}{c}.$

which accounts for the total number of virons resulted from a single viron through the virus-to-cell infection mod. $1-f$ is the fraction of productive infection that leads to viral production, and $fK$ represents the contribution to productively infected cells from activation of latently infected cells. $\Re_0$ will serves as threshold value for (1), which completely determines the global behaviors of equilibria of (1).

Direct calculation yields that if $\Re_0>1$, then (1) admits a unique infection equilibrium $P^* = (T^*, e^*(a), i^*(b), V^*)$ with

$\begin{align}\label{*} &\ T^* = \frac{T^0}{\Re_0}, \ e^*(a) = fh\left(1-\frac{1}{\Re_0}\right)\Omega(a), \nonumber\\ &\ i^*(b) = (1-f+ fK)h\left(1-\frac{1}{\Re_0}\right)\Gamma(b), \ V^* = \frac{1}{c}\int_0^\infty p(b)i^*(b) db. \end{align}$

(12)

In summary, we have shown the following result.

Proposition 4. (ⅰ) System (1) always has an infection-free equilibrium $P^0$.

(ⅱ) If $\Re_0>1$, then (1) admits a unique infection equilibrium $P^*$, which is defined by (12).

---

3. Asymptotic smoothness of $\Phi(t,X_0)$

By Proposition 2 and 3, the semiflow is point-dissipative and $\Phi(\mathbb{R}_+\times B)$ is bounded for every bounded subset $B$ of $\mathcal{Y}$. By Theorem 3.4.6 in [5], the semiflow has a compact attractor of bounded sets if it is asymptotically smooth. To give the existence of compact attractor, we follow the approach in [43,Theorem 4.2 of Chapter IV].

Definition 3.1. [28] A set $A$ in $\mathcal{Y}$ is called a compact attractor of a set $B\subseteq X$ if $A$ is compact, invariant, and non-empty and $\Phi_t(B)\rightarrow A$ as $t\rightarrow\infty$. The last means that, for every open subset $U$ of $\mathcal{Y}$ with $A\subseteq U$, there is some $r > 0$ such that $\Phi_t(B)\subseteq U$ for all $t\geq r$ (i.e. $\Phi([r, \infty)\times B)\subseteq U$).

Recall that $M(t)$ and $N(t)$ are defined by (3). The following Proposition is devoted to prove basic properties of the functions $M(t)$ and $N(t)$ using Proposition 2, Assumption 1.2 and [38,Proposition 4.1].

Proposition 5. For any solution of (1), the associated functions $M(t)$ and $N(t)$ are Lipschitz continuous on $\mathbb{R}_{+}$.

Proof. Let $t\geq 0$ and $h> 0$. We can check that

$\begin{align}\label{ML} &M(t+h)-M(t) =\ \int_0^{\infty} \xi(a)e(a, t+h)da-\int_0^{\infty} \xi(a)e(a, t)da\nonumber\\ \leq &\ \int_0^h \xi(a)e(a, t+h)da+\int_h^{\infty} \xi(a)e(a, t+h)da-\int_0^{\infty} \xi(a)e(a, t)da\nonumber\\ \leq &\ \int_0^h \xi(a)e(0, t+h-a)\Omega(a)da\nonumber\\ &\ +\int_h^{\infty} \xi(a)e(a, t+h)da-\int_0^{\infty} \xi(a)e(a, t)da. \end{align}$

(13)

By applying $\xi(a)\leq \bar{\xi}, \ e(0,t)\leq f\beta A^2$ and $\Omega(a)\leq 1$ for the first integral, and making the substitution $\sigma = a-h$ for the second integral to (13), we get

$M(t + h) - M(t) \le f\beta {A^2}\bar \xi h + \int_0^\infty \xi (\sigma + h)e(\sigma + h,t + h)d\sigma - \int_0^\infty \xi (a)e (a,t)da$

It follows from (4) that

$e(\sigma+h,t+h) = e(\sigma,t)\frac{\Omega(\sigma+h)}{\Omega(\sigma)}.$

Thus,

$\begin{align*} M(t+h)-M(t) \leq&\ f\beta A^2\bar{\xi} h+\int_0^{\infty} \left(\xi(a+h)\frac{\Omega(a+h)}{\Omega(a)}-\xi(a)\right)e(a,t)da\nonumber\\ =&\ f\beta A^2\bar{\xi} h+\int_0^{\infty}\left(\xi(a+h)e^{-\int_{a}^{a+h}\theta_1(s)ds}-\xi(a)\right)e(a,t)da\nonumber\\ =&\ f\beta A^2\bar{\xi} h+\int_0^{\infty} \xi(a+h)\left(e^{-\int_{a}^{a+h}\theta_1(s)ds}-1\right)e(a,t)da\nonumber\\ &\ + \int_0^{\infty} (\xi(a+h)-\xi(a))e(a,t)da. \end{align*}$

From (ⅱ) of Assumption 1.2, we obtain $-\theta_1h\leq -\int_a^{a+h}\theta_1(s)ds \leq 0$. It follows that $1\geq e^{-\int_a^{a+h}\theta_1(s)ds}\geq e^{-\bar{\theta_1}h}\geq1-\bar{\theta_1}h$. Therefore,

$0\leq \xi(a+h)\left|e^{-\int_a^{a+h}\theta_1(s)ds}-1\right|\leq \bar{\xi}\bar{\theta_1}h.$

Recall that $\int_0^{\infty} e(a,t)da\leq \left\| \Phi_t (X_0) \right\|_{\mathcal{Y}}\leq A$. From (ⅲ) of Assumption 1.2, we obtain the following estimate,

$\begin{align*} M(t+h)-M(t)\leq f\beta A^2\bar{\xi}h+\bar{\xi}\bar{\theta_1}Ah+M_\xi Ah. \end{align*}$

Hence, $M(t)$ is Lipschitz continuous with Lipschitz coefficients $L_M = (\bar{\xi}f\beta A+\bar{\xi}\bar{\theta_1}+M_{\xi}) A$. Similarly, it is easy to check that $N(t)$ is Lipschitz continuous with Lipschitz coefficients $L_N = [\bar{p}(1-f)\beta A+\bar{p}\bar{\xi}+\bar{p}\bar{\theta_2}+M_{p}] A$.

Next we divide $\Phi:\ \mathbb{R}_+ \times \mathcal{Y} \to \mathcal{Y}$ into the following two operators $\Theta$, $\Psi:\ \mathbb{R}_+ \times \mathcal{Y} \to \mathcal{Y}$:

$\begin{eqnarray*} \Theta ( t,X_0 ) & := & (0,\tilde{\varphi}_e(\cdot,t),\tilde{\varphi}_i(\cdot,t),0 ), \\ \Psi ( t,X_0 ) & := & (T(t),\tilde{e}(\cdot,t),\tilde{i}(\cdot,t),V(t) ), \end{eqnarray*}$

where

$\begin{align*} \tilde{\varphi}_e ( a,t ) =\left\{ \begin{array}{ll} 0, \ \ &\ \mbox{if}\ \ t>a \geq 0,\\ e( a,t ),\ \ &\ \mbox{if}\ \ a \geq t \geq 0; \end{array} \right. &\ \tilde{\varphi}_i ( b,t ) =\left\{ \begin{array}{ll} 0,\ \ &\mbox{if}\ \ t>b \geq 0,\\ i( b,t ),\ \ &\mbox{if}\ \ b \geq t \geq 0; \end{array} \right. \end{align*}$

$\tilde e(a,t) = \left\{ {\begin{array}{*{20}{l}} {e(a,t),\;}&{{\rm{if}}\;\;t > a \ge 0,}\\ {0,\;}&{{\rm{if}}\;\;a \ge t \ge 0;} \end{array}} \right.\;\tilde i(b,t) = \left\{ {\begin{array}{*{20}{l}} {i(b,t),\;}&{{\rm{if}}\;\;t > b \ge 0,}\\ {0,}&{{\rm{if}}\;\;b \ge t \ge 0.} \end{array}} \right.$

Then $\Phi(t, X_0) =\Theta(t, X_0) +\Psi(t,X_0)$ for $t\geq 0$. Following the proof of [42,Proposition 3.13], we can arrive at the following main result of this section.

Theorem 3.2. For $X_0\in \Xi$, the orbit $\left\{ \Phi ( t, X_0 ) \ | \ t\geq 0 \right\}$ has a compact closure in $\mathcal{Y}$ if the following two conditions hold,

(ⅰ) There exists a function $\Delta:\ \mathbb{R}_+ \times \mathbb{R}_+ \to \mathbb{R}_+$ such that, for any $r>0$, $\lim_{t\to \infty}\Delta \left( t,r\right) = 0$ and if $X_0 \in \Omega$ with $\left\| X_0 \right\|_{\mathcal{Y}} \leq r$ then $\left\| \Theta \left( t,X_0 \right) \right\|_{\mathcal{Y}} \leq \Delta \left( t,r\right)$ for $t\geq 0$;

(ⅱ) For $t\geq 0$, $\Psi \left( t,\cdot \right)$ maps any bounded sets of $\Xi$ into sets with compact closure in $\mathcal{Y}$.

proof Proof of (ⅰ) of Theorem 3.2. Let $\Delta \left(t, r \right) = e^{-\mu_0 t} r$, then $\lim_{t\to \infty} \Delta( t,r) = 0$. By (4) and (5),

${\tilde \varphi _e}(a,t) = \left\{ {\begin{array}{*{20}{l}} {0,\;\;}&{{\rm{if}}\;\;t > a \ge 0,}\\ {{e_0}(a - t)\frac{{\Omega (a)}}{{\Omega (a - t)}},\;}&{{\rm{if}}\;\;a \ge t \ge 0;} \end{array}} \right.$

and

${\tilde \varphi _i}(b,t) = \left\{ {\begin{array}{*{20}{l}} {0,\;\;}&{{\rm{if}}\;\;t > b \ge 0,}\\ {{i_0}(b - t)\frac{{\Gamma (b)}}{{\Gamma (b - t)}},\;}&{{\rm{if}}\;\;b \ge t \ge 0.} \end{array}} \right.$

Then, for $X_0 \in \Xi$ satisfying $\left\| X_0\right\|_{\mathcal{Y}} \leq r$ and for $t\ge 0$, we have

$\begin{align*} \left\| \Theta \left( t,X_0 \right) \right\|_{\mathcal{Y}} = &\ \left| 0\right| +\| \tilde{\varphi}_e(\cdot,t)\|_{L^1} +\| \tilde{\varphi}_i ( \cdot,t ) \|_{L^1} + \left| 0\right| \\ = &\ \int_t^\infty\left| e_0 (a-t)\frac{\Omega(a)}{\Omega(a-t)}\right| da +\int_t^\infty\left| i_0 (b-t)\frac{\Gamma(b)}{\Gamma(b-t)}\right| db \\ =&\ \int_0^\infty\left| e_0 (\sigma)\frac{\Omega(\sigma+t)}{\Omega(\sigma)}\right| d\sigma + \int_0^\infty\left| i_0 (\sigma)\frac{\Gamma(\sigma+t)}{\Gamma(\sigma)}\right| d\sigma \\ = &\ \int_0^\infty\left| e_0 (\sigma)e^{-\int_\sigma^{\sigma+t}\theta_1(\tau) d\tau}\right| d\sigma + \int_0^\infty\left| i_0 (\sigma)e^{-\int_\sigma^{\sigma+t}\theta_2(\tau) d\tau}\right| d\sigma \\ \leq &\ e^{-\mu_0t}\| e_0\|_{L^1} + e^{-\mu_0t}\| i_0\|_{L^1} \\ \leq &\ e^{-\mu_0t}\left\| X_0 \right\|_{\mathcal{Y}}. \end{align*}$

Proof of (ⅱ) of Theorem 3.2. It is sufficient to show that $\Psi \left( t,\cdot \right)$ maps any bounded sets of $\Xi$ into sets with compact closure in $\mathcal{Y}$. From Proposition 2, $T(t)$ and $V(t)$ remains in the compact set $\left[ 0, h/\tilde{\mu}_0 \right]\subset \left[ 0, A \right]$. Thus it remains unknown that whether $\tilde{e}\left( a,t\right)$ and $\tilde{i}\left( b,t\right)$ remain in a precompact subset of $L_+^1 \left( 0,\infty \right)$, which is independent of $X_0 \in \Xi$. To this end, we next to verify the following conditions for $\tilde{e}\left(a, t\right)$ and similar ones for $\tilde{i}\left( b,t\right)$ (see, for example, [25,Theorem B.2]).

(ⅰ) The supremum of $\|\tilde{e}\left( \cdot,t\right)\|_{L^1}$ with respect to $X_0 \in \Xi$ is finite;

(ⅱ) $\lim_{h\to \infty} \int_{h}^{\infty} \tilde{e}\left( a,t\right) da = 0$ uniformly with respect to $X_0 \in \Xi$;

(ⅲ) $\lim_{h\to 0+} \int_{0}^{\infty} \left| \tilde{e} \left( a+h,t \right) - \tilde{e} \left( a,t\right) \right| da = 0$ uniformly with respect to $X_0 \in \Xi$;

(ⅳ) $\lim_{h\to 0+} \int_{0}^{h} \tilde{e}\left( a,t\right) da = 0$ uniformly with respect to $X_0 \in \Xi$.

It follows from (4), (5), Proposition 3 and (2) that $\tilde{e}(a,t)\leq f\beta A^2e^{-\mu_0 a}$, $\tilde{i}(b,t)\leq \left[(1-f)\beta A^2+\bar{\xi}A\right]e^{-\mu_0 b}$. Thus, (ⅰ), (ⅱ) and (ⅳ) are directly satisfied.

Next we verify condition (ⅲ). For sufficiently small $h\in \left( 0,t\right)$, we have

$\begin{array}{l} \int_0^\infty {\left| {\tilde e(a + h,t) - \tilde e(a,t)} \right|} da = \int_0^{t - h} | e(a + h,t) - e(a,t)|da + \int_{t - h}^t | 0 - e(a,t)|da\\ \quad \quad \quad \quad \quad \quad \quad \quad \quad \quad = \int_0^{t - h} | e(0,t - a - h)\Omega (a + h) - e(0,t - a)\Omega (a)|da\\ \quad \quad \quad \quad \quad \quad \quad \quad \quad \quad \quad + \int_{t - h}^t | e(0,t - a)\Omega (a)|da\\ \quad \quad \quad \quad \quad \quad \quad \quad \quad \quad \le \;{\Delta _1} + {\Delta _2} + f\beta {A^2}h, \end{array}$

where

$\Delta_1 = \int_0^{t-h} e(0,t-a-h)\Big|\Omega(a+h)-\Omega(a)\Big| da$

and

$\Delta_2 =\int_0^{t-h} \Big| e(0,t-a-h)-e(0,t-a)\Big|\Omega(a)da.$

We first get an estimate of $\Delta_1$. Since

$\begin{align*} \int_0^{t-h}|\Omega(a+h)-\Omega(a)| da =&\ \int_0^{t-h}\Big(\Omega(a)-\Omega(a+h)\Big) da\\ = &\ \int_0^{t-h}\Omega(a) da-\int_h^t\Omega(a) da\\ = &\ \int_0^{t-h}\Omega(a) da-\int_h^{t-h}\Omega(a) da -\int_{t-h}^t\Omega(a) da\\ = &\ \int_0^h\Omega(a) da-\int_{t-h}^t\Omega(a) da \\ \leq &\ h, \end{align*}$

it follows from Proposition 3 that

$\Delta_1\le f\beta A^2h.$

Next we estimate $\Delta_2$. We rewrite $\Delta_2$ as

$\begin{align*} \Delta_2 =&\ \int_0^{t-h} \Big|f\beta T(t-a-h)V(t-a-h)-f\beta T(t-a)V(t-a)\Big|\Omega(a) da. \end{align*}$

It is easy to see that $T(t)$ and $V(t)$ are both Lipschitz continuous on $\mathbb{R}_+$ with Lipschitz constants $M_T = h+dA+\beta A^2$ and $M_V = (\bar{p}+c)A$, respectively. According to [13,Proposition 6], we conclude that $T(t)V(t)$is Lipschitz continuous with Lipschitz constants $M_{TV} = AM_V+A M_T$. Denote that $G = f\beta M_{TV}$. This estimate immediately yields

$\Delta_2\le Gh \int_0^{t-h} e^{-\mu_0a}da \le \frac{Gh}{\mu_0}.$

Hence

$\int_0^\infty \left| \tilde{e}(a+h,t)-\tilde{e}(a,t)\right| da\leq\ \left(2f\beta A^2+\frac{G}{\mu_0}\right) h,$

and condition (ⅲ) directly follows.

As to $\tilde{i}\left( b,t\right)$, we have

$\begin{array}{l} \int_0^\infty {\left| {\tilde i(b + h,t) - \tilde i(b,t)} \right|} db = \;\int_0^{t - h} | i(b + h,t) - i(b,t)|db + \int_{t - h}^t | 0 - i(b,t)|db\\ \quad \quad \quad \quad \quad \quad \quad \quad \quad \;\; = \;\int_0^{t - h} | i(0,t - b - h)\Gamma (b + h) - i(0,t - b)\Gamma (b)|db\\ \quad \quad \quad \quad \quad \quad \quad \quad \quad \;\;\quad \; + \int_{t - h}^t | i(0,t - b)\Gamma (b)|db\\ \quad \quad \quad \quad \quad \quad \quad \quad \quad \;\; \le {\Upsilon _1} + {\Upsilon _2} + \left[ {(1 - f)\beta {A^2} + \bar \xi A} \right]h, \end{array}$

where

$\Upsilon_1 = \int_0^{t-h} i(0,t-b-h)\Big|\Gamma(b+h)-\Gamma(b)\Big| db$

and

$\Upsilon_2 =\int_0^{t-h} \Big| i(0,t-b-h)-i(0,t-b)\Big|\Gamma(b)db.$

Similarly, we have $\int_0^{t-h}\Big|\Gamma(b+h)-\Gamma(b)\Big|db\leq h$. Hence from Proposition 3, we can conclude that

$\Upsilon_1\le \left[(1-f)\beta A^2+\bar{\xi}A\right] h.$

Next we estimate $\Upsilon_2$. Firstly, we have

$\begin{array}{l} {\Upsilon _2} = \;\int_0^{t - h} | (1 - f)\beta T(t - a - h)V(t - a - h) + M(t - a - h)\\ \quad \quad \quad \quad \quad - (1 - f)\beta T(t - a)V(t - a) - M(t - a)|\Gamma (b)da\\ \quad \quad \le (1 - f)\beta \int_0^{t - h} | T(t - a - h)V(t - a - h) - T(t - a)V(t - a)|\Gamma (b)da\\ \quad \quad \quad \quad + \int_0^{t - h} | M(t - a - h) - M(t - a)|\Gamma (b)da. \end{array}$

As before, $M_{TV} = AM_V+A M_T$. Recall that $M(t)$ is Lipschitz continuous on $\mathbb{R}_{+}$ with Lipschitz constants $L_M = (\bar{\xi}f\beta A+\bar{\xi}\bar{\theta_1}+M_{\xi}) A$. Set $H = (1-f)\beta M_{TV}+L_M$, By a zero-trick, then we have

$\Upsilon_2\leq M_1h\int_0^{t-h}e^{-\mu_0b}db\leq \frac{Hh}{\mu_0}.$

Finally, we have

$\int_0^\infty \left| \tilde{i}(b+h,t)-\tilde{i}(b,t)\right| db \leq \left\{2\left[(1-f)\beta A^2+\bar{\xi}A\right] +\frac{H}{\mu_0}\right\}h,$

thus condition (ⅲ) directly follows. This completes the proof.

Consequently, we have the following theorem for the semi-flow $\{\Phi(t)\}_{t\geq 0}$, which establish the existence of global attractors by Smith and Thieme [25].

Theorem 3.3. The semi-flow $\{\Phi(t)\}_{t\geq 0}$ has a global attractor $\mathcal{A}$ in $\mathcal{Y}$, which attracts any bounded subset of $\mathcal{Y}$.

---

4. The uniform persistence

This section is spent on proving that (1) is uniformly persistent under the condition $\Re_0>1$, which indicates that $\Re_0>1$ is a threshold index for infection persistence.

Let $\hat{e}(t): = e(0,t)$ and $\hat{i}(t): = i(0,t)$. We rewrite the first three equations of (1) as

$\left\{ {\begin{array}{*{20}{l}} {\frac{{dT(t)}}{{dt}} = h - dT(t) - \frac{1}{f}\hat e(t),}\\ {e(a,t) = \left\{ {\begin{array}{*{20}{l}} {\hat e(t - a)\Omega (a),\;\;}&{{\rm{if}}\;\;t \ge a \ge 0,}\\ {{e_0}(a - t)\frac{{\Omega (a)}}{{\Omega (a - t)}},\;}&{{\rm{if}}\;\;a \ge t \ge 0;} \end{array}} \right.}\\ {i(b,t) = \left\{ {\begin{array}{*{20}{l}} {\hat i(t - b)\Gamma (b),\;}&{{\rm{if}}\;\;t \ge b \ge 0,}\\ {{i_0}(b - t)\frac{{\Gamma (b)}}{{\Gamma (b - t)}},\;}&{{\rm{if}}\;\;b \ge t \ge 0,} \end{array}} \right.} \end{array}} \right.$

(14)

where

$\begin{align}\label{e_int_1} \hat{e}(t) =&\ f\beta T(t) V(t) \end{align}$

(15)

and

$\begin{array}{l} \begin{array}{*{20}{l}} {\hat i(t) = }&{\;(1 - f)\beta T(t)V(t) + \int_0^t \xi (a)\Omega (a)\hat e(t - a)da} \end{array}\\ \quad \quad \quad \quad \quad \quad + \int_t^\infty \xi (a)\frac{{\Omega (a)}}{{\Omega (a - t)}}{e_0}(a - t)da. \end{array}$

(16)

Lemma 4.1. If $\Re_0 > 1$, then there exists a positive constant $\epsilon_0 > 0$ such that

$\limsup\limits_{t \to \infty} \hat{e}(t) > \epsilon_0.$

(17)

Proof. We first get an estimate on $\hat{i}(t)$ as follows. By (16), we have

$\hat{i}(t) \geq (1-f)\beta T(t) V(t)+\int_{0}^{t} \xi(a) \Omega(a) \hat{e}(t-a) da.$

(18)

Solving the fourth equation of (1) with initial condition $V(0) = V_0$, we have that

$V(t) = V_0e^{-ct}+\int_0^t\int_0^\infty p(b)i(b,\tau)db\cdot e^{-c(t-\tau)}d\tau$

Then

$V(t) \geq \int_0^t e^{-c(t-\tau)}\int_0^\tau p(b)i(b,\tau) db d\tau =\int_0^t e^{-c(t-\tau)}\int_0^\tau p(b)\Gamma(b)\hat{i}(\tau-b) db d\tau.$

This, combined with (18), gives us

$\begin{array}{*{20}{l}} {\hat i(t) \ge }&{\;(1 - f)\beta T(t)\int_0^t {{e^{ - c(t - \tau )}}} \int_0^\tau p (b)\Gamma (b)\hat i(\tau - b)dbd\tau }\\ {}&{ + f\beta \int_0^t \xi (a)\Omega (a)T(t - a)\int_0^{t - a} {{e^{ - c(t - a - \tau )}}} \int_0^\tau p (b)\Gamma (b)\hat i(\tau - b)dbd\tau da.} \end{array}$

(19)

Since $\Re_0 > 1$, there exists a sufficiently small $\epsilon_1> 0(\epsilon_1 = \frac{1}{f}\epsilon_0)$ such that

$\frac{(1-f)\beta}{c} \frac{h-\epsilon_1}{d} \int_{0}^{\infty} p(b) \Gamma(b) db +\frac{f\beta}{c} \frac{h-\epsilon_1}{d} \int_{0}^{\infty} p(b) \Gamma(b) db \int_0^\infty \xi(a)\Omega(a) da> \ 1. \label{R0_e}$

We claim that (17) holds for this $\varepsilon_0$.Suppose that there exists a $T>0$ such that

$\hat{e}(t) \leq \epsilon_0\ \ \mbox{for all}\ \ t\geq T.$

Then it follows from (14) that $\frac{ dT(t)}{ dt} \geq h -d T(t)-\epsilon_1$ for $t\geq T$. This implies that $\liminf\limits_{t\to\infty} T(t)\ge \frac{h-\varepsilon_1}{d}$. Thus there exists $\hat{T}>T$ such that $T(t) \geq \frac{h-\epsilon_1}{d}$ for all $t\ge \hat{T}$ and hence (19) becomes

$\begin{align}\label{i_int_3} \hat{i}(t) \geq &\ (1-f)\beta \frac{h-\epsilon_1}{d} \int_{0}^{t} e^{-c(t-\tau)} \int_{0}^{\tau} p(b) \Gamma(b) \hat{i}(\tau-b) db d\tau \nonumber \\ &\ + f\beta \frac{h-\epsilon_1}{d}\int_{0}^{t} \xi(a) \Omega(a)\int_{0}^{t-a}e^{-c(t-a-\tau)}\int_0^\tau p(b)\Gamma(b)\hat{i}(\tau-b) db d\tau da \end{align}$

(20)

for all $t\geq \hat{T}$. Without loss of generality, we can assume that (20) holds for all $t\ge 0$ (just replace $X_0$ by $\Phi(\hat{T}, X_0)$). Then taking the Laplace transforms of both sides of (20), we obtain

$\begin{array}{l} {\cal L}[\hat i] \ge \;(1 - f)\beta \frac{{h - {_1}}}{d}\int_0^\infty {{e^{ - \lambda t}}} \int_0^t {{e^{ - c(t - \tau )}}} \int_0^\tau p (b)\Gamma (b)\hat i(\tau - b)dbd\tau dt + \\ \quad \quad f\beta \frac{{h - {_1}}}{d}\int_0^\infty {{e^{ - \lambda t}}} \int_0^t \xi (a)\Omega (a)\int_0^{t - a} {{e^{ - c(t - a - \tau )}}} \int_0^\tau p (b)\Gamma (b)\hat i(\tau - b)dbd\tau dadt \end{array}$

$\begin{align*} = &\ (1-f)\beta \frac{h-\epsilon_1}{d}\frac{1}{c+\lambda} \int_{0}^{\infty} e^{-\lambda b} p(b) \Gamma(b) db\mathcal{L}[\hat{i}] \\ &\ +f\beta \frac{h-\epsilon_1}{d}\frac{1}{c+\lambda} \int_{0}^{\infty} e^{-\lambda b} p(b) \Gamma(b) db \int_0^\infty e^{-\lambda a} \xi(a) \Omega(a) da \mathcal{L}[\hat{i}]. \end{align*}$

Here $\mathcal{L}[\hat{i}]$ denotes the Laplace transform of $\hat{i}$, which is strictly positive because of (15) and Assumption 1.2. Dividing both sides of the above inequality by $\mathcal{L}[\hat{i}]$ and letting $\lambda \to 0$ give us

$1\ge \frac{(1-f)\beta }{c}\frac{h-{{\epsilon }_{1}}}{d}\int_{0}^{\infty }{p}(b)\Gamma (b)db+\frac{f\beta }{c}\frac{h-{{\epsilon }_{1}}}{d}\int_{0}^{\infty }{p}(b)\Gamma (b)db\int_{0}^{\infty }{\xi }(a)\Omega (a)da,$

which yields a contradiction.

In order to apply a technique used by Smith and Thieme [25,Chapter 9] (see also McCluskey [13,Section 8]), we consider a total $\Phi$-trajectory of (1) in space $\mathcal{Y}$. A total trajectory of $\Phi$ is a function $X:\mathbb{R}\rightarrow\mathcal{Y}$ such that $\Phi_s(X(t)) = X(t+s)$ for all $t\in \mathbb{R}$ and all $s\geq0$. For a non-empty compact set $\tilde{A}$, it is said to be a compact attractor of a class $\mathcal{C}$ of set if $\tilde{A}$ is invariant and $\textrm{d}\left(\Phi_t(C),\tilde{A}\right)\rightarrow0$ for each $C\in \mathcal{C}$. For each $X_0\in \tilde{A}$, there exists a total trajectory $X$ such that $X(0) = X_0$ and $X(t)\in\tilde{A}$ for all $t\in \mathbb{R}$.

Let $\phi: \mathbb{R}\to \mathcal{Y}$ be a total $\Phi$-trajectory such that $\phi(r) = \left( T(r),e(\cdot,r),i(\cdot,r),V(r) \right)$, $r\in \mathbb{R}$. Then

$\begin{align*} \phi(r+t) =&\ \Phi(t,\phi(r))\ \ \mbox{for}\ \ t\geq 0\ \ \mbox{and}\ \ r\in \mathbb{R}, \\ e(a,r) =&\ e\left(0,r-a \right) \Omega (a) =\hat{e}(r-a) \Omega (a)\ \ \mbox{for}\ \ r\in \mathbb{R} \ \mbox {and}\ \ a \geq 0, \\ i(b,r) =&\ i\left(0,r-b \right) \Gamma (b) =\hat{i}(r-b) \Gamma (b)\ \ \mbox{for} \ \ r\in \mathbb{R}\ \ \mbox {and}\ \ b \geq 0. \end{align*}$

So it follows from (14)-(15) that

$\left\{ \begin {array}{lll} \displaystyle \frac{ dT(r)}{ dr} =h -d T(r)-\frac{1}{f}\hat{e}(r), \\ \displaystyle \hat{e}(r) =f\beta T(r) V(r), \\ \displaystyle \hat{i}(r) =(1-f)\beta T(r) V(r)+\int_{0}^{\infty} \xi(a) \Omega(a) \hat{e}(r-a) da, \\ \displaystyle \frac{ dV(r)}{ dr} =\int_{0}^{\infty} p(b) \Gamma(b) \hat{i}(r-b) db-c V(r),\ \ \mbox{for}\ \ r\in \mathbb{R}. \end {array}\right.$

By the similar arguments as in McCluskey [13,Section 8] and Wang et al. [38,Section 5] and a slight modification of the proof in [17, Lemma 4.1], actually, a total $\Phi$-trajectory $\phi$ enjoys the following nice properties.

Thus Lemma 4.1 tells us that if $\Re_0>1$ then the semi-flow $\Phi$ is uniformly weakly $\rho$-persistent. Moreover, with the help of Theorem 3.3 and the Lipschitz continuity of $\hat{i}$ (which immediately follows from Proposition 5), we can apply Theorem 5.2 of Smith and Thieme [25] to conclude that the uniform weak $\rho$-persistence of the semi-flow $\Phi$ implies its uniform (strong) $\rho$-persistence, that is, we have obtained the following result.

Theorem 4.2. If $\Re_0 >1$, then the semi-flow $\Phi$ is uniformly (strongly) $\rho$-persistent.

When $\Re_0>1$, the uniform persistence of (1) immediately follows from Theorem 4.2. In fact, it follows from (14) that $\left\| e(\cdot,t) \right\|_{L^1} \geq \int_{0}^{t} \hat{e}(t-a) \Omega (a) da$ and hence from a variation of the Lebesgue-Fatou lemma [26,Section B.2] we get

$\liminf\limits_{t\to \infty} \left\| e(\cdot ,t) \right\|_{L^1} \geq \hat{e}^\infty \int_{0}^{\infty} \Omega (a) da,$

where $\hat{e}^\infty =\liminf_{t\to \infty} \hat{e}(t)$. Under Theorem 4.2, there exists a positive constant $\epsilon > 0$ such that $\hat{e}^\infty > \epsilon$ if $\Re_0 > 1$ and hence the persistence of $e(a,t)$ with respect to $\left\| \cdot \right\|_{L^1}$ follows. By a similar argument, we can prove that $T(t)$ and $V(t)$ are persistent with respect to $|\cdot|$ and $i(a,t)$ is persistent with respect to $\|\cdot\|_{L^1}$. In a summary, we get the following result.

Theorem 4.3. If $\Re_0 > 1$, then the semiflow $\{\Phi(t)\}_{t\geq0}$ is uniformly persistent in $\mathcal{Y}$, that is, there exists a constant $\epsilon > 0$ such that, for each $X_0 \in \mathcal{Y}$,

$\liminf\limits_{t\to \infty} T(t)\geq \epsilon, \liminf\limits_{t\to \infty}\|e(\cdot,t)\|_{L^1} \geq \epsilon, \liminf\limits_{t\to \infty}\|i(\cdot,t)\|_{L^1} \geq \epsilon, \liminf\limits_{t\to \infty} V(t)\geq \epsilon.$

---

5. The local stability of equilibria

This section is devoted to investigate the local stability of equilibria of (1).

Theorem 5.1. (ⅰ) If $\Re_0 < 1$, the infection-free equilibrium $P^0$ of (1) is locally asymptotically stable while it is unstable if $\Re_0 > 1$.

(ⅱ) If $\Re_0 > 1$, the infection equilibrium $P^*$ of (1) is locally asymptotically stable.

Proof. Proof. of (ⅰ) of Theorem 5.1. Linearizing (1) around the infection-free equilibrium $P^0$ by using

$\begin{align*} x_1(t) = T(t)-\frac{h}{d},\ x_2(a,t) = e(a, t),\ x_3(b, t) = i(b,t),\ x_4(t) = V(t), \end{align*}$

we get

$\left\{ \begin{array}{lll} \displaystyle \frac{ dx_1(t)}{ dt} =-d x_1(t)-\frac{\beta h}{d}x_4(t), \\ \displaystyle \left(\frac{\partial}{\partial t}+\frac{\partial}{\partial a}\right)x_2(a,t) =-\theta_1(a)x_2(a,t),\\ \displaystyle \left(\frac{\partial}{\partial t}+\frac{\partial}{\partial b}\right)x_3(b,t) =-\theta_2(b)x_3(b,t),\\ \displaystyle \frac{ dx_4(t)}{ dt} =\int_0^\infty p(b)x_3(b,t) db-cx_4(t),\\ \displaystyle x_2(0,t) = f\beta\frac{h}{d}x_4(t),\\ \displaystyle x_3(0,t) = (1-f)\beta\frac{h}{d}x_4(t)+\int_0^\infty \xi(a)x_2(a,t)da. \end{array}\right.$

(21)

Set

$\begin{align}\label{z1} x_1(t) = x_1^0e^{\lambda t}, ~~x_2(a, t) = x_2^0(a)e^{\lambda t}, ~~x_3(b, t) = x_3^0(b)e^{\lambda t}, ~~x_4(t) = x_4^0 e^{\lambda t}, \end{align}$

(22)

where $x_1^0, x_2^0(a), x_3^0(b), x_4^0$ are to be determined later. Substituting (22) into (21), we have

$\begin{align*} \lambda x_1^0 = -dx_1^0-\frac{h\beta}{d}x_4^0, \end{align*}$

$\left\{ \begin{array}{ll} \displaystyle \lambda x_2^0(a)+\frac{dx_2^0(a)}{d a} = -\theta_1(a)x_2^0(a),\\ \displaystyle x_2^0(0) = f \beta\frac{ h}{d}x_4^0, \end{array}\right.$

(23)

$\label{Lsx_3} \left\{ \begin{array}{ll} \displaystyle \lambda x_3^0(b)+\frac{dx_3^0(b)}{d b} = -\theta_2(b)x_3^0(b),\\ \displaystyle x_3^0(0) = (1-f)\frac{h \beta}{d}x_4^0+\int_0^\infty \xi(a)x_2^0(a)da, \end{array}\right.$

(24)

$\begin{align}\label{Lsx_4} \lambda x_4^0 = \int_0^\infty p(b)x_3^0(b)db-cx_4^0. \end{align}$

(25)

We integrate the first equation of Equ (23) and Equ (24) from $0$ to $a$,

$\begin{align} \label{Ls2} x_2^0(a) =&\ x_2^0(0)e^{-\lambda a-\int_0^a\theta_1(s)ds}, \end{align}$

(26)

and

$\begin{align} \label{Ls3} x_3^0(b) =&\ x_3^0(0)e^{-\lambda b-\int_0^b\theta_2(s)ds} \notag\\ = &\ \left[(1-f)\frac{h \beta}{d}x_4^0+\int_0^\infty \xi(a)x_2^0(a)da\right]e^{-\lambda b-\int_0^b\theta_2(s)ds}. \end{align}$

(27)

From (25), (26) and (27),

$\begin{align} \label{Ls4} x_4^0 =&\ \frac{\int_0^\infty p(b)x_3^0(b)db}{\lambda+c} \notag\\ = &\ \frac{1-f}{\lambda+c}\frac{h \beta}{d}x_4^0\int_0^\infty p(b)e^{-\lambda b-\int_0^b\theta_2(s)ds}db\notag\\ &\ +\frac{x_2^0(0)}{\lambda+c}\frac{h}{d}\int_0^\infty \xi(a)e^{-\lambda a-\int_0^a\theta_1(s)ds}da\int_0^\infty p(b)e^{-\lambda b-\int_0^b\theta_2(s)ds} db. \end{align}$

(28)

Combining (23) into (28), it follows the equation that

$\mathcal{W}(\lambda) = 1,\label{CE}$

(29)

where

$\begin{align*} \mathcal{W}(\lambda) = &\ \frac{f\beta}{\lambda+c}\frac{h}{d}\int_0^\infty \xi(a) e^{-\lambda a-\int_0^a\theta_1(s) ds} da \int_0^\infty p(b)e^{-\lambda b-\int_0^b\theta_2(s) ds} db \\ &\ +\frac{1-f}{\lambda+c}\frac{h\beta}{d} \int_0^\infty p(b)e^{-\lambda b-\int_0^b\theta_2(s) ds} db. \end{align*}$

Since $\mathcal{W}$ have the properties that

$\lim\limits_{\lambda\to\infty}\mathcal{W}(\lambda) = 0,\ \lim\limits_{\lambda\to-\infty}\mathcal{W}(\lambda) = \infty,\ \mathcal{W}'(\lambda)<0,$

Thus (29) admits a unique real root, $\lambda^*$. Recall that $\mathcal{W}(0) = \Re_0$, it follows that, $\lambda^*<0$ if $\Re_0<1$ and $\lambda^*>0$ if $\Re_0>1$, that is, $P^0$ is unstable if $\Re_0>1$. Suppose that $\Re_0<1$. Let $\lambda = \mu+\nu i$ be an arbitrary complex root of (29). It is easy to see that

$1 = |\mathcal{W}(\lambda)| = |\mathcal{W}(\mu+\nu i)|\leq \mathcal{W}(\mu),$

which implies that $0>\lambda^*\ge \mu$. Hence all roots of (29) have negative real parts, that is $P^0$ is locally asymptotically stable if $\Re_0<1$.

Proof of (ⅱ) of Theorem 5.1. Linearizing the system (1) at $P^*$ by using

$\begin{align*} &\ y_1(t) = T(t)-T^*,\ y_2(a,t) = e(a, t)-e^*(a),\\ &\ y_3(b, t) = i(b,t)-i^*(b),\ y_4(t) = V(t)-V^*, \end{align*}$

we get

$\left\{ \begin{array}{lll} \displaystyle \frac{ dy_1(t)}{ dt} =-d \Re_0 y_1(t)-\beta T^*y_4(t), \\ \displaystyle \left(\frac{\partial}{\partial t}+\frac{\partial}{\partial a}\right)y_2(a,t) =-\theta_1(a)y_2(a,t),\\ \displaystyle \left(\frac{\partial}{\partial t}+\frac{\partial}{\partial b}\right)y_3(b,t) =-\theta_2(b)y_3(b,t),\\ \displaystyle \frac{ dy_4(t)}{ dt} =\int_0^\infty p(b)y_3(b,t) db-cy_4(t),\\ \displaystyle y_2(0,t) = fd(\Re_0-1)y_1(t)+f\beta T^*y_4(t),\\ \displaystyle y_3(0,t) = (1-f)d(\Re_0-1)y_1(t)+(1-f)\beta T^*y_4(t)+\int_0^\infty \xi(a)y_2(a,t)da, \end{array}\right.$

(30)

Set

$\begin{align}\label{z2} y_1(t) = y_1^0e^{\lambda t}, ~~y_2(a, t) = y_2^0(a)e^{\lambda t}, ~~y_3(b, t) = y_3^0(b)e^{\lambda t}, ~~y_4(t) = y_4^0 e^{\lambda t}, \end{align}$

(31)

where $y_1^0, y_2^0(a), y_3^0(b), y_4^0$ are to be determined. Substituting (31) into (30) yields

$\begin{align}\label{Lsx_5} \lambda y_1^0 = -d\Re_0y_1^0-\beta T^*y_4^0, \end{align}$

(32)

$\left\{ \begin{array}{ll} \displaystyle \lambda y_2^0(a)+\frac{dy_2^0(a)}{d a} = -\theta_1(a)y_2^0(a),\\ \displaystyle y_2^0(0) = fd(\Re_0-1)y_1^0+f\beta T^*y_4^0, \end{array}\right.$

(33)

$\left\{ {\begin{array}{*{20}{l}} {\lambda y_3^0(b) + \frac{{dy_3^0(b)}}{{db}} = - {\theta _2}(b)y_3^0(b),}\\ {y_3^0(0) = (1 - f)d({\Re _0} - 1)y_1^0 + (1 - f)\beta {T^*}y_4^0 + y_2^0(0)\int\limits_0^\infty \xi (a){e^{ - \lambda a - \int\limits_0^a {{\theta _1}} (s)ds}}da,} \end{array}} \right.$

(34)

and

$\begin{align}\label{Lsx_8} \lambda y_4^0 = \int_0^\infty p(b)y_3^0(b)db-cy_4^0. \end{align}$

(35)

We integrate the first equation of (33), (34) from $0$ to $a$,

$\begin{align*} y_2^0(a) =&\ y_2^0(0)e^{-\lambda a-\int_0^a \theta_1(s)ds}, \end{align*}$

and

$\begin{align*} y_3^0(b) =&\ y_3^0(0)e^{-\lambda b-\int_0^b\theta_2(s)ds} \notag\\ = &\ \left[(1-f)d(\Re_0-1)y_1^0+(1-f)\beta T^*y_4^0\right]e^{-\lambda b-\int_0^b\theta_2(s)ds}\\ &\ +y_2^0(0)\int_0^\infty \xi(a)e^{-\lambda a-\int_0^a \theta_1(s)ds}da\cdot e^{-\lambda b-\int_0^b\theta_2(s)ds}. \end{align*}$

and from (35), we have

$\begin{align} \label{Ls8} y_4^0 =&\ \frac{\int_0^\infty p(b)y_3^0(b)db}{\lambda+c} \notag\\ = &\ \frac{1-f}{\lambda+c}\left(d(\Re_0-1)y_1^0+\beta T^*y_4^0\right)\int_0^\infty p(b)e^{-\lambda b-\int_0^b\theta_2(s)ds}db\nonumber\\ & +\frac{y_2^0(0)}{\lambda+c}\int_0^\infty \xi(a)e^{-\lambda a-\int_0^a\theta_1(s)ds}da\int_0^\infty p(b)e^{-\lambda b-\int_0^b\theta_2(s)ds} db. \end{align}$

(36)

Combining (32), (33) into (36), yields the characteristic equation at $P^*$ that

$G(\lambda) = (\lambda+d)\mathcal{W}_1(\lambda)-\lambda-d\Re_0 = 0,$

(37)

where

$\begin {eqnarray*} \mathcal{W}_1(\lambda) & = &\ \frac{(1-f)\beta T^*}{\lambda+c} \int_0^\infty p(b)e^{-\lambda b-\int_0^b\theta_2(s) ds} db \\ && +\frac{f\beta T^*}{\lambda+c}\int_0^\infty \xi(a)e^{-\lambda a-\int_0^a\theta_1(s) ds} da \int_0^\infty p(b)e^{-\lambda b-\int_0^b\theta_2(s) ds} db. \end{eqnarray*}$

It is sufficient to show that (37) has no roots with non-negative real parts. Suppose that it has a root $\lambda = \mu+\nu i$ with $\mu\ge 0$. Then we have

$(\mu+\nu i+d)\mathcal{W}_1(\mu+\nu i)-\mu -\nu i -d\Re_0 = 0.$

Separating the real part of the above equality gives

$\mathrm{Re}\ \mathcal{W}_1(\mu+\nu i) = \frac{(\mu+d\Re_0)(\mu+d)+\nu^2}{(\mu+d)^2+\nu^2}>1.$

(38)

Noticing that $\mathcal{W}_1(0) =T^* \frac{\Re_0}{T_0} =1$ and $\mathcal{W}_1$ is a decreasing function, we have

$\mathrm{Re}\ \mathcal{W}_1(\mu+\nu i)\leq |\mathcal{W}_1(\mu)| = \mathcal{W}_1(\mu)\leq \mathcal{W}_1(0) = 1,$

which yields a contradiction. This completes the proof.

---

6. Global stability of equilibria

This section is devoted to investigate the global stability of the equilibria by using Lyapunov functionals under the threshold value. In what follows, we introduce an important function $g$ on $(0,\infty)$ defined by $g(x) = x-1-\ln x$ for $x\in(0,\infty)$. This function is continuous and concave up with $g(1) = 0$. By Theorem 5.1, it is suffice to show that equilibria of (1) are globally attractive in $\mathcal{Y}$.

Theorem 6.1. The infection-free equilibrium $P^0$ of (1) is globally attractive if $\Re_0 \leq 1$.

Proof. Considering the candidate Lyapunov function as follows,

$\mathcal{L}_{IFE}(t) = \mathcal{L}_1(t)+\mathcal{L}_2(t)+\mathcal{L}_3(t)+\mathcal{L}_4(t),$

where $\mathcal{L}_1(t) = T^0g\left(\frac{T(t)}{T^0}\right)$, $\mathcal{L}_2(t) = \int_0^\infty \phi(a)e(a,t) da$, $\mathcal{L}_3(t) =\int_0^\infty \psi(b)i(b,t) db$, and $\mathcal{L}_4(t) = \frac{\beta T^0}{c}V(t)$. Here the nonnegative kernel functions $\phi(a)$ and $\psi(b)$ will be determined later. Firstly, we calculate the derivative of $\mathcal{L}_i$, $i = 1,2,3,4$, respectively,

$\begin{align*} \frac{ d\mathcal{L}_1(t)}{ dt} = - dT_0\bigg(\frac{T^0}{T}+\frac{T}{T_0}-2\bigg)-\beta TV+\beta T^0V. \end{align*}$

By integration by parts, we calculate the derivative of $\mathcal{L}_2$,

$\begin{align*} \frac{ d\mathcal{L}_2(t)}{ dt} = &\ \int_0^\infty \phi(a)\frac{\partial e(a,t)}{\partial t} da =-\int_0^\infty \phi(a) \bigg[\theta_1(a)e(a,t)+\frac{\partial e(a,t)}{\partial a}\bigg] da \\ = &\ -\phi(a)e(a,t) \bigg|_0^\infty+\int_0^\infty \phi'(a)e(a,t) da-\int_0^\infty \phi(a)\theta_1(a)e(a,t) da \\ = &\ \phi(0)e(0,t)+\int_0^\infty \bigg(\phi'(a)-\phi(a)\theta_1(a)\bigg)e(a,t) da. \end{align*}$

An argument similar to the one used in calculating the derivative of $\mathcal{L}_2$, we get

$\frac{ d\mathcal{L}_3(t)}{ dt} = \psi(0)i(0,t)+\int_0^\infty \bigg(\psi'(b)-\psi(b)\theta_2(b)\bigg)i(b,t) db.$

We calculate the derivative of $\mathcal{L}_4$,

$\frac{ d\mathcal{L}_4(t)}{ dt} =\frac{\beta T^0}{c}\int_0^\infty p(b)i(b,t) db-\beta T^0V.$

Secondly, we have

$\begin{array}{l} \frac{{d{{\cal L}_{IFE}}(t)}}{{dt}} = \; - d{T^0}(\frac{{{T^0}}}{T} + \frac{T}{{{T_0}}} - 2) - \beta TV + \phi (0)f\beta TV + \psi (0)(1 - f)\beta TV\\ \quad \quad \quad \quad \;\; + \int_0^\infty ( \phi '(a) - \phi (a){\theta _1}(a) + \psi (0)\xi (a))e(a,t)da\\ \quad \quad \quad \quad \;\; + \int_0^\infty ( \psi '(b) - \psi (b){\theta _2}(b) + \frac{{\beta {T^0}}}{c}p(b))i(b,t)db. \end{array}$

Choosing

$\left\{ \begin{array}{lll} \displaystyle \psi(b) = \int_b^\infty \frac{\beta T^0}{c}p(u) e^{-\int_b^u \theta_2(\omega) d\omega} du,\\ \displaystyle \phi(a) = \int_a^\infty \psi(0)\xi(u) e^{-\int_a^u \theta_1(\omega) d\omega} du.\\ \end{array}\right.$

Then it is easy to see that

$\left\{ \begin{array}{lll} \displaystyle \psi(0) = \frac {\beta T^0J}{c},\ \phi(0) = \frac {\beta T^0JK}{c},\\ \displaystyle \psi'(b)-\psi(b)\theta_2(b)+\frac{\beta T^0}{c}p(b) = 0,\\ \displaystyle \phi'(a)-\phi(a)\theta_1(a)+\psi(0)\xi(a) = 0. \end{array}\right.$

Consequently, $\mathcal{L}_{IFE}$ satisfies

$\frac{ d\mathcal{L}_{IFE}(t)}{ dt} =- dT_0\bigg(\frac{T^0}{T}+\frac{T}{T_0}-2\bigg)+\big(\Re_0-1\big)\beta TV.$

Notice that $\frac{ d\mathcal{L}_{IFE}(t)}{ dt} =0$ implies that $T =T^0$. It can be verified that the largest invariant set where $\frac{ d\mathcal{L}_{IFE}(t)}{ dt} =0$ is the singleton $\{P^0\}$. Therefore, by the invariance principle, $P^0$ is globally attractive when $\Re_0\leq 1$.

To establish the global stability of the infection equilibrium, we introduce the following Lemma.

Lemma 6.2. Suppose that $\Re_0>1$. Then, for any solution $(T(t), e(a,t), i(b,t), V(t))$ of (1), the following equalities hold,

$\begin {eqnarray} \label{L1} (1-f)\beta T^*V^*\bigg[1-\frac{e(0, t)i^*(0)}{e^*(0) i(0, t)}\bigg]+\int_0^\infty \xi(a) e^*(a)\bigg[1-\frac{e(a, t)i^*(0)}{e^*(a) i(0, t)}\bigg] d a & =& 0, \end{eqnarray}$

(39)

Proof. We give the proof for (39). In fact,

$\begin{align*} &\ (1-f)\beta T^*V^*+\int_0^\infty \xi(a) e^*(a)da\\ &\ -(1-f)\beta T^*V^*\frac{e(0,t)i^*(0)}{e^*(0) i(0,t)}-\int_0^\infty \xi(a) e^*(a)\frac{e(a,t)i^*(0)}{e^*(a) i(0,t)} d a \end{align*}$

$\begin{align*} = &\ i^*(0)- \left((1-f)\beta TV+\int_0^\infty \xi(a) e(a,t)da\right)\frac{i^*(0)}{i(0,t)}\\ = &\ 0 \end{align*}$

This immediately gives (39).

Theorem 6.3. If $\Re_0>1$, then the infection equilibrium $P^* = (T^*,e^*(a),i^*(a),V^*)$ of (1) is globally attractive.

Proof. Let

$G[x,y] = x-y-y \ln \frac{x}{y},\ \mbox{for}\ x, y>0.$

It is easy to see that $G$ is non-negative on $(0,\infty)\times (0,\infty)$ with the minimum value $0$ only when $x = y$. Furthermore, it is easy to verify that $xG_x[x,y]+yG_y[x,y] = G[x,y]$.

Considering the following candidate Lyapunov function,

$\mathcal{L}_{EE}(t) = \mathcal{H}_1(t)+\mathcal{H}_2(t)+\mathcal{H}_3(t)+\mathcal{H}_4(t),$

where

$\begin{align*} &\ \mathcal{H}_1(t) =G[T, T^*],\ \mathcal{H}_2(t) =\int_0^\infty \phi_1(a)G \big[ e(a,t),e^*(a)\big]da,\\ &\ \mathcal{H}_3(t) =\int_0^\infty \psi_1(b) G \big[i(b,t), i^*(b) \big] db,\ \mathcal{H}_4(t) =\frac{\beta T^*}{c}G [V,V^*]. \end{align*}$

We define $\phi_1(a)$ and $\psi_1 (b)$ as

$\psi_1 (b) = \int_b^\infty \frac{\beta T^*}{c}p(u) e^{-\int_b^u \theta_2(\omega) d\omega} du,$

and

$\phi_1 (a) = \int_a^\infty \psi_1(0)\xi(u) e^{-\int_a^u \theta_1(\omega) d\omega} du,$

it follows that $\psi_1(0) = \frac{\beta T^*J}{c}$, $\phi_1(0) = \frac{\beta T^*KJ}{c}$ and

$\begin{align*} \psi_1'(b)-\psi_1(b)\theta_2(b) =&\ -\frac{\beta T^*}{c}p(b). \end{align*}$

$\begin{align*} \phi_1'(a)-\phi_1(a)\theta_1(a) =&\ -\psi_1(0)\xi(a). \end{align*}$

Firstly, we calculate the derivative of $\mathcal{H}_i$, $i = 1,2,3,4$, respectively,

$\frac{ d\mathcal{H}_1(t)}{ dt} =-dT^*\bigg(\frac{T}{T^*}+\frac{T^*}{T}-2\bigg)+\frac{1}{f}\left(1-\frac{T^*}{T}\right)\left(e^*(0)-e(0,t)\right)$

By using (4),

$\begin{align*} \mathcal{H}_2(t) =&\ \int_0^t \phi_1(a) G[e(0,t-a)\Omega(a), e^*(a)] da\\ &\ +\int_t^\infty \phi_1(a) G[e_0(a-t)e^{-\int_{a-t}^a \theta_1(\omega) d\omega}, e^*(a)] da\\ = &\ \int_0^t \phi_1(t-r) G[e(0,r)\Omega(t-r), e^*(t-r)] dr\\ &\ +\int_0^\infty \phi_1(t+r) G[e_0(r)e^{-\int_r^{t+r} \theta_1(\omega) d\omega}, e^*(t+r)] dr\\ = &\ \mathcal{B_1}(t)+\mathcal{B_2}(t). \end{align*}$

The derivative of $\mathcal{B_1}$ and $\mathcal{B_2}$ take the following form,

$\begin{align*} &\ \frac{ d\mathcal{B_1}(t)}{ dt}\\& =\phi_1(0) G[e(0,t),e^*(0)] +\int_0^t \phi_1'(t-r) G\left[e(0,r)e^{-\int_0^{t-r} \theta_1(\omega) d\omega},e^*(t-r)\right] dr\\ &\ -\int_0^t \phi_1(t-r) \theta_1(t-r)\bigg[e(0,r)e^{-\int_0^{t-r} \theta_1(\omega) d\omega} G_x\left[e(0,r)e^{-\int_0^{t-r} \theta_1(\omega) d\omega},e^*(t-r)\right]\\ &\ + e^*(t-r) G_y\left[e(0,r)e^{-\int_0^{t-r} \theta_1(\omega) d\omega},e^*(t-r)\right] \bigg] dr, \end{align*}$

and

$\begin{align*} &\ \frac{ d\mathcal{B_2}(t)}{ dt} =\int_0^\infty \phi_1'(t+r) G\left[e_0(r)e^{-\int_r^{t+r} \theta_1(\omega) d\omega},e^*(t+r)\right] dr\\ &\ -\int_0^\infty \phi_1(t+r) \theta_1(t+r) \bigg[e_0(r)e^{-\int_r^{t+r} \theta_1(\omega) d\omega} G_x\left[e_0(r)e^{-\int_r^{t+r} \theta_1(\omega) d\omega}, e^*(t+r)\right]\\ &\ +e^*(t+r) G_y\left[e_0(r)e^{-\int_r^{t+r} \theta_1(\omega) d\omega}, e^*(t+r)\right]\bigg] dr. \end{align*}$

We obtain the derivative of $\mathcal{H}_2(t)$,

$\begin{array}{l} \frac{{d{{\cal H}_2}(t)}}{{dt}} = \;{\phi _1}(0)G[e(0,t),{e^*}(0)] + \int_0^\infty [ {\phi _{1'}}(a) - {\phi _1}(a){\theta _1}(a)]G[e(a,t),{e^*}(a)]da\\ \quad \quad \quad \; = {\phi _1}(0)G[e(0,t),{e^*}(0)] - \int_0^\infty {{\psi _1}} (0)\xi (a)G[e(a,t),{e^*}(a)]da. \end{array}$

A similar argument as in the derivative of $\mathcal{H}_2$, we calculate the derivative of $\mathcal{H}_3$,

$\begin{align*} \frac{ d\mathcal{H}_3(t)}{ dt} =&\ \psi_1(0) G[i(0,t),i^*(0)]+\int_0^\infty \big[\psi'(b)-\psi(b)\theta_2(b) \big] G[i(b,t),i^*(b)] db \\ =&\ \psi_1(0) G[i(0,t),i^*(0)]-\int_0^\infty \frac{\beta T^*}{c}p(b) G[i(b,t),i^*(b)] db. \end{align*}$

We calculate the derivative of $\mathcal{H}_4$,

$\frac{{d{{\cal H}_4}(t)}}{{dt}} = \frac{{\beta {T^*}}}{c}\int_0^\infty p (b)i(b,t)db - \beta {T^*}V + \beta {T^*}{V^*} - \frac{{\beta {T^*}{V^*}}}{{cV}}\int_0^\infty p (b)i(b,t)db.$

If follows from $\psi_1(0) = \frac{\beta T^*J}{c}$ and $\phi_1(0) = \frac{\beta T^*KJ}{c}$ that

$\begin{array}{*{20}{l}} {\frac{{d{{\cal L}_{EE}}}}{{dt}} = }&{\; - d{T^*}(\frac{T}{{{T^*}}} + \frac{{{T^*}}}{T} - 2) + \frac{1}{f}\left( {1 - \frac{{{T^*}}}{T}} \right)\left( {{e^*}(0) - e(0,t)} \right)}\\ {}&{ + {\phi _1}(0)G[e(0,t),{e^*}(0)] - \int_0^\infty {{\psi _1}} (0)\xi (a)G[e(a,t),{e^*}(a)]da}\\ {}&{ + {\psi _1}(0)G[i(0,t),{i^*}(0)] - \int_0^\infty {\frac{{\beta {T^*}}}{c}} p(b)G[i(b,t),{i^*}(b)]db}\\ {}&{ + \int_0^\infty {\frac{{\beta {T^*}}}{c}} p(b)i(b,t)db + \beta {T^*}{V^*} - \beta {T^*}V - \frac{{{V^*}}}{V}\int_0^\infty {\frac{{\beta {T^*}}}{c}} p(b)i(b,t)db.{\rm{ }}} \end{array}$

(40)

Recall that

$(1-f)(\beta T^*V^*-\beta TV)+\int_0^\infty \xi(a)\bigg(e^*(a)-e(a,t)\bigg) da = i^*(0)-i(0,t),$

and

$\frac{f\beta T^* K J}{c}+\frac{(1-f)\beta T^* J}{c} = \left(\frac{f\beta K J}{c}+\frac{(1-f)\beta J}{c}\right)\frac{T^0}{\Re_0} = 1.$

Thus (40) becomes

$\begin{array}{l} \frac{{d{{\cal L}_{EE}}(t)}}{{dt}} = \; - d{T^*}(\frac{T}{{{T^*}}} + \frac{{{T^*}}}{T} - 2) + \frac{1}{f}\left( {1 - \frac{{{T^*}}}{T}} \right)\left( {{e^*}(0) - e(0,t)} \right)\\ \quad \quad \quad \quad + \frac{1}{f}G[e(0,t),{e^*}(0)] - \int_0^\infty {\frac{{\beta {T^*}}}{c}} p(b)G[i(b,t),{i^*}(b)]db\\ \quad \quad \quad \quad + \frac{{\beta {T^*}J}}{c}\left[ {(1 - f)\beta {T^*}{V^*}\ln \frac{{e(0,t){i^*}(0)}}{{{e^*}(0)i(0,t)}} + \int_0^\infty \xi (a){e^*}(a)\ln \frac{{e(a,t){i^*}(0)}}{{{e^*}(a)i(0,t)}}da} \right]\;\\ \quad \quad \quad \quad + \int_0^\infty {\frac{{\beta {T^*}}}{c}} p(b)i(b,t)db + \beta {T^*}{V^*} - \beta {T^*}V - \frac{{{V^*}}}{V}\int_0^\infty {\frac{{\beta {T^*}}}{c}} p(b)i(b,t)db. \end{array}$

It follows that,

$\begin{array}{*{20}{l}} {\frac{{d{{\cal L}_{EE}}(t)}}{{dt}}}&{ = \; - d{T^*}(\frac{T}{{{T^*}}} + \frac{{{T^*}}}{T} - 2) - \frac{1}{f}{e^*}(0)\left( {\frac{{{T^*}}}{T} - \ln \frac{{e(0,t)}}{{{e^*}(0)}}} \right)}\\ {}&{ - \int_0^\infty {\frac{{\beta {T^*}}}{c}} p(b)G[i(b,t),{i^*}(b)]db}\\ {}&{ + \frac{{\beta {T^*}J}}{c}\left[ {(1 - f)\beta {T^*}{V^*}\ln \frac{{e(0,t){i^*}(0)}}{{{e^*}(0)i(0,t)}} + \int_0^\infty \xi (a){e^*}(a)\ln \frac{{e(a,t){i^*}(0)}}{{{e^*}(a)i(0,t)}}da} \right]}\\ {}&{ + \int_0^\infty {\frac{{\beta {T^*}}}{c}} p(b)i(b,t)db + \beta {T^*}{V^*} - \frac{{{V^*}}}{V}\int_0^\infty {\frac{{\beta {T^*}}}{c}} p(b)i(b,t)db.} \end{array}$

(41)

Recall that $e^*(0) = f\beta T^*V^*$ and $\int_0^\infty p(b)i^*(b)db = cV^*$ in (11). Collecting the terms of (41) yields

$\begin{align*} &\ \frac{d\mathcal{L}_{EE}(t)}{dt} =-dT^*\bigg(\frac{T}{T^*}+\frac{T^*}{T}-2\bigg)\\ &\ +\frac{\beta T^* J}{c}\left[(1-f)\beta T^*V^*\ln\frac{e(0,t)i^*(0)}{e^*(0) i(0,t)}+\int_0^\infty \xi(a) e^*(a)\ln\frac{e(a,t)i^*(0)}{e^*(a) i(0,t)} d a\right]\\ &\ +\int_0^\infty\frac{\beta T^*}{c} p(b)i^*(b)\bigg(2+\ln\frac{i(b,t)}{i^*(b)}-\frac{T^*}{T}-\ln \frac{e(0,t)}{e^*(0)}-\frac{V^* i(b,t)}{V i^*(b)}\bigg) db. \end{align*}$

Further, we have

$\begin{array}{*{20}{l}} {\frac{{d{{\cal L}_{EE}}(t)}}{{dt}} = }&{\; - d{T^*}(\frac{T}{{{T^*}}} + \frac{{{T^*}}}{T} - 2)}\\ {}&{ + \frac{{\beta {T^*}J}}{c}(1 - f)\beta {T^*}{V^*}\left( {1 - \frac{{e(0,t){i^*}(0)}}{{{e^*}(0)i(0,t)}} + \ln \frac{{e(0,t){i^*}(0)}}{{{e^*}(0)i(0,t)}}} \right)}\\ {}&{ + \frac{{\beta {T^*}J}}{c}\int_0^\infty \xi (a){e^*}(a)\left( {1 - \frac{{e(a,t){i^*}(0)}}{{{e^*}(a)i(0,t)}} + \ln \frac{{e(a,t){i^*}(0)}}{{{e^*}(a)i(0,t)}}} \right)da}\\ {}&{ + \int_0^\infty {\frac{{\beta {T^*}p(b)}}{c}} {i^*}(b)(2 - \frac{{{T^*}}}{T} + \ln \frac{{{T^*}}}{T} - \frac{{{V^*}i(b,t)}}{{V{i^*}(b)}} + \ln \frac{{{V^*}i(b,t)}}{{V{i^*}(b)}})db} \end{array}$

(42)

$\begin{align*} &\ -\frac{\beta T^* J}{c}\bigg\{(1-f)\beta T^*V^*\bigg[1-\frac{e(0,t)i^*(0)}{e^*(0) i(0,t)}\bigg]\nonumber\\ &\ +\int_0^\infty \xi(a) e^*(a)\bigg[1-\frac{e(a,t)i^*(0)}{e^*(a) i(0,t)}\bigg] d a \bigg\}.\nonumber \end{align*}$

Recall that Lemma 6.2 holds. Putting (39) into (42), we have

$\begin{align*} \frac{d\mathcal{L}_{EE}(t)}{dt} =&\ -dT^*\bigg(\frac{T}{T^*}+\frac{T^*}{T}-2\bigg)\\ &\ -\frac{\beta T^*J}{c}\bigg[\int_0^\infty \xi(a) e^*(a)g\bigg(\frac{e(a,t)i^*(0)}{e^*(a) i(0,t)}\bigg)da\\ &\ +(1-f)\int_0^\infty\frac{\beta T^* }{c} p(b)i^*(b)g\bigg(\frac{e(0,t)i^*(0)}{e^*(0) i(0,t)}\bigg)db\bigg]\\ &\ -\int_0^\infty\frac{\beta T^* }{c} p(b)i^*(b)\bigg[g\bigg(\frac{T^*}{T}\bigg) +g\bigg(\frac{V^* i(b,t)}{Vi^*(b)}\bigg)\bigg] db\\ \leq &\ 0 \end{align*}$

and $\frac{ d\mathcal{L}_{EE}(t)}{ dt} =0$ implies that $T =T^*$ and

$\frac{i(b,t)}{i^*(b)}=\frac{i(0,t)}{i^*(0)}=\frac{V}{V^*}=\frac{e(0,t)}{e^*(0)}=\frac{e(a,t)}{e^*(a)}, \ \ \mbox{for all $a$, $b\geq 0$.}$

It is not difficult to check that the largest invariant subset $\{\frac{d\mathcal{L}_{EE}(t)}{ dt} =0\}$ is the singleton $\{P^*\}$. By the invariance principle, $P^*$ is globally attractive and hence the proof is complete.

---

7. Discussion

This paper is devoted to the global dynamics of an HIV infection model subject to latency age and infection age, where the model formulation, basic reproduction number computation, and rigorous mathematical analysis, such as relative compactness and persistence of the solution semi-flows, and existence of a global attractor are involved. We have shown that the existence of a compact attractor of all compact sets of nonnegative initial data and used the Lyapunov functional to show that this attractor is the singleton set containing the equilibrium. Given that the model is so complex, the proof does require some rigorous calculation. The dynamics (at least the long-term dynamics) of the model do not appear to have been altered by adding the $e(a,t)$ component. We hope the model studied here have a contribution to improve the broader contexts of investigating viral infection subject to age structure.

---

Acknowledgments

The authors would like to thank the anonymous referees and editor for very helpful suggestions and comments which led to improvements of our original manuscript. J. Wang is supported by National Natural Science Foundation of China (No. 11226255 and No. 11201128), the Science and Technology Innovation Team in Higher Education Institutions of Heilongjiang Province (No. 2014TD005). X. Dong is supported by Graduate Students Innovation Research Program of Heilongjiang University (No. YJSCX2017-177HLJU).

---