\documentclass[reqno]{amsart} \usepackage{hyperref} \AtBeginDocument{{\noindent\small \emph{Electronic Journal of Differential Equations}, Vol. 2013 (2013), No. 94, pp. 1--11.\newline ISSN: 1072-6691. URL: http://ejde.math.txstate.edu or http://ejde.math.unt.edu \newline ftp ejde.math.txstate.edu} \thanks{\copyright 2013 Texas State University - San Marcos.} \vspace{9mm}} \begin{document} \title[\hfilneg EJDE-2013/94\hfil Global stability] {Global stability of a delay differential equation of hepatitis B virus infection with immune response} \author[J. Wang, X. Tian \hfil EJDE-2013/94\hfilneg] {Jinliang Wang, Xinxin Tian} \address{Jinliang Wang \newline School of Mathematical Science, Heilongjiang University, Harbin 150080, China} \email{jinliangwang@yahoo.cn} \address{Xinxin Tian \newline School of Mathematical Science, Heilongjiang University, Harbin 150080, China} \email{xinxintian@yahoo.cn} \thanks{Submitted March 4, 2013. Published April 11, 2013.} \subjclass[2000]{34K20, 92D25} \keywords{HBV infection model; delay; CTLs; global stability} \begin{abstract} The global stability for a delayed HBV infection model with CTL immune response is investigated. We show that the global dynamics is determined by two sharp thresholds, basic reproduction number $\Re_0$ and CTL immune-response reproduction number $\Re_1$. When $\Re_0 \leq 1$, the infection-free equilibrium is globally asymptotically stable, which means that the viruses are cleared and immune is not active; when $\Re_1 \leq 1 < \Re_0$, the CTL-inactivated infection equilibrium exists and is globally asymptotically stable, which means that CTLs immune response would not be activated and viral infection becomes chronic; and when $\Re_1 > 1$, the CTL-activated infection equilibrium exists and is globally asymptotically stable, in this case the infection causes a persistent CTLs immune response. Our model is formulated by incorporating a Cytotoxic T lymphocytes (CTLs) immune response to recent work [Gourley, Kuang, Nagy, J. Bio. Dyn., 2(2008), 140-153] to model the role in antiviral by attacking virus infected cells. Our analysis provides a quantitative understandings of HBV replication dynamics in vivo and has implications for the optimal timing of drug treatment and immunotherapy in chronic HBV infection. \end{abstract} \maketitle \numberwithin{equation}{section} \newtheorem{theorem}{Theorem}[section] \newtheorem{lemma}[theorem]{Lemma} \allowdisplaybreaks \section{Introduction} Approximately more than 350 million people worldwide live with chronic hepatitis B virus (HBV) infection\cite{WHO}, and 25-40 percent of these chronic infection carrier will at risk of developing chronic liver disease, cirrhosis and hepatocellular carcinoma \cite{Nowak2}. HBV infection therefore represents a significant global public health problem. A basic within-host viral infection model introduced by Nowak et al \cite{Nowak1, Nowak2} and Perelson et al \cite{Perelson3} have been widely used in the studies of HBV and HIV infection dynamics and its treatment with the reverse transcriptase inhibitor lamivudine. After then several mathematical models have been modified to study of anti-HBV infection treatment and its dynamics. Most of these models focus on cell-free viral spread in a compartment such as the bloodstream, see, for example, In \cite{Song}, saturated mass action incidence rates $\beta xv/(1 +\alpha v)$ was proposed under the assumption that a less than linear response in $v$ could occur due to saturation at high virus concentration. Min et al \cite{Min} and Zheng et al \cite{Zheng} employed a standard incidence function instead of the mass action incidence to describe the hepatitis B virus infection as follows \begin{equation}\label{M1} \begin{gathered} \dot{x}(t)= \lambda-dx(t)-\frac{\beta x(t)v(t)}{x(t)+y(t)}, \\ \dot{y}(t)= \frac{\beta x(t)v(t)}{x(t)+y(t)}-ay(t), \\ \dot{v}(t)= ky(t)-\mu y(t), \end{gathered} \end{equation} where $x$, $y$ and $v$ are numbers of uninfected (susceptible) liver cells, infected liver cells and free virions, respectively. Uninfected liver cells are assumed to be produced at a constant rate, $\lambda$, to maintain tissue homeostasis in the face of hepatocyte turnover, described by the linear term $dx$, where $d$ is the per-capita death rate. Infected liver cells are killed by immune cells at rate $ay$. Free virions are generated from infected cells at the rate of $ky$ and decay by lymphatic and other mechanisms at the rate of $\mu v$, where $k$ is the so-called ``burst'' constant. Upon infection with HIV-1, there is a short intracellular ``eclipse phase'' (often referred to as ``latency'' in the literature), during which the cell is infected but has not yet begun producing virus. There are two methods to model this eclipse phase, by a time delay or by an explicit class of latently infected cells. Recently, Gourley et al \cite{Gourley} proposed the following model (As an extension of this model \eqref{M1}) under some biologically motivated modifications: \begin{equation}\label{M2} \begin{gathered} \dot{x}(t)= \lambda-dx(t)-\frac{\beta kx(t)y(t)}{\mu(x(t)+y(t))}, \\ \dot{e}(t)= -de(t)+\frac{\beta kx(t)y(t)}{\mu(x(t)+y(t))} -\frac{\beta ke^{-d\tau}x(t-\tau)y(t-\tau)}{\mu(x(t-\tau)+y(t-\tau))}, \\ \dot{y}(t)= \frac{\beta e^{-d\tau} kx(t-\tau)y(t-\tau)} {\mu(x(t-\tau)+y(t-\tau))}-ay(t), \end{gathered} \end{equation} where $e(t)$ represents the number of exposed cells (i.e., cells that have acquired the virus but are not yet producing new virions). Exposed cells begin shedding virions after $\tau$ units of time, representing the time required to construct, transcribe and translate the episomal viral genome, construct and then release mature virions. Other parameters are the same as in the basic virus model \eqref{M1}. Model \eqref{M2} is obtained from the following three observations: (1) A typical chronically infected HBV patient has a total serum load of about $2\times 10^{11}$ to $3\times 10^{12}$ virions \cite{Nowak1}. The average human liver has about an equal number of cells (assuming a liver mass of about 1.5 kg). These large numbers suggest that a more plausible HBV model should employ a standard incidence function, instead of the mass action incidence. On the other hand, the time delay associated with virus production is on the order of a day or two \cite{Murray}, much shorter than the life expectancy of a typical hepatocyte, which is 6-12 months or more \cite{Seeger}. This makes $e$ much smaller than $x$ and $y$. Hence, $e$ can be omitted from the denominators of the infection term. (2) The HBV incubation period, which varies from 45 to 180 days, and the delay in viral shedding mentioned previously both suggest that viral production delay may significantly impact infection dynamics and, hence, should be explicitly modeled. (3) Variable $v$ (virus particles) evolves on much faster time scale than the liver cells do, so a quasi-steady state assumption can be applied to $v$; i.e., to a good approximation, $v = ky/\mu$. In fact, it has been reported (Dimitrov et al. \cite{Dimitrov} and Sato et al. \cite{Sato}) that cell-to-cell spread of virus is favored over infections with cell-free virus inocula. For example, HTLV-I infection in vivo is achieved through cell-to-cell contact among healthy and infected CD4+ T cells \cite{Bangham}. It is evidently that cell-to-cell infection is the predominant route of viral spread since viral replication in a system with rapid cell turnover kinetics depends on cell-to-cell transfer of virus(see e.g. Gummuluru et al. \cite{Gummuluru}, Haase et al \cite{Haase1, Haase2}, Philips et al \cite{Philips}). Then the above HBV infection model cab be termed as cell to cell infection model. Note that in above simpler model \eqref{M2}, the $x$ and $y$ equations do not involve variable $e$ and form a closed subsystem of two equations. Guo and Cai \cite{GuoCai} resolved the global stability of infection equalibrium of model \eqref{M2}, without other additional conditions, which is left as an open problem in \cite{Gourley}. They showed that the infected equilibrium of system \eqref{M2} is always globally asymptotically stable as long as it exists by constructing suitable Lyapunov functional and LaSalle invariance principle. In most viral infections, cytotoxic T lymphocyte cells (CTLs), which attack infected cells, and antibody cells, which attack viruses, play a key role in antiviral defense. Chronic HBV infection is often the result of exposure early in life, leading to viral persistence in the absence of strong antibody or cellular immune responses \cite{Ferrari}. Therapy of HBV carriers can aim to either inhibit viral replication or enhance immunological responses against the virus, or both \cite{Regenstein}. It is currently widely accepted that HBV infection is non-cytopathic. Infected hepatocytes are killed not by the virus but by HBV-specific cytotoxic T lymphocytes (CTLs) \cite{Guidotti, Murray}. This mortality is somehow amplified by inflammation responses within the liver, but CTLs appear to be the major mediator of hepatitis B pathogenesis \cite{Iannacone}. Therefore, one of the dynamics of viral infection model with CTLs response have recently drawn much attention of researchers in the related areas and the interaction between infected cells and CTLs response in vivo has been studied by ordinary differential equations (ODEs) or delay differential equations (DDEs) (see e.g.\cite{Arnaout, Culshaw, WWL}). In this article, letting $z(t)$ be the density of CTLs, we propose the model \begin{equation}\label{M3} \begin{gathered} \dot{x}(t)= \lambda-d_1x(t)-\frac{\beta x(t)y(t)}{x(t)+y(t)}, \\ \dot{y}(t)= \frac{\beta e^{-d_1\tau} x(t-\tau)y(t-\tau)}{x(t-\tau) +y(t-\tau)}-d_2y(t)-ay(t)z(t), \\ \dot{z}(t)= py(t)z(t)-d_3z(t), \end{gathered} \end{equation} where the infected cells $y(t)$ are removed at a rate $ayz$ by the CTL immune response and the virus-specific CTL cells proliferate at a rate $pyz$ by contact with the infected cells, and die at a rate $d_3z$. The aim of the present paper is to carry out a complete mathematical analysis of model \eqref{M3}, we will show that the global properties of model \eqref{M3} for $\Re_1\leq 1<\Re_0$ and $\Re_1>1$ without any further conditions on the parameters. More precisely, we show that if $\Re_1\leq 1<\Re_0$, CTL-inactivated infection equilibrium $E_1$ is globally asymptotically stable; if $\Re_1>1$, CTL-activated infection equilibrium $E_2$ is globally asymptotically stable. The global stabilities of these models are established by constructing Lyapunov functionals and Lyapunov-LaSalle invariance principle (see e.g. \cite{Haddock}). Similar methods and techniques had been engaged by motivated by the works by Huang et al \cite{Huang}, Korobeinikov \cite{Korobeinikov}, McCluskey \cite{McCluskey} and Wang et al \cite{Wang1, Wang2}. The organization of this paper is as follows. In Section 2, we give the preliminaries of model \eqref{M3} including basic reproduction number, CTL immune-response reproduction number and equilibria. In Section 3, The global stability results is proved by Lyapunov functionals. A brief discussion is given in Section 4 to conclude this work. \section{Preliminaries} We denote by $\mathcal{C}$ the Banach space of continuous real-valued functions $\mathcal{C} = \mathcal{C}([-\tau, 0], \mathbb{R}^3)$ with the sup-norm \begin{equation} \|\varphi\|= \max\Big\{\sup_{-\tau\leq\theta\leq 0}|\varphi_1(\theta)|, \sup_{-\tau\leq\theta\leq 0}|\varphi_2(\theta)|, \sup_{-\tau\leq\theta\leq 0}|\varphi_3(\theta)|\Big\} \end{equation} for $\varphi=(\varphi_1, \varphi_2, \varphi_3)\in \mathcal{C}$. Further, the nonnegative cone of $\mathcal{C}$ is defined as $\mathcal{C}^+ = \mathcal{C}([-\tau, 0], \mathbb{R}^3_+)$. The initial conditions of system \eqref{M3} at $t=0$ are given as $x(\theta) =\varphi_1(\theta)$, $y(\theta) = \varphi_2(\theta)$, $z(\theta) = \varphi_3(\theta)$, $\theta\in [-\tau, 0]$. where \begin{equation}\label{Ic} \varphi=(\varphi_1, \varphi_2, \varphi_3)\in \mathcal{C}^+, \quad \varphi(0)>0. \end{equation} The following theorem establishes the positivity and boundedness of solutions for system \eqref{M3} with initial conditions \eqref{Ic}. \begin{theorem} Under the preceding initial conditions \eqref{Ic}, then $x(t), y(t)$ and $z(t)$ are all nonnegative and bounded for all $t$ at which the solution exists. \end{theorem} \begin{proof} By the existence and uniqueness theorem \cite[Theorem 2.1]{Kuang} of delay differential equations, there exists a $t_0 > 0$ such that there exists a solution $(x(t), y(t), z(t))$ of system \eqref{M3} for $0 0$ be the first time such that $x(t_1) = 0$. If $x(t)$ ere to lose its non-negativity, there would have to be $x' (t_1)\leq 0$, by the first equation of system \eqref{M1}, this is clearly impossible given the equation for $x(t)$ in system \eqref{M3}. It follows that $x(t) > 0$ for $t > 0$ as long as $x(t)$ exists. By the second equation of system \eqref{M3}, we have \begin{align*}\label{E11} y(t)&=\ y(0)\exp\Big(-d_2t-a\int_0^tz(\theta)d\theta\Big)\\ &\quad +\int_0^t\frac{\beta e^{-d\tau}x(\theta-\tau)y(\theta-\tau)} {x(\theta-\tau)+y(\theta-\tau)}e^{d_2(\theta-t)} \exp\Big(-a\int_{\theta}^tz(\sigma)d\sigma\Big)d\theta. \end{align*} It follows that $y(t) > 0$ for $t > 0$. From the third equation of system \eqref{M3}, we have $z(t)=z(0)\exp[(py-d_3)t]$. This shows that $z(t)\geq 0$ for $ 0 \leq t < t_1$. Next we show that positive solutions of \eqref{M3} are ultimately uniformly bounded for $t \geq 0$. Let $$ G(t)=e^{-d_1\tau}x(t)+y(t+\tau)+\frac{a}{p}z(t). $$ Adding all the equations of \eqref{M3} we obtain \begin{align*} G'(t)&= \lambda e^{-d_1\tau}-d_1e^{-d_1\tau}x(t) -d_2y(t+\tau)-\frac{d_3}{a}z(t) \\ &\leq \lambda e^{-d_1\tau}-dG(t), \end{align*} where $d=\min\{d_1, d_2, d_3\}$. Then $G(t)\leq M_1$ for some $M_1 > 0$ for sufficiently large $t$. For example, we can take as $M_1 = \frac{2\lambda e^{-d\tau}}{d}$, which implies that $G(t)$ is ultimately bounded, and so are $x(t), y(t)$ and $z(t)$. This proof is complete. \end{proof} System \eqref{M3} always exists an infection-free equilibrium $E_0=(x_0,0,0)$, where $x_0=\frac{\lambda}{d_1}$, which represents the state that the viruses are absent. The basic reproduction number of system \eqref{M3} is given by $$ \Re_0=\frac{\beta e^{-d\tau}}{d_2}. $$ If $\Re_0\leq 1$, an infection-free equilibrium $E_0$ is the unique equilibrium, corresponding to the extinction of free viruses. If $\Re_0 > 1$, in addition to $E_0$, there exists an CTL-inactivated infection equilibrium $E_1(x_1, y_1,0)$, where $$ x_1=\frac{\lambda}{d_1+d_2e^{d_1\tau}(\Re_0-1)},\quad y_1=\frac{\lambda (\Re_0-1)}{d_1+d_2e^{d_1\tau}(\Re_0-1)}, $$ which represents the state that the viruses are present whereas the CTLs are absent. We introduce a CTL immune-response reproduction number $$ \Re_1=\frac{d_1+d_2 e^{d_1\tau}(\Re_0-1)}{p\lambda}(py_1-d_3)+1. $$ Given $\Re_1>1$, then system \eqref{M3} has an CTL-activated infection equilibrium $E_2(x_2, y_2, z_2)$, where \begin{gather*} x_2=\frac{(\lambda p-d_1d_3-\beta d_3) +\sqrt{(\lambda p-d_1d_3-\beta d_3)^2+4d_1d_3\lambda p}}{2d_1p}, \\ y_2=\frac{d_3}{p},\quad z_2=\frac{\beta pe^{-d_1\tau}x_2}{ax_2p+d_3}-\frac{d_2}{a}. \end{gather*} Clearly, the endemic equilibrium represents the state that both the viruses and CTL response are present. \section{Main results} Throughout the article, we let $g(x)=x-1-\ln x$, to simplify many of the expressions which follow. Note that $g: \mathbb{R}_{+}\rightarrow \mathbb{R}_{+}$ has strict global minimum $g(1)=0$. \begin{theorem}\label{T1} If $\Re_0 \leq 1$, then the disease free equilibrium $E_0$ is globally asymptotically stable. \end{theorem} \begin{proof} Define a Lyapunov functional \begin{equation} L_0=e^{-d_1\tau}\Big[x-x_0-\int_{x_0}^{x} \frac{x_0(\theta+y)}{\theta(x_0+y)}d\theta\Big]+y +\frac{a}{p}z+\beta e^{-d_1\tau}\int_{t-\tau}^t\frac{x(\theta) y(\theta)}{x(\theta)+y(\theta)}d\theta. \end{equation} Calculating the time derivative of $L_0$ along the solution of \eqref{M3}, it follows that \begin{align*} &\frac{dL_0}{dt}\Big|_{\eqref{M3}}\\ &= e^{-d_1\tau}\big[1-\frac{x_0(x+y)}{x(x_0+y)}\big]\dot{x}+\dot{y}+\frac{a}{p}\dot{z}+\frac{\beta e^{-d_1\tau}xy}{x+y}-\frac{\beta e^{-d_1\tau}x(t-\tau)y(t-\tau)}{x(t-\tau)+y(t-\tau)}\\ &= e^{-d_1\tau}\big[1-\frac{x_0(x+y)}{x(x_0+y)}\big] \big[\lambda-d_1x-\frac{\beta xy}{x+y}\big]+\frac{\beta e^{-d_1\tau} x(t-\tau)y(t-\tau)}{x(t-\tau)+y(t-\tau)}\\ &\quad -d_2y-ayz+\frac{a}{p}(pyz-d_3z) +\frac{\beta e^{-d_1\tau}xy}{x+y} -\frac{\beta e^{-d_1\tau}x(t-\tau)y(t-\tau)}{x(t-\tau)+y(t-\tau)}\\ &= e^{-d_1\tau}d_1(x_0-x)\big[1-\frac{x_0(x+y)}{x(x_0+y)}\big] +\frac{\beta e^{-d_1\tau}x_0y}{x_0+y}-d_2y-\frac{ad_3z}{p}\\ &= e^{-d_1\tau}d_1(x_0-x)\big[1-\frac{x_0(x+y)}{x(x_0+y)}\big] +\frac{d_2x_0y(\Re_0-1)}{x_0+y}-\frac{d_2y^2}{x_0+y}-\frac{ad_3z}{p}, \end{align*} where $$ e^{-d_1\tau}d_1(x_0-x)\big[1-\frac{x_0(x+y)}{x(x_0+y)}\big] =-\frac{e^{-d_1\tau}d_1y(x_0-x)^2}{x(x_0+y)}\leq 0. $$ Therefore, $\Re_0\leq 1$ ensures that $dL_0/dt \leq 0$ for all $x >0$, $y \geq 0$, $z \geq 0$, and $dL_0/dt = 0$ holds if and only if $x = x_0, y = 0$, and $z(t) = 0$ for $\Re_0 \leq 1$. If follows that the largest invariant set in $\{(x_t, y_t, v_t, z_t)|dL_0/dt = 0\}$ is $E_0$. The classical Lyapunov-LaSalle invariance principle \cite[Theorem 2.5.3]{Kuang} shows that $E_0$ is globally asymptotically stable when $\Re_0\leq 1$. \end{proof} \begin{theorem} \label{T2} If $\Re_1\leq 1<\Re_0$, then CTL-inactivated infection equilibrium $E_1$ is globally asymptotically stable. \end{theorem} \begin{proof} Define a Lyapunov functional \begin{align*} L_1&= x-x_1-\int_{x_1}^{x}\frac{x_1(\theta+y_1)}{\theta(x_1+y_1)}d\theta +e^{d_1\tau}y_1g\big(\frac{y}{y_1}\big)+\frac{ae^{d_1\tau}}{p}z\\ &\quad + e^{d_1\tau}d_2y_1 \int_{t-\tau}^tg \Big(\frac{\beta x(\theta)y(\theta)}{e^{d_1\tau}d_2y_1(x(\theta)+y(\theta))}\Big) d\theta. \end{align*} Calculating the time derivative of $L_1$ along the solution of \eqref{M3}, it follows that \begin{equation} \label{E1} \begin{aligned} \frac{dL_1}{dt}\Big|_{\eqref{M1}} &= \big[1-\frac{x_1(x+y_1)}{x(x_1+y_1)}\big]\dot{x} +e^{d_1\tau}\big(1-\frac{y_1}{y}\big)\dot{y}+\frac{ae^{d_1\tau}}{p}\dot{z} \\ &\quad +\frac{\beta xy}{x+y}-\frac{\beta x(t-\tau)y(t-\tau)}{x(t-\tau) +y(t-\tau)}-e^{d_1\tau}d_2y_1\ln\frac{\beta xy}{e^{d_1\tau}d_2y_1(x+y)} \\ &\quad + e^{d_1\tau}d_2y_1\ln\frac{\beta x(t-\tau)y(t-\tau)}{e^{d_1\tau} d_2y_1(x(t-\tau)+y(t-\tau))} \\ &= \big[1-\frac{x_1(x+y_1)}{x(x_1+y_1)}\big] \Big[d_1(x_1-x)+\Big(\frac{\beta x_1y_1}{x_1+y_1}-\frac{\beta xy}{x+y}\Big)\Big] \\ &\quad + \big[\frac{\beta x(t-\tau)y(t-\tau)}{x(t-\tau)+y(t-\tau)} -e^{d_1\tau} (d_2y(t)-ay(t)z(t))\big]\big(1-\frac{y_1}{y}\big) \\ &\quad +\frac{ae^{d_1\tau}}{p}(py-d_3)z +\frac{\beta xy}{x+y} -\frac{\beta x(t-\tau)y(t-\tau)}{x(t-\tau)+y(t-\tau)} \\ &\quad -e^{d_1\tau}d_2y_1\ln\frac{\beta xy}{e^{d_1\tau}d_2y_1(x+y)} \\ &\quad + e^{d_1\tau}d_2y_1\ln\frac{\beta x(t-\tau)y(t-\tau)}{e^{d_1 \tau}d_2y_1(x(t-\tau)+y(t-\tau))}. \end{aligned} \end{equation} Here we used that \begin{equation}\label{E2} \lambda=d_1x_1+\frac{\beta x_1y_1}{x_1+y_1},\quad d_2y_1=\frac{\beta e^{-d_1\tau x_1y_1}}{x_1+y_1}. \end{equation} Combining the \eqref{E1} and \eqref{E2} we obtain \begin{align*} &\frac{dL_1}{dt}\Big|_{\eqref{M1}}\\ &= d_1(x_1-x)\big[1-\frac{x_1(x+y_1)}{x(x_1+y_1)}\big] +\frac{\beta x_1y_1}{x_1+y_1}\big[1-\frac{x_1(x+y_1)}{x(x_1+y_1)} +\frac{y(x+y_1)}{y_1(x+y)}-\frac{y}{y_1}\big]\\ &\quad -\frac{y_1}{y}\frac{\beta x(t-\tau)y(t-\tau)}{x(t-\tau)+y(t-\tau)} +e^{d_1\tau}ay_1z-\frac{e^{d_1\tau} ad_3z}{p}+\frac{\beta x_1y_1}{x_1+y_1}\\ &\quad -e^{d_1\tau}d_2y_1\ln\frac{\beta xy}{e^{d_1\tau}d_2y_1(x+y)}+ e^{d_1\tau}d_2y_1\ln\frac{\beta x(t-\tau)y(t-\tau)} {e^{d_1\tau}d_2y_1(x(t-\tau)+y(t-\tau))}\\ &= d_1(x_1-x)\big[1-\frac{x_1(x+y_1)}{x(x_1+y_1)}\big] +\frac{\beta x_1y_1}{x_1+y_1}\Big[\Big(1-\frac{y(x+y_1)}{y_1(x+y)}\Big) \Big(\frac{x+y}{x+y_1}-1\Big) \\ &\quad -\Big(\frac{x+y}{x+y_1}-1-\ln\frac{x+y}{x+y_1}\Big) -\Big(\frac{x_1(x+y_1)}{x(x_1+y_1)}-1-\ln\frac{x_1(x+y_1)}{x(x_1+y_1)}\Big)\\ &\quad -\ln\frac{x+y}{x+y_1}-\ln\frac{x_1(x+y_1)}{x(x_1+y_1)}\Big] +\frac{ae^{d_1\tau}}{p}(py_1-d_3) \\ &\quad -\frac{\beta x_1y_1}{x_1+y_1} \Big[\frac{(x_1+y_1)x(t-\tau)y(t-\tau)}{x_1y(x(t-\tau)+y(t-\tau))}-1 -\ln\frac{(x_1+y_1)x(t-\tau)y(t-\tau)}{x_1y(x(t-\tau)+y(t-\tau))}\Big] \\ &\quad -\frac{\beta x_1y_1}{x_1+y_1} \Big[\ln\frac{(x_1+y_1)x(t-\tau)y(t-\tau)}{x_1y(x(t-\tau)+y(t-\tau))} +\ln\frac{\beta xy}{e^{d_1\tau}d_2y_1(x+y)} \\ &\quad -\ln\frac{\beta x(t-\tau)y(t-\tau)}{e^{d_1\tau}d_1y_1(x(t-\tau) +y(t-\tau))}\Big]\\ &= d_1(x_1-x)\Big[1-\frac{x_1(x+y_1)}{x(x_1+y_1)}\Big] +\frac{\beta x_1y_1}{x_1+y_1}\Big(1-\frac{y(x+y_1)}{y_1(x+y)}\Big) \Big(\frac{x+y}{x+y_1}-1\Big)\\ &\quad -\frac{\beta x_1y_1}{x_1+y_1}g\Big(\frac{x+y}{x+y_1}\Big) -\frac{\beta x_1y_1}{x_1+y_1}g\Big(\frac{x_1(x+y_1)}{x(x_1+y_1)}\Big)\\ &\quad -\frac{\beta x_1y_1}{x_1+y_1}g \Big(\frac{(x_1+y_1)x(t-\tau)y(t-\tau)}{x_1y(x(t-\tau)+y(t-\tau))}\Big) +\frac{ae^{d_1\tau}}{p}(py_1-d_3), \end{align*} where $$ d_1(x_1-x)\big[1-\frac{x_1(x+y_1)}{x(x_1+y_1)}\big] =-\frac{d_1y(x_1-x)^2}{x(x_1+y_1)}\leq 0. $$ Note that \[ \Re_1=\frac{d_1+d_2 e^{d_1\tau}(\Re_0-1)}{p\lambda}(py_1-d_3)+1, \] which implies that $\Re_1 \leq 1$ is equivalent to $py_1/d_3 \leq 1$. The latter $py_1/d_3$ is seen to be the immune reproductive number, which expresses the average number of activated CTLs generated from one CTL during its life time $1/d_3$ through the stimulation of the infected cells $y_1$. It is reasonable that immune is activated in the case where $\Re_1 > 1$. Hence $\frac{dL_1}{dt}$ is always non-positive under the condition $\Re_1\leq 1<\Re_0$, and it can be verified that $\frac{dL_1}{dt}=0$ if and only if $x=x_1$ and $\frac{x+y}{x+y_1}=\frac{x_1(x+y_1)}{x(x_1+y_1)} =\frac{(x_1+y_1)x(t-\tau)y(t-\tau)}{x_1y(x(t-\tau)+y(t-\tau))}=1$. Using the first two equations of system \eqref{M3}, we have \begin{gather*} 0= \dot{x}(t)= \lambda-d_1x_1-\frac{\beta x_1y(t)}{x_1+y(t)}, \\ 0= \dot{y}(t)= \frac{\beta e^{-d_1\tau} x_1y(t-\tau)}{x_1+y(t-\tau)} -d_2y(t)-ay(t)z(t). \end{gather*} This gives $y=y_1$, $z=0$. So, the global asymptotic stability of $E_1$ follows from the LaSalle's invariant principle. \end{proof} \begin{theorem} \label{T3} If $\Re_1>1$, then CTL-activated infection equilibrium $E_2$ is globally asymptotically stable; i.e., $E_2$ is globally asymptotically stable whenever it exists. \end{theorem} \begin{proof} Define a Lyapunov functional \begin{align*} L_2&= x(t)-x_2-\int_{x_2}^{x(t)}\frac{x_2(\theta+y_2)}{\theta(x_2+y_2)}d\theta +e^{d_1\tau}y_2g\Big(\frac{y(t)}{y_2}\Big) +\frac{ae^{d_1\tau}}{p}z_2g\left(\frac{z(t)}{z_2}\right)\\ &\quad + e^{d_1\tau}(d_2y_2+ay_2z_2) \int_{t-\tau}^tg \Big(\frac{\beta x(\theta)y(\theta)}{e^{d_1\tau}(d_2y_2+ay_2z_2)(x(\theta) +y(\theta))}\Big)d\theta. \end{align*} Calculating the time derivative of $L_2$ along the solution of \eqref{M1}, we obtain \begin{equation} \label{E3} \begin{aligned} \frac{dL_2}{dt}\Big|_{\eqref{M1}} &= \Big[1-\frac{x_2(x+y_2)}{x(x_2+y_2)}\Big] \Big[d_1(x_2-x)+\Big(\frac{\beta x_2y_2}{x_2+y_2}-\frac{\beta xy}{x+y}\Big)\Big] \\ &\quad + \Big[\frac{\beta x(t-\tau)y(t-\tau)}{x(t-\tau)+y(t-\tau)} -e^{d_1\tau} (d_2y(t)-ay(t)z(t))\Big]\big(1-\frac{y_2}{y}\big) \\ &\quad +\frac{ae^{d_1\tau}}{p}(pyz-d_3z)\big(1-\frac{z_2}{z}\big) +\frac{\beta xy}{x+y}-\frac{\beta x(t-\tau)y(t-\tau)}{x(t-\tau)+y(t-\tau)} \\ &\quad -\frac{\beta x_2y_2}{x_2+y_2}\ln\frac{\beta xy}{\frac{\beta x_2y_2}{x_2+y_2}(x+y)}+ \frac{\beta x_2y_2}{x_2+y_2}\ln\frac{\beta x(t-\tau)y(t-\tau)} {\frac{\beta x_2y_2}{x_2+y_2}(x(t-\tau)+y(t-\tau))}. \end{aligned} \end{equation} Here we used that \begin{equation}\label{E4} \lambda=d_1x_2+\frac{\beta x_2y_2}{x_2+y_2},\quad d_2y_2+ay_2z_2=\frac{\beta e^{-d_1\tau}x_2y_2}{x_2+y_2}, \quad py_2=d_3. \end{equation} Combining \eqref{E3} and \eqref{E4} we obtain \begin{align*} &\frac{dL_2}{dt}\Big|_{\eqref{M1}}\\ &= d_1(x_2-x)\Big[1-\frac{x_2(x+y_2)}{x(x_2+y_2)}\Big] +\frac{\beta x_2y_2}{x_2+y_2} \Big[1-\frac{x_2(x+y_2)}{x(x_2+y_2)}+\frac{y(x+y_2)}{y_2(x+y)}-\frac{y}{y_2} \Big]\\ &\quad +\frac{\beta x_2y_2}{x_2+y_2}-\frac{y_2}{y} \frac{x(t-\tau)y(t-\tau)}{x(t-\tau)+y(t-\tau)}\\ &\quad -\frac{\beta x_2y_2}{x_2+y_2}\ln\frac{\beta xy}{\frac{\beta x_2y_2} {x_2+y_2}(x+y)}+ \frac{\beta x_2y_2}{x_2+y_2}\ln\frac{\beta x(t-\tau)y(t-\tau)} {\frac{\beta x_2y_2}{x_2+y_2}(x(t-\tau)+y(t-\tau))} \\ &= d_1(x_2-x)\Big[1-\frac{x_2(x+y_2)}{x(x_2+y_2)}\Big] +\frac{\beta x_2y_2}{x_2+y_2}\Big[\Big(1-\frac{y(x+y_2)}{y_2(x+y)}\Big) \Big(\frac{x+y}{x+y_2}-1\Big)\\ &\quad -\Big(\frac{x+y}{x+y_2}-1-\ln\frac{x+y}{x+y_2}\Big) -\Big(\frac{x_2(x+y_2)}{x(x_2+y_2)}-1-\ln\frac{x_2(x+y_2)}{x(x_2+y_2)}\Big)\\ &\quad -\ln\frac{x+y}{x+y_2}-\ln\frac{x_2(x+y_2)}{x(x_2+y_2)}\Big]\\ &\quad -\frac{\beta x_2y_2}{x_2+y_2} \Big[\frac{(x_2+y_2)x(t-\tau)y(t-\tau)}{x_2y(x(t-\tau) +y(t-\tau)}-1-\ln\frac{(x_2+y_2)x(t-\tau)y(t-\tau)}{x_2y(x(t-\tau) +y(t-\tau)}\Big]\\ &\quad -\frac{\beta x_2y_2}{x_2+y_2} \Big[\ln\frac{(x_2+y_2)x(t-\tau)y(t-\tau)}{x_2y(x(t-\tau)+y(t-\tau)} +\ln\frac{\beta xy}{\frac{\beta x_2y_2}{x_2+y_2}(x+y)} \\ &\quad -\ln\frac{\beta x(t-\tau)y(t-\tau)} {\frac{\beta x_2y_2}{x_2+y_2}(x(t-\tau)+y(t-\tau))}\Big]\\ &= d_1(x_2-x)\Big[1-\frac{x_2(x+y_2)}{x(x_2+y_2)}\Big] +\frac{\beta x_2y_2}{x_2+y_2}\Big(1-\frac{y(x+y_2)}{y_2(x+y)}\Big) \Big(\frac{x+y}{x+y_2}-1\Big)\\ &\quad -\frac{\beta x_2y_2}{x_2+y_2}g\Big(\frac{x+y}{x+y_2}\Big) -\frac{\beta x_2y_2}{x_2+y_2}g\Big(\frac{x_2(x+y_2)}{x(x_2+y_2)}\Big)\\ &\quad -\frac{\beta x_2y_2}{x_2+y_2}g \Big(\frac{(x_2+y_2)x(t-\tau)y(t-\tau)}{x_2y(x(t-\tau)+y(t-\tau))}\Big). \end{align*} Similar to the proof of Theorem \ref{T2}, the terms of $dL_2/dt$ always are non-positive. Hence $dL_2/dt$ for all $x > 0, y > 0 $and $z > 0$, and $dL_2/dt=0$ if and only if $x = x_2$ and $y=y_2, z=z_2$. The largest invariant set in $\{(x_t, y_t,z_t)\mid dL_2/dt=0\}$ is $E_2$. From the Lyapunov-LaSalle invariance principle, it shows that equilibrium $E_2(x_1, y_2, z_2)$ is globally asymptotically stable. \end{proof} \section{Summary and Discussion} In this article, we have modified the delay differential equation model for cell-to-cell infection of HBV in tissue cultures proposed by Gourley et al \cite{Gourley} by incorporating a Cytotoxic T lymphocytes (CTLs) immune response to model the role in antiviral by attacking virus infected cells. Since immune response after viral infection is universal and necessary to eliminate or control the disease. Our analysis provides a quantitative understandings of HBV replication dynamics in vivo and has implications for the optimal timing of drug treatment and immunotherapy in chronic HBV infection. By constructing Lyapunov functionals, we obtain the global stability of the equilibria of \eqref{M3} that depends only on the basic reproductive number $\Re_0$ and the basic immune reproductive number $\Re_1$. For delay differential equations model \eqref{M3}, the basic reproductive number is given by $\Re_0=\frac{\beta e^{-d\tau}}{d_2}$ and it is a decreasing function on intracellular delay $\tau$ such that $\Re_0(\infty) = 0$. Theorems \ref{T1}-\ref{T3} show that when $\Re_0 \leq 1$, the infection-free equilibrium is globally asymptotically stable, which means that the viruses are cleared and immune is not active; when $\Re_1 \leq 1 < \Re_0$, the CTL-inactivated infection equilibrium exists and is globally asymptotically stable, which means that CTLs immune response would not be activated and viral infection becomes chronic but with a low level of proviral load; and when $\Re_1 > 1$, the CTL-activated infection equilibrium exists and is globally asymptotically stable, in this case the infection causes a persistent CTLs immune response and is chronic with a high level of proviral load. We can see that under the condition $\Re_1(\tau) > 1$, as delay $\tau$ increases, the number of CTLs immune response does not change in this situation. However, when the delay $\tau$ is sufficiently large, and brings $\Re_1(\tau)$) to a level lower than unity, the CTL-inactivated infection equilibrium $E_1$ becomes globally asymptotically stable. It has been repointed in Nowak \cite{Nowak1} that \lq\lq Treatment of chronic HBV infections with lamivudine leads to a rapid and sustained decline of plasma virus levels, but clinical benefit with a reduced risk of cirrhosis and development of liver cancer will greatly depend on the decline of infected cells. Immunotherapy without antiviral treatment could be problematic because of the very large number of infected liver cells in the typical HBV carrier. Therefore, the drugs which can prolong the latent period, and/or decrease the needed time of immune response activation and/or inhibit infection can slow down the virus production process. This gives us a good guidance on the development of treatment strategies. On the other hand, cell-to-cell models may be applicable to study the within-host dynamics of other types of viral infections such as human T-cell leukaemia virus type 1 (HTLV-1), hepatitis C, etc. We leave the modeling and study of the cell-to-cell HTLV-1 infection for future consideration. Other realistic modifications can be made. For example, we can modify target-cell dynamics to be a mitosis component given by a logistic term and the loss/gain term as nonlinear incidence function. Another possible modification would be to incorporate diffusion term into the delayed model to more accurately reflect the realistic situation in tissue cultures. \subsection*{Acknowledgments} The authors want to thank the anonymous referees and the editor for their valuable suggestions and comments. J. Wang is supported by National Natural Science Foundation of China (TianYuan No.11226255). \begin{thebibliography}{99} \bibitem{Arnaout} R. Arnaout, M. A. Nowak, D. Wodarz; HIV-1 dynamics revisited: biphasic decay by cytotoxic lymphocyte killing?, \emph{Proceedings of the Royal Society B: Biological Sciences}, \textbf{265} (2000), 1347--1354. \bibitem{Bangham} C. R. M. Bangham; The immune control and cell-to-cell spread of human T-lymphotropic virus type 1, \emph{J. Gen. Virol.}, \textbf{84} (2003), 3177--3189. \bibitem{Culshaw} R. Culshaw, S. Ruan, R. Spiteri; Optimal HIV treatment by maximising immune response, \emph{J. Math. Biol.}, \textbf{48} (2004), 545--562. \bibitem{Dimitrov} D. S. Dimitrov, R. L. Willey, H. Sato, L. J. Chang, R. Blumenthal, M. A. Martin; Quantitation of human immunodeficiency virus type 1 infection kinetics, \emph{J. Virol.}, \textbf{67} (1993), 2182--2190. \bibitem{Ferrari} C. Ferrari, A. Penna, et al.; Cellular immune response to hepatitis B virus-encoded antigens in acute and chronic hepatitis B virus infection, \emph{J. Immunol.}, \textbf{145} (1990), 3442--3449. \bibitem{Gourley} S. A. Gourley, Y. Kuang, J. D. Nagy; Dynamics of a delay differential equation o hepatitis B virus inection, \emph{J. Bio. Dyn.}, \textbf{2} (2008), 140--153. \bibitem{Gummuluru} S. Gummuluru, C. M. Kinsey, M. Emerman; An in vitro rapid-turnover assay for human immunodeficiency virus type 1 replication selects for cell-to-cell spread of virua, \emph{J. Virol.}, \textbf{74} (2000), 10882--10891. \bibitem{GuoCai} B. Z. Guo, L. M. Cai; A note for the global stability of a delay differential equation of hepatitis B virus, \emph{Math. Bios. Eng.}, \textbf{8} (2011), 689--694. \bibitem{Guidotti} L. G. Guidotti, F. V. Chisari; Immunobiology and pathogenesis of viral hepatitis, \emph{Ann. Rev.: Pathology Mech. Disease}, \textbf{1} (2006), 23--61. \bibitem{Haase1} A. T. Haase; Population biology of HIV-1 infection: viral and CD4+ T cell demographics and dynamics in lymphatic tissues, \emph{Annu. Rev. Immunol.}, \textbf{17} (1999), 625--656. \bibitem{Haase2} A. T. Haase, K. Henry, M. Zupancic, et al.; Quantitative image analysis of HIV-1 infection in lymphoid tissue, \emph{Science}, \textbf{274} (1996), 985--989. \bibitem{Huang} G. Huang, Y. Takeuchi, W. MA, D. Wei; Global tability for delay SIR and SEIR epidemic models with nonlinear incidence rate, \emph{Bull. Math. Biol.}, \textbf{72} (2010), 1192--1207. \bibitem{Haddock} J. R. Haddock, J. Terj\'{a}ki; Liapunov-Razumikhin functions and an invariance principle for functionaldifferential equations, \emph{J. Differ. Equ.}, \textbf{48} (1983), 95--122. \bibitem{Iannacone} M. Iannacone; HBV pathogenesis in animal models: Recent advances on the role of platelets, \emph{J. Hepatology}, \textbf{46} (2007), 719--726. \bibitem{Kuang} Y. Kuang; \emph{Delay Differential Equations with Applications in Population Dynamics}, Boston: Academics Press, 1993. \bibitem{Korobeinikov} A. Korobeinikov; Global properties of basic virus dynamics models, \emph{Bull. Math. Biol.}, \textbf{66} (2004), 879--883. \bibitem{Min} L.Q. Min, Y.M. Su, Y. Kuang; Mathematical analysis of a basic virus infection model with application to HBV infection, \emph{Rocky Mountain J. Math.}, \textbf{38} (2008), 1573--1585. \bibitem{McCluskey} C. C. McCluskey; Global stability for an SEIR epidemiological model with varying infectivity and infinite delay, \emph{Math. Biosci. Eng.}, \textbf{6} (2009), 603--610. \bibitem{Murray} J. M. Murray, R. H. Purcell, S. F. Wieland; The half-life of hepatitis B virions, \emph{Hepatology}, \textbf{44} (2006), 1117--1121. \bibitem{Nowak1} M. A. Nowak, et al.; Viral dynamics in hepatitis B virus infection, \emph{Proc. Nat. Acad. Sci., USA}, \textbf{93} (1996), 4398--4402. \bibitem{Nowak2} M. A. Nowak, R. M. May; \emph{Viral Dynamics}, Oxford University Press, Oxford, 2000. \bibitem{Nowak3} M. A. Nowak, C. Bangham; Population dynamics of immune response to persistent viruses, \emph{Science}, \textbf{272} (1996), 74--79. \bibitem{Nowak4} M. A. Nowak, R. M. May; Mathematical biology of HIV infections: Antigenic varia- tion and diversity threshold, \emph{Math. Biosci.}, \textbf{106} (1991), 1--21. \bibitem{Philips} D. M. Philips; The role of cell-to-cell transmission in HIV infection, \emph{AIDS}, \textbf{8} (1994), 719--731. \bibitem{Perelson3} A. S. Perelson, P. W. Nelson; Mathematical analysis of HIV-1 dynamics in vivo, \emph{SIAM Rev.}, \textbf{41} (1999), 3--44. \bibitem{Regenstein} F. Regenstein; New approaches to the treatment of chronic viral hepatitis B and C, \emph{Am. J. Med.}, \textbf{96} (1994), 47--51. \bibitem{Seeger} C. Seeger, W. S. Mason; Hepatitis B virus biology, \emph{Microbiol. Mol. Biol. Rev.}, \textbf{64} (2000), 51--68. \bibitem{Sato} H. Sato, J. Orenstein, D. S. Dimitrov, M. A. Martin; Cell-to-cell spread of HIV-1 occurs with minutes and may not involve the participation of virus particles, \emph{Virology}, \textbf{186} (1992), 712--724. \bibitem{Song} X. Q. Song, A. U. Neuman; Global stability and periodic solution of the viral dynamics, \emph{J. Math. Anal. Appl.}, \textbf{329} (2007), 281--297. \bibitem{WHO} WHO, 2008. Fact sheet N204 Hepatitis B, available at http://www.who.int/media centre/factsheets/fs204/en/index.html. \bibitem{Wang1} J. Wang, G. Huang, Y. Takeuch, S. Liu; SVEIR epidemiological model with varying infectivity and distributed delays, \emph{Math. Biosci. Eng.}, \textbf{8} (2011), 875--888. \bibitem{Wang2} J. Wang, G. Huang, Y. Takeuchi; Global asymptotic stability for HIV-1 dynamics with two distributed delays, \emph{Math. Medic. Bio.}, doi:10.1093/imammb/dqr009. \bibitem{WWL} K. F. Wang, W. D. Wang, X. N. Liu; Global stability in a viral infection model with lytic and nonlytic immune response, \emph{Computers and Mathematics with Applications}, \textbf{51} (2006), 1593--1610. \bibitem{Zheng} Y. Zheng, L. Q. Min, Y. Ji, Y. M. Su, Y. Kuang; Global stability of endemic equilibrium point of basic virus infection model with application to HBV infection, \emph{J. Syst. Sci. Complex.}, \textbf{23} (2010), 1221--1230. \end{thebibliography} \end{document}