A Non-Standard Finite Difference Scheme for a Diffusive HIV-1 Infection Model with Immune Response and Intracellular Delay

Abstract

In this paper, we propose and study a diffusive HIV infection model with infected cells delay, virus mature delay, abstract function incidence rate and a virus diffusion term. By introducing the reproductive numbers for viral infection $R_0$ and for CTL immune response number $R_1$, we show that $R_0$ and $R_1$ act as threshold parameter for the existence and stability of equilibria. If $R_0\le 1$, the infection-free equilibrium $E_0$ is globally asymptotically stable, and the viruses are cleared; if $R_1\le 1<R_0$, the CTL-inactivated equilibrium $E_1$ is globally asymptotically stable, and the infection becomes chronic but without persistent CTL response; if $R_1>1$, the CTL-activated equilibrium $E_2$ is globally asymptotically stable, and the infection is chronic with persistent CTL response. Next, we study the dynamic of the discreted system of our model by using non-standard finite difference scheme. We find that the global stability of the equilibria of the continuous model and the discrete model is not always consistent. That is, if $R_0\le 1$, or $R_1\le 1<R_0$, the global stability of the two kinds model is consistent. However, if $R_1>1$, the global stability of the two kinds model is not consistent. Finally, numerical simulations are carried out to illustrate the theoretical results and show the effects of diffusion factors on the time-delay virus model.

1. Introduction

In the past few years, host-virus dynamics models have been developed to explain the interactions between virus and target T cells, much attention has been given to the role of the immune response to human immunodeficiency virus (HIV) infection. Many different mechanisms of immune system, defenses against viral infections are of interest because lots of the diseases caused by them, e.g., hepatitis B and AIDS, are chronic and incurable [1,2]. With the new coronavirus epidemic rages around the world [3,4,5], virus dynamics has become a hot spot again. In the immune response mechanism in vivo for viral infections, the cytotoxic T lymphocyte ($CTL$) plays a particularly important role, therefore many authors have examined various $CTL$ dynamics.

A virus must take over host cells and use them to replicate because it can not replicate on its own. HIV targets the $CD{4}^{+}T$ cells, often referred to as “helper” T cells, when it invades the body. These cells can be considered “messengers”, or the command centres of the immune system. They send signal to other immune cells that an invader is to be fought. Once invaded by the viruses, these infected cells will cause a cytotoxic T-lymphocyte ($CTL$) response from the immune system. The immune response cells, or cytotoxic lymphocytes, respond to this message and set out to eliminate infection by killing infected cells. Through the lysis of the infected cells, the viruses are prevented from further replication [2]. The $CTL$ response is also striking in that it sometime does damage to the body when it tries to clear the virus. Over half the tissue damage caused by hepatitis is actually caused by the $CTL$ response [1,6].

If the immune system is functioning normally, these components work together efficiently and an infection is eliminated quickly, causing only temporary discomfort to the host. However, over time HIV is able to deplete the population of $CD{4}^{+}T$ cells. What remains unknown is the exact mechanism by which this occurs, but several models have been suggested. For a variety of different hypotheses of how this occurs, we refer the reader to papers [7,8,9]. The natural killer cells may be fit to eliminate infection, but they are never deployed, which is the the impact of the depletion of $CD{4}^{+}T$ cells on the host. This then culminates in a clinical problem wherein the patient becomes vulnerable to infections that a healthy immune system would normally handle.

Quite a lot of mathematical models of HIV have been set up. The classical model is a system with three ordinary differential equations [10,11]. To better understanding the dynamics of these infections, many mathematical models have been proposed by using different kinds of differential equations, see [12,13,14,15,16] and references therein. For example, Yang et al. [15] studied the following model

