1. Introduction

In recent years, virus dynamics attracts more and more attentions of researchers and plays a crucial role in many diseases research, including AIDS, hepatitis and influenza. Many mathematical models have provided insights into virus infection and dynamics, as well as on how an infection can be managed, reduced or even eradicated ([3], [4], [7], [15], [17], [27], [38], [43], [44]). Since the basic three-dimensional viral infection model was proposed by Nowak et al. [21], Perelson et al. [26], Perelson and Nelson [25], Nowak and May [20], many people have established different within-host infection model, which help us to better understand virus infection and various drug therapy strategies by mathematical analysis, numerical simulations and clinical data ([13], [19], [22], [28], [29]). Note that immune responses play a critical part in the process of viral infections. Concretely, cytotoxic T lymphocyte (CTL) cells can attack infected cells, and antibody cells can neutralize viruses. To better understand the role of the immune function during virus infection, Wodarz proposed the following model with both CTL and antibody immune responses [41],

$\label{model1} \left\{ \begin{array}{l} \dot{T}(t)=\lambda-d_{1}T(t)-\beta T(t)V(t), \\ \dot{I}(t)=\beta T(t)V(t)- d_{2}I(t)-pI(t)C(t), \\ \dot{V}(t)=rd_{2}I(t)-d_{3}V(t)-qA(t)V(t), \\ \dot{C}(t)=k_{1}I(t)C(t)-d_{4}C(t), \\ \dot{A}(t)=k_{2}A(t)V(t)-d_{5}A(t), \end{array} \right.$

(1)

where a dot denotes the differentiation with respect to time $t$, $T(t), I(t), V(t), C(t)$ and $A(t)$ are the concentrations of healthy cells, infected cells, free virus, CTL cells and antibody cells at time $t$, respectively. $d_{1}, d_2, d_3, d_4, d_5$ are the death rates of healthy cell, infected cell, free virus, CTL cells and antibody cells, respectively. $\lambda$ represents a constant production of the healthy cells. The term $\beta TV$ represents the rate for the heathy $T$ cells to be infected by virus. Furthermore, infected cells are killed by CTL cells at a rate $pIC$. Free virus are produced by infected cells at a rate $rd_{2}I$, and are killed by antibody responses at a rate $qAV$. CTL immune responses are activated at a rate proportional to the abundance of CTL cells and infected cells, $k_{1}IC$. Antibody responses are produced at a rate proportional to the abundance of antibodies and free virus, $k_{2}AV$. Biologically, all parameters are positive.

After that, some researchers have taken into account the effect of immune responses including CTL responses or antibody responses ([24], [35], [36], [37], [39]). Some other researchers have incorporated the effect of CTL responses and intracellular delays ([11], [16], [18], [32], [45]). Concretely, the global dynamics of (1) with and without intracellular time delay is given in [24] and [42], respectively. Note that model (1) assumes that CTL and antibody responses are produced at bilinear rates. However, De Boer [5] pointed out that the bilinear rates cannot model several immune responses that are together controlling a chronic infection. In [5], De Boer has proposed an immune response function with the saturation. Incorporating the saturation effects of immune responses and the delay, [12] also obtained the global stability of the model, which is totally determined by the corresponding reproductive numbers. These results preclude the complicated behaviors such as the backward bifurcations and Hopf bifurcations which may be induced by saturation factors and time delay.

Note also that most of models assume CTL responses are activated by infected cells/antigenic stimulation, and antibody responses are activated by virus in these studies. However, as pointed out by Nowak and May [20], CTL responses have another function of self-regulating, i.e., the CTL responses are triggered by encountering foreign antigen and then adopts a constant level which is independent of the concentration of virions or infected cell. Bocharvor et al. have provided evidence the export of precursor CTL cells from the thymus [2]. Pang and Cui et al. have studied the export of specific precursor CTL cells from the thymus in [23], but they didn't considered intracellular time delay and antibody responses. Similarly, Wang and Wang have considered that neutralizing antibodies are produced at a constant rate after the injection [37], but they didn't take into account the effect of CTL responses and intracellular time delay.

Motivated by the above studies, we will formulate and analyze a virus dynamics model with the recruitment of immune responses, saturation effects of immune responses and an intracellular time delay, which can be described by the following functional differential equations:

$\label{model} \left\{ \begin{array}{l} \dot{T}(t)=\lambda-d_{1}T(t)-\beta T(t)V(t), \\ \dot{I}(t)=\beta T(t-\tau)V(t-\tau)e^{-s\tau}- d_{2}I(t)-pI(t)C(t), \\ \dot{V}(t)=rd_{2}I(t)-d_{3}V(t)-qA(t)V(t), \\ \dot{C}(t)=\lambda_1+\frac{k_{1}I(t)C(t)}{h_{1}+C(t)}-d_{4}C(t), \\ \dot{A}(t)=\lambda_2+\frac{k_{2}A(t)V(t)}{h_{2}+A(t)}-d_{5}A(t). \end{array} \right.$

(2)

Here, we use $\lambda_1$ to describe the export of specific precursor CTL cells from the thymus, $\lambda_2$ to describe the recruitment rate of antibody responses. The production of infected cells is delayed in such a way that the production of new virus at time $t$ depends on the population of virus and infected cells at a previous time $t-\tau$, and only a fraction of $e^{-s\tau}$ can survive after the interval $\tau$, where $1/s$ is the average lifetime of infected cells without reproduction. $h_1, h_2$ are saturation constants. Other parameters are same as that in model (1).

The main aim of the present paper is to explore the effects of the recruitment of immune responses on virus infection. The organization of this paper is as follows. In the next section, some preliminary analyzes of the model (2) will be given. Stability of all equilibria are given in Section 3. In Section 4, some numerical simulations are given to explain the effects of $\lambda_1, \lambda_2$, i.e., the term of recruitment of immune responses. Lastly, some brief conclusions are given in Section 5.

---

2. Preliminary analyses of the model

In this section, we will first prove the positivity and boundedness of solutions, and then derive the expression of the basic reproduction number for model (2).

---

2.1. Positivity and boundedness of solutions

Let ${X}:=C([-\tau, 0], R^{5})$ be the Banach space of continuous functions from $[-\tau, 0]$ to $R^{5}$. For $\phi=(\phi_1, \phi_2, \phi_3, \phi_4,$ $\phi_5) \in {X}$, define $\parallel\phi\parallel=\sum_{i=1}^{5} \parallel\phi_i\parallel_\infty$, in which $\parallel\phi_i\parallel_\infty=\max\limits_{\theta\in [-\tau, 0]} |\phi_i(\theta)|$. The initial functions for model (2) are provided with $\phi\in X^+=C([-\tau, 0], R^{5}_+)$.

Proposition 1. Under the above initial conditions, all solutions of model (2) are nonnegative. In particular, the solution $(T(t), I(t), V(t), C(t), A(t))$ of model (2) is positive for $t>0$ in its existence interval if $T(0)>0, I(0)>0, V(0)>0, C(0)>0, A(0)>0$.

Proof. We first verify that $T(t)$ is positive in the existence interval of the solution. Suppose not. Then there is $t_{1}>0$ such that $T(t_{1})=0$ and $T(t)>0, t\in[0, t_{1})$. Indeed, if $T(0)=0$, we have $\dot{T}(0)=\lambda>0$. Thus, $T(t)>0$ for small positive $t$. Evidently, this remains valid if $T(0)>0$. As a result, the existence of $t_1$ follows if the claim is not true. Note that $\dot{T}(t_{1})=\lambda>0$. Thus, there is a sufficiently small $\varepsilon>0$, such that $T(t) < 0$ for all $t\in (t_{1}-\varepsilon, t_{1})$. We get a contradiction because $T(t)>0, t\in[0, t_{1})$. Hence, $T(t)$ is positive in the existence interval of the solution. In the same way, we obtain $C(t)$ and $A(t)$ are positive in the existence interval of the solution. Then, we verify that $I(t)$ and $V(t)$ are positive. In the same way, we assume that there is a first time $t_{2}$ such that $V(t_{2})=0$, from the third equation of model (2) we have

$\dot{V}(t_{2})=rd_{2}I(t_{2}).$

By solving the second equation of model (2), we obtain

$I(t_2)=e^{\int_{0}^{t_2}-(d_{2}+pC(\xi))d\xi}[I(0)+\int_{0}^{t_2}\beta T(\theta-\tau)V(\theta-\tau)e^{-s\tau}e^{\int_{0}^{\theta}(d_{2}+pC(\xi))d\xi}d\theta]>0.$

It follows that $\dot{V}(t_{2})>0$, hence $V(t)>0$. Furthermore,

$I(t)=e^{\int_{0}^{t}-(d_{2}+pC(\xi))d\xi}[I(0)+\int_{0}^{t}\beta T(\theta-\tau)V(\theta-\tau)e^{-s\tau}e^{\int_{0}^{\theta}(d_{2}+pC(\xi))d\xi}d\theta].$

From the above expression of $I(t)$ solution, we get $I(t)>0$.

It follows easily that $I(t)\ge 0, V(t)\ge 0$ in the existence interval of the solution if the initial functions are in $X^+$, and $I(t)>0, V(t)>0, C(t)>0, A(t)>0$ in the existence interval of the solution if $I(0)>0, V(0)>0, C(0)>0, A(0)>0$.

Proposition 2. All solutions of model (2) in $X^+$ are ultimately bounded.

Proof. Set

$L(t)=T(t)+I(t+\tau)+\frac{1}{3r}V(t+\tau)+\frac{d_2}{3k_1}C(t+\tau) +\frac{d_3}{4k_{2}r}A(t+\tau).$

Calculating the derivative of $L$ along the solution of (2), we get

$\begin{align} &\dot{L}(t)=\lambda -{{d}_{1}}T(t)-\beta T(t)V(t)+\beta T(t)V(t)-{{d}_{2}}I(t+\tau )-pI(t+\tau )C(t+\tau ) \\ &+\frac{{{d}_{2}}}{3}I(t+\tau )-\frac{{{d}_{3}}}{3r}V(t+\tau )-\frac{q}{3r}A(t+\tau )V(t+\tau )+\frac{{{d}_{2}}}{3{{k}_{1}}}{{\lambda }_{1}} \\ &+\frac{{{d}_{2}}}{3}\frac{I(t+\tau )C(t+\tau )}{{{h}_{1}}+C(t+\tau )}-\frac{{{d}_{2}}}{3{{k}_{1}}}{{d}_{4}}C(t+\tau )+\frac{{{d}_{3}}}{4{{k}_{2}}r}{{\lambda }_{2}} \\ &+\frac{{{d}_{3}}}{4r}\frac{A(t+\tau )V(t+\tau )}{{{h}_{2}}+A(t+\tau )}-\frac{{{d}_{3}}}{4{{k}_{2}}r}{{d}_{5}}A(t+\tau ). \\ \end{align}$

(3)

Since

$\frac{C(t+\tau)}{h_1+C(t+\tau)}\leq1, \;\; \frac{A(t+\tau)}{h_2+A(t+\tau)}\leq1,$

we obtain

$\begin{align} &\dot{L}(t)\le \lambda -{{d}_{1}}T(t)-{{d}_{2}}I(t+\tau )+\frac{{{d}_{2}}}{3}I(t+\tau )-\frac{{{d}_{3}}}{3r}V(t+\tau )+\frac{{{d}_{2}}}{3{{k}_{1}}}{{\lambda }_{1}}+\frac{{{d}_{2}}}{3}I(t+\tau ) \\ &-\frac{{{d}_{2}}}{3{{k}_{1}}}{{d}_{4}}C(t+\tau )+\frac{{{d}_{3}}}{4{{k}_{2}}r}{{\lambda }_{2}}+\frac{{{d}_{3}}}{4r}V(t+\tau )-\frac{{{d}_{3}}}{4{{k}_{2}}r}{{d}_{5}}A(t+\tau ) \\ &\le \lambda +\frac{{{d}_{2}}}{3{{k}_{1}}}{{\lambda }_{1}}+\frac{{{d}_{3}}}{4{{k}_{2}}r}{{\lambda }_{2}}-{{d}_{1}}T(t)-\frac{{{d}_{2}}}{3}I(t+\tau )-\frac{{{d}_{3}}}{4}\frac{1}{3r}V(t+\tau ) \\ &-{{d}_{4}}\frac{{{d}_{2}}}{3{{k}_{1}}}C(t+\tau )-{{d}_{5}}\frac{{{d}_{3}}}{4{{k}_{2}}r}A(t+\tau ) \\ &\le \lambda +\frac{{{d}_{2}}}{3{{k}_{1}}}{{\lambda }_{1}}+\frac{{{d}_{3}}}{4{{k}_{2}}r}{{\lambda }_{2}}-mL(t), \\ \end{align}$

where $m=\min{\{d_{1}, d_{2}/3, d_{3}/4, d_{4}, d_{5}\}}$. It follows that the nonnegative solutions of (2) exist on $[0, \infty)$ and are ultimately bounded. Moreover,

$\label{Lupper2} \limsup\limits_{t\longrightarrow\infty}L(t)\leq\frac{\lambda}{m}+\frac{d_2\lambda_1}{3k{_1}m}+\frac{d_3\lambda_2}{4k_{2}rm}.$

From the first equation of model (2), we get

$\dot{T}(t)\leq\lambda-d_{1}T(t).$

It follows that

$\limsup\limits_{t\longrightarrow\infty}T(t)\leq\frac{\lambda}{d_{1}}.$

Set $F(t)=T(t)+I(t+\tau)$, then

$\dot{F}(t)=\dot{T}(t)+\dot{I}(t+\tau)\leq \lambda-nF(t), n=\min{\{d_1, d_2\}},$

thus

$\limsup\limits_{t\longrightarrow\infty}F(t)\leq\frac{\lambda}{n}.$

Then,

$\label{eq1} \limsup\limits_{t\longrightarrow\infty}(T(t)+I(t+\tau))\leq\frac{\lambda}{n}.$

(4)

From the third equation of model (2) and (4), we have

$\limsup\limits_{t\longrightarrow\infty}V(t)\leq\frac{rd_{2}\lambda}{d_{3}n}.$

Further, let

$M=\max{\{\frac{\lambda}{d_{1}}, \frac{\lambda}{n}, \frac{rd_{2}\lambda}{d_{3}n}, \frac{\lambda}{m}+\frac{d_2\lambda_1}{3k{_1}m}+\frac{d_3\lambda_2}{4k_{2}rm}\}}.$

The dynamics of model (2) can be analyzed in the following bounded feasible region

$\begin{align} &\Gamma =\{(T,I,V,C,A)\mid 0\le T\le M,0\le T+I\le M,0\le V\le M, \\ &0\le C\le M,0\le A\le M\}. \\ \end{align}$

---

2.2. The basic reproductive number

Based on the concept of the basic reproductive number for an epidemic disease presented in [6, 35], we know the basic reproductive number $R_0$ of virus is the expected number of viruses that one virion gives rise to in an totally uninfected cell population during its lifetime.

From model (2), it is clear that healthy cells, CTL cells and antibody cells will stabilize to $\lambda/d_1, \lambda_1/d_4$ and $\lambda_2/d_5$ if there is not infection, i.e., $I(t)=V(t)=0$. In this case, suppose that one virion is introduced, it can produce a maximum amount of $P_1(\varphi_1) = \beta\frac{\lambda}{d_1} e^{-s\tau}\varphi_{1}$ infections during its mean lifetime of $\varphi_1=1/ (d_3+q\frac{\lambda_2}{d_5})$. In addition, an infected cell has an average lifetime of $\varphi_2=1/(d_2+p\frac{\lambda_1}{d_4})$, and hence, it can averagely generate $P_2(\varphi_2) = r d_2 \varphi_2$ virus. Therefore, the basic reproduction number of virus for model (2) can be defined as

$R_0 = P_1(\varphi_1) P_2(\varphi_2) = \beta\frac{\lambda}{d_1}\frac{rd_2}{(d_2+p\frac{\lambda_1}{d_4})} \frac{1}{(d_3+q\frac{\lambda_2}{d_5})}e^{-s\tau}.$

Based on the above expression, we know that there are inverse proportional relationship between the basic reproduction number of virus ($R_0$) and the recruitment rate of immune responses ($\lambda_1$ and $\lambda_2$). Thus, $R_0$ will decreases along with $\lambda_1, \lambda_2$ increasing, which means that ignoring the effects of recruitment rate of immune responses will overestimate the basic reproduction number of virus.

---

3. Stability of the equilibria

In this section, we first discuss the existence of infection-free equilibrium, and then analyze its stability. Besides, using the uniform persistence theory, we obtain the existence of an endemic equilibrium. After that, the stability of an endemic equilibrium was proved by constructing Lyapunov functional.

---

3.1. Infection-free equilibrium

Apparently, there is always an infection-free equilibrium in system (2): $E_{0}=(T_0, 0, 0, C_0, A_0)$, where

$T_0=\frac{\lambda}{d_{1}}, \; C_0=\frac{\lambda_1}{d_4}, \; A_0=\frac{\lambda_2}{d_5}.$

Next, we discuss the stability of the infection-free equilibrium $E_0$.

Theorem 3.1. When $R_{0} < 1$, the infection-free equilibrium $E_{0}$ is globally asymptotically stable in region $\Gamma$.

Proof. First we define a Lyapunov functional $L_{0}$ by

$\begin{align*} L_0=&\int_{T_0}^{T(t)}\frac{(S-T_0)}{S}dS+e^{s\tau}I(t)+(\frac{1}{r}+\frac{pC_0}{rd_2})e^{s\tau}V(t) \\ &+\frac{p e^{s\tau}}{k_{1}}\int_{C_0}^{C(t)}\frac{(h_1+S)(S-C_0)}{S}dS\\ &+\frac{(d_2+pC_0)} {rd_2}\frac{q e^{s\tau}}{k_2}\int_{A_0}^{A(t)}\frac{(h_2+S)(S-A_0)}{S}dS\\ &+\int_{-\tau}^{0}\beta T(t+\theta)V(t+\theta)d\theta. \end{align*}$

Calculating the time derivative of $L_{0}$ along the solution of system (2), we obtain

$\begin{align*} \dot{L_{0}}=&\lambda-d_{1}T(t)-\beta T(t)V(t)-\frac{T_{0}}{T(t)}(\lambda-d_{1}T(t)-\beta T(t)V(t))\\ &+\beta T(t-\tau)V(t-\tau)-d_{2}I(t)e^{s\tau}-pI(t)C(t)e^{s\tau}+\frac{1}{r}(rd_{2}I(t)e^{s\tau}\\ &-d_{3}V(t)e^{s\tau}-qA(t)V(t)e^{s\tau})+\frac{pC_0}{rd_2}(rd_{2}I(t)e^{s\tau}-d_{3}V(t)e^{s\tau}\\ &-qA(t)V(t)e^{s\tau}) +\frac{p e^{s\tau}}{k_{1}}(h_1+C(t))\{\lambda_1-d_{4}C(t)+\frac{k_{1}I(t)C(t)}{h_{1}+C(t)}\\ & -\frac{C_0}{C(t)}(\lambda_1+\frac{k_{1}I(t)C(t)}{h_{1}+C(t)}-d_{4}C(t))\} +\frac{(d_2+pC_0)}{rd_2}\frac{q e^{s\tau}}{k_2}(h_2\\ &+A(t))\{\lambda_2 +\frac{k_{2}A(t)V(t)}{h_{2}+A(t)}-d_{5}A(t) -\frac{A_0}{A(t)}(\lambda_2+\frac{k_{2}A(t)V(t)}{h_{2}+A(t)}\\ &-d_{5}A(t))\}+\beta T(t)V(t)-\beta T(t-\tau)V(t-\tau). \end{align*}$

Since $\lambda=d_{1}T_{0}, \lambda_1=d_{4}C_0, \lambda_2=d_{5}A_0$, it follows that

$\begin{align} &\overset{.}{\mathop{{{L}_{0}}}}\,=2{{d}_{1}}{{T}_{0}}-{{d}_{1}}T(t)-\frac{{{T}_{0}}}{T(t)}{{d}_{1}}{{T}_{0}}+\beta {{T}_{0}}V(t)-pI(t)C(t){{e}^{s\tau }} \\ &-\frac{{{d}_{3}}}{r}V(t){{e}^{s\tau }}-\frac{q}{r}A(t)V(t){{e}^{s\tau }}+pI(t){{C}_{0}}{{e}^{s\tau }}-\frac{p{{C}_{0}}}{r{{d}_{2}}}{{d}_{3}}V(t){{e}^{s\tau }} \\ &-\frac{p{{C}_{0}}}{r{{d}_{2}}}qA(t)V(t){{e}^{s\tau }}+\frac{p{{e}^{s\tau }}}{{{k}_{1}}}{{\lambda }_{1}}({{h}_{1}}+C(t))+pI(t)C(t){{e}^{s\tau }}-\frac{p{{e}^{s\tau }}}{{{k}_{1}}}{{d}_{4}}C(t)({{h}_{1}}+C(t)) \\ &-\frac{p{{e}^{s\tau }}}{{{k}_{1}}}{{\lambda }_{1}}({{h}_{1}}+C(t))\frac{{{C}_{0}}}{C(t)}-pI(t){{C}_{0}}{{e}^{s\tau }}+\frac{p{{e}^{s\tau }}}{{{k}_{1}}}{{d}_{4}}{{C}_{0}}({{h}_{1}}+C(t)) \\ &+\frac{({{d}_{2}}+p{{C}_{0}})}{r{{d}_{2}}}\frac{q{{e}^{s\tau }}}{{{k}_{2}}}{{\lambda }_{2}}({{h}_{2}}+A(t))+\frac{({{d}_{2}}+p{{C}_{0}})}{r{{d}_{2}}}{{e}^{s\tau }}qA(t)V(t) \\ &-\frac{({{d}_{2}}+p{{C}_{0}})}{r{{d}_{2}}}\frac{q{{e}^{s\tau }}}{{{k}_{2}}}{{d}_{5}}A(t)({{h}_{2}}+A(t))-\frac{({{d}_{2}}+p{{C}_{0}})}{r{{d}_{2}}}\frac{q{{e}^{s\tau }}}{{{k}_{2}}}{{\lambda }_{2}}\frac{{{A}_{0}}}{A(t)}({{h}_{2}} \\ &+A(t))-\frac{({{d}_{2}}+p{{C}_{0}}){{e}^{s\tau }}}{r{{d}_{2}}}q{{A}_{0}}V(t)+\frac{({{d}_{2}}+p{{C}_{0}})}{r{{d}_{2}}}\frac{q}{{{k}_{2}}}{{e}^{s\tau }}{{d}_{5}}{{A}_{0}}({{h}_{2}}+A(t)) \\ &={{d}_{1}}{{T}_{0}}(2-\frac{T(t)}{{{T}_{0}}}-\frac{{{T}_{0}}}{T(t)})+\frac{({{d}_{2}}+p{{C}_{0}})({{d}_{3}}+q{{A}_{0}}){{e}^{s\tau }}}{r{{d}_{2}}}({{R}_{0}}-1)V(t) \\ &-\frac{p{{d}_{4}}{{e}^{s\tau }}}{{{k}_{1}}}{{(C(t)-{{C}_{0}})}^{2}}+\frac{p}{{{k}_{1}}}{{\lambda }_{1}}{{h}_{1}}{{e}^{s\tau }}(2-\frac{{{C}_{0}}}{C(t)}-\frac{C(t)}{{{C}_{0}}}) \\ \end{align}$

$-\frac{({{d}_{2}}+p{{C}_{0}})q{{e}^{s\tau }}}{r{{d}_{2}}{{k}_{2}}}{{d}_{5}}{{(A(t)-{{A}_{0}})}^{2}}+\frac{({{d}_{2}}+p{{C}_{0}})q{{e}^{s\tau }}}{r{{d}_{2}}{{k}_{2}}}{{\lambda }_{2}}{{h}_{2}}(2-\frac{{{A}_{0}}}{A(t)}-\frac{A(t)}{{{A}_{0}}}).$

Since the geometric mean is less than or equal to the arithmetical mean, it follows from $R_{0} < 1$ that $\dot{L_{0}}\leq0$. Let

$D_{0}=\{(T(t), I(t), V(t), C(t), A(t))|\dot{L_{0}}=0\}.$

It is easy to show that $E_{0}=(T_{0}, 0, 0, C_0, A_0)$ is the largest invariant set in $D_{0}$. By the Lyapunov-LaSalle invariance principle [8], $E_{0}$ is globally asymptotically stable.

---

3.2. Uniform persistence

In order to obtain the the existence of an endemic equilibrium, in this subsection, we investigate the uniform persistence of (2). We first introduce a preliminary theory. Let $X$ be a metric space and $\Phi$ be a semiflow on $X$. Suppose that $X^0$ is an open set in $X$, $X^0\subset X, X^0\cap X_0=\emptyset$, and $X^0\cup X_0=X$. Define $M_\partial=\{x\in X_0: \Phi_t(x)\in X_0, t\geq 0\}$, which may be empty. A continuous $\mathfrak{p}: X\rightarrow[0, \infty)$ satisfying condition: $\mathfrak{p}(\Phi_t(x))>0$ for $t>0$ if either $\mathfrak{p}(x)=0$ and $x\in X^0$ or if $\mathfrak{p}(x)>0$, will be called a generalized distance function for $\Phi$.

Lemma 3.2. ([31], Theorem 3) Let $\mathfrak{p}$ be a generalized distance function for semiflow $\Phi$. Assume that

(H1) $\Phi$ has a global attractor $\widetilde{A}$.

(H2) There exists a finite sequence $\widetilde{M}=\{M_1, ..., M_k\}$ of pairwise disjoint, compact and isolated sets in $X_0$ with the following properties:

(ⅰ) $\bigcup_{x\in \widetilde{M}_{\partial}}\omega(x) \subset \bigcup_{i=1}^{k} M_i$,

(ⅱ) no subset of $\widetilde{M}$ forms a cycle in $X_0$,

(ⅲ) $M_i$ is isolated in $X$,

(ⅳ) $W^{s}(M_i) \bigcap \mathfrak{p}^{-1}(0, \infty)=\emptyset$, $i=1, ...k$.

Then there exists $\delta>0$ such that for any compact chain transitive set $L$ with $L \not \subset M_i$ for $i=1, ..k$, there holds $\min_{x\in L} \mathfrak{p}(x)>\delta$.

By applying Lemma 3.2 to (2), we can obtain the following result for the uniform persistence of (2).

Theorem 3.3. If $R_0>1$, then system (2) is uniformly persistent, i.e., there exists $\varepsilon>0$ (independent of initial conditions), such that, $\liminf\limits_{t\rightarrow+\infty}T(t)\geq \varepsilon, \liminf\limits_{t\rightarrow+\infty} I(t)\geq \varepsilon, \liminf\limits_{t\rightarrow+\infty} V(t)\geq \varepsilon, \liminf\limits_{t\rightarrow+\infty} C(t)\geq \varepsilon, \liminf\limits_{t\rightarrow+\infty } A(t)\geq \varepsilon$, for all solutions of (2) with initial condition.

Proof. Let

$\begin{align} &{{X}_{0}}=\{\tilde{\phi }\in {{X}^{+}}:{{{\tilde{\phi }}}_{2}}(\theta )\equiv 0,{{{\tilde{\phi }}}_{3}}(\theta )\equiv 0\ \ for\ \ \theta \in [-\tau,0]\}, \\ &{{X}^{0}}={{X}^{+}}\setminus {{X}_{0}},{{M}_{\partial }}=\{\psi \in {{X}^{+}}:{{\Phi }_{t}}(\psi )\in {{X}_{0}},t\ge 0\}. \\ \end{align}$

Basic analysis of (2) implies that $X_0$ is a positive invariant set for (2). The positive invariance of $X^0$ follows from Proposition 1. For any initial value condition $\phi_0 \in X^+$, define $\Phi_{t}(\phi_0)$ for $t\geq 0$ as $\Phi_{t}(\phi_0):=(T_{t}(\theta), I_{t}(\theta), V_{t}(\theta) , C_{t}(\theta), A_{t}(\theta))$ for $\theta \in[-\tau, 0]$, where $(T_{t}(\theta), I_{t}(\theta), V_{t}(\theta), C_{t}(\theta), A_{t}(\theta))$ is the solution of (2) with initial condition $\phi_0$. By Proposition 1 and Proposition 1 we have $\Phi_{t}(\phi_0)$ is dissipative in $X^+$, and hence by Arzel$\grave{a}$-Ascoli theorem, condition (H1) of Lemma 3.2 is satisfied.

Let $\omega(\psi)$ be the omega limit set of the orbit $\Phi(t)$ through $\psi \in X^+$. Note that system (2) has an unique boundary equilibrium $E_0=(\lambda / d_1, 0, 0, \lambda_1 /d_4, \lambda_2/d_5).$ For any ${\psi}\in M_\partial$, i.e., $\Phi_{t}({\psi})\in X_0$, we have $I_{t}({\psi})\equiv 0, V_{t}({\psi})\equiv 0$ for all $t\geq0$ and

$\label{equC} \left\{ \begin{array}{l} \dot{T}(t)=\lambda-d_{1}T(t), \\ \dot{C}(t)=\lambda_1-d_4C(t), \\ \dot{A}(t)=\lambda_2-d_5A(t). \end{array} \right.$

(5)

It then follows from the result in [14] that $(T_{t}(\psi), C_{t}(\psi), A_{t}(\psi))\rightarrow(T_0, C_0, A_0)$ as $t\rightarrow+\infty$. Thus, $\bigcup_{\psi\in M_{\partial} }\omega(\psi)= E_0$. Furthermore, by Theorem 3.1, $E_0$ is unstable if $R_0>1$. Then $\{E_0\}$ is isolated and acyclic covering and the conditions (ⅰ), (ⅱ) and (ⅲ) of Lemma 3.2 are satisfied.

Since

$R_0=\beta \frac{\lambda}{d_1}\frac{ rd_2}{d_2+p\frac{\lambda_1}{d_4}}\frac{1}{d_3+q\frac{\lambda_2}{d_5}}e^{-s\tau}>1,$

we have

$\label{eq13} (d_2+p\frac{\lambda_1}{d_4})(d_3+q\frac{\lambda_2}{d_5}) < \beta \frac{\lambda}{d_1}rd_2e^{-s\tau}.$

(6)

Thus, there is sufficiently small $\sigma$ such that

$(d_2+p(\frac{\lambda_1}{d_4}+\sigma))(d_3+q(\frac{\lambda_2}{d_5}+\sigma)) < \beta (\frac{\lambda}{d_1}-\sigma)rd_2e^{-s\tau}.$

Suppose $W^{s}(E_0)\cap X^0\neq \emptyset$. There exists a positive solution $(T^*(t), I^*(t), V^*(t),$ $C^*(t), A^*(t))$ such that

$(T^*(t), I^*(t), V^*(t), C^*(t), A^*(t))\rightarrow (\lambda/d_1, 0, 0, \lambda_1/d_4, \lambda_2/d_5)\;\ as \;\ t\rightarrow+\infty.$

For sufficiently large $t^*_1$, when $t\geq t^*_1$, we have

$\begin{align} &\frac{\lambda }{{{d}_{1}}}-\sigma < {{T}^{*}}(t) < \frac{\lambda }{{{d}_{1}}}+\sigma, \\ &\frac{{{\lambda }_{1}}}{{{d}_{4}}}-\sigma < {{C}^{*}}(t) < \frac{{{\lambda }_{1}}}{{{d}_{4}}}+\sigma, \\ &\frac{{{\lambda }_{2}}}{{{d}_{5}}}-\sigma < {{A}^{*}}(t) < \frac{{{\lambda }_{2}}}{{{d}_{5}}}+\sigma, \\ \end{align}$

if $t>t^*_1+\tau$, it follows that

$\label{modelC} \left\{ \begin{array}{l} \dot{I}^{*}(t)\geq -d_{2}I^*(t)+\beta V^*(t)(\frac{\lambda}{d_1}-\sigma)e^{-s\tau}-pI^*(t)(\frac{\lambda_1}{d_4}+\sigma), \\ \dot{V}^{*}(t)\geq -d_{3}V^*(t)+r d_2 I^*(t)-q(\frac{\lambda_2 }{d_5}+\sigma)V^*(t).\\ \end{array} \right.$

(7)

Since

$\label{eq33} \ A_\sigma= \left( \begin{array}{cc} -d_2-p(\frac{\lambda_1}{d_4}+\sigma)&\beta (\frac{\lambda}{d_1}-\sigma)e^{-s\tau} \\ rd_2& -d_3-q(\frac{\lambda_2}{d_5}+\sigma) \end{array} \right),$

(8)

the non-diagonal elements of (8) are positive, and from (6), we obtain $|A_\sigma| < 0$. By Perron-Frobrniuss Theorem, we can obtain the maximum eigenvalue $\alpha>0$ of $A_\sigma$, and it has an eigenvector $u=(u_1, u_2), u_1>0, u_2>0$, then choose sufficiently small $l$ such that $I^*(t_0)> l u_1, V^*(t_0)> l u_2$.

Now consider the following auxiliary system

$\label{modelA1} \left\{ \begin{array}{l} \dot{I}^{*}(t)=-d_{2}I^*(t)+\beta V^*(t)(\frac{\lambda}{d_1}-\sigma)-pI^*(t)(\frac{\lambda_1}{d_4}+\sigma), \\ \dot{V}^{*}(t)=-d_{3}V^*(t)+r d_2 I^*(t)-q(\frac{\lambda_2 }{d_5}+\sigma)V^*(t).\\ \end{array} \right.$

(9)

Note $I(t), V(t)$ are solutions of (9) with initial condition $I(t_0)=lu_1, V(t_0)=lu_2$. Since the model (9) is monotone, and $A_{\sigma}u>0$, by [[30] Corollary $5.2.2$], we have $I(t)\rightarrow +\infty, V(t)\rightarrow +\infty$ as $t\rightarrow +\infty$. Using the comparison theorem, we have $I^*(t)\rightarrow +\infty, V^*(t)\rightarrow +\infty$ for $t\rightarrow +\infty$, which is a contradiction. Hence, $W^{s}(E_0)\cap X^0=\emptyset$.

Define a continuous function $\mathfrak{p}: X^+\rightarrow R^+$ by

$\mathfrak{p}(\overline{\phi})=\min \{\overline{\phi}_2(0), \overline{\phi}_3(0)\}, \forall \overline{\phi} \in X^+.$

It is clear that $\mathfrak{p}^{-1}(0, \infty)\subset X^0$ and if $\mathfrak{p}(\overline{\phi})>0$ then $\mathfrak{p}(\Phi_{t}(\overline{\phi}))>0$ for all $t>0$ (see e.g. [1]). Based on the above proof, and by Lemma 3.2, it then follows that there exists $\delta>0$ such that $\liminf \limits_{t\rightarrow\infty}\mathfrak{p}(\Phi_{t}(\overline{\phi}))\geq\delta$ for all $\overline{\phi}\in X^0$, which implies that there is a positive $\eta_0$ such that $\liminf \limits_{t\rightarrow\infty}(I(t), V(t))\geq(\eta_0, \eta_0)$ i.e., $\liminf \limits_{t\rightarrow\infty}I(t)\geq\eta_0$, $\liminf \limits_{t\rightarrow\infty}V(t)\geq\eta_0$.

Furthermore, from the first equation of (2), Proposition 1 and the above results, we have

$\dot{T}(t)=\lambda-d_1T(t)-\beta T(t)V(t)>\lambda-d_1T(t)-\beta MT(t)=\lambda-(d_1+\beta M)T(t),$

Thus,

$\liminf\limits_{t\rightarrow+\infty} T(t)> \frac{\lambda}{d_1+\beta M}.$

From the fourth equation of (2),

$\dot{C}(t)=\lambda_1+\frac{k_1 I(t)C(t)}{h_1+C(t)}-d_4C(t)\geq \lambda_1-d_4C(t),$

we have

$\liminf\limits_{t\rightarrow+\infty} C(t) \geq \frac{\lambda_1}{d_4}.$

From the fifth equation of (2),

$\dot{A}(t)=\lambda_2+\frac{k_2 A(t)V(t)}{h_1+A(t)}-d_5A(t)\geq \lambda_2-d_5A(t),$

Therefore, taking $\varepsilon=\min{\{\eta_0, \lambda_1/d_4, \lambda_2/d_5, \lambda/(d_1+\beta M)\}}$, we can conclude that

$\liminf\limits_{t\rightarrow+\infty } T(t)\geq \varepsilon, \liminf\limits_{t\rightarrow+\infty } I(t)\geq \varepsilon, \liminf\limits_{t\rightarrow+\infty } V(t)\geq \varepsilon, \liminf\limits_{t\rightarrow+\infty } C(t)\geq \varepsilon, \liminf\limits_{t\rightarrow+\infty } A(t)\geq \varepsilon$

are valid for any solution of system (2) with initial condition in $X^0$. This completes the proof.

From the Theorem 3.1, we are easy to get that $E_0$ is unstable if $R_0>1$, and by Proposition 1 and the uniformly persistent of system (2), we can obtain that if $R_{0}>1$, system (2) exists at least one endemic equilibrium $E_{1}=(T_{1}, I_{1}, V_{1}, C_1, A_1).$

---

3.3. The endemic equilibrium

Now, we discuss the stability of the endemic equilibrium $E_1$.

Theorem 3.4. When $R_{0}>1$, the endemic equilibrium $E_{1}$ is globally asymptotically stable in region $\Gamma$.

Proof. Set

$m_1=\frac{\beta T_{1}V_{1}}{rd_{2}I_{1}}.$

Define a Lyapunov functional $L_{1}$ by

$\begin{align} &{{L}_{1}}=\int_{{{T}_{1}}}^{T(t)}{\frac{(S-{{T}_{1}})}{S}}dS+{{e}^{s\tau }}\int_{{{I}_{1}}}^{I(t)}{\frac{(S-{{I}_{1}})}{S}}dS+{{m}_{1}}\int_{{{V}_{1}}}^{V(t)}{\frac{(S-{{V}_{1}})}{S}}dS \\ &+\frac{p{{e}^{s\tau }}}{{{k}_{1}}}\int_{{{C}_{1}}}^{C(t)}{\frac{({{h}_{1}}+C(t))(S-{{C}_{1}})}{S}}dS+\frac{{{m}_{1}}q}{{{k}_{2}}}\int_{{{A}_{1}}}^{A(t)}{\frac{({{h}_{2}}+A(t))(S-{{A}_{1}})}{S}}dS \\ &+\beta {{T}_{1}}{{V}_{1}}\int_{-\tau }^{0}{(\frac{T(t+\theta )V(t+\theta )}{{{T}_{1}}{{V}_{1}}}-1-\ln \frac{T(t+\theta )V(t+\theta )}{{{T}_{1}}{{V}_{1}}})}d\theta . \\ \end{align}$

Calculating the time derivative of $L_{1}$ along the solution of system (2), we obtain

$\begin{align} &\overset{.}{\mathop{{{L}_{1}}}}\,=\lambda -{{d}_{1}}T(t)-\beta T(t)V(t)-\frac{{{T}_{1}}}{T(t)}(\lambda -{{d}_{1}}T(t)-\beta T(t)V(t))-pI(t)C(t){{e}^{s\tau }} \\ &+\beta T(t-\tau )V(t-\tau )-{{d}_{2}}I(t){{e}^{s\tau }}-\frac{{{I}_{1}}}{I(t)}(\beta T(t-\tau )V(t-\tau )-{{d}_{2}}I(t){{e}^{s\tau }} \\ &-pI(t)C(t){{e}^{s\tau }})+{{m}_{1}}(r{{d}_{2}}I(t)-{{d}_{3}}V(t)-qA(t)V(t))-{{m}_{1}}\frac{{{V}_{1}}}{V(t)}(r{{d}_{2}}I(t) \\ &-{{d}_{3}}V(t)-qA(t)V(t))+\frac{p{{e}^{s\tau }}}{{{k}_{1}}}({{h}_{1}}+C(t))\{{{\lambda }_{1}}+\frac{{{k}_{1}}I(t)C(t)}{{{h}_{1}}+C(t)}-{{d}_{4}}C(t) \\ &-\frac{{{C}_{1}}}{C(t)}({{\lambda }_{1}}+\frac{{{k}_{1}}I(t)C(t)}{{{h}_{1}}+C(t)}-{{d}_{4}}C(t))\}+\frac{q{{m}_{1}}}{{{k}_{2}}}({{h}_{2}}+A(t))\{{{\lambda }_{2}}+\frac{{{k}_{2}}A(t)V(t)}{{{h}_{2}}+A(t)} \\ &-{{d}_{5}}A(t)-\frac{{{A}_{1}}}{A(t)}({{\lambda }_{2}}+\frac{{{k}_{2}}A(t)V(t)}{{{h}_{2}}+A(t)}-{{d}_{5}}A(t))\}+\beta T(t)V(t) \\ &-\beta T(t-\tau )V(t-\tau )+\beta {{T}_{1}}{{V}_{1}}\ln \frac{T(t-\tau )V(t-\tau )}{T(t)V(t)}. \\ \end{align}$

Since

$\begin{align} &\lambda ={{d}_{1}}{{T}_{1}}+\beta {{T}_{1}}{{V}_{1}},\ \beta {{T}_{1}}{{V}_{1}}=({{d}_{2}}{{I}_{1}}+p{{I}_{1}}{{C}_{1}}){{e}^{s\tau }},\ \ r{{d}_{2}}{{I}_{1}}={{d}_{3}}{{V}_{1}}+q{{A}_{1}}{{V}_{1}}, \\ &{{\lambda }_{1}}+\frac{{{k}_{1}}{{I}_{1}}{{C}_{1}}}{{{h}_{1}}+{{C}_{1}}}={{d}_{4}}{{C}_{1}},\ {{\lambda }_{2}}+\frac{{{k}_{2}}{{A}_{1}}{{V}_{1}}}{{{h}_{2}}+{{A}_{1}}}={{d}_{5}}{{A}_{1}}, \\ \end{align}$

we have

$\begin{align*} \dot{L_{1}}&=d_{1}T_{1}(2-\frac{T(t)}{T_{1}}-\frac{T_{1}}{T(t)})+\beta T_{1}V_{1}-\beta T_{1}V_{1}\frac{T_{1}}{T(t)}+\beta T_1 V(t)\\ &-\frac{I_{1}}{I(t)}\beta T(t-\tau)V(t-\tau)+d_2I_1e^{s\tau}+pI_1C(t)e^{s\tau}-m_1d_3V(t) \\ &-m_1qA(t)V(t)-\beta T_1V_1\frac{V_1I(t)}{V(t)I_1}+m_1d_3V_1+m_1qA(t)V_1 \\ &+\frac{pe^{s\tau}}{k_1}\lambda_1(h_1+C(t))-\frac{p e^{s\tau}}{k_1}d_4C(t)(h_1+C(t)) -\frac{p e^{s\tau}}{k_1}\frac{C_1}{C(t)} \lambda_1(h_1\\ &+C(t))+\frac{pe^{s\tau}}{k_1}d_4C_1(h_1+C(t))+m_1\frac{q}{k_2}\lambda_2(h_2+A(t)) +m_1qA(t)V(t) \\ &-m_1\frac{q}{k_2}d_5A(t)(h_2+A(t)) -m_1\frac{q}{k_2}\frac{A_1}{A(t)}\lambda_2(h_2+A(t))-m_1q A_1 V(t) \\ &+m_1\frac{q}{k_2}d_5A_1(h_2+A(t))+\beta T_1V_1\ln \frac{T(t-\tau)V(t-\tau)}{T(t)V(t)}\\ =&d_{1}T_{1}(2-\frac{T(t)}{T_{1}}-\frac{T_{1}}{T(t)}) +\beta T_1V_1(1-\frac{T_1}{T(t)}+\ln\frac{T_1}{T(t)})+\beta T_1V_1(1 \end{align*}$

$\begin{align} &-\frac{{{I}_{1}}T(t-\tau )V(t-\tau )}{I(t){{T}_{1}}{{V}_{1}}}+\ln \frac{{{I}_{1}}T(t-\tau )V(t-\tau )}{I(t){{T}_{1}}{{V}_{1}}})+\beta {{T}_{1}}{{V}_{1}}(1-\frac{{{V}_{1}}I(t)}{{{I}_{1}}V(t)} \\ &+\ln \frac{{{V}_{1}}I(t)}{{{I}_{1}}V(t)})+\frac{p}{{{k}_{1}}}{{\lambda }_{1}}{{h}_{1}}{{e}^{s\tau }}(2-\frac{C(t)}{{{C}_{1}}}-\frac{{{C}_{1}}}{C(t)})-\frac{p}{{{k}_{1}}}{{d}_{4}}{{e}^{s\tau }}{{(C(t)-{{C}_{1}})}^{2}} \\ &+{{m}_{1}}\frac{q}{{{k}_{2}}}{{\lambda }_{2}}{{h}_{2}}(2-\frac{{{A}_{1}}}{A(t)}-\frac{A(t)}{{{A}_{1}}})-{{m}_{1}}\frac{q}{{{k}_{2}}}{{d}_{5}}{{(A(t)-{{A}_{1}})}^{2}}. \\ \end{align}$

Since the geometric mean is less than or equal to the arithmetical mean and $1-x+\ln x\leq0$ for any $x>0$, it follows that $\dot{L_{1}}\leq0$. Let

$D_{1}=\{(T(t), I(t), V(t), C(t), A(t))|\dot{L_{1}}=0\}.$

It is easy to verify that $\dot{L_{1}}(t)=0$ if and only if

$\frac{T_{1}}{T(t)}=\frac{I_1T(t-\tau)V(t-\tau)}{I(t)T_1V_1}=\frac{V_1I(t)}{I_1V(t)}=1.$

Thus, $T(t)=T_{1}$ and

$\dot{T}(t)=\lambda-d_{1}T_{1}-\beta T_{1}V(t)=0.$

As a result, we have $V(t)=V_{1}$, and then $I(t)=I_{1}$. From the second equation and the third equation of model (2), we have

$\left\{ \begin{align} &\dot{I}(t)=\beta {{T}_{1}}{{V}_{1}}{{e}^{-s\tau }}-{{d}_{2}}{{I}_{1}}-p{{I}_{1}}{{C}_{1}}=0, \\ &\dot{V}(t)=r{{d}_{2}}{{I}_{1}}-{{d}_{3}}{{V}_{1}}-q{{A}_{1}}{{V}_{1}}=0, \\ \end{align} \right.$

which implies $C(t)=C_1, A(t)=A_1$. Therefore, the largest invariant set in $D_{1}$ is $E_{1}$. Thus, when $R_{0}>1$, all positive solutions converge to $E_{1}$ by the LaSalle invariance principle [8].

---

4. Numerical simulations

In this section, we implement numerical simulations to explore the effects of the recruitment of immune responses ($\lambda_1$ and $\lambda_2$) on the infected cells ($I_1$) and virus load ($V_1$) at the endemic equilibrium $E_{1}$. The parameter values are chosen from literatures ([2], [26], [33], [32], [41], [38], [40], [46]). Especially, according to [9], we choose the range of $\lambda_1, \lambda_2$ from 0 to 1.

The all parameter values are shown in Table 1.

Figure 1 illustrates that $I_1, V_1$ decrease along with $\lambda_1, \lambda_2$ increasing in which implies that ignoring the recruitment of immune responses will overestimate the severity of the infection.

---

5. Conclusions

In this paper, the global dynamics of a within-host model with immune responses and intracellular time delay has been studied. By the method of Lyapunov functional and persistence theory, we obtain the global stability of the model (2) are completely determined by the values of the reproductive number. The results imply that the complicated behaviors such as backward bifurcations and Hopf bifurcations do not exist in the model with both immune responses and time delay.

Considering the basic reproductive number of virus

$R_0=R(\tau)=\frac{\lambda \beta r d_2 e^{-s\tau}}{d_1 (d_2+p\frac{\lambda_1}{d_4})(d_3+q\frac{\lambda_2}{d_5})}$

as a function of $\tau$, we can find that it is decreasing in $\tau$ and it tends to $0$ if the time delay tends to $\infty$. Furthermore, comparing with the previous studies ([11], [12], [16], [23], [24], [39], [42], [45]), we find that the expression of $R_0$ for model (2) is different, i.e., it includes the parameters $\lambda_1, \lambda_2$ which reflect the recruitment of immune responses. This implies that ignoring the recruitment of immune responses will result in overestimation of the reproductive number. Numerical simulations also show that the part of $V_1, I_1$ at steady state will decrease along with $\lambda_1, \lambda_2$ increasing (see Fig. 1), which mean that the recruitment of immune responses play a significant role in eradication of diseases. To sum up, we can conclude that ignoring the recruitment of immune responses will overestimate the infection degree and the severity of disease. These may provide a new insight for developing antiviral drug therapy strategies, which is to increase the recruitment rate of immune responses.

---

Acknowledgments

The authors are very grateful to the anonymous referees for their valuable comments and suggestions. This research is supported by the National Natural Science Fund of P. R. China (No. 11271369).

---