(1)$\left\{\begin{array}{c}\frac{\partial T(x,t)}{\partial t}=\lambda -d_{1}T(x,t)-{\beta }_{1}T(x,t)V(x,t),\hfill \\ \frac{\partial I(x,t)}{\partial t}={\beta }_{1}T(x,t)V(x,t)-d_{2}I(x,t),\hfill \\ \frac{\partial V(x,t)}{\partial t}=d\Delta V(x,t)+\gamma I(x,t)-d_{3}V(x,t),\hfill \end{array}\right.$

where $T(x,t),I(x,t)$ and $V(x,t)$ denote the densities of uninfected cells, infected cells and free virus cells at position x at time t, respectively. $\lambda$ stands for the recruitment rate of the uninfected cells; ${\beta }_1$ is the virus-to-cell infection rate; $d_1,d_2$ and $d_3$ represent death rates of uninfected cells, infected cells and free viruses; $\gamma$ stands for the recruitment rate for free viruses; d stands for the diffusion coefficient and $\Delta$ is the Laplacian operator.

To help the body heal, cytotoxic T-lymphocyte effectors ($CTLe$) of the immune system will remove the infected cells to prevent further viral replications. To model these extra dynamics, researchers have studied the model of viral interaction with $CTL$ response [10,17]

(2)$\left\{\begin{array}{c}\dot{x}=\lambda -dx-\beta xy,\hfill \\ \dot{y}=\beta xy-ay-pyz,\hfill \\ \dot{z}=cyz-bz,\hfill \end{array}\right.$

where variables $x,y$ and z denote the density of the healthy cells, the infected cells, and the $CTLs$ populations, respectively. Healthy cells are produced at rate $\lambda$ and their natural mortality is $dx$; these cells may come into contact with the virus and become infected cells at rate $\beta xy$, infected cells’s natural mortality is $ay$, and they are removed by $CTLs$ at rate $pyz$; the $CTL$ population increases at the rate $cyz$ and they are removed at the rate $bz$.

In [18,19], researchers studied a mathematical model for HIV-1 infection with both intracellular delay and cell-mediated immune response:

(3)$\left\{\begin{array}{c}\frac{dx\left(t\right)}{dt}=\lambda -dx\left(t\right)-\beta xv,\hfill \\ \frac{dy\left(t\right)}{dt}=e^{-a\tau }\beta x(t-\tau )v(t-\tau )-ay\left(t\right)-py\left(t\right)z\left(t\right),\hfill \\ \frac{dv\left(t\right)}{dt}=ky\left(t\right)-\mu v\left(t\right),\hfill \\ \frac{dz\left(t\right)}{dt}=cy\left(t\right)z\left(t\right)-bz\left(t\right).\hfill \end{array}\right.$

Researchers obtain the global stability of the infection-free equilibrium and give many conditions for the local stability of the two infection equilibria: one without $CTLs$ being activated and the other with. There are many references in the dynamics of HIV-1 infection with $CTLs$ response (see, e.g., [17,20,21,22] and the references therein).

However, there is no diffusion term and only one delay in (3). As we know, the virus is not stationary in space, the movement of the virus in space leads to the spatial spread of the disease, and mostly with general nonlinear incidence rate. Fickian diffusion can reasonably describe the spread of this virus in space and this diffusion process is often represented by the Laplace operator. Inspired by [16,23], in this paper, we extend the classic model of virus dynamics to a diffusive infection model with intracellular delay and cell-mediated immune response, with two delays and general nonlinear incidence rate, as follows

(4)$\left\{\begin{array}{c}\frac{\partial T(x,t)}{\partial t}=\lambda -d_{1}T(x,t)-{\beta }_{1}T(x,t)f\left(V(x,t)\right)-{\beta }_{2}T(x,t)g\left(I(x,t)\right),\hfill \\ \frac{\partial I(x,t)}{\partial t}=e^{-{\mu }_1{\tau }_1}\left({\beta }_{1}T(x,t-{\tau }_1)f\left(V(x,t-{\tau }_1)\right)+{\beta }_{2}T(x,t-{\tau }_1)g\left(I(x,t-{\tau }_1)\right)\right)\hfill \\ -d_{2}I(x,t)-p_{1}I(x,t)Z(x,t),\hfill \\ \frac{\partial V(x,t)}{\partial t}=D\Delta V(x,t)+p_2{e}^{-{\mu }_2{\tau }_2}I(x,t-{\tau }_2)-d_{3}V(x,t),\hfill \\ \frac{\partial Z(x,t)}{\partial t}=qI(x,t)Z(x,t)-d_{4}Z(x,t),\hfill \end{array}\right.$

here $T(x,t),I(x,t),V(x,t)$ and $Z(x,t)$ stand for the densities of uninfected cells, infected cells, virus cells and $CTLs$ at position x at time t, respectively. $\lambda$ and $d_1$ denote the natural produce and mortality rate of uninfected cells, and uninfected cells are infected with a rate ${\beta }_2$; and ${\beta }_1$ is the virus-to-cell infection rate; and ${\beta }_1$ is the virus-to-cell infection rate; the natural mortality rate of the infected cells are $d_2$ and are killed by $CTL$ with a rate $p_1$ (Note that $d_2$ reflects the combined effects of natural death rate of uninfected cells, $d_1$, and any additional cytotoxic effects the virus may have); ${\mu }_1$ represents the death rate for infected but not yet virus-producing cells, ${\tau }_1$ represents the latent delay, i.e., the time period from being infected to becoming productive infected cells. Therefore, the probability of surviving from time $t-{\tau }_1$ to time t is $e^{-{\mu }_1{\tau }_1}$; the probability of survival of immature virus is denoted by $e^{-{\mu }_2{\tau }_2}$ and the average life time of an immature virus is $\frac{1}{{\mu }_2}$; where ${\tau }_2$ represents the time necessary for the newly produced virus to become mature; D is the diffusion coefficient and $\Delta$ is the Laplacian operator; $p_2$ is the recruitment rate for free viruses. Virus particles are removed from the system at rate $d_3$; q stands for the $CTL$ responsiveness and $d_4$ denotes decay rate for $CTLs$ in the absence of stimulation.

Here, the incidences are assumed to be the nonlinear responses to the concentrations of virus particles and infected cells, using the forms ${\beta }_{1}Tf\left(V\right)$ and ${\beta }_{2}Tg\left(I\right)$, where $f\left(V\right)$ and $g\left(I\right)$ are the force of infection by virus particles and infected cells and satisfy the following properties [24]:

$f\left(0\right)=g\left(0\right)=0,f^{\prime }\left(V\right)>0,g^{\prime }\left(I\right)>0,f^{″}\left(V\right)\le 0,g^{″}\left(I\right)\le 0.\left(A_1\right)$

It follows from $\left(A_1\right)$ and the Mean Value Theorem that

$f^{\prime }\left(V\right)V\le f\left(V\right)\le f^{\prime }\left(0\right)V,g^{\prime }\left(I\right)I\le g\left(I\right)\le g^{\prime }\left(0\right)I,forI,V\ge 0.\left(A_2\right)$

Epidemiologically, condition $\left(A_1\right)$ implies that: $\left(1\right)$ the disease cannot spread if there is no infection; $\left(2\right)$ the incidences ${\beta }_{1}Tf\left(V\right)$ and ${\beta }_{2}Tg\left(I\right)$ become faster when the densities of the virus particles and infected cells increase; $\left(3\right)$ the per capita infection rates by virus particles and infected cells will slow down as certain inhibiting effect since $\left(A_2\right)$ implies that ${\left(\frac{f\left(V\right)}{V}\right)}^{\prime }\le 0$ and ${\left(\frac{g\left(I\right)}{I}\right)}^{\prime }\le 0$. The incidence rate with condition $\left(A_1\right)$ contains the bilinear and the saturation incidences.

In this paper, we will consider the system (4) with initial conditions

(5)$\begin{array}{c}T(x,s)={\varphi }_1(x,s)\ge 0,I(x,s)={\varphi }_2(x,s)\ge 0,\hfill \\ V(x,s)={\varphi }_3(x,s)\ge 0,Z(x,s)={\varphi }_4(x,s)\ge 0,(x,s)\in \overline{\Omega }\times [-\tau ,0]\hfill \end{array}$

and homogeneous Neumann boundary conditions

(6)$\begin{array}{c}\frac{\partial V}{\partial n}=0,t>0,x\in \partial \Omega .\hfill \end{array}$

where $\tau =max\{{\tau }_1,{\tau }_2\}$ and $\Omega$ is a bounded domain in $R^4$ with smooth boundary $\partial \Omega$, and $\frac{\partial V}{\partial n}$ stands for the outward normal derivative on $\partial \Omega$.

Usually, the exact solution for a system as (1) is difficult or even impossible to be determined. Hence, researchers seek numerical ones instead. However, how to choose the proper discrete scheme so that the global dynamics of solutions of the corresponding continuous models can be efficiently preserved is still an open problem [25]. Mickens has made an attempt in this connection, by presenting a robust non-standard finite difference (NSFD)scheme [26], which has been widely employed in the study of different epidemic models [23,25,26,27,28,29,30,31,32]. For example, Yang et al. [30] applied the NSFD scheme to discretize system (1) and found that the dynamical behaviors of the discrete model are consistent with the original system. Motivated by the work of [23,25,26,27,28,29,30,31,32], we apply the NSFD scheme to discretize system (4) and obtain:

(7)$\left\{\begin{array}{c}\frac{T_{n+1}^m-T_n^m}{\Delta t}=\lambda -d_1{T}_{n+1}^m-{\beta }_1{T}_{n+1}^{m}f\left(V_n^m\right)-{\beta }_2{T}_{n+1}^{m}g\left(I_n^m\right),\hfill \\ \frac{I_{n+1}^m-I_n^m}{\Delta t}=e^{-{\mu }_1{\tau }_1}\left({\beta }_1{T}_{n-m_1+1}^{m}f\left(V_{n-m_1}^m\right)+{\beta }_2{T}_{n-m_1+1}^{m}g\left(I_{n-m_1}^m\right)\right)\hfill \\ -d_2{I}_{n+1}^m-p_1{I}_{n+1}^m{Z}_n^m,\hfill \\ \frac{V_{n+1}^m-V_n^m}{\Delta t}=D\frac{V_{n+1}^{m+1}-2V_{n+1}^m+V_{n+1}^{m-1}}{{(\Delta x)}^2}+p_2{e}^{-{\mu }_2{\tau }_2}{I}_{n-m_2+1}^m-d_3{V}_{n+1}^m,\hfill \\ \frac{Z_{n+1}^m-Z_n^m}{\Delta t}=q{I}_{n+1}^m{Z}_n^m-d_4{Z}_{n+1}^m.\hfill \end{array}\right.$

Here, we assume that $x\in \Omega =[a,b]$, let $\Delta t>0$ be the time step size and $\Delta x=\frac{b-a}{N}$ be the space step size with N a positive integer. Suppose that there exist two integers $m_1,m_2\in \aleph$ with ${\tau }_1=m_1\Delta t,{\tau }_2=m_2\Delta t$. Denote the mesh grid point as $\{(x_m,t_n),m=0,1,2,\cdots ,N,n\in N\}$ with $x_m=a+m\Delta x$ and $t_n=n\Delta t$. At each point, we use approximations of $\left(T(x_m,t_n),I(x_m,t_n),V(x_m,t_n),Z(x_m,t_n)\right)$ by $(T_n^m,I_n^m,V_n^m,Z_n^m)$. We set all the approximation solutions at the time $t_n$ by the $N+1$-dimensional vector $U_n={(U_n^0,U_n^1,\cdots ,U_n^N)}^T$, where $U_n^{(·)}\in \left\{(T_n,I_n,V_n,Z_n)\right\}$ and the notation ${(·)}^T$ is the transposition of a vector. $U_n\ge 0$ means that all components of a vector $U_n$ are nonnegative. The discrete initial conditions of system (7) are given as

(8)$\begin{array}{c}{T}_s^m={\varphi }_1(x_m,t_s)\ge 0,I_s^m={\varphi }_2(x_m,t_s)\ge 0,\hfill \\ V_s^m={\varphi }_3(x_m,t_s)\ge 0,Z_s^m={\varphi }_4(x_m,t_s)\ge 0,\hfill \end{array}$

for all $s=-l,-l+1,\cdots ,0,l=max\{m_1,m_2\}$, and the discrete boundary conditions are

$V_n^{-1}=V_n^0,V_n^N=V_n^{N+1},forn\in \aleph .$

The main purpose of this paper is to investigate the asymptotic stability of system (4) and (7). Another purpose of this paper is to discuss, whether the discretized system (7) that derived by using NSFD scheme can efficiently preserves the global asymptotic stability of the equilibria to the original system (4) or not.

The paper is organized as follows. In Section 2.1, the model is introduced, and, under some assumptions, positivity and boundedness properties of the solutions are proved by using nonlinear functional analysis methods. In Section 2.2, we consider the existence of infection-free equilibrium, $CTL$-inactivated equilibrium and infection equilibrium with immunity. In Section 2.3, by introducing the reproductive numbers for viral infection $R_0$ and for $CTL$ immune response number $R_1$, we show that $R_0$ and $R_1$ act as threshold parameter for the existence and stability of equilibria. If $R_0\le 1$, the infection-free equilibrium $E_0$ is globally asymptotically stable, and the viruses are cleared; if $R_1\le 1<R_0$, the $CTL$-inactivated equilibrium $E_1$ is globally asymptotically stable, and the infection becomes chronic but without persistent $CTL$ response; if $R_1>1$, the $CTL$-activated equilibrium $E_2$ is globally asymptotically stable, and the infection is chronic with persistent $CTL$ response. In Section 3, we investigate the global dynamics of discrete system (7) correponding to the continuous system (4), by using nonstandard finite difference scheme. We find that the global stability of the equilibria of the continuous model and the discrete model is not always consistent. That is, if $R_0\le 1$, or $R_1\le 1<R_0$, the global stability of the two kinds model is consistent. However, if $R_1>1$, the global stability of the two kinds model is not consistent. In Section 4, some numerical simulations are given to illustrate the theoretical results and show the effects of diffusion factors on the time-delay virus model. The paper ends with a discussion in Section 5.

2. Dynamical Behaviors of Continuous System

2.1. Positivity and Boundedness of Solutions

In order to study positivity and boundedness of solutions to system (4), we first introduce some notations.

Assume $X=C(\overline{\Omega },{\mathbb{R}}^4)$ be the space of continuous functions from the topological space $\overline{\Omega }$ into the space ${\mathbb{R}}^4$. Let $C=C\left(\right[-\tau ,0],X)$ be the Banach space of continuous functions from $[-\tau ,0]$ into X with the usual supremum normal. $\varphi \in C$ is defined by

$\varphi (x,s)=\varphi \left(s\right)\left(x\right).$

Define $x_t\left(s\right)=x(t+s),s\in [-\tau ,0]$, where $x(·):[-\tau ,\sigma )\to X$ is a continuous function from $[0,\sigma )$ to C.

Theorem 1.

For any $\varphi \in C$,

$\left(a\right)$ system (4)–(6) has a unique solution defined on $[0,+\infty )$; and

$\left(b\right)$ the solution of (4)–(6) is nonnegative and bounded for all $t\ge 0$.

Proof. 

For any $\varphi ={({\varphi }_1,{\varphi }_2,{\varphi }_3,{\varphi }_4)}^T\in C$ and $x\in \overline{\Omega }$, assume

$F=(F_1,F_2,F_3,F_4):C\to X$

by

$\left\{\begin{array}{c}{F}_1\left(\varphi \right)\left(x\right)=\lambda -d_1{\varphi }_1(x,0)-{\beta }_1{\varphi }_1(x,0)f\left({\varphi }_3(x,0)\right)-{\beta }_2{\varphi }_1(x,0)g\left({\varphi }_2(x,0)\right),\hfill \\ F_2\left(\varphi \right)\left(x\right)=e^{-{\mu }_1{\tau }_1}\left[{\beta }_1{\varphi }_1(x,-{\tau }_1)f\left({\varphi }_3(x,-{\tau }_1)\right)+{\beta }_2{\varphi }_1(x,-{\tau }_1)g\left({\varphi }_2(x,-{\tau }_1)\right)\right]\hfill \\ -d_2{\varphi }_2(x,0)-p_1{\varphi }_2(x,0){\varphi }_4(x,0),\hfill \\ F_3\left(\varphi \right)\left(x\right)=p_2{e}^{-{\mu }_2{\tau }_2}k{\varphi }_2(x,-{\tau }_2)-d_3{\varphi }_3(x,0),\hfill \\ F_4\left(\varphi \right)\left(x\right)=q{\varphi }_2(x,0){\varphi }_4(x,0)-d_4{\varphi }_4(x,0).\hfill \end{array}\right.$

Then system (4)–(6) can be rewritten as following form

(9)$\left\{\begin{array}{c}{U}^{\prime }\left(t\right)=AU+F\left(U_t\right),t>0,\hfill \\ U\left(0\right)=\varphi \in X,\hfill \end{array}\right.$

where $U={(T,I,V,Z)}^T$, $\varphi ={({\varphi }_1,{\varphi }_2,{\varphi }_3,{\varphi }_4)}^T$ and $AU={(0,0,d\Delta v,0)}^T$. It is clear that the operator F is locally Lipschitz in space X. From [27,32,33,34,35,36], we conclude that system (9) has a unique local solution on $t\in [0,T_{max})$, where $T_{max}$ is the maximal existence time for solution of system (4). In addition, it follows from 0 is a sub-solution of each equation of system (4) that $T(x,t)\ge 0,I(x,t)\ge 0,V(x,t)\ge 0,Z(x,t)\ge 0$.

Next, we prove the boundedness of solutions. Let

$G_1(x,t)=e^{-{\mu }_1{\tau }_1}T(x,t-{\tau }_1)+I(x,t)+\frac{p_1}{q}Z(x,t),$

then

$\begin{array}{ccc}\hfill \frac{\partial G_1(x,t)}{\partial t}& =& \lambda e^{-{\mu }_1{\tau }_1}-d_1{e}^{-{\mu }_1{\tau }_1}T(x,t-{\tau }_1)-d_{2}I(x,t)-\frac{p_1{d}_4}{q}Z(x,t)\hfill \\ & \le & \lambda -\tilde{d}{G}_1(x,t),\hfill \end{array}$

where $\tilde{d}=min\{d_1,d_2,d_4\}$, then

$G_1(x,t)\le max\left\{\frac{\lambda }{\tilde{d}},ma{x}_{x\in \overline{\Omega }}\left\{e^{-{\mu }_1{\tau }_1}{\varphi }_1(x,-{\tau }_1)+{\varphi }_2(x,0)+\frac{p_1}{q}{\varphi }_3(x,0)\right\}\right\}={\xi }_1,$

so $T(x,t)$,$I(x,t)$ and $Z(x,t)$ are bounded.

From the boundedness of $I(x,t)$ and system (4)–(6), $V(x,t)$ satisfies the following system

(10)$\left\{\begin{array}{c}\frac{\partial V}{\partial t}-D\Delta V\le p_2{e}^{-{\mu }_2{\tau }_2}{\xi }_1-d_{3}V,\hfill \\ \frac{\partial V}{\partial n}=0,\hfill \\ V(x,0)={\varphi }_3(x,0)\ge 0.\hfill \end{array}\right.$

Assume $V_1\left(t\right)$ be a solution to the ordinary differential equation

(11)$\left\{\begin{array}{c}\frac{d{V}_1}{dt}=p_2{e}^{-{\mu }_2{\tau }_2}{\xi }_1-d_3{V}_1,\hfill \\ V_1\left(0\right)=ma{x}_{x\in \overline{\Omega }}{\varphi }_3(x,0),\hfill \end{array}\right.$

then

$V_1\left(t\right)\le max\left\{\frac{p_2{e}^{-{\mu }_2{\tau }_2}{\xi }_1}{d_3},ma{x}_{x\in \overline{\Omega }}{\varphi }_3(x,0)\right\},\forall t\in [0,T_{max}).$

It follows from the comparison principle [37] that $V(x,t)\le V_1\left(t\right)$. Therefore

$V(x,t)\le max\left\{\frac{p_2{e}^{-{\mu }_2{\tau }_2}{\xi }_1}{d_3},ma{x}_{x\in \overline{\Omega }}{\varphi }_3(x,0)\right\}={\xi }_2,\forall (x,t)\in \overline{\Omega }\times [0,T_{max}).$

From the above, $T(x,t),I(x,t),V(x,t)$ and $Z(x,t)$ are bounded in $\overline{\Omega }\times [0,T_{max})$. Furthermore, it follows from the standard theory for semilinear parabolic systems [38] that $T_{max}=+\infty$. □

2.2. Existence of Equilibria

It is clear that system (4) always has an infection-free equilibrium

$E_0=(T_0,0,0,0),$

where $T_0=\frac{\lambda }{d}$, corresponding to the maximal level of healthy $C{D}_4^{+}T$ cells. It is the only biologically meaningful equilibrium if

$R_0=\frac{\lambda e^{-{\mu }_1{\tau }_1}({\beta }_1{p}_2{f}^{\prime -{\mu }_2{\tau }_2}+{\beta }_2{d}_3{g}^{\prime }\left(0\right))}{d_1{d}_2{d}_3}<1,$

where $R_0$ is basic reproduction number.

At an equilibrium of model (4), we have

(12)$\left\{\begin{array}{c}\lambda =d_{1}T+{\beta }_{1}Tf\left(V\right)+{\beta }_{2}Tg\left(I\right),\hfill \\ e^{-{\mu }_1{\tau }_1}\left({\beta }_{1}Tf\left(V\right)+{\beta }_{2}Tg\left(I\right)\right)=d_{2}I+p_{1}IZ,\hfill \\ p_2{e}^{-{\mu }_2{\tau }_2}I=d_{3}V,\hfill \\ qIZ=d_{4}Z,\hfill \end{array}\right.$

if $Z=0$, then a short calculation

$\lambda -d_{1}T=\frac{d_2{d}_3{e}^{{\mu }_1{\tau }_1+{\mu }_2{\tau }_2}}{p_2}V,I=\frac{d_3{e}^{{\mu }_2{\tau }_2}}{p_2}V,$

which implies that in order to have $T\ge 0$ and $V>0$ at an equilibrium, then $V\in (0,\frac{\lambda p_2}{d_2{d}_3{e}^{{\mu }_1{\tau }_1+{\mu }_2{\tau }_2}}]$. From the second equation of (12), we have

$T=\frac{d_2{d}_3{e}^{{\mu }_1{\tau }_1+{\mu }_2{\tau }_2}}{p_2\left({\beta }_{1}f\left(V\right)+{\beta }_{2}g\left(\frac{d_3{e}^{{\mu }_2{\tau }_2}}{p_2}V\right)\right)}V,$

then substituting T into the first equation of (12)

$\lambda =\frac{d_1{d}_2{d}_3{e}^{{\mu }_1{\tau }_1+{\mu }_2{\tau }_2}}{p_2\left({\beta }_{1}f\left(V\right)+{\beta }_{2}g\left(\frac{d_3{e}^{{\mu }_2{\tau }_2}}{p_2}V\right)\right)}V+\frac{d_2{d}_3\mu e^{{\mu }_1{\tau }_1+{\mu }_2{\tau }_2}}{p_2}V=H\left(V\right).$

According to $\left(A_2\right)$, for all $V>0$, we have

$\begin{array}{ccc}\hfill H^{\prime }\left(V\right)& =& \frac{d_1{d}_2{d}_3{e}^{{\mu }_1{\tau }_1+{\mu }_2{\tau }_2}\left({\beta }_1(f\left(V\right)-V{f}^{\prime }\left(V\right))+{\beta }_2\left(g\left(\frac{d_3{e}^{{\mu }_2{\tau }_2}}{p_2}V\right)-\frac{d_3{e}^{{\mu }_2{\tau }_2}}{p_2}V{g}^{\prime }\left(\frac{d_3{e}^{{\mu }_2{\tau }_2}}{p_2}V\right)\right)\right)}{p_2{({\beta }_{1}f\left(V\right)+{\beta }_{2}g\left(\frac{d_3{e}^{{\mu }_2{\tau }_2}}{p_2}V\right))}^2}\hfill \\ & +& \frac{d_2{d}_3{e}^{{\mu }_1{\tau }_1+{\mu }_2{\tau }_2}}{p_2}>0,\hfill \end{array}$

further, from $\left(A_1\right)$

$\underset{V\to 0^{+}}{lim}H\left(V\right)=\frac{d_1{d}_2{d}_3{e}^{{\mu }_1{\tau }_1+{\mu }_2{\tau }_2}}{p_2{\beta }_1{f}^{\prime }\left(0\right)+d_3{\beta }_2{e}^{{\mu }_2{\tau }_2}{g}^{\prime }\left(0\right)}=\frac{\lambda }{R_0},$

and

$H\left(\frac{\lambda p_2}{d_2{d}_3{e}^{{\mu }_1{\tau }_1+{\mu }_2{\tau }_2}}\right)=\lambda +\frac{\lambda d_1}{{\beta }_{1}f\left(\frac{\lambda p_2}{d_2{d}_3{e}^{{\mu }_1{\tau }_1+{\mu }_2{\tau }_2}}\right)+{\beta }_{2}g\left(\frac{\lambda }{d_2{e}^{{\mu }_1{\tau }_1}}\right)}>\lambda ,$

this implies that there exists a CTL-inactivated equilibrium $E_1=(T_1,I_1,V_1,0)$ when $R_0>1.$

Define

$R_1=\frac{\lambda q{e}^{-{\mu }_1{\tau }_1}\left({\beta }_{1}f\left(\frac{d_4{p}_2{e}^{-{\mu }_2{\tau }_2}}{d_{3}q}\right)+{\beta }_{2}g\left(\frac{d_4}{q}\right)\right)}{d_2{d}_4\left(d_1+{\beta }_{1}f\left(\frac{d_4{p}_2{e}^{-{\mu }_2{\tau }_2}}{d_{3}q}\right)+{\beta }_{2}g\left(\frac{d_4}{q}\right)\right)},$

which stands for the immune response activation number and determines whether a persistent immune response can be established or not. If $Z\ne 0$, then from (12), we have

$T_2=\frac{\lambda }{d_1+{\beta }_{1}f\left(V_2\right)+{\beta }_{2}g\left(I_2\right)},I_2=\frac{d_4}{q},$

$V_2=\frac{d_4{p}_2{e}^{-{\mu }_2{\tau }_2}}{d_{3}q},Z_2=\frac{d_2}{p_1}(R_1-1),$

then, the infection equilibrium with immunity $E_2=(T_2,I_2,V_2,Z_2)$ exists if $R_1>1.$ From the above, we have the following result.

Lemma 1.

For system (4),

• (1) 

  if $R_0<1$, then there exists a unique infection-free equilibrium $E_0$.

• (2) 

  if $R_1\le 1<R_0$, then there exists a unique infection equilibrium without immunity $E_1$ besides $E_0$.

• (3) 

  if $R_1>1$, then there exists a unique infection equilibrium with immunity $E_2$ besides $E_0$ and $E_1$.

2.3. Global Asymptotic Stability

In this section, we will investigate the global asymptotic stability of the system (4). Assume $\phi \left(u\right)=u-1-lnu$ for $u\in (0,+\infty )$, then $\phi \left(x\right)\ge \phi \left(1\right)=0$.

Theorem 2.

For system (4), if $R_0\le 1$, the infection-free equilibrium $E_0$ is globally asymptotically stable.

Proof. 

Define the Lyapunov function as follows

$\begin{array}{ccc}\hfill L_1& =& {\int }_{\Omega }\{T_0\phi \left(\frac{T}{T_0}\right)+e^{{\mu }_1{\tau }_1}I+\frac{{\beta }_1{p}_2{T}_0{e}^{-{\mu }_2{\tau }_2}{f}^{\prime }\left(0\right)}{R_0{d}_3}{\int }_{t-{\tau }_2}^{t}I\left(s\right)ds+\frac{{\beta }_1{T}_0{f}^{\prime }\left(0\right)}{R_0{d}_3}V\hfill \\ & +& \frac{p_1{e}^{{\mu }_1{\tau }_1}}{q}Z+{\int }_{t-{\tau }_1}^t\left[{\beta }_{1}T\left(s\right)f\left(V\left(s\right)\right)+{\beta }_{2}T\left(s\right)g\left(I\left(s\right)\right)\right]ds\}dx,\hfill \end{array}$

then $L_1\ge 0$, calculating $\frac{d{L}_1}{dt}$ along the solutions of system (4) and using $\lambda =d_1{T}_0$, we have

$\begin{array}{ccc}\frac{d{L}_1}{dt}\hfill & =\hfill & {\int }_{\Omega }\{\left(1-\frac{T_0}{T(x,t)}\right)\left(d_1{T}_0-d_{1}T(x,t)-{\beta }_{1}T(x,t)f\left(V(x,t)\right)\right)\hfill \\ \hfill & -\hfill & {\beta }_{2}T(x,t)g\left(I(x,t)\right)+{\beta }_{1}T(x,t-{\tau }_1)f\left(V(x,t-{\tau }_1)\right)\hfill \\ \hfill & +\hfill & {\beta }_{2}T(x,t-{\tau }_1)g\left(I(x,t-{\tau }_1)\right)-d_2{e}^{{\mu }_1{\tau }_1}I(x,t)-p_1{e}^{{\mu }_1{\tau }_1}I(x,t)Z(x,t)\hfill \\ \hfill & +\hfill & \frac{{\beta }_1{T}_0{f}^{\prime }\left(0\right)}{R_0{d}_3}\left[D\Delta V(x,t)+p_2{e}^{-{\mu }_2{\tau }_2}I(x,t-{\tau }_2)-d_{3}V(x,t)\right]\hfill \\ \hfill & +\hfill & \frac{p_1{e}^{{\mu }_1{\tau }_1}}{q}\left[qI(x,t)Z(x,t)-d_{4}Z(x,t)\right]+{\beta }_{1}T(x,t)f\left(V(x,t)\right)\hfill \\ \hfill & +\hfill & {\beta }_{2}T(x,t)g\left(I(x,t)\right)-{\beta }_{1}T(x,t-{\tau }_1)f\left(V(x,t-{\tau }_1)\right)\hfill \\ \hfill & -\hfill & {\beta }_{2}T(x,t-{\tau }_1)g\left(I(x,t-{\tau }_1)\right)+\frac{{\beta }_1{p}_2{T}_0{e}^{-{\mu }_2{\tau }_2}{f}^{\prime }\left(0\right)}{R_0{d}_3}\left[I(x,t)-I(x,t-{\tau }_2)\right]\}dx,\hfill \\ \hfill & =\hfill & {\int }_{\Omega }\{d_1{T}_0\left(2-\frac{T_0}{T(x,t)}-\frac{T(x,t)}{T_0}\right)+\left(\frac{T_0}{T(x,t)}-1\right)[{\beta }_{1}T(x,t)f\left(V(x,t)\right)\hfill \\ \hfill & +\hfill & {\beta }_{2}T(x,t)g\left(I(x,t)\right)]-d_2{e}^{{\mu }_1{\tau }_1}I(x,t)+\frac{{\beta }_1{T}_0{f}^{\prime }\left(0\right)}{R_0{d}_3}D\Delta V(x,t)\hfill \\ \hfill & +\hfill & \frac{{\beta }_1{T}_0{f}^{\prime }\left(0\right)p_2{e}^{-{\mu }_2{\tau }_2}}{R_0{d}_3}I(x,t-{\tau }_2)-\frac{{\beta }_1{T}_0{f}^{\prime }\left(0\right)}{R_0}V(x,t)-\frac{p_1{d}_4{e}^{{\mu }_1{\tau }_1}}{q}Z(x,t)\hfill \\ \hfill & +\hfill & {\beta }_{1}T(x,t)f\left(V(x,t)\right)+{\beta }_{2}T(x,t)g\left(I(x,t)\right)+\frac{{\beta }_1{T}_0{f}^{\prime }\left(0\right)p_2{e}^{-{\mu }_2{\tau }_2}}{R_0{d}_3}I(x,t)\hfill \\ \hfill & -\hfill & \frac{{\beta }_1{T}_0{f}^{\prime }\left(0\right)p_2{e}^{-{\mu }_2{\tau }_2}}{R_0{d}_3}I(x,t-{\tau }_2)\}dx,\hfill \\ \hfill & \hfill \end{array}$

from ${\int }_{\Omega }\Delta V(x,t)dx=0$ and condition $\left(A_2\right)$, we obtain

$d_2{e}^{{\mu }_1{\tau }_1}-\frac{{\beta }_1{p}_2{T}_0{f}^{\prime -{\mu }_2{\tau }_2}}{R_0{d}_3}=\frac{{\beta }_2{T}_0{g}^{\prime }\left(0\right)}{R_0},$

therefore

$\begin{array}{ccc}\frac{d{L}_1}{dt}\hfill & =\hfill & {\int }_{\Omega }\{d_1{T}_0\left[2-\frac{T_0}{T(x,t)}-\frac{T(x,t)}{T_0}\right]+{\beta }_1{T}_{0}f\left(V(x,t)\right)+{\beta }_2{T}_{0}g\left(I(x,t)\right)\hfill \\ \hfill & -\hfill & d_2{e}^{{\mu }_1{\tau }_1}I(x,t)-\frac{{\beta }_1{T}_0{f}^{\prime }\left(0\right)}{R_0}V-\frac{p_1{d}_4{e}^{{\mu }_1{\tau }_1}}{q}Z(x,t)\hfill \\ \hfill & +\hfill & \frac{{\beta }_1{T}_0{f}^{\prime }\left(0\right)p_2{e}^{-{\mu }_2{\tau }_2}}{R_0{d}_3}I(x,t)\}dx.\hfill \\ \hfill & \le \hfill & {\int }_{\Omega }\{d_1{T}_0\left[2-\frac{T_0}{T(x,t)}-\frac{T(x,t)}{T_0}\right]+\frac{{\beta }_1{T}_0{f}^{\prime }\left(0\right)}{R_0}V(x,t)(R_0-1)\hfill \\ \hfill & +\hfill & \frac{{\beta }_2{T}_0{g}^{\prime }\left(0\right)}{R_0}I(x,t)(R_0-1)-\frac{d_4{p}_1{e}^{{\mu }_1{\tau }_1}}{q}Z(x,t)\}dx.\hfill \end{array}$

It is follows from $R_0\le 1$ that $\frac{d{L}_1}{dt}\le 0$. Furthermore, the largest invariant set of $\{\frac{d{L}_1}{dt}=0\}$ is the singleton $\left\{E_0\right\}$. Then, the classical LaSalle’s invariance principle implies that $E_0$ is globally asymptotically stable. This completes the proof. □

Theorem 3.

For system (4), if $R_1\le 1<R_0$, the $CTL$-inactivated infection equilibrium $E_1$ is globally asymptotically stable.

Proof. 

Define the Lyapunov function as follows

$\begin{array}{ccc}{L}_2\hfill & =\hfill & {\int }_{\Omega }\{T_1\phi \left(\frac{T}{T_1}\right)+e^{{\mu }_1{\tau }_1}{I}_1\phi \left(\frac{I}{I_1}\right)+\frac{{\beta }_1{T}_{1}f\left(V_1\right)}{p_2{e}^{-{\mu }_2{\tau }_2}{I}_1}{V}_1\phi \left(\frac{V}{V_1}\right)+\frac{p_1{e}^{{\mu }_1{\tau }_1}}{q}Z\hfill \\ \hfill & +\hfill & {\beta }_1{T}_{1}f\left(V_1\right){\int }_{t-{\tau }_1}^t\phi \left(\frac{T\left(\theta \right)f\left(V\left(\theta \right)\right)}{T_{1}f\left(V_1\right)}\right)d\theta +{\beta }_2{T}_{1}g\left(I_1\right){\int }_{t-{\tau }_1}^t\phi \left(\frac{T\left(\theta \right)g\left(I\left(\theta \right)\right)}{T_{1}g\left(I_1\right)}\right)d\theta \hfill \\ \hfill & +\hfill & {\beta }_1{T}_{1}f\left(V_1\right){\int }_{t-{\tau }_2}^t\phi \left(\frac{I\left(\theta \right)}{I_1}\right)d\theta \}dx.\hfill \end{array}$

The Lyapunov derivative along system (4) is

$\begin{array}{ccc}\frac{d{L}_2}{dt}\hfill & =\hfill & {\int }_{\Omega }\left\{\left(1-\frac{T_1}{T(x,t)}\right)\right[\lambda -d_{1}T(x,t)-{\beta }_{1}T(x,t)f\left(V(x,t)\right)-{\beta }_{2}T(x,t)g\left(I(x,t)\right)\hfill \\ \hfill & +\hfill & \left(1-\frac{I_1}{I(x,t)}\right)\left[{\beta }_{1}T(x,t-{\tau }_1)f\left(V(x,t-{\tau }_1)\right)+{\beta }_{2}T(x,t-{\tau }_1)g\left(I(x,t-{\tau }_1)\right)\right]\hfill \\ \hfill & -\hfill & d_2{e}^{{\mu }_1{\tau }_1}I(x,t)-p_1{e}^{{\mu }_1{\tau }_1}I(x,t)Z(x,t)]\hfill \\ \hfill & +\hfill & \frac{{\beta }_1{T}_{1}f\left(V_1\right)}{p_2{e}^{-{\mu }_2{\tau }_2}{I}_1}\left(1-\frac{V_1}{V(x,t)}\right)\left[D\Delta V(x,t)+p_2{e}^{-{\mu }_2{\tau }_2}I(x,t-{\tau }_2)-d_{3}V(x,t)\right]\hfill \\ \hfill & +\hfill & \frac{p_1{e}^{{\mu }_1{\tau }_1}}{q}\left[qI(x,t)Z(x,t)-d_{4}Z(x,t)\right]\hfill \\ \hfill & +\hfill & {\beta }_1{T}_{1}f\left(V_1\right)[\frac{T(x,t)f\left(V(x,t)\right)}{T_{1}f\left(V_1\right)}-\frac{T(x,t-{\tau }_1)f\left(V(x,t-{\tau }_1)\right)}{T_{1}f\left(V_1\right)}\hfill \\ \hfill & +\hfill & ln\frac{T(x,t-{\tau }_1)f\left(V(x,t-{\tau }_1)\right)}{T(x,t)f\left(V(x,t)\right)}]+{\beta }_2{T}_{1}g\left(I_1\right)[\frac{T(x,t)g\left(I(x,t)\right)}{T_{1}g\left(I_1\right)}\hfill \\ \hfill & -\hfill & \frac{T(x,t-{\tau }_1)g\left(I(x,t-{\tau }_1)\right)}{T_{1}g\left(I_1\right)}+ln\frac{T(x,t-{\tau }_1)g\left(I(x,t-{\tau }_1)\right)}{T(x,t)g\left(I(x,t)\right)}]\hfill \\ \hfill & +\hfill & {\beta }_1{T}_{1}f\left(V_1\right)\left[\frac{I(x,t)}{I_1}-\frac{I(x,t-{\tau }_2)}{I_1}+ln\frac{I(x,t-{\tau }_2)}{I_1}\right]\}dx.\hfill \\ \hfill & \hfill \end{array}$

According to the equilibrium conditions of $E_1$, that

$\lambda =d_1{T}_1+{\beta }_1{T}_{1}f\left(V_1\right)+{\beta }_2{T}_{1}g\left(I_1\right),$

${\beta }_1{T}_{1}f\left(V_1\right)+{\beta }_2{T}_{1}g\left(I_1\right)=d_2{e}^{{\mu }_1{\tau }_1}{I}_1,p_2{e}^{-{\mu }_2{\tau }_2}{I}_1=d_3{V}_1,$

also recall ${\int }_{\Omega }\Delta V(x,t)dx=0$ and ${\int }_{\Omega }\frac{\Delta V(x,t)}{V(x,t)}dx={\int }_{\Omega }\frac{{\parallel \nabla V(x,t)\parallel }^2}{V^2(x,t)}dx$, we have

$\begin{array}{ccc}\frac{d{L}_2}{dt}\hfill & =\hfill & {\int }_{\Omega }\{d_1{T}_1\left(1-\frac{T_1}{T(x,t)}\right)\left(1-\frac{T(x,t)}{T_1}\right)+\left(1-\frac{T_1}{T(x,t)}\right)[{\beta }_1{T}_{1}f\left(V_1\right)+{\beta }_2{T}_{1}g\left(I_1\right)\hfill \\ \hfill & -\hfill & {\beta }_{1}T(x,t)f\left(V(x,t)\right)-{\beta }_{2}T(x,t)g\left(I(x,t)\right)]\hfill \\ \hfill & +\hfill & \left(1-\frac{I_1}{I(x,t)}\right)\left[{\beta }_{1}T(x,t-{\tau }_1)f\left(V(x,t-{\tau }_1)\right)+{\beta }_{2}T(x,t-{\tau }_1)g\left(I(x,t-{\tau }_1)\right)\right]\hfill \\ \hfill & -\hfill & \left(1-\frac{I_1}{I(x,t)}\right)\frac{I(x,t)}{I_1}\left({\beta }_1{T}_{1}f\left(V_1\right)+{\beta }_2{T}_{1}g\left(I_1\right)\right)\hfill \\ \hfill & -\hfill & p_1{e}^{{\mu }_1{\tau }_1}I(x,t)Z(x,t)+p_1{e}^{{\mu }_1{\tau }_1}{I}_{1}Z(x,t)+{\beta }_1{T}_{1}f\left(V_1\right)\left(1-\frac{V_1}{V(x,t)}\right)(\frac{I(x,t-{\tau }_2)}{I_1}\hfill \\ \hfill & -\hfill & \frac{V(x,t)}{V_1})+p_1{e}^{{\mu }_1{\tau }_1}I(x,t)Z(x,t)-\frac{p_1{e}^{{\mu }_1{\tau }_1}{d}_4}{q}Z(x,t)\hfill \\ \hfill & +\hfill & {\beta }_1{T}_{1}f\left(V_1\right)\left[\frac{T(x,t)f\left(V(x,t)\right)}{T_{1}f\left(V_1\right)}-\frac{T(x,t-{\tau }_1)f\left(V(x,t-{\tau }_1)\right)}{T_{1}f\left(V_1\right)}+ln\frac{T(x,t-{\tau }_1)f\left(V(x,t-{\tau }_1)\right)}{T(x,t)f\left(V(x,t)\right)}\right]\hfill \\ \hfill & +\hfill & {\beta }_2{T}_{1}g\left(I_1\right)\left[\frac{T(x,t)g\left(I(x,t)\right)}{T_{1}g\left(I_1\right)}-\frac{T(x,t-{\tau }_1)g\left(I(x,t-{\tau }_1)\right)}{T_{1}g\left(I_1\right)}+ln\frac{T(x,t-{\tau }_1)g\left(I(x,t-{\tau }_1)\right)}{T(x,t)g\left(I(x,t)\right)}\right]\hfill \\ \hfill & +\hfill & {\beta }_1{T}_{1}f\left(V_1\right)\left[\frac{I(x,t)}{I_1}-\frac{I(x,t-{\tau }_1)}{I_1}+ln\frac{I(x,t-{\tau }_2)}{I(x,t)}\right]\}dx\hfill \\ \hfill & -\hfill & \frac{{\beta }_1{T}_{1}f\left(V_1\right)D{V}_1}{p_2{I}_1{e}^{-{\mu }_2{\tau }_2}}{\int }_{\Omega }\frac{{\parallel \nabla V(x,t)\parallel }^2}{V^2(x,t)}dx\hfill \\ \hfill & =\hfill & {\int }_{\Omega }\{d_1{T}_1\left(1-\frac{T_1}{T(x,t)}\right)\left(1-\frac{T(x,t)}{T_1}\right)+{\beta }_1{T}_{1}f\left(V_1\right)[3-\frac{T_1}{T(x,t)}-\frac{V_{1}I(x,t-{\tau }_2)}{I_{1}V(x,t)}\hfill \\ \hfill & -\hfill & \frac{T(x,t-{\tau }_1)f\left(V(x,t-{\tau }_1)\right)I_1}{T_{1}f\left(V_1\right)I(x,t)}+\frac{f\left(V\right)}{f\left(V_1\right)}-\frac{V}{V_1}+ln\frac{T(x,t-{\tau }_1)f\left(V(x,t-{\tau }_1)I(x,t-{\tau }_2)\right)}{T(x,t)f\left(V(x,t)\right)I(x,t)}]\hfill \\ \hfill & +\hfill & {\beta }_2{T}_{1}g\left(I_1\right)[2-\frac{T_1}{T(x,t)}-\frac{T(x,t-{\tau }_1)g\left(I(x,t-{\tau }_1)\right)I_1}{T_{1}g\left(I_1\right)I(x,t)}+ln\frac{T(x,t-{\tau }_1)g\left(I(x,t-{\tau }_1)\right)}{T(x,t)g\left(I(x,t)\right)}\hfill \\ \hfill & +\hfill & \frac{g\left(I\right)}{g\left(I_1\right)}-\frac{I(x,t)}{I_1}]+p_1{I}_1{e}^{{\mu }_1{\tau }_1}Z(x,t)-\frac{p_1{e}^{{\mu }_1{\tau }_1}{d}_4}{q}Z(x,t)\}dx\hfill \\ \hfill & -\hfill & \frac{{\beta }_1{T}_{1}f\left(V_1\right)D{V}_1}{p_2{I}_1{e}^{-{\mu }_2{\tau }_2}}{\int }_{\Omega }\frac{{\parallel \nabla V(x,t)\parallel }^2}{V^2(x,t)}dx\hfill \\ \hfill & =\hfill & {\int }_{\Omega }\{d_1{T}_1\left(1-\frac{T_1}{T(x,t)}\right)\left(1-\frac{T(x,t)}{T_1}\right)+{\beta }_1{T}_{1}f\left(V_1\right)[-\phi \left(\frac{T_1}{T(x,t)}\right)\hfill \\ \hfill & -\hfill & \phi \left(\frac{T(x,t-{\tau }_1)f\left(V(x,t-{\tau }_1)\right)I_1}{T_{1}f\left(V_1\right)I(x,t)}\right)-\phi \left(\frac{V_{1}I(x,t-{\tau }_2)}{I_{1}V(x,t)}\right)+\frac{f\left(V\right(x,t\left)\right)}{f\left(V_1\right)}-\frac{V(x,t)}{V_1}+ln\frac{f\left(V_1\right)V(x,t)}{f\left(V(x,t)\right)V_1}]\hfill \\ \hfill & +\hfill & {\beta }_2{T}_{1}g\left(I_1\right)[-\phi \left(\frac{T_1}{T(x,t)}\right)-\phi (\frac{T(x,t-{\tau }_1)g\left(I(x,t-{\tau }_1)\right)I_1}{T_{1}g\left(I_1\right)I(x,t))}+\frac{g\left(I(x,t)\right)}{g\left(I_1\right)}-\frac{I(x,t)}{I_1}\hfill \\ \hfill & +\hfill & ln\frac{g\left(I_1\right)I(x,t)}{g\left(I(x,t)\right)I_1}]+p_1{e}^{{\mu }_1{\tau }_1}(I_1-I_2)Z(x,t)\}dx\hfill \\ \hfill & -\hfill & \frac{{\beta }_1{T}_{1}f\left(V_1\right)D{V}_1}{p_2{I}_1{e}^{-{\mu }_2{\tau }_2}}{\int }_{\Omega }\frac{{\parallel \nabla V(x,t)\parallel }^2}{V^2(x,t)}dx\hfill \\ \hfill & =\hfill & {\int }_{\Omega }\{d_1{T}_1\left(1-\frac{T_1}{T(x,t)}\right)(1-\frac{T(x,t)}{T_1})+{\beta }_1{T}_{1}f\left(V_1\right)[-\phi \left(\frac{T_1}{T(x,t)}\right)-\phi \left(\frac{V_{1}I(x,t-{\tau }_2)}{I_{1}V(x,t)}\right)\hfill \\ \hfill & -\hfill & \phi \left(\frac{T(x,t-{\tau }_1)f\left(V(x,t-{\tau }_1)\right)I_1}{T_{1}f\left(V_1\right)I(x,t)}\right)-\phi \left(\frac{f\left(V_1\right)V(x,t)}{f\left(V(x,t)\right)V_1}\right)+\left(\frac{f\left(V(x,t)\right)}{f\left(V_1\right)}-\frac{V(x,t)}{V_1}\right)\left(1-\frac{f\left(V_1\right)}{f\left(V\right(x,t\left)\right)}\right)]\hfill \\ \hfill & +\hfill & {\beta }_2{T}_{1}g\left(I_1\right)[-\phi \left(\frac{T_1}{T(x,t)}\right)-\phi (\frac{T(x,t-{\tau }_1)g\left(I(x,t-{\tau }_1)\right)I_1}{T_{1}g\left(I_{1}I(x,t)\right)}-\phi \left(\frac{g\left(I_1\right)I(x,t)}{g\left(I(x,t)\right)I_1}\right)\hfill \\ \hfill & +\hfill & \left(\frac{g\left(I(x,t)\right)}{g\left(I_1\right)}-\frac{I(x,t)}{I_1}\right)\left(1-\frac{g\left(I_1\right)}{g\left(I\right(x,t\left)\right)}\right)]+p_1{e}^{{\mu }_1{\tau }_1}(I_1-I_2)Z(x,t)\}dx\hfill \\ \hfill & -\hfill & \frac{{\beta }_1{T}_{1}f\left(V_1\right)D{V}_1}{p_2{I}_1{e}^{-{\mu }_2{\tau }_2}}{\int }_{\Omega }\frac{{\parallel \nabla V(x,t)\parallel }^2}{V^2(x,t)}dx.\hfill \\ \hfill & \hfill \end{array}$

It follows from $\left(A_1\right)$ that

$\left(\frac{f\left(V(x,t)\right)}{f\left(V_1\right)}-\frac{V(x,t)}{V_1}\right)\left(1-\frac{f\left(V_1\right)}{f\left(V(x,t)\right)}\right)\le 0,$

$\left(\frac{g\left(I(x,t)\right)}{g\left(I_1\right)}-\frac{I(x,t)}{I_1}\right)\left(1-\frac{g\left(I_1\right)}{g\left(I(x,t)\right)}\right)\le 0.$

As $\phi \left(u\right)\ge 0$ for $u>0$, similar to [23], $sgn(I_1-I_2)=sgn(R_1-1)$, then $\frac{d{L}_2}{dt}\le 0$, therefore $E_1$ is stable, and $\frac{d{L}_2}{dt}=0$ holds if and only if $T(x,t)=T_1,I(x,t)=I_1,V(x,t)=V_1$ and $Z(x,t)=0$ when $R_1<1$, or $T(x,t)=T_1,I(x,t)=I_1,V(x,t)=V_1$ when $R_1=1$. The largest invariance set of $\{\frac{d{L}_2}{dt}=0\}$ is the singleton $\left\{E_1\right\}$. It follows from the classical LaSalle’s invariance principle that $E_1$ is globally asymptotically stable when $R_1\le 1<R_0$. This completes the proof. □

Theorem 4.

For system (4), if $R_1>1$, the interior equilibrium $E_2$ is globally asymptotically stable.

Proof. 

Define the Lyapunov function as follows

$\begin{array}{ccc}\hfill L_3& =& {\int }_{\Omega }\{T_2\phi \left(\frac{T}{T_2}\right)+e^{{\mu }_1{\tau }_1}{I}_2\phi \left(\frac{I}{I_2}\right)+\frac{{\beta }_1{T}_{2}f\left(V_2\right)}{p_2{I}_2{e}^{-{\mu }_2{\tau }_2}}{V}_2\phi \left(\frac{V}{V_2}\right)\hfill \\ & +& \frac{p_1{e}^{{\mu }_1{\tau }_1}}{q}{Z}_2\phi \left(\frac{Z}{Z_2}\right)+{\beta }_1{T}_{2}f\left(V_2\right){\int }_{t-{\tau }_1}^t\phi \left(\frac{T\left(\theta \right)f\left(V\right(\theta \left)\right)}{T_{2}f\left(V_2\right)}\right)d\theta \hfill \\ & +& {\beta }_2{T}_{2}g\left(I_2\right){\int }_{t-{\tau }_1}^t\phi \left(\frac{T\left(\theta \right)g\left(I\left(\theta \right)\right)}{T_{2}g\left(I_2\right)}\right)d\theta +{\beta }_1{T}_{2}f\left(V_2\right){\int }_{t-{\tau }_2}^t\phi \left(\frac{I\left(\theta \right)}{I_2}\right)d\theta \}dx,\hfill \end{array}$

calculating $\frac{d{L}_3}{dt}$ along the solutions of system (4), we have

$\begin{array}{ccc}\frac{d{L}_3}{dt}\hfill & =\hfill & {\int }_{\Omega }\{\left(1-\frac{T_2}{T(x,t)}\right)\left(\lambda -d_{1}T(x,t)-{\beta }_{1}T(x,t)f\left(V(x,t)\right)-{\beta }_{2}T(x,t)g\left(I(x,t)\right)\right)\hfill \\ \hfill & +\hfill & \left(1-\frac{I_2}{I(x,t)}\right)\left[{\beta }_{1}T(x,t-{\tau }_1)f\left(V(x,t-{\tau }_1)\right)+{\beta }_{2}T(x,t-{\tau }_1)g\left(I(x,t-{\tau }_1)\right)\right]\hfill \\ \hfill & -\hfill & e^{{\mu }_1{\tau }_1}\left(1-\frac{I_2}{I(x,t)}\right)\left[d_{2}I(x,t)-p_{1}I(x,t)Z(x,t)\right]\hfill \\ \hfill & +\hfill & \frac{{\beta }_1{T}_{2}f\left(V_2\right)}{p_2{I}_2{e}^{-{\mu }_2{\tau }_2}}\left(1-\frac{V_2}{V(x,t)}\right)\left[\Delta V(x,t)+p_2{e}^{-{\mu }_2{\tau }_2}I(x,t-{\tau }_2)-d_{3}V(x,t)\right]\hfill \\ \hfill & +\hfill & \frac{p_1{e}^{{\mu }_1{\tau }_1}}{q}\left(1-\frac{Z_2}{Z(x,t)}\right)\left[qI(x,t)Z(x,t)-d_{4}Z(x,t)\right]\hfill \\ \hfill & +\hfill & {\beta }_1{T}_{2}f\left(V_2\right)\left[\frac{T(x,t)f\left(V(x,t)\right)}{T_{2}f\left(V_2\right)}-\frac{T(x,t-{\tau }_1)f\left(V(x,t-{\tau }_1)\right)}{T_{2}f\left(V_2\right)}+ln\frac{T(x,t-{\tau }_1)f\left(V(x,t-{\tau }_1)\right)}{T(x,t)f\left(V(x,t)\right)}\right]\hfill \\ \hfill & +\hfill & {\beta }_2{T}_{2}g\left(I_2\right)\left[\frac{T(x,t)g\left(I(x,t)\right)}{T_{2}g\left(I_2\right)}-\frac{T(x,t-{\tau }_1)g\left(I(x,t-{\tau }_1)\right)}{T_{2}g\left(I_2\right)}+ln\frac{T(x,t-{\tau }_1)g\left(I(x,t-{\tau }_1)\right)}{T(x,t)g\left(I(x,t)\right)}\right]\hfill \\ \hfill & +\hfill & {\beta }_1{T}_{2}f\left(V_2\right)\left[\frac{I(x,t)}{I_2}-\frac{I(x,t-{\tau }_2)}{I_2}+ln\frac{I(x,t-{\tau }_2)}{I(x,t)}\right]\}dx,\hfill \\ \hfill & \hfill \end{array}$

using the equilibrium conditions of $E_2$, then

$\lambda =d_1{T}_2+{\beta }_1{T}_{2}f\left(V_2\right)+{\beta }_2{T}_{2}g\left(I_2\right),p_2{e}^{-{\mu }_2{\tau }_2}{I}_2=d_3{V}_2,$

${\beta }_1{T}_{2}f\left(V_2\right)+{\beta }_2{T}_{2}g\left(I_2\right)=e^{{\mu }_1{\tau }_1}(d_2{I}_2+p_1{I}_2{Z}_2),I_2=\frac{d_4}{q},$

also recall ${\int }_{\Omega }\Delta V(x,t)dx=0$ and ${\int }_{\Omega }\frac{\Delta V(x,t)}{V(x,t)}dx={\int }_{\Omega }\frac{{\parallel \nabla V(x,t)\parallel }^2}{V^2(x,t)}dx$, we have

$\begin{array}{ccc}\frac{d{L}_3}{dt}\hfill & =\hfill & {\int }_{\Omega }\{d_1{T}_2\left(1-\frac{T_2}{T(x,t)}\right)\left(1-\frac{T(x,t)}{T_2}\right)+{\beta }_1{T}_{2}f\left(V_2\right)[3-\frac{T_2}{T(x,t)}-\frac{V_{2}I(x,t-{\tau }_2)}{I_{2}V(x,t)}\hfill \\ \hfill & -\hfill & \frac{T(x,t-{\tau }_1)f\left(V(x,t-{\tau }_1)\right)I_2}{T_{2}f\left(V_2\right)I(x,t)}+\frac{f\left(V(x,t)\right)}{f\left(V_2\right)}-\frac{V(x,t)}{V_2}+ln\frac{T(x,t-{\tau }_1)f\left(V(x,t-{\tau }_1)\right)I(t-{\tau }_2)}{T(x,t)f\left(V\right(x,t\left)\right)I(x,t)}]\hfill \\ \hfill & +\hfill & {\beta }_2{T}_{2}g\left(I_2\right)[2-\frac{T_2}{T(x,t)}-\frac{T(x,t-{\tau }_1)g\left(I(x,t-{\tau }_1)\right)I_2}{T_{2}g\left(I_2\right)I(x,t)}+\frac{g\left(I(x,t)\right)}{g\left(I_2\right)}-\frac{I(x,t)}{I_2}\hfill \\ \hfill & +\hfill & ln\frac{T(x,t-{\tau }_1)g\left(I(x,t-{\tau }_1)\right)}{T(x,t)g\left(I(x,t)\right)}\left]\right\}dx-\frac{{\beta }_1{T}_{2}f\left(V_2\right)D{V}_2}{p_2{I}_2{e}^{-{\mu }_2{\tau }_2}}{\int }_{\Omega }\frac{{\parallel \nabla V(x,t)\parallel }^2}{V^2(x,t)}dx\hfill \\ \hfill & =\hfill & {\int }_{\Omega }\{d_1{T}_2\left(1-\frac{T_2}{T(x,t)}\right)(1-\frac{T(x,t)}{T_2})+{\beta }_1{T}_{2}f\left(V_2\right)[-\phi \left(\frac{T_2}{T(x,t)}\right)-\phi \left(\frac{V_{2}I(x,t-{\tau }_2)}{I_{2}V(x,t)}\right)\hfill \\ \hfill & -\hfill & \phi \left(\frac{T(x,t-{\tau }_1)f\left(V(x,t-{\tau }_1)\right)I_2}{T_{2}f\left(V_2\right)I(x,t)}\right)+\frac{f\left(V(x,t)\right)}{f\left(V_2\right)}-\frac{V(x,t)}{V_2}+ln\frac{f\left(V_2\right)V(x,t)}{V_{2}f\left(V(x,t)\right)}]\hfill \\ \hfill & +\hfill & {\beta }_2{T}_{2}g\left(I_2\right)[-\phi \left(\frac{T_2}{T(x,t)}\right)-\phi \left(\frac{T(x,t-{\tau }_1)g\left(I(x,t-{\tau }_1)\right)I_2}{T_{2}g\left(I_2\right)I(x,t)}\right)\hfill \\ \hfill & +\hfill & \frac{g\left(I(x,t)\right)}{g\left(I_2\right)}-\frac{I(x,t)}{I_2}+ln\frac{g\left(I_2\right)I(x,t)}{I_{2}g\left(I(x,t)\right)}\left]\right\}dx-\frac{{\beta }_1{T}_{2}f\left(V_2\right)D{V}_2}{p_2{I}_2{e}^{-{\mu }_2{\tau }_2}}{\int }_{\Omega }\frac{{\parallel \nabla V(x,t)\parallel }^2}{V^2(x,t)}dx\hfill \\ \hfill & =\hfill & {\int }_{\Omega }\{d_1{T}_2\left(1-\frac{T_2}{T(x,t)}\right)\left(1-\frac{T(x,t)}{T_2}\right)+{\beta }_1{T}_{2}f\left(V_2\right)[-\phi \left(\frac{T_2}{T(x,t)}\right)-\phi \left(\frac{V_{2}I(x,t-{\tau }_2)}{I_{2}V(x,t)}\right)-\hfill \\ \hfill & \phi \hfill & \left(\frac{T(x,t-{\tau }_1)f\left(V(x,t-{\tau }_1)\right)I_2}{T_{2}f\left(V_2\right)I(x,t)}\right)-\phi \left(\frac{f\left(V_2\right)V(x,t)}{V_{2}f\left(V(x,t)\right)}\right)+\left(\frac{f\left(V\right(x,t\left)\right)}{f\left(V_2\right)}-\frac{V(x,t)}{V_2}\right)\left(1-\frac{f\left(V_2\right)}{f\left(V(x,t)\right)}\right)]\hfill \\ \hfill & +\hfill & {\beta }_2{T}_{2}g\left(I_2\right)[-\phi \left(\frac{T_2}{T(x,t)}\right)-\phi \left(\frac{T(x,t-{\tau }_1)g\left(I(x,t-{\tau }_1)\right)I_2}{T_{2}g\left(I_2\right)I(x,t)}\right)-\phi \left(\frac{g\left(I_2\right)I(x,t)}{I_{2}g\left(I(x,t)\right)}\right)\hfill \\ \hfill & +\hfill & \left(\frac{g\left(I(x,t)\right)}{g\left(I_2\right)}-\frac{I(x,t)}{I_2}\right)\left(1-\frac{g\left(I_2\right)}{g\left(I\right(x,t\left)\right)}\right)\left]\right\}dx-\frac{{\beta }_1{T}_{2}f\left(V_2\right)D{V}_2}{k{I}_2}{\int }_{\Omega }\frac{{\parallel \nabla V(x,t)\parallel }^2}{V^2(x,t)}dx,\hfill \\ \hfill & \hfill \end{array}$

from $\left(A_1\right)$, it is easy to see that

$\left(\frac{f\left(V\right(x,t\left)\right)}{f\left(V_1\right)}-\frac{V(x,t)}{V_1}\right)\left(1-\frac{f\left(V_1\right)}{f\left(V\right(x,t\left)\right)}\right)\le 0,$

$\left(\frac{g\left(I\right(x,t\left)\right)}{g\left(I_1\right)}-\frac{I(x,t)}{I_1}\right)\left(1-\frac{g\left(I_1\right)}{g\left(I\right(x,t\left)\right)}\right)\le 0.$

As $\phi \left(u\right)\ge 0$ for $u>0$, then $\frac{d{L}_3}{dt}\le 0$. The largest invariant set of $\{\frac{d{L}_3}{dt}=0\}$ is the single point $\left\{E_2\right\}$, similar to the proof of Theorem 3, $E_2$ is globally asymptotically stable. This completes the proof. □

3. Dynamical Behaviors of Discrete System

In preceding section, by introducing Lyapunov functions, we have shown by using continuous Lyapunov functionals that the global asymptotic stability of the equilibria of the continuous system (4) is completely determined by the basic reproduction number. $R_0$ and $R_1$ act as threshold parameter for the existence and stability of equilibria. This arises a natural question that whether the global asymptotic stability of the equilibria of the discrete system (7) can be preserved. In this section, we will discuss this problem.

Obviously, the discrete system (7) has the same equilibria as system (4). Similarly, $E_0=(T_0,0,0,0)$ is the infection-free equilibrium, $E_1=(T_1,I_1,V_1,0)$ stands for the $CTL$-inactivated equilibrium and $E_2=(T_2,I_2,V_2,Z_2)$ is the $CTL$-activated equilibrium.

Rewriting the discrete system (7) yields

(13)$\left\{\begin{array}{c}{T}_{n+1}^m=\frac{\lambda \Delta t+T_n^m}{1+\Delta t(d_1+{\beta }_{1}f(V_n^m+{\beta }_{2}g\left(I_n^m\right))},\hfill \\ I_{n+1}^m=\frac{I_n^m+\Delta t{e}^{-{\mu }_1{\tau }_1}({\beta }_1{T}_{n-m_1+1}^{m}f\left(V_{n-m_1}^m\right)+{\beta }_2{T}_{n-m_1+1}^{m}g\left(I_{n-m_1}^m\right)}{1+\Delta t(d_2+p_1{Z}_n^m)},\hfill \\ A{V}_{n+1}=V_n+\Delta t{p}_2{e}^{-{\mu }_2{\tau }_2}{I}_{n-m_2+1},\hfill \\ Z_{n+1}^m=\frac{(1+\Delta tq{I}_{n+1}^m)}{1+\Delta t{d}_4}{Z}_n^m,\hfill \end{array}\right.,$

where the square matrix A of dimension $(N+1)\times (N+1)$ is given by

$\left(\begin{array}{ccccccc}{c}_1& c_2& 0& \cdots & 0& 0& 0\\ c_2& c_3& c_2& \cdots & 0& 0& 0\\ 0& c_2& c_3& \cdots & 0& 0& 0\\ ⋮& ⋮& ⋮& \ddots & ⋮& ⋮& ⋮\\ 0& 0& 0& \cdots & c_3& c_2& 0\\ 0& 0& 0& \cdots & c_2& c_3& c_2\\ 0& 0& 0& \cdots & 0& c_2& c_1\end{array}\right)$

with $c_1=1+D\Delta t/{(\Delta x)}^2+d_3\Delta t,c_2=-D\Delta t/{(\Delta x)}^2,c_3=1+2D\Delta t/{(\Delta x)}^2+d_3\Delta t$. It is clear that A is strictly diagonally dominant matrix, therefore A is non-singular. From the third equation of the above system, we have

$V_{n+1}=A^{-1}(V_n+\Delta t{p}_2{e}^{-{\mu }_2{\tau }_2}{I}_{n-m_2+1}).$

Theorem 5.

For any $\Delta t>0,\Delta x>0$, the solutions of the system (7) remain nonnegative and bounded for all $n\in \mathbb{N}$.

Proof. 

Since all parameters in (7) are positive, then using the induction, it is easy to deduce from (13) that all solutions of system (7) remain nonnegative provided that the initial value are nonnegative, for all $n\in \aleph$.

Next, we establish the boundedness of solutions. Define a sequence $G_n$ as follows

$G_n^m=T_n^m+I_n^m+\frac{p_1}{q}{Z}_n^m+\Delta t\sum _{j=n-m_1}^{n-1}\left[{\beta }_{1}f\left(V_j^m\right)+{\beta }_{2}g\left(I_j^m\right)\right]T_{j+1}^m{e}^{-\Delta t{\mu }_1(n-j)},$

then

$\begin{array}{ccc}\hfill G_{n+1}^m-G_n^m& =& \Delta t\left(\lambda -d_1{T}_{n+1}^m-{\beta }_1{T}_{n+1}^{m}f\left(V_n^m\right)-{\beta }_2{T}_{n+1}^{m}g\left(I_n^m\right)\right)\hfill \\ & +& \Delta t[e^{-{\mu }_1{\tau }_1}\left({\beta }_1{T}_{n-m_1+1}^{m}f\left(V_{n-m_1}^m\right)+{\beta }_2{T}_{n-m_1+1}^{m}g\left(I_{n-m_1}^m\right)\right)\hfill \\ & -& d_2{I}_{n+1}^m-p_1{I}_{n+1}^m{Z}_n^m]+\frac{p_1}{q}\Delta t(q{I}_{n+1}^m{Z}_n^m-d_4{Z}_{n+1}^m)\hfill \\ & +& \Delta t\sum _{j=n-m_1+1}^n\left[{\beta }_{1}f\left(V_j^m\right)+{\beta }_{2}g\left(I_j^m\right)\right]T_{j+1}^m{e}^{-\Delta t{\mu }_1(n-j+1)}\hfill \\ & -& \Delta t\sum _{j=n-m_1}^{n-1}\left[{\beta }_{1}f\left(V_j^m\right)+{\beta }_{2}g\left(I_j^m\right)\right]T_{j+1}^m{e}^{-\Delta t{\mu }_1(n-j)}\hfill \\ & =& \Delta t\{\lambda -d_1{T}_{n+1}^m-d_2{I}_{n+1}^m-\frac{p_1{d}_4}{q}{Z}_{n+1}^m\hfill \\ & +& \left(1-e^{\Delta t{\mu }_1}\right)\sum _{j=n-m_1+1}^n\left[{\beta }_{1}f\left(V_j^m\right)+{\beta }_{2}g\left(I_j^m\right)\right]T_{j+1}^m{e}^{-\Delta t{\mu }_1(n-j+1)}\}\hfill \\ & \le & \Delta t(\lambda -\xi G_{n+1}^m),\hfill \end{array}$

where $\xi =min\{d_1,d_2,d_4,\frac{e^{\Delta t{\mu }_1}-1}{\Delta t}\}$, then we have

$G_{n+1}^m\le \frac{1}{1+\Delta t\xi }{G}_n^m+\frac{\Delta t\lambda }{1+\Delta t\xi },$

it follows from the induction that

$G_n^m\le {\left(\frac{1}{1+\Delta t\xi }\right)}^n{G}_0^m+\frac{\lambda }{\xi }\left[1-{\left(\frac{1}{1+\Delta t\xi }\right)}^n\right],$

therefore

$\underset{n\to \infty }{lim\; sup}{G}_n^m\le \frac{\lambda }{\xi },forallm\in \{0,1,\cdots ,N\},$

this implies that $\left\{G_n\right\}$ is bounded, then $\left\{T_n\right\}$, $\left\{I_n\right\}$ and $\left\{Z_n\right\}$ are bounded.

From the third equation of system (7)

$\sum _{m=0}^N{V}_{n+1}^m=\frac{1}{1+d_3\Delta t}\left(\sum _{m=0}^N{V}_n^m+\Delta t{p}_2{e}^{-{\mu }_2{\tau }_2}\sum _{m=0}^N{I}_{n-m_2+1}^m\right),$

since $\left\{I_n\right\}$ is bounded, then there exists $\eta >0$ such that $I_n^m\le \eta$ for all $n\in \{-m_2,-m_2+1,\cdots ,0,1,\cdots \},m\in \{0,1,\cdots ,\aleph \}$, then

$\sum _{m=0}^N{V}_{n+1}^m\le \frac{1}{1+d_3\Delta t}\left[\sum _{m=0}^N{V}_n^m+\Delta t{p}_2{e}^{-{\mu }_2{\tau }_2}\eta (N+1)\right],$

by induction, we have

$\begin{array}{ccc}\hfill \sum _{m=0}^N{V}_n^m& \le & \frac{1}{{(1+d_3\Delta t)}^n}\sum _{m=0}^N{V}_0^m+\frac{p_2{e}^{-{\mu }_2{\tau }_2}\eta (N+1)}{d_3}\left[1-\frac{1}{{(1+d_3\Delta t)}^n}\right]\hfill \\ & \le & \sum _{m=0}^N{V}_0^m+\frac{p_2{e}^{-{\mu }_2{\tau }_2}\eta (N+1)}{d_3},\hfill \end{array}$

therefore $\left\{V_n\right\}$ is bounded. This completes the proof. □

Global Stability

In this section, we will study the global stability of the equilibria of system (7).

Theorem 6.

For system (7), if $R_0\le 1$, the infection-free equilibrium $E_0$ is globally asymptotically stable.

Proof. 

Define the discrete Lyapunov function as follows

$\begin{array}{ccc}\hfill W_n& =& \sum _{m=0}^N\left\{\frac{1}{\Delta t}\right[T_0\phi \left(\frac{T_n^m}{T_0}\right)+e^{{\mu }_1{\tau }_1}\left(1+\Delta t\frac{{\beta }_2{T}_0{g}^{\prime }\left(0\right)}{R_0{e}^{{\mu }_1{\tau }_1}}\right)I_n^m+\frac{{\beta }_1{T}_0{f}^{\prime }\left(0\right)}{d_3{R}_0}\left(1+\Delta t{d}_3\right)V_n^m\hfill \\ & +& \frac{p_1{e}^{{\mu }_1{\tau }_1}}{q}\left(1+\Delta t{d}_4\right)Z_n^m]+\sum _{j=n-m_1}^{n-1}\left({\beta }_{1}f\left(V_j^m\right)+{\beta }_{2}g\left(I_j^m\right)\right)T_{j+1}^m\hfill \\ & +& \frac{{\beta }_1{p}_2{T}_0{f}^{\prime -{\mu }_2{\tau }_2}}{R_0{d}_3}\sum _{j=n-m_2}^{n-1}{I}_{j+1}^m\}.\hfill \end{array}$

It follows from $u-1\ge lnu$ for all $u>0$, that $W_n\ge 0$ for all $n\in \aleph$. Then, along the trajectory of (7)

$\begin{array}{ccc}\hfill W_{n+1}-W_n& =& \sum _{m=0}^N\left\{\frac{1}{\Delta t}\right[T_{n+1}^m-T_n^m+T_{0}ln\frac{T_n^m}{T_{n+1}^m}+e^{{\mu }_1{\tau }_1}\left(1+\Delta t\frac{{\beta }_2{T}_0{g}^{\prime }\left(0\right)}{R_0{e}^{{\mu }_1{\tau }_1}}\right)(I_{n+1}^m-I_n^m)\hfill \\ & +& \frac{{\beta }_1{T}_0{f}^{\prime }\left(0\right)}{d_3{R}_0}\left(1+\Delta t{d}_3\right)\left(V_{n+1}^m-V_n^m\right)+\frac{p_1{e}^{{\mu }_1{\tau }_1}}{q}\left(1+\Delta t{d}_4\right)\left(Z_{n+1}^m-Z_n^m\right)]\hfill \\ & +& \sum _{j=n-m_1+1}^n\left({\beta }_{1}f\left(V_j^m\right)+{\beta }_{2}g\left(I_j^m\right)\right)T_{j+1}^m\hfill \\ & -& \sum _{j=n-m_1}^{n-1}\left({\beta }_{1}f\left(V_j^m\right)+{\beta }_{2}g\left(I_j^m\right)\right)T_{j+1}^m\hfill \\ & +& \frac{{\beta }_1{p}_2{T}_0{f}^{\prime -{\mu }_2{\tau }_2}}{R_0{d}_3}\left[\sum _{j=n-m_2+1}^n{I}_{j+1}^m-\sum _{j=n-m_2}^{n-1}{I}_{j+1}^m\right]\},\hfill \end{array}$

using the equilibrium condition of $E_0$, we have

$\begin{array}{ccc}\hfill W_{n+1}-W_n& \le & \sum _{m=0}^N\{\left(1-\frac{T_0}{T_{n+1}^m}\right)\left(d_1{T}_0-d_1{T}_{n+1}^m-\left({\beta }_{1}f\left(V_n^m\right)+{\beta }_{2}g\left(I_n^m\right)\right)T_{n+1}^m\right)\hfill \\ & +& \left(1+\Delta t\frac{{\beta }_2{T}_0{g}^{\prime }\left(0\right)}{R_0{e}^{{\mu }_1{\tau }_1}}\right)\left({\beta }_1{T}_{n-m_1+1}^{m}f\left(V_{n-m_1}^m\right)+{\beta }_2{T}_{n-m_1+1}^{m}g\left(I_{n-m_1}^m\right)\right)\hfill \\ & +& e^{{\mu }_1{\tau }_1}\left(1+\Delta t\frac{{\beta }_2{T}_0{g}^{\prime }\left(0\right)}{R_0{e}^{{\mu }_1{\tau }_1}}\right)\left(-d_2{I}_{n+1}^m-p_1{I}_{n+1}^m{Z}_n^m\right)\hfill \\ & +& \frac{{\beta }_1{T}_0{f}^{\prime }\left(0\right)}{d_3{R}_0}\left(1+\Delta t{d}_3\right)D\frac{V_{n+1}^{m+1}-2V_{n+1}^m+V_{n+1}^{m-1}}{{(\Delta x)}^2}\hfill \\ & +& \frac{{\beta }_1{T}_0{f}^{\prime }\left(0\right)}{d_3{R}_0}\left(1+\Delta t{d}_3\right)\left(p_2{e}^{{\mu }_2{\tau }_2}{I}_{n-m_2+1}^m-d_3{V}_{n+1}^m\right)\hfill \\ & +& \frac{p_1{e}^{{\mu }_1{\tau }_1}}{q}\left(1+\Delta t{d}_4\right)\left(q{I}_{n+1}^m{Z}_n^m-d_4{Z}_{n+1}^m\right)+\left({\beta }_{1}f\left(V_n^m\right)+{\beta }_{2}g\left(I_n^m\right)\right)T_{n+1}^m\hfill \\ & -& \left({\beta }_{1}f\left(V_{n-m_1}^m\right)+{\beta }_{2}g\left(I_{n-m_1}^m\right)\right)T_{n-m_1+1}^m\hfill \\ & -& \frac{{\beta }_1{p}_2{T}_0{f}^{\prime -{\mu }_2{\tau }_2}}{R_0{d}_3}\left(I_{n+1}^m-I_{n-m_2+1}^m\right)\}\hfill \\ & =& \sum _{m=0}^N\{d_1{T}_0\left(2-\frac{T_0}{T_{n+1}^m}-\frac{T_{n+1}^m}{T_0}\right)+{\beta }_1{T}_{0}f\left(V_n^m\right)+{\beta }_2{T}_{0}g\left(I_n^m\right)\hfill \\ & -& \frac{{\beta }_2{T}_0{g}^{\prime }\left(0\right)}{R_0}{I}_n^m-\frac{{\beta }_1{T}_0{f}^{\prime }\left(0\right)}{R_0}{V}_n^m-\frac{d_4{p}_1{e}^{{\mu }_1{\tau }_1}}{q}{Z}_n^m\}\hfill \\ & +& \frac{{\beta }_1{T}_0{f}^{\prime }\left(0\right)D}{d_3{R}_0{(\Delta x)}^2}\left(V_{n+1}^{N+1}-V_{n+1}^N+V_{n+1}^{-1}-V_{n+1}^0\right)\hfill \\ & \le & \sum _{m=0}^N\{d_1{T}_0\left(2-\frac{T_0}{T_{n+1}^m}-\frac{T_{n+1}^m}{T_0}\right)+\frac{{\beta }_1{T}_0{f}^{\prime }\left(0\right)}{R_0}{V}_n^m\left(R_0-1\right)\hfill \\ & +& \frac{{\beta }_2{g}^{\prime }\left(0\right)T_0}{R_0}{I}_n^m\left(R_0-1\right)-\frac{d_4{p}_1{e}^{{\mu }_1{\tau }_1}}{q}{Z}_n^m\},\hfill \end{array}$

the last inequality is deduced from condition $\left(A_2\right)$, if $R_0\le 1$, then $W_{n+1}-W_n\le 0$, for all $n\in \aleph$, therefore, $\left\{W_n\right\}$ is monotone decreasing sequence. It follows from $W_n\ge 0$ that ${lim}_{n\to \infty }{W}_n\ge 0,$ then

$\underset{n\to \infty }{lim}(W_{n+1}-W_n)=0,$

therefore

$\left(1\right)$ If $R_0<1$, then ${lim}_{n\to \infty }(W_{n+1}-W_n)=0$ implies that ${lim}_{n\to \infty }{T}_n^m=T_0$, ${lim}_{n\to \infty }{V}_n^m=0$, ${lim}_{n\to \infty }{Z}_n^m=0$, ${lim}_{n\to \infty }{I}_n^m=0$.

$\left(2\right)$ If $R_0=1$, then ${lim}_{n\to \infty }(W_{n+1}-W_n)=0$ implies that ${lim}_{n\to \infty }{T}_n^m=T_0$, ${lim}_{n\to \infty }{Z}_n^m=0$, from system (7), we obtain ${lim}_{n\to \infty }{I}_n^m=0$, ${lim}_{n\to \infty }{V}_n^m=0$.

Hence, $E_0$ is globally asymptotically stable when $R_0\le 1$. This completes the proof. □

Theorem 7.

For system (7), if $R_1<1<R_0$, the $CTL$-inactivated infection equilibrium $E_1$ of is globally asymptotically stable.

Proof. 

Define the discrete Lyapunov function as follows

$\begin{array}{ccc}\hfill {\tilde{W}}_n& =& \sum _{m=0}^N\{\frac{1}{\Delta t}\left[T_1\phi \left(\frac{T_n^m}{T_1}\right)+e^{{\mu }_1{\tau }_1}{I}_1\phi \left(\frac{I_n^m}{I_1}\right)+\frac{{\beta }_1{T}_{1}f\left(V_1\right)}{p_2{e}^{-{\mu }_2{\tau }_2}{I}_1}{V}_1\phi \left(\frac{V_n^m}{V_1}\right)+\frac{p_1{e}^{{\mu }_1{\tau }_1}}{q}{Z}_n^m\right]\hfill \\ & +& {\beta }_1{T}_{1}f\left(V_1\right)\sum _{j=n-m_1}^{n-1}\phi \left(\frac{T_{j+1}^{m}f\left(V_j^m\right)}{T_{1}f\left(V_1\right)}\right)+{\beta }_2{T}_{1}g\left(I_1\right)\sum _{j=n-m_1}^{n-1}\phi \left(\frac{T_{j+1}^{m}g\left(I_j^m\right)}{T_{1}g\left(I_1\right)}\right)\hfill \\ & +& {\beta }_1{T}_{1}f\left(V_1\right)\sum _{j=n-m_2}^{n-1}\phi \left(\frac{I_{j+1}^m}{I_1}\right)+{\beta }_1{T}_{1}f\left(V_1\right)\phi \left(\frac{f\left(V_n^m\right)}{f\left(V_1\right)}\right)+{\beta }_2{T}_{1}g\left(I_1\right)\phi \left(\frac{g\left(I_n^m\right)}{g\left(I_1\right)}\right)\}.\hfill \end{array}$

Since $u-1\ge lnu$ for all $u>0$, then ${\tilde{W}}_n\ge 0$ for all $n\in \aleph$. The Lyapunov derivative along (7) is

$\begin{array}{ccc}\hfill {\tilde{W}}_{n+1}-{\tilde{W}}_n& =& \sum _{m=0}^N\left\{\frac{1}{\Delta t}\right[T_{n+1}^m-T_n^m+T_{1}ln\frac{T_n^m}{T_{n+1}^m}+e^{{\mu }_1{\tau }_1}\left(I_{n+1}^m-I_n^m+I_{1}ln\frac{I_n^m}{I_{n+1}^m}\right)\hfill \\ & +& \frac{{\beta }_1{T}_1{V}_1}{p_2{e}^{-{\mu }_2{\tau }_2}{I}_1}\left(V_{n+1}^m-V_n^m+V_{1}ln\frac{V_n^m}{V_{n+1}^m}\right)+\frac{p_1{e}^{{\mu }_1{\tau }_1}}{q}\left(Z_{n+1}^m-Z_n^m\right)]\hfill \\ & +& {\beta }_1{T}_{1}f\left(V_1\right)\left[\sum _{j=n-m_1+1}^n\phi \left(\frac{T_{j+1}^{m}f\left(V_j^m\right)}{T_{1}f\left(V_1\right)}\right)-\sum _{j=n-m_1}^{n-1}\phi \left(\frac{T_{j+1}^{m}f\left(V_j^m\right)}{T_{1}f\left(V_1\right)}\right)\right]\hfill \\ & +& {\beta }_2{T}_{1}g\left(I_1\right)\left[\sum _{j=n-m_1+1}^n\phi \left(\frac{T_{j+1}^{m}g\left(I_j^m\right)}{T_{1}g\left(I_1\right)}\right)-\sum _{j=n-m_1}^{n-1}\phi \left(\frac{T_{j+1}^{m}g\left(I_j^m\right)}{T_{1}g\left(I_1\right)}\right)\right]\hfill \\ & +& {\beta }_1{T}_{1}f\left(V_1\right)\left[\sum _{j=n-m_2+1}^n\phi \left(\frac{I_{j+1}^m}{I_1}\right)-\sum _{j=n-m_2}^{n-1}\phi \left(\frac{I_{j+1}^m}{I_1}\right)\right]\hfill \\ & +& {\beta }_1{T}_{1}f\left(V_1\right)\left(\frac{f\left(V_{n+1}^m\right)}{f\left(V_1\right)}-\frac{f\left(V_n^m\right)}{f\left(V_1\right)}+ln\frac{f\left(V_n^m\right)}{f\left(V_{n+1}^m\right)}\right)\hfill \\ & +& {\beta }_2{T}_{1}g\left(I_1\right)\left(\frac{g\left(I_{n+1}^m\right)}{g\left(I_1\right)}-\frac{g\left(I_n^m\right)}{g\left(I_1\right)}+ln\frac{g\left(I_n^m\right)}{g\left(I_{n+1}^m\right)}\right)\}\hfill \\ & \le & \sum _{m=0}^N\{\frac{1}{\Delta t}[\left(1-\frac{T_1}{T_{n+1}^m}\right)\left(T_{n+1}^m-T_n^m\right)+e^{{\mu }_1{\tau }_1}\left(1-\frac{I_1}{I_{n+1}^m}\right)\left(I_{n+1}^m-I_n^m\right)\hfill \\ & +& \frac{{\beta }_1{T}_{1}f\left(V_1\right)}{p_2{e}^{-{\mu }_2{\tau }_2}{I}_1}\left(1-\frac{V_1}{V_{n+1}^m}\right)\left(V_{n+1}^m-V_n^m\right)+\frac{p_1{e}^{{\mu }_1{\tau }_1}}{q}(Z_{n+1}^m-Z_n^m)]\hfill \\ & +& {\beta }_1{T}_{1}f\left(V_1\right)\left[\phi \left(\frac{T_{n+1}^{m}f\left(V_n^m\right)}{T_{1}f\left(V_1\right)}\right)-\phi \left(\frac{T_{n-m_1+1}^{m}f\left(V_{n-m_1}^m\right)}{T_{1}f\left(V_1\right)}\right)\right]\hfill \\ & +& {\beta }_2{T}_{1}g\left(I_1\right)\left[\phi \left(\frac{T_{n+1}^{m}g\left(I_n^m\right)}{T_{1}g\left(I_1\right)}\right)-\phi \left(\frac{T_{n-m_1+1}^{m}g\left(I_{n-m_1}^m\right)}{T_{1}g\left(I_1\right)}\right)\right]\hfill \\ & +& {\beta }_1{T}_{1}f\left(V_1\right)\left[\phi \left(\frac{I_{n+1}^m}{I_1}\right)-\phi \left(\frac{I_{n-m_2+1}^m}{I_1}\right)\right]\hfill \\ & +& {\beta }_1{T}_{1}f\left(V_1\right)\left[\frac{f\left(V_{n+1}^m\right)}{f\left(V_1\right)}-\frac{f\left(V_n^m\right)}{f\left(V_1\right)}+ln\frac{f\left(V_n^m\right)}{f\left(V_{n+1}^m\right)}\right]\hfill \\ & +& {\beta }_2{T}_{1}g\left(I_1\right)\left[\frac{g\left(I_{n+1}^m\right)}{g\left(I_1\right)}-\frac{g\left(I_n^m\right)}{g\left(I_1\right)}+ln\frac{g\left(I_n^m\right)}{g\left(I_{n+1}^m\right)}\right]\}.\hfill \end{array}$

As $E_1$ satisfies

$\lambda =d_1{T}_1+{\beta }_1{T}_{1}f\left(V_1\right)+{\beta }_2{T}_{1}g\left(I_1\right),$

${\beta }_1{T}_{1}f\left(V_1\right)+{\beta }_2{T}_{1}g\left(I_1\right)=e^{-{\mu }_1{\tau }_1}{d}_2{I}_1,p_2{e}^{-{\mu }_1{\tau }_1}{I}_1=d_3{V}_1,$

then

$\begin{array}{ccc}\hfill {\tilde{W}}_{n+1}-{\tilde{W}}_n& \le & \sum _{m=0}^N\{d_1{T}_1\left(1-\frac{T_1}{T_{n+1}^m}\right)\left(1-\frac{T_{n+1}^m}{T_1}\right)+\left(1-\frac{T_1}{T_{n+1}^m}\right)\left({\beta }_1{T}_{1}f\left(V_1\right)+{\beta }_2{T}_{1}g\left(I_1\right)\right)\hfill \\ & -& \left(1-\frac{T_1}{T_{n+1}^m}\right)\left({\beta }_1{T}_{n+1}^{m}f\left(V_n^m\right)+{\beta }_2{T}_{n+1}^{m}g\left(I_n^m\right)\right)\hfill \\ & +& \left(1-\frac{I_1}{I_{n+1}^m}\right)\left({\beta }_1{T}_{n-m_1+1}^{m}f\left(V_{n-m_1}^m\right)+{\beta }_2{T}_{n-m_1+1}^{m}g\left(I_{n-m_1}^m\right)\right)\hfill \\ & -& {\beta }_1{T}_{1}f\left(V_1\right)\frac{I_{n+1}^m}{I_1}\left(1-\frac{I_1}{I_{n+1}^m}\right)-{\beta }_2{T}_{1}g\left(I_1\right)\frac{I_{n+1}^m}{I_1}\left(1-\frac{I_1}{I_{n+1}^m}\right)\hfill \\ & +& {\beta }_1{T}_{1}f\left(V_1\right)\left(1-\frac{V_1}{V_{n+1}^m}\right)\left(\frac{I_{n-m_2+1}^m}{I_1}-\frac{V_{n+1}^m}{V_1}\right)-\frac{p_1{d}_4{e}^{{\mu }_1{\tau }_1}}{q}{Z}_{n+1}^m\hfill \\ & +& {\beta }_1{T}_{1}f\left(V_1\right)\left[\phi \left(\frac{T_{n+1}^{m}f\left(V_n^m\right)}{T_{1}f\left(V_1\right)}\right)-\phi \left(\frac{T_{n-m_1+1}^{m}f\left(V_{n-m_1}^m\right)}{T_{1}f\left(V_1\right)}\right)\right]\hfill \\ & +& {\beta }_2{T}_{1}g\left(I_1\right)\left[\phi \left(\frac{T_{n+1}^{m}g\left(I_n^m\right)}{T_{1}g\left(I_1\right)}\right)-\phi \left(\frac{T_{n-m_1+1}^{m}g\left(I_{n-m_1}^m\right)}{T_{1}g\left(I_1\right)}\right)\right]\hfill \\ & +& {\beta }_1{T}_{1}f\left(V_1\right)\left[\phi \left(\frac{I_{n+1}^m}{I_1}\right)-\phi \left(\frac{I_{n-m_2+1}^m}{I_1}\right)\right]\hfill \\ & +& {\beta }_1{T}_{1}f\left(V_1\right)\left(\frac{f\left(V_{n+1}^m\right)}{f\left(V_1\right)}-\frac{f\left(V_n^m\right)}{f\left(V_1\right)}+ln\frac{f\left(V_n^m\right)}{f\left(V_{n+1}^m\right)}\right)\hfill \\ & +& {\beta }_2{T}_{1}g\left(I_1\right)\left(\frac{g\left(I_{n+1}^m\right)}{g\left(I_1\right)}-\frac{g\left(I_n^m\right)}{g\left(I_1\right)}+ln\frac{g\left(I_n^m\right)}{g\left(I_{n+1}^m\right)}\right)\}\hfill \\ & -& \sum _{m=0}^N\frac{{\beta }_1{T}_{1}f\left(V_1\right)D}{p_2{I}_1{e}^{-{\mu }_2{\tau }_2}{(\Delta x)}^2}\left(1-\frac{V_1}{V_{n+1}^m}\right)\left(V_{n+1}^{m+1}-2V_{n+1}^m+V_{n+1}^{m-1}\right)\hfill \\ & =& \sum _{m=0}^N\{d_1{T}_1\left(1-\frac{T_1}{T_{n+1}^m}\right)\left(1-\frac{T_{n+1}^m}{T_1}\right)+{\beta }_1{T}_{1}f\left(V_1\right)[3-\frac{T_1}{T_{n+1}^m}-\frac{V_1{I}_{n-m_2+1}^m}{V_{n+1}^m{I}_1}\hfill \\ & -& \frac{T_{n-m_1+1}^{m}f\left(V_{n-m_1}^m\right)I_1}{T_{1}f\left(V_1\right)I_{n+1}^m}+\frac{f\left(V_{n+1}^m\right)}{f\left(V_1\right)}-\frac{V_{n+1}^m}{V_1}+ln\frac{T_{n-m_1+1}f\left(V_{n-m_1}^m\right)I_{n-m_2+1}^m}{T_{n+1}^{m}f\left(V_{n+1}^m\right)I_{n+1}^m}]\hfill \\ & +& {\beta }_2{T}_{1}g\left(I_1\right)[2-\frac{T_1}{T_{n+1}^m}-\frac{T_{n-m_1+1}^{m}g\left(I_{n-m_1}^m\right)I_1}{T_{1}g\left(I_1\right)I_{n+1}^m}+\frac{g\left(I_{n+1}^m\right)}{g\left(I_1\right)}\hfill \\ & -& \frac{I_{n+1}^m}{I_1}+ln\frac{T_{n-m_1+1}^{m}g\left(I_{n-m_1}^m\right)}{T_{n+1}^{m}g\left(I_{n+1}^m\right)}]\hfill \\ & -& \frac{p_1{d}_4{e}^{{\mu }_1{\tau }_1}}{q}{Z}_{n+1}^m\}-\frac{{\beta }_1{T}_{1}f\left(V_1\right)D}{p_2{I}_1{e}^{-{\mu }_2{\tau }_2}{(\Delta x)}^2}{V}_1\sum _{m=0}^{N-1}\frac{{(V_{n+1}^{m+1}-V_{n+1}^m)}^2}{V_{n+1}^{m+1}{V}_{n+1}^m}\hfill \\ & =& \sum _{m=0}^N\{d_1{T}_1\left(1-\frac{T_1}{T_{n+1}^m}\right)\left(1-\frac{T_{n+1}^m}{T_1}\right)+{\beta }_1{T}_{1}f\left(V_1\right)[-\phi \left(\frac{T_1}{T_{n+1}^m}\right)\hfill \\ & -& \phi \left(\frac{T_{n-m_1+1}^{m}f\left(V_{n-m_1}^m\right)I_1}{T_{1}f\left(V_1\right)I_{n+1}^m}\right)-\phi \left(\frac{V_1{I}_{n-m_2+1}^m}{V_{n+1}^m{I}_1}\right)+\frac{f\left(V_{n+1}^m\right)}{f\left(V_1\right)}-\frac{V_{n+1}^m}{V_1}\hfill \\ & +& ln\frac{f\left(V_1^m\right)V_{n+1}^m}{f\left(V_{n+1}^m\right)V_1}]+{\beta }_2{T}_{1}g\left(I_1\right)[-\phi \left(\frac{T_1}{T_{n+1}^m}\right)-\phi \left(\frac{T_{n-m_1+1}^{m}g\left(I_{n-m_1}^m\right)I_1}{T_{1}g\left(I_1\right)I_{n+1}^m}\right)\hfill \\ & +& \frac{g\left(I_{n+1}^m\right)}{g\left(I_1\right)}-\frac{I_{n_1}^m}{I_1}+ln\frac{g\left(I_1^m\right)I_{n+1}^m}{g\left(I_{n+1}^m\right)I_1}]\hfill \\ & -& \frac{p_1{d}_4{e}^{{\mu }_1{\tau }_1}}{q}{Z}_{n+1}^m\}-\frac{{\beta }_1{T}_{1}f\left(V_1\right)D}{p_2{I}_1{e}^{-{\mu }_2{\tau }_2}{(\Delta x)}^2}{V}_1\sum _{m=0}^{N-1}\frac{{(V_{n+1}^{m+1}-V_{n+1}^m)}^2}{V_{n+1}^{m+1}{V}_{n+1}^m}\hfill \\ & =& \sum _{m=0}^N\{d_1{T}_1\left(1-\frac{T_1}{T_{n+1}^m}\right)\left(1-\frac{T_{n+1}^m}{T_1}\right)+{\beta }_1{T}_{1}f\left(V_1\right)[-\phi \left(\frac{T_1}{T_{n+1}^m}\right)\hfill \\ & -& \phi \left(\frac{T_{n-m_1+1}^{m}f\left(V_{n-m_1}^m\right)I_1}{T_{1}f\left(V_1\right)I_{n+1}^m}\right)-\phi \left(\frac{V_1{I}_{n-m_2+1}^m}{V_{n+1}^m{I}_1}\right)\hfill \\ & -& \phi \left(\frac{f\left(V_1^m\right)V_{n+1}^m}{f\left(V_{n+1}^m\right)V_1}\right)+\left(\frac{f\left(V_{n+1}^m\right)}{f\left(V_1\right)}-\frac{V_{n+1}^m}{V_1}\right)\left(1-\frac{f\left(V_{n+1}^m\right)V_1}{f\left(V_1^m\right)V_{n+1}^m}\right)]\hfill \\ & +& {\beta }_2{T}_{1}g\left(I_1\right)[-\phi \left(\frac{T_1}{T_{n+1}^m}\right)-\phi \left(\frac{T_{n-m_1+1}^{m}g\left(I_{n-m_1}^m\right)I_1}{T_{1}g\left(I_1\right)I_{n+1}^m}\right)-\phi \left(\frac{g\left(I_1^m\right)I_{n+1}^m}{g\left(I_{n+1}^m\right)I_1}\right)\hfill \\ & +& \left(\frac{g\left(I_{n+1}^m\right)}{g\left(I_1\right)}-\frac{I_{n+1}^m}{I_1}\right)\left(1-\frac{g\left(I_1\right)}{g\left(I_{n+1}^m\right)}\right)]-\frac{p_1{d}_4{e}^{{\mu }_1{\tau }_1}}{q}{Z}_{n+1}^m\}\hfill \\ & -& -\frac{{\beta }_1{T}_{1}f\left(V_1\right)D}{p_2{I}_1{e}^{-{\mu }_2{\tau }_2}{(\Delta x)}^2}{V}_1\sum _{m=0}^{N-1}\frac{{(V_{n+1}^{m+1}-V_{n+1}^m)}^2}{V_{n+1}^{m+1}{V}_{n+1}^m}.\hfill \end{array}$

Similar to the proof of Theorem 3, we have

$\left(\frac{f\left(V\right(x,t\left)\right)}{f\left(V_1\right)}-\frac{V(x,t)}{V_1}\right)\left(1-\frac{f\left(V_1\right)}{f\left(V\right(x,t\left)\right)}\right)\le 0,$

$\left(\frac{g\left(I\right(x,t\left)\right)}{g\left(I_1\right)}-\frac{I(x,t)}{I_1}\right)\left(1-\frac{g\left(I_1\right)}{g\left(I\right(x,t\left)\right)}\right)\le 0.$

It follows from $\phi \left(u\right)\ge 0$ that $({\tilde{W}}_{n+1}-{\tilde{W}}_n)\le 0$, for all $n\in \aleph$, this implies that $\left\{{\tilde{W}}_n\right\}$ is monotone decreasing sequence. As ${\tilde{W}}_n\ge 0$, then ${lim}_{n\to \infty }{\tilde{W}}_n\ge 0$, ${lim}_{n\to \infty }({\tilde{W}}_{n+1}-{\tilde{W}}_n)=0$, so that ${lim}_{n\to \infty }{T}_n^m=T_1$. Combined with system (7), we obtain ${lim}_{n\to \infty }{I}_n^m=I_1$, ${lim}_{n\to \infty }{V}_n^m=V_1$ and ${lim}_{n\to \infty }{Z}_n^m=0$, for all $m\in \{0,1,\cdots ,N\}$, then $E_1$ of system (7) is globally asymptotically stable. This completes the proof. □

Theorem 8.

For system (7), if $R_1>1$, the interior equilibrium $E_2$ is not globally asymptotically stable.

Proof. 

Define the discrete Lyapunov function as follows

$\begin{array}{ccc}\hfill {\overline{W}}_n& =& \sum _{m=0}^N\left\{\frac{1}{\Delta t}\right[T_2\phi \left(\frac{T_n^m}{T_2}\right)+e^{{\mu }_1{\tau }_1}{I}_2\phi \left(\frac{I_n^m}{I_2}\right)+\frac{{\beta }_1{T}_{2}f\left(V_2\right)}{p_2{e}^{-{\mu }_2{\tau }_2}{I}_2}{V}_2\phi \left(\frac{V_n^m}{V_2}\right)\hfill \\ & +& \frac{p_1{e}^{{\mu }_1{\tau }_1}}{q}{Z}_2\phi \left(\frac{Z_n^m}{Z_2}\right)+\Delta t{p}_1{e}^{{\mu }_1{\tau }_1}{I}_2{Z}_n^m]+{\beta }_1{T}_{2}f\left(V_2\right)\sum _{j=n-m_1}^{n-1}\phi \left(\frac{T_{j+1}^{m}f\left(V_j^m\right)}{T_{2}f\left(V_2\right)}\right)\hfill \\ & +& {\beta }_2{T}_{2}g\left(I_2\right)\sum _{j=n-m_1}^{n-1}\phi \left(\frac{T_{j+1}^{m}g\left(I_j^m\right)}{T_{2}g\left(I_2\right)}\right)+{\beta }_1{T}_{2}f\left(V_2\right)\sum _{j=n-m_2}^{n-1}\phi \left(\frac{I_{j+1}^m}{I_2}\right)\hfill \\ & +& {\beta }_1{T}_{2}f\left(V_2\right)\phi \left(\frac{f\left(V_n^m\right)}{f\left(V_2\right)}\right)+{\beta }_2{T}_{2}g\left(I_2\right)\phi \left(\frac{g\left(I_n^m\right)}{g\left(I_2\right)}\right)\},\hfill \end{array}$

it follows from $u-1\ge lnu$ that ${\overline{W}}_n\ge 0$ for all $n\in \aleph$. Then, along the trajectory of (7)

$\begin{array}{ccc}\hfill {\overline{W}}_{n+1}-{\overline{W}}_n& =& \sum _{m=0}^N\left\{\frac{1}{\Delta t}\right[T_{n+1}^m-T_n^m+T_{2}ln\frac{T_n^m}{T_{n+1}^m}+e^{{\mu }_1{\tau }_1}\left(I_{n+1}^m-I_n^m+I_{2}ln\frac{I_n^m}{I_{n+1}^m}\right)\hfill \\ & +& \frac{{\beta }_1{T}_{2}f\left(V_2\right)}{p_2{e}^{-{\mu }_2{\tau }_2}{I}_2}\left(V_{n+1}^m-V_n^m+V_{2}ln\frac{V_n^m}{V_{n+1}^m}\right)+\frac{p_1{e}^{{\mu }_1{\tau }_1}}{q}(Z_{n+1}^m-Z_n^m\hfill \\ & +& Z_{2}ln\frac{Z_n^m}{Z_{n+1}^m})+\Delta t{p}_1{e}^{{\mu }_1{\tau }_1}{I}_2\left(Z_{n+1}^m-Z_n^m\right)]\hfill \\ & +& {\beta }_1{T}_{2}f\left(V_2\right)\left[\sum _{j=n-m_1+1}^n\phi \left(\frac{T_{j+1}^{m}f\left(V_j^m\right)}{T_{2}f\left(V_2\right)}\right)-\sum _{j=n-m_1}^{n-1}\phi \left(\frac{T_{j+1}^{m}f\left(V_j^m\right)}{T_{2}f\left(V_2\right)}\right)\right]\hfill \\ & +& {\beta }_2{T}_{2}g\left(I_2\right)\left[\sum _{j=n-m_1+1}^n\phi \left(\frac{T_{j+1}^{m}g\left(I_j^m\right)}{T_{2}g\left(I_2\right)}\right)-\sum _{j=n-m_1}^{n-1}\phi \left(\frac{T_{j+1}^{m}g\left(I_j^m\right)}{T_{2}g\left(I_2\right)}\right)\right]\hfill \\ & +& {\beta }_1{T}_{2}f\left(V_2\right)\left[\sum _{j=n-m_2+1}^n\phi \left(\frac{I_{j+1}^m}{I_2}\right)-\sum _{j=n-m_2}^{n-1}\phi \left(\frac{I_{j+1}^m}{I_2}\right)\right]\hfill \\ & +& {\beta }_1{T}_{2}f\left(V_2\right)\left(\frac{f\left(V_{n+1}^m\right)}{f\left(V_2\right)}-\frac{f\left(V_n^m\right)}{f\left(V_2\right)}+ln\frac{f\left(V_n^m\right)}{f\left(V_{n+1}^m\right)}\right)\hfill \\ & +& {\beta }_2{T}_{2}g\left(I_2\right)\left(\frac{g\left(I_{n+1}^m\right)}{g\left(I_2\right)}-\frac{g\left(I_n^m\right)}{g\left(I_2\right)}+ln\frac{g\left(I_n^m\right)}{g\left(I_{n+1}^m\right)}\right)\}\hfill \\ & \le & \sum _{m=0}^N\left\{\frac{1}{\Delta t}\right[\left(1-\frac{T_2}{T_{n+1}^m}\right)\left(T_{n+1}^m-T_n^m\right)+e^{{\mu }_1{\tau }_1}\left(1-\frac{I_2}{I_{n+1}^m}\right)\left(I_{n+1}^m-I_n^m\right)\hfill \\ & +& \frac{{\beta }_1{T}_{2}f\left(V_2\right)}{p_2{e}^{-{\mu }_2{\tau }_2}{I}_2}\left(1-\frac{V_2}{V_{n+1}^m}\right)\left(V_{n+1}^m-V_n^m\right)\hfill \\ & +& \frac{p_1{e}^{{\mu }_1{\tau }_1}}{q}\left(1-\frac{Z_2}{Z_{n+1}^m}\right)\left(Z_{n+1}^m-Z_n^m\right)+\Delta t{p}_1{e}^{{\mu }_1{\tau }_1}{I}_2\left(Z_{n+1}^m-Z_n^m\right)]\hfill \\ & +& {\beta }_1{T}_{2}f\left(V_2\right)\left[\phi \left(\frac{T_{n+1}^{m}f\left(V_n^m\right)}{T_{2}f\left(V_2\right)}\right)-\phi \left(\frac{T_{n-m_1+1}^{m}f\left(V_{n-m_1}^m\right)}{T_{2}f\left(V_2\right)}\right)\right]\hfill \\ & +& {\beta }_2{T}_{2}g\left(I_2\right)\left[\phi \left(\frac{T_{n+1}^{m}g\left(I_n^m\right)}{T_{2}g\left(I_2\right)}\right)-\phi \left(\frac{T_{n-m_1+1}^{m}g\left(I_{n-m_1}^m\right)}{T_{2}g\left(I_2\right)}\right)\right]\hfill \\ & +& {\beta }_1{T}_{2}f\left(V_2\right)\left[\phi \left(\frac{I_{n+1}^m}{I_2}\right)-\phi \left(\frac{I_{n-m_2+1}^m}{I_2}\right)\right]\hfill \\ & +& {\beta }_1{T}_{2}f\left(V_2\right)\left(\frac{f\left(V_{n+1}^m\right)}{f\left(V_2\right)}-\frac{f\left(V_n^m\right)}{f\left(V_2\right)}+ln\frac{f\left(V_n^m\right)}{f\left(V_{n+1}^m\right)}\right)\hfill \\ & +& {\beta }_2{T}_{2}g\left(I_2\right)\left(\frac{g\left(I_{n+1}^m\right)}{g\left(I_2\right)}-\frac{g\left(I_n^m\right)}{g\left(I_2\right)}+ln\frac{g\left(I_n^m\right)}{g\left(I_{n+1}^m\right)}\right)\}.\hfill \end{array}$

From the equilibrium condition of $E_2$, we have

$\lambda =d_1{T}_2+{\beta }_1{T}_{2}f\left(V_2\right)+{\beta }_2{T}_{2}g\left(I_2\right),p_2{e}^{-{\mu }_2{\tau }_2}{I}_2=d_3{V}_2,$

${\beta }_1{T}_{2}f\left(V_2\right)+{\beta }_2{T}_{2}g\left(I_2\right)=e^{{\mu }_1{\tau }_1}(d_2{I}_2+p_1{I}_2{Z}_2),I_2=\frac{d_4}{q},$

then

$\begin{array}{ccc}\hfill {\overline{W}}_{n+1}-{\overline{W}}_n& \le & \sum _{m=0}^N\{d_2{T}_2\left(1-\frac{T_2}{T_{n+1}^m}\right)\left(1-\frac{T_{n+1}^m}{T_2}\right)+\left(1-\frac{T_2}{T_{n+1}^m}\right)({\beta }_1{T}_{2}f\left(V_2\right)\hfill \\ & +& {\beta }_2{T}_{2}g\left(I_2\right))-\left(1-\frac{T_2}{T_{n+1}^m}\right)\left({\beta }_1{T}_{n+1}^{m}f\left(V_n^m\right)+{\beta }_2{T}_{n+1}^{m}g\left(I_n^m\right)\right)\hfill \\ & +& \left(1-\frac{I_2}{I_{n+1}^m}\right)\left({\beta }_1{T}_{n-m_1+1}^{m}f\left(V_{n-m_1}^m\right)+{\beta }_2{T}_{n-m_1+1}^{m}g\left(I_{n-m_1}^m\right)\right)\hfill \\ & -& d_2{e}^{{\mu }_1{\tau }_1}{I}_{n+1}^m-p_1{e}^{{\mu }_1{\tau }_1}{I}_{n+1}^m{Z}_n^m+{\beta }_1{T}_{2}f\left(V_2\right)+{\beta }_2{T}_{2}g\left(I_2\right)\hfill \\ & -& p_1{e}^{{\mu }_1{\tau }_1}{I}_2{Z}_2+p_1{e}^{{\mu }_1{\tau }_1}{I}_2{Z}_n^m+p_1{e}^{{\mu }_1{\tau }_1}{I}_2(Z_{n+1}^m-Z_n^m)\hfill \\ & +& {\beta }_1{T}_{2}f\left(V_2\right)\left(1-\frac{V_2}{V_{n+1}^m}\right)\left(\frac{I_{n-m_2+1}^m}{I_2}-\frac{V_{n+1}^m}{V_2}\right)\hfill \\ & +& p_1{e}^{{\mu }_1{\tau }_1}{I}_{n+1}^m{Z}_n^m-\frac{p_1{d}_4{e}^{{\mu }_1{\tau }_1}}{q}{Z}_{n+1}^m-p_1{e}^{{\mu }_1{\tau }_1}{Z}_2\frac{I_{n+1}^m{Z}_n^m}{Z_{n+1}^m}\hfill \\ & +& p_1{e}^{{\mu }_1{\tau }_1}{I}_2{Z}_2+{\beta }_1{T}_{2}f\left(V_2\right)\left[\phi \left(\frac{T_{n+1}^{m}f\left(V_n^m\right)}{T_{2}f\left(V_2\right)}\right)-\phi \left(\frac{T_{n-m_1+1}^{m}f\left(V_{n-m_1}^m\right)}{T_{2}f\left(V_2\right)}\right)\right]\hfill \\ & +& {\beta }_2{T}_{2}g\left(I_2\right)\left[\phi \left(\frac{T_{n+1}^{m}g\left(I_n^m\right)}{T_{2}g\left(I_2\right)}\right)-\phi \left(\frac{T_{n-m_1+1}^{m}g\left(I_{n-m_1}^m\right)}{T_{2}g\left(I_2\right)}\right)\right]\hfill \\ & +& {\beta }_1{T}_{2}f\left(V_2\right)\left[\phi \left(\frac{I_{n+1}^m}{I_2}\right)-\phi \left(\frac{I_{n-m_2+1}^m}{I_2}\right)\right]\hfill \\ & +& {\beta }_1{T}_{2}f\left(V_2\right)\left(\frac{f\left(V_{n+1}^m\right)}{f\left(V_2\right)}-\frac{f\left(V_n^m\right)}{f\left(V_2\right)}+ln\frac{f\left(V_n^m\right)}{f\left(V_{n+1}^m\right)}\right)\hfill \\ & +& {\beta }_2{T}_{2}g\left(I_2\right)\left(\frac{g\left(I_{n+1}^m\right)}{g\left(I_2\right)}-\frac{g\left(I_n^m\right)}{g\left(I_2\right)}+ln\frac{g\left(I_n^m\right)}{g\left(I_{n+1}^m\right)}\right)\}\hfill \\ & -& \frac{D{V}_2{\beta }_1{T}_{2}f\left(V_2\right)}{p_2{e}^{-{\mu }_2{\tau }_2}{I}_2{(\Delta x)}^2}\sum _{m=0}^{N-1}\frac{{(V_{n+1}^{m+1}-V_{n+1}^m)}^2}{V_{n+1}^{m+1}{V}_{n+1}^m}\hfill \\ & =& \sum _{m=0}^N\{d_1{T}_2\left(1-\frac{T_2}{T_{n+1}^m}\right)\left(1-\frac{T_{n+1}^m}{T_2}\right)+{\beta }_1{T}_{2}f\left(V_2\right)[3-\frac{T_2}{T_{n+1}^m}-\frac{V_2{I}_{n-m_2+1}^m}{V_{n+1}^m{I}_2}\hfill \\ & -& \frac{T_{n-m_1+1}^{m}f\left(V_{n-m_1}^m\right)I_2}{T_{2}f\left(V_2\right)I_{n+1}^m}+\frac{f\left(V_{n+1}^m\right)}{f\left(V_2\right)}-\frac{V_{n+1}^m}{V_2}+ln\frac{T_{n-m_1+1}^{m}f\left(V_{n-m_1}^m\right)I_{n-m_2+1}}{T_{n+1}^{m}f\left(V_{n+1}^m\right)I_{n+1}}]\hfill \\ & +& {\beta }_2{T}_{2}g\left(I_2\right)[2-\frac{T_2}{T_{n+1}^m}-\frac{T_{n-m_1+1}^{m}g\left(I_{n-m_1}^m\right)I_2}{T_{2}g\left(I_2\right)I_{n+1}^m}+\frac{g\left(I_{n+1}^m\right)}{g\left(I_2\right)}\hfill \\ & -& \frac{I_{n+1}^m}{I_2}+ln\frac{T_{n-m_1+1}^{m}g\left(I_{n-m_1}^m\right)}{T_{n+1}^{m}g\left(I_{n+1}^m\right)}]-d_2{e}^{{\mu }_1{\tau }_1}{I}_{n+1}^m-p_1{e}^{{\mu }_1{\tau }_1}{Z}_2\frac{I_{n+1}^m{Z}_n^m}{Z_{n+1}^m}\}\hfill \\ & -& \frac{D{V}_2{\beta }_1{T}_{2}f\left(V_2\right)}{p_2{e}^{-{\mu }_2{\tau }_2}{I}_2{(\Delta x)}^2}\sum _{m=0}^{N-1}\frac{{(V_{n+1}^{m+1}-V_{n+1}^m)}^2}{V_{n+1}^{m+1}{V}_{n+1}^m}\hfill \\ & =& \sum _{m=0}^N\{d_1{T}_2\left(1-\frac{T_2}{T_{n+1}^m}\right)\left(1-\frac{T_{n+1}^m}{T_2}\right)+{\beta }_1{T}_{2}f\left(V_2\right)[-\phi \left(\frac{T_2}{T_{n+1}^m}\right)\hfill \\ & -& \phi \left(\frac{V_2{I}_{n-m_2+1}^m}{V_{n+1}^m{I}_2}\right)-\phi \left(\frac{T_{n-m_1+1}^{m}f\left(V_{n-m_1}^m\right)I_2}{T_{2}f\left(V_2\right)I_{n+1}^m}\right)\hfill \\ & +& \left(\frac{f\left(V_{n+1}\right)}{f\left(V_2\right)}-\frac{V_{n+1}}{V_2}\right)\left(1-\frac{f\left(V_2\right)}{f\left(V_{n+1}\right)}\right)]\hfill \\ & +& {\beta }_2{T}_{2}g\left(I_2\right)[-\phi \left(\frac{T_2}{T_{n+1}^m}\right)-\phi \left(\frac{T_{n-m_1+1}^{m}g\left(I_{n-m_1}^m\right)I_2}{T_{2}g\left(I_2\right)I_{n+1}^m}\right)+(\frac{g\left(I_{n+1}\right)}{g\left(I_2\right)}\hfill \\ & -& \frac{I_{n+1}^m}{I_2}\left)\left(1-\frac{g\left(I_2\right)}{g\left(I_{n+1}\right)}\right)\right]-d_2{e}^{{\mu }_1{\tau }_1}{I}_{n+1}^m-p_1{e}^{{\mu }_1{\tau }_1}{Z}_2\frac{I_{n+1}^m{Z}_n^m}{Z_{n+1}^m}\}\hfill \\ & -& \frac{D{V}_2{\beta }_1{T}_{2}f\left(V_2\right)}{p_2{e}^{-{\mu }_2{\tau }_2}{I}_2{(\Delta x)}^2}\sum _{m=0}^{N-1}\frac{{(V_{n+1}^{m+1}-V_{n+1}^m)}^2}{V_{n+1}^{m+1}{V}_{n+1}^m}.\hfill \end{array}$

Similar to the proof of Theorem 3, we have

$\left(\frac{f\left(V\right(x,t\left)\right)}{f\left(V_1\right)}-\frac{V}{V_1}\right)\left(1-\frac{f\left(V_1\right)}{f\left(V\right(x,t\left)\right)}\right)\le 0,$

$\left(\frac{g\left(I\right(x,t\left)\right)}{g\left(I_1\right)}-\frac{I(x,t)}{I_1}\right)\left(1-\frac{g\left(I_1\right)}{g\left(I\right(x,t\left)\right)}\right)\le 0,$

this implies that $\left\{{\overline{W}}_n\right\}$ is a monotone decreasing sequence, then ${\overline{W}}_n\ge 0$, ${lim}_{n\to \infty }{\overline{W}}_n\ge 0$, therefore

$\underset{n\to \infty }{lim}({\overline{W}}_{n+1}-{\overline{W}}_n)=0.$

According to the system (7), we claim that the CTL-activated equilibrium $E_2$ is not globally asymptotically stable. In fact, if the CTL-activated equilibrium $E_2$ is globally asymptotically stable, from the above inequality, we have

$0\le -d_2{e}^{{\mu }_1{\tau }_1}{I}_2-p_1{e}^{{\mu }_1{\tau }_1}{Z}_2{I}_2<0,$

this is a contradiction. This completes the proof. □

4. Numerical Simulation

In this section, we choose $f\left(V\right)=V,g\left(I\right)=I$, some numerical results of system (4) are presented for supporting our analytic results. Based on biological meanings of virus dynamics model from papers [39,40], we have estimated the values of our model parameters as follows:

If we choose $D=3$, then we can give a numerical simulation of the stability of system (4). Using the data in Table 1, we first show in a simulation that the interior equilibrium is stable (see Figure 1).

From Figure 1, we can see that the population has gradually stabilized after a sharp fluctuation.

If we choose ${\beta }_1=0.0003$ and ${\beta }_2=0.004$, then $R_0<1$. We can simulate that the infection-free equilibrium is globally asymptotically stable (see Figure 2).

From Figure 2, we can see that the number of infected cells, virus and CTLs tends to zero, except uninfected cells.

If we choose $q=0.000004$ and $p_2=0.9$, then $R_1\le 1<R_0$. This moment the CTL-inactivated equilibrium is globally asymptotically stable (see Figure 3).

From Figure 3, we can see that the population in the compartment $CTLs$ tends to 0. In addition, except for $CTLs$, the number of uninfected cells, infected cells, virus tends to certain constants.

The novelty of this paper is that we consider the effects of diffusion, time delay, and abstract functions on the spread of viruses. In order to see the impact of proliferation on the spread of the virus more intuitively, we first choose $q=0.04$. Next, we select $D=0$ and $D=300$ decibels to simulate the image of I while other parameters keep the values in the Table 1.

The left image in Figure 4 is an image without time delay, and the right image is an image with time delay equal to 300. Since we are simulating long-term dynamic behavior, from the overall image of the two figures, there is no obvious difference in either the stable position or the growth rate. So where is the effect of diffusion reflected? We believe that the effect of diffusion should be reflected in the growth of I. Therefore, we project the two graphs in Figure 4 on the time-quantity axis (Figure 5).

From the left image of Figure 5, we can clearly see that when there is no time delay, the image rises smoothly and the curve is smooth. When the time lag is equal to 300, the image is not a smooth curve, which shows that the proliferation brings about the proliferation of infected cells and the uneven fluctuation.

5. Conclusions and Discussion

It is necessary to understand the dynamics model for HIV infection since these infected cells usually cause a $CTL$ response from the immune system. In this paper, we first developed a diffusive infection model (4) with general nonlinear incidence rate and two delays on the base of model (3), we show that the global stability of equilibria is completely determined by the reproductive numbers for viral infection $R_0$ and for $CTL$ immune response $R_1$. Second, we considered the corresponding discretization of the continuous model by using nonstandard finite difference scheme, and then studied the global stability of the discrete system. Some numerical simulations were also presented to support our analytic results. In general, systems of PDE cannot be solved explicitly, and numerical solutions have to be studied instead. By using the NSFD scheme, we showed that the proposed discrete model partly preserves the global stability of equilibria of the corresponding continuous model. We plan to address how other diffusive terms (for infected and uninfected cells) affect the model in future work.

Author Contributions

X.-L.L.: writing original draft preparation, methodology, formal analysis, editing, validation, funding acquisition; C.-C.Z.: validation, investigation, data curation, visualization. All authors have read and agreed to the published version of the manuscript.

Funding

This project is supported by the Natural Science Foundation of Jiangsu Province, China (Grant No. BK20190578)and the science and technology foundation of Suqian(S201818, Z2021131), China.

Institutional Review Board Statement

Not applicable.

Informed Consent Statement

Not applicable.

Data Availability Statement

Not applicable.

Acknowledgments

The authors are very grateful to the reviewers for their invaluable and expert coments.

Conflicts of Interest

The authors declare no conflict of interest.

Publisher’s Note: MDPI stays neutral with regard to jurisdictional claims in published maps and institutional affiliations.

Figures and Table

Figure 1. When [Forumla omitted. See PDF.], the interior equilibrium [Forumla omitted. See PDF.] is globally asymptotically stable.

Figure 2. When [Forumla omitted. See PDF.], the infection-free equilibrium [Forumla omitted. See PDF.] is globally asymptotically stable.

Figure 3. When [Forumla omitted. See PDF.], the CTL-inactivated equilibrium [Forumla omitted. See PDF.] is globally asymptotically stable.

Figure 4. Comparison of compartment I at [Forumla omitted. See PDF.] and [Forumla omitted. See PDF.].

Figure 5. Comparison projection of compartment I when [Forumla omitted. See PDF.] and [Forumla omitted. See PDF.].

Table 1

State variables and parameters of HIV-1 infection model.

Parameter	Description
$\lambda$	$0.9$	Reference [40]
$d_1$	$0.03$	Reference [39]
$d_2$	$0.5$	Reference [39]
$d_3$	$0.1$	Reference [40]
$d_4$	$0.3$	Reference [40]
${\beta }_1$	$0.3$	Reference [40]
${\beta }_2$	$0.4$	Reference [40]
$p_1$	$0.08{\mathrm{day}}^{-1}$	Estimate
$p_2$	$0.5{\mathrm{day}}^{-1}$	Reference [40]
q	$0.4$	Estimate

References

1. Male, D.; Brostoff, J.; Roth, D.; Roitt, I. Immunology; 7th ed. Elsevier: Amsterdam, The Netherlands, 2006.

2. Nowak, M.; May, R. Vitus Dynamics; Oxford University Press: Oxford, UK, 2000.

3. Zhu, C.C.; Zhu, J. Spread trend of COVID-19 epidemic outbreak in China: Using exponential attractor method in a spatial heterogeneous SEIQR model. Math. Biosci. Eng.; 2020; 17, pp. 3062-3087. [DOI: https://dx.doi.org/10.3934/mbe.2020174] [PubMed: https://www.ncbi.nlm.nih.gov/pubmed/32987517]

4. Zhu, C.C.; Zhu, J. Dynamic analysis of a delayed COVID-19 epidemic with home quarantine in temporal-spatial heterogeneous via global exponential attractor method. Chaos Solitons Fractals; 2021; 143, 110546. [DOI: https://dx.doi.org/10.1016/j.chaos.2020.110546] [PubMed: https://www.ncbi.nlm.nih.gov/pubmed/33519115]

5. Zhu, C.C.; Zhu, J. The effect of self-limiting on the prevention and control of the diffuse COVID-19 epidemic with delayed and temporal-spatial heterogeneous. BMC Infect. Dis.; 2021; 21, 1145. [DOI: https://dx.doi.org/10.1186/s12879-021-06670-y] [PubMed: https://www.ncbi.nlm.nih.gov/pubmed/34753451]

6. Wodarz, D. Killer Cell Dynamics: Mathematical and Computational Approaches to Immunology; Springer: Berlin/Heidelberg, Germany, 2007.

7. Weber, J.N.; Weiss, R.A. HIV infection: The cellular picture. Sci. Am.; 1988; 259, pp. 101-109. [DOI: https://dx.doi.org/10.1038/scientificamerican1088-100] [PubMed: https://www.ncbi.nlm.nih.gov/pubmed/2854299]

8. Kirschner, D.E.; Webb, G.; Cloyd, M. Model of HIV-1 disease progression based on virus-induced lymph node homing and homing-induced apotosis of CD4⁺ lymphocytes. J. AIDS; 2000; 24, pp. 352-362.

9. Anderson, R.; May, R. Infections Diseases of Humans: Dynamics and Control; Oxford University Press: Oxford, UK, 1991.

10. Boer, R.D.; Perelson, A. Target cell limited and immune control models of HIV infection: A comparison. J. Theoret. Biol.; 1998; 190, pp. 201-214. [DOI: https://dx.doi.org/10.1006/jtbi.1997.0548] [PubMed: https://www.ncbi.nlm.nih.gov/pubmed/9514649]

11. Nowak, M.A.; Bonhoeffer, S.; Hill, A.M.; Boehme, R.; Thomas, H.C.; McDade, H. Viral dynamicss in hepatitis B virus infection. Proc. Natl. Acad. Sci. USA; 1996; 93, pp. 4398-4402. [DOI: https://dx.doi.org/10.1073/pnas.93.9.4398]

12. Lai, X.; Zou, X. Modeling HIV-1 virus dynamics with both virus-to-cell infection and cell-to-cell transimission. SIAM J. Appl. Math.; 2014; 74, pp. 898-917. [DOI: https://dx.doi.org/10.1137/130930145]

13. Xu, J.; Zhou, Y. Bifurcation analysis of HIV-1 infection model with cell-to-cell transimission and immune response delay. Math. Biosci. Eng.; 2016; 13, pp. 343-367. [DOI: https://dx.doi.org/10.3934/mbe.2015006]

14. Yang, Y.; Zou, L.; Ruan, S. Global dynamics of a delayed within-host viral infection model with both virus-to-cell and cell-to-cell transimissions. Math. Biosci.; 2015; 270, pp. 183-191. [DOI: https://dx.doi.org/10.1016/j.mbs.2015.05.001] [PubMed: https://www.ncbi.nlm.nih.gov/pubmed/25998145]

15. Wang, K.; Wang, W. Propagation of HBV with spatial dependence. Math. Biosci.; 2007; 210, pp. 78-95. [DOI: https://dx.doi.org/10.1016/j.mbs.2007.05.004]

16. Yang, Y.; Zhou, J.l.; Ma, X.S.; Zhang, T.H. Nonstandard finite difference scheme for a diffusive within-host virus dynamics model with both virus-to-cell and cell-to-cell transmissions. Comput. Math. Appl.; 2016; 72, pp. 1013-1020. [DOI: https://dx.doi.org/10.1016/j.camwa.2016.06.015]

17. Nowak, M.; Bangham, C. Population dynamics of immune responses to persistent viruses. Proc. Natl. Acad. Sci. USA; 1996; 272, pp. 74-79. [DOI: https://dx.doi.org/10.1126/science.272.5258.74] [PubMed: https://www.ncbi.nlm.nih.gov/pubmed/8600540]

18. Zhu, H.; Zou, X. Dynamics of a HIV-1 infection model with cell-mediated immune response and intracellular delay. Discrete Contin. Dyn. Syst. Ser. B; 2009; 12, pp. 511-524. [DOI: https://dx.doi.org/10.3934/dcdsb.2009.12.511]

19. Wang, J.L.; Guan, L.J. Global stability for a HIV-1 infection model with cell-mediated immune response and intrecellular delay. Discret. Contin. Dyn. Syst. Ser. B; 2012; 17, pp. 297-302.

20. Arnaout, R.; Nowak, M.; Wodarz, D. HIV-1 dynamics revisited: Biphasic decay by cytotoxic lymphocyte killing. Proc. R. Soc. Lond. B; 2000; 265, pp. 1347-1354. [DOI: https://dx.doi.org/10.1098/rspb.2000.1149] [PubMed: https://www.ncbi.nlm.nih.gov/pubmed/10972131]

21. Culshaw, R.; Ruan, S.; Spiteri, R. Optimal HIV treatment by maxinising immune response. J. Math. Biol.; 2004; 48, pp. 545-562. [DOI: https://dx.doi.org/10.1007/s00285-003-0245-3] [PubMed: https://www.ncbi.nlm.nih.gov/pubmed/15133623]

22. Wang, K.; Wang, W.; Liu, X. Global stability in a viral infection model with lytic and nonlytic immune response. Comput. Math. Appl.; 2006; 51, pp. 1593-1610. [DOI: https://dx.doi.org/10.1016/j.camwa.2005.07.020]

23. Xu, J.H.; Geng, Y. Dynamic consistent NSFD Scheme for a Delayed Viral Infection Model with Immune Response and Nonlinear Incidence. Discret. Dyn. Nat. Soc.; 2017; 2017, 3141736. [DOI: https://dx.doi.org/10.1155/2017/3141736]

24. Sigdel, R.P.; Mcluskey, C.C. Global stability for an SEI model of infection disease with immigration. Appl. Math. Comput.; 2014; 243, pp. 684-689.

25. Enatsu, Y.; Nakata, Y.; Muroya, Y.; Izzo, G.; Vecchio, A. Global dynamics of difference equations for SIR epidemic models with a class of nonlinear incidence rates. J. Differ. Equ. Appl.; 2012; 18, pp. 1163-1181. [DOI: https://dx.doi.org/10.1080/10236198.2011.555405]

26. Mickens, R.E. Nonstandard Finite Difference Models of Differential Equations; World Scientific: Singapore, 1994.

27. Travis, C.C.; Webb, G.F. Existence and stability for partial functional differential equations. Trans. Am. Math. Soc.; 1974; 200, pp. 395-418. [DOI: https://dx.doi.org/10.1090/S0002-9947-1974-0382808-3]

28. Shi, P.; Dong, L. Dynamical behaviors of a discrete HIV-1 virus model with bilinear infective rate. Math. Methods Appl. Sci.; 2014; 37, pp. 2271-2280. [DOI: https://dx.doi.org/10.1002/mma.2974]

29. Hattaf, K.; Yousfi, N. Global properties of a discrete viral infection model with general incidence rate. Math. Methods Appl. Sci.; 2016; 39, pp. 998-1004. [DOI: https://dx.doi.org/10.1002/mma.3536]

30. Yang, Y.; Zhou, J. Global dynamics of a PDE in-host viral model. Appl. Anal.; 2014; 93, pp. 2312-2329.

31. Wang, J.P.; Teng, Z.D.; Miao, H. Global dynamics for discrete-time analog of viral infection model with nonlinear incidence and CTL immune response. Adv. Differ. Equ.; 2016; 1, 143. [DOI: https://dx.doi.org/10.1186/s13662-016-0862-y] [PubMed: https://www.ncbi.nlm.nih.gov/pubmed/32226448]

32. Zhu, C.C.; Zhu, J. Stability of a reaction-diffusion alcohol model with the impact of tax policy. Comput. Math. Appl.; 2017; 74, pp. 613-633. [DOI: https://dx.doi.org/10.1016/j.camwa.2017.05.005]

33. Martin, R.H.; Smith, H.L. Abstract functional differential equations and reaction-diffusion systems. Trans. Am. Math. Soc.; 1990; 321, pp. 1-44.

34. Martin, R.H.; Smith, H.L. Reaction-diffusion systems with time delays: Monotonicity, invariance, comparison and convergence. J. Reine Angew. Math.; 1991; 413, pp. 1-35.

35. Wu, J. Theory and Applications of Parfunctional Differential Equations; Springer: New York, NY, USA, 1996.

36. Fitzgibbon, W.E. Semilinear functional differential equations in Banach space. J. Differ. Equ.; 1978; 29, pp. 1-14. [DOI: https://dx.doi.org/10.1016/0022-0396(78)90037-2]

37. Protter, M.H.; Weinberger, H.F. Maximum Principles in Differential Equations; Prentice Hall: Engle-Wood Cliffs, NJ, USA, 1967.

38. Henry, D. Gerometric Theory of Semilinear Parabolic Equation; Lecture Notes in Mathematics Springer: Berlin, Germany, 1993; 840.

39. Nelson, P.; Murray, J.; Perelson, A. A model of HIV-1 pathogenesis that includes an intracellular delay. Math. Biosci.; 2000; 163, pp. 201-215. [DOI: https://dx.doi.org/10.1016/S0025-5564(99)00055-3]

40. Xiang, H.; Tang, Y.L.; Huo, H.F. A viral model with intracellular delay and humoral immunity. Bull. Malays. Math. Sci. Soc.; 2017; 40, pp. 1011-1023. [DOI: https://dx.doi.org/10.1007/s40840-016-0326-2]