Global Stability of a MERS-CoV Infection Model with CTL Immune Response and Intracellular Delay

Abstract

In this paper, we propose and study a Middle East respiratory syndrome coronavirus (MERS-CoV) infection model with cytotoxic T lymphocyte (CTL) immune response and intracellular delay. This model includes five compartments: uninfected cells, infected cells, viruses, dipeptidyl peptidase 4 (DPP4), and CTL immune cells. We obtained an immunity-inactivated reproduction number $R_0$ and an immunity-activated reproduction number $R_1$. By analyzing the distributions of roots of the corresponding characteristic equations, the local stability results of the infection-free equilibrium, the immunity-inactivated equilibrium, and the immunity-activated equilibrium were obtained. Moreover, by constructing suitable Lyapunov functionals and combining LaSalle’s invariance principle and Barbalat’s lemma, some sufficient conditions for the global stability of the three types of equilibria were obtained. It was found that the infection-free equilibrium is globally asymptotically stable if $R_0\le 1$ and unstable if $R_0>1$; the immunity-inactivated equilibrium is locally asymptotically stable if $R_0>1>R_1$ and globally asymptotically stable if $R_0>1>R_1$ and condition (H1) holds, but unstable if $R_1>1$; and the immunity-activated equilibrium is locally asymptotically stable if $R_1>1$ and is globally asymptotically stable if $R_1>1$ and condition (H1) holds.

1. Introduction

Middle East respiratory syndrome (MERS) is a viral respiratory disease caused by the Middle East respiratory syndrome coronavirus (MERS-CoV), which has a high mortality rate (approximately 35%) and has become an important public health problem in many countries since it was first reported in Saudi Arabia in 2012 [1,2]. There is no vaccine or specific treatment for MERS, and treatment is mainly supportive based on the clinical status of MERS patients [2]. It is important to study MERS-CoV dynamics in the host to provide insights into the pathogenesis and treatment of MERS-CoV. In [3], the authors provide a systematic review and meta-analysis of severe acute respiratory syndrome coronavirus 2 (SARS-CoV-2), SARS-CoV, and MERS-CoV host virus dynamics. In [4], the authors used mathematical models combined with published viral load data to compare in detail the similarities and differences in the within-host viral dynamics of SARS-CoV-2, MERS-CoV, and SARS-CoV.

Dipeptidyl peptidase-4 (DPP4, also known as CD26) is the functional receptor for MERS-CoV [5]. The engagement of MERS-CoV spike proteins with DPP4 mediates viral attachment to host cells and virus–cell membrane fusion, thus playing a key role in MERS-CoV infection [6,7]. DPP4 receptors are present on the epithelial surfaces of various human organs (such as the kidney, intestine, liver, thymocytes), and their systematic distribution facilitates the transmission of viruses in the human body [8,9]. In the last decades, some classical virus dynamics models have been proposed to explore the relationship between uninfected cells (target cells), infected cells, and viral load (see, e.g., [10,11,12,13]). These studies have been of great help to understanding the virus dynamics in the host and to devising effective virus control strategies. Based on some classical works related to the modeling of virus dynamics in [14], to describe the effect of DPP4 on MERS-CoV infection in the host, the authors constructed the following four-dimensional ordinary differential equation model:

$$\left\{\begin{array}{cc}\hfill \dot{T}\left(t\right)=& \lambda -\beta D\left(t\right)v\left(t\right)T\left(t\right)-dT\left(t\right),\hfill \\ \hfill \dot{I}\left(t\right)=& \beta D\left(t\right)v\left(t\right)T\left(t\right)-d_{1}I\left(t\right),\hfill \\ \hfill \dot{v}\left(t\right)=& d_{1}MI\left(t\right)-cv\left(t\right),\hfill \\ \hfill \dot{D}\left(t\right)=& {\lambda }_1-{\beta }_1\beta D\left(t\right)v\left(t\right)T\left(t\right)-\gamma D\left(t\right).\hfill \end{array}\right.$$

(1)

Here, $T\left(t\right)$, $I\left(t\right)$, $v\left(t\right)$, and $D\left(t\right)$ denote the concentrations of uninfected cells, infected cells, free virus, and DPP4 on the surface of uninfected cells at time t, respectively. The constant $\lambda >0$ is the rate at which uninfected cells are produced. The constant $\beta >0$ is the rate at which uninfected cells are infected by the free virus (i.e., infected cells are increased at a mount of $\beta D\left(t\right)v\left(t\right)T\left(t\right)$ because uninfected cells are infected by the free virus). The constants $d>0$ and $d_1>0$ denote the death rates of uninfected cells and infected cells, respectively. The constant $M>0$ denotes the number of the free viruses released by the lysis of each infected cell after death. The constant $c>0$ denotes the death rate of the free viruses. The constant ${\lambda }_1>0$ denotes the rate at which DPP4 is produced on the surface of uninfected cells. The constant ${\beta }_1>0$ denotes the rate at which DPP4 is decreased (i.e., DPP4 is decreased at a mount of ${\beta }_1\beta D\left(t\right)v\left(t\right)T\left(t\right)$, because uninfected cells are infected by the free virus). The constant $\gamma >0$ denotes the hydrolysis rate of DPP4. In [14], the authors mainly studied the local and global stability of the infection-free equilibrium and the infected equilibrium of model (1). In addition, in [15], the authors further extend model (1) to its periodic case and obtain sufficient conditions for the existence of positive periodic solutions of the model by using the continuation theorem of the coincidence degree theory and constructing the appropriate auxiliary function.

Based on the important role of the CTL immune response in the control and clearance of MERS-CoV (CTL immune cells can attack virus-infected cells [16]), in [17], the authors further considered the existence of positive periodic solutions of the periodic model with CTL immune response. Several studies of viral infection models have also shown that the CTL immune response plays a positive role in reducing viral load, such as in the HIV model [10] and HCV model [18]. To incorporate the intracellular phase of the viral life cycle, in [19], the authors assumed that virus production lags by a delay $\tau$ behind the infection of a cell. The intracellular delay refers to the time between the entry of the virus into the target cell and the production of new viral particles. Viral infection models with intracellular delay have been studied by many scholars and have yielded many excellent results (e.g., [20,21,22,23,24,25,26,27]). Motivated by the above studies, in this paper, we will consider the following MERS-CoV infection model with CTL immune response and intracellular delay (see Figure 1):

$$\left\{\begin{array}{c}\dot{T}\left(t\right)=\lambda -\beta D\left(t\right)v\left(t\right)T\left(t\right)-dT\left(t\right),\hfill \\ \dot{I}\left(t\right)=e^{-d_1\tau }\beta D(t-\tau )v(t-\tau )T(t-\tau )-d_{1}I\left(t\right)-pI\left(t\right)Z\left(t\right),\hfill \\ \dot{v}\left(t\right)=d_{1}MI\left(t\right)-cv\left(t\right),\hfill \\ \dot{D}\left(t\right)={\lambda }_1-{\beta }_1\beta D\left(t\right)v\left(t\right)T\left(t\right)-\gamma D\left(t\right),\hfill \\ \dot{Z}\left(t\right)=qI\left(t\right)Z\left(t\right)-bZ\left(t\right).\hfill \end{array}\right.$$

(2)

In model (2), the state variables $T\left(t\right)$, $I\left(t\right)$, $v\left(t\right)$, and $D\left(t\right)$ as well as the parameters $\lambda$, $\beta$, d, $d_1$, M, c, ${\lambda }_1$, ${\beta }_1$, and $\gamma$ have the same biological meanings as in model (1). $Z\left(t\right)$ denotes the concentration of CTL immune cells at time t. CTL immune cells increase at a rate of $qI\left(t\right)Z\left(t\right)$ by the viral antigen of the infected cells and decay at rate $bZ\left(t\right)$, and infected cells are killed by the CTL immune response at rate $pI\left(t\right)Z\left(t\right)$. The constant $\tau \ge 0$ denotes the time between viral entry into an uninfected cell and the production of new virions. The term $e^{-d_1\tau }$ denotes the surviving rate of infected cells before they become productively infected [24]. All the parameters of model (2) are assumed to be positive constants except for the delay $\tau$. The main purpose of this paper is to study the local and global dynamics of model (2).

The rest of the paper is organized as follows. The boundedness of the solutions of model (2), and the classification of equilibria of model (2) are described in Section 2. Section 3 shows how by analyzing the distributions of roots of the corresponding characteristic equations, the local stability of the infection-free equilibrium, and the immunity-inactivated equilibrium, the immunity-activated equilibrium can be obtained. Section 4 mainly includes the following: (i) The global stability result of the infection-free equilibrium is obtained by constructing suitable Lyapunov functional and using LaSalle’s invariance principle. (ii) By some ingenious analytical techniques, explicit estimates of the ultimate lower bounds for the concentrations of viruses and infected cells are obtained. (iii) The global stability results of the immunity-inactivated equilibrium and the immunity-activated equilibrium are obtained by constructing suitable Lyapunov functionals and using Barbalat’s lemma. The last section presents a few numerical simulations and the conclusions of this paper.

2. Preliminaries

Let $X=C([-\tau ,0],{\mathbb{R}}^5)$ be the Banach space of continuous functions mapping from $[-\tau ,0]$ to ${\mathbb{R}}_{+}^5$ equipped with the sup-norm $\parallel \varphi \parallel ={sup}_{-\tau \le \theta \le 0}\left|\varphi \left(\theta \right)\right|$. The initial condition of (2) is given as follows:

$$\begin{array}{cc}\hfill & T\left(\theta \right)={\varphi }_1\left(\theta \right),I\left(\theta \right)={\varphi }_2\left(\theta \right),v\left(\theta \right)={\varphi }_3\left(\theta \right),D\left(\theta \right)={\varphi }_4\left(\theta \right),Z\left(\theta \right)={\varphi }_5\left(\theta \right),\theta \in [-\tau ,0],\hfill \end{array}$$

(3)

where $\varphi ={({\varphi }_1,{\varphi }_2,{\varphi }_3,{\varphi }_4,{\varphi }_5)}^T\in X^{+}:=\{\varphi \in X:\varphi \ge 0\}$.

2.1. The Well-Posedness and Dissipativeness

By the standard theory of functional differential equations (see [28]), it is easy to ascertain that the solution ${(T\left(t\right),I\left(t\right),v\left(t\right),D\left(t\right),Z\left(t\right))}^T$ of model (2) with the initial condition (3) is existent, unique, and nonnegative.

Theorem 1.

The solution ${(T\left(t\right),I\left(t\right),v\left(t\right),D\left(t\right),Z\left(t\right))}^T$ of model (2) with the initial condition (3) satisfies

$$\begin{array}{cc}\hfill & \underset{t\to +\infty }{lim\; sup}T\left(t\right)\le \frac{\lambda }{d}:=T_0,\hfill \\ \hfill & \underset{t\to +\infty }{lim\; sup}I\left(t\right)\le \frac{\lambda }{{\mu }_1}{e}^{-d_1\tau }:=I_{max},\hfill \\ \hfill & \underset{t\to +\infty }{lim\; sup}v\left(t\right)\le \frac{d_{1}M{I}_{max}}{c}:=v_{max},\hfill \\ \hfill & \underset{t\to +\infty }{lim\; sup}D\left(t\right)\le \frac{{\lambda }_1}{\gamma }:=D_0,\hfill \\ \hfill & \underset{t\to +\infty }{lim\; sup}Z\left(t\right)\le \frac{q\beta e^{-d_1\tau }{D}_0{T}_0{v}_{max}}{p{\mu }_2}:=Z_{max},\hfill \end{array}$$

(4)

where ${\mu }_1=min\{d,d_1\}$ and ${\mu }_2=min\{d_1,b\}$.

Proof. 

From the first and fourth equations of model (2), it follows that for $t\ge 0$,

$$\dot{T}\left(t\right)\le \lambda -dT\left(t\right),\dot{D}\left(t\right)\le {\lambda }_1-\gamma D\left(t\right),$$

which implies that

$$\underset{t\to +\infty }{lim\; sup}T\left(t\right)\le \frac{\lambda }{d}=T_0,\underset{t\to +\infty }{lim\; sup}D\left(t\right)\le \frac{{\lambda }_1}{\gamma }=D_0.$$

(5)

Define the following:

$$S_1\left(t\right)=e^{-d_1\tau }T(t-\tau )+I\left(t\right).$$

From the first and second equations of model (2), it follows that for $t\ge \tau$,

$$\begin{array}{cc}\hfill {\dot{S}}_1\left(t\right)& =e^{-d_1\tau }\lambda -e^{-d_1\tau }dT(t-\tau )-d_{1}I\left(t\right)-pI\left(t\right)Z\left(t\right)\hfill \\ & \le e^{-d_1\tau }\lambda -{\mu }_1{S}_1\left(t\right),\hfill \end{array}$$

which implies that

$$\underset{t\to +\infty }{lim\; sup}{S}_1\left(t\right)\le \frac{\lambda e^{-d_1\tau }}{{\mu }_1}.$$

Thus, we have

$$\underset{t\to +\infty }{lim\; sup}I\left(t\right)\le \frac{\lambda e^{-d_1\tau }}{{\mu }_1}=I_{max}.$$

(6)

From (6) and the third equation of model (2), it is not difficult to ascertain the following:

$$\underset{t\to +\infty }{lim\; sup}v\left(t\right)\le \frac{d_{1}M{I}_{max}}{c}=v_{max}.$$

(7)

From (5), (6), and (7), it follows that for any positive constant $\epsilon <min\{T_0,v_{max},D_0\}$, there exists a ${\varsigma }_0\left(\epsilon \right)>0$ such that, for $t\ge {\varsigma }_0\left(\epsilon \right)$,

$$T\left(t\right)<T_0+\epsilon ,v\left(t\right)<v_{max}+\epsilon ,D\left(t\right)<D_0+\epsilon .$$

We define the following:

$$S_2\left(t\right)=qI\left(t\right)+pZ\left(t\right).$$

Then, from the second and last equation of model (2), it follows that for $t\ge {\varsigma }_0\left(\epsilon \right)+\tau$,

$$\begin{array}{cc}\hfill S_2\left(t\right)& =q{e}^{-d_1\tau }\beta D(t-\tau )v(t-\tau )T(t-\tau )-q{d}_{1}I\left(t\right)-pbZ\left(t\right)\hfill \\ & \le q{e}^{-d_1\tau }\beta (D_0+\epsilon )(v_{max}+\epsilon )(T_0+\epsilon )-{\mu }_2{S}_2\left(t\right),\hfill \end{array}$$

which implies that

$$\underset{t\to +\infty }{lim\; sup}{S}_2\left(t\right)\le \frac{q{e}^{-d_1\tau }\beta (D_0+\epsilon )(v_{max}+\epsilon )(T_0+\epsilon )}{{\mu }_2}.$$

Since, the above inequality holds for any sufficiently small $\epsilon >0$, we have

$$\underset{t\to +\infty }{lim\; sup}{S}_2\left(t\right)\le \frac{q{e}^{-d_1\tau }\beta D_0{v}_{max}{T}_0}{{\mu }_2},$$

which leads to ${lim\; sup}_{t\to +\infty }Z\left(t\right)\le Z_{max}$. □

Theorem 1 shows that the solution of model (2) is ultimately bounded. Biologically, it indicates that the viral load in the host varies within a finite range with the evolution of time t.

2.2. The Equilibria

It is clear that model (2) always has an infection-free equilibrium $E_0=(T_0,0,0,D_0,0)$.

By using the next generation method in [29], we can obtain the immunity-inactivated reproduction number $R_0$ of model (2) as follows. First, we define the following matrices:

$$F=\left(\begin{array}{ccc}0& e^{-d_1\tau \beta D_0{T}_0}& 0\\ 0& 0& 0\\ 0& 0& 0\end{array}\right),V=\left(\begin{array}{ccc}{d}_1& 0& 0\\ -d_{1}M& c& 0\\ 0& 0& b\end{array}\right).$$

Then, we define the following:

$$R_0=\rho \left(F{V}^{-1}\right)=\frac{e^{-d_1\tau }\beta D_0{T}_{0}M}{c}=\frac{e^{-d_1\tau }\beta \lambda {\lambda }_{1}M}{cd\gamma },$$

where $\rho \left(F{V}^{-1}\right)$ is the spectrum radius of $F{V}^{-1}$.

Suppose $(T,I,v,D,0)$ ($T>0$, $I>0$, $v>0$, $D>0$) is an immunity-inactivated equilibrium of model (2). From model (2), it is not difficult to obtain the following relationships:

$$\begin{array}{cc}\hfill & DT=\frac{e^{d_1\tau }c}{M\beta },v=\frac{d_{1}M}{c}I,\hfill \\ \hfill & T=\frac{1}{d}\left(\lambda -\beta \frac{e^{d_1\tau }c}{M\beta }v\right)=\frac{1}{d}\left(\lambda -e^{d_1\tau }{d}_{1}I\right),\hfill \\ \hfill & D=\frac{1}{\gamma }\left({\lambda }_1-{\beta }_1\beta \frac{e^{d_1\tau }c}{M\beta }v\right)=\frac{1}{\gamma }\left({\lambda }_1-e^{d_1\tau }{\beta }_1{d}_{1}I\right).\hfill \end{array}$$

(8)

From (8) and the second equation of model (2), we can ascertain that I satisfies the following equation:

$$\begin{array}{cc}\hfill F_1\left(I\right)\equiv & R_0\left(1-\frac{e^{d_1\tau }{d}_{1}I}{\lambda }\right)\left(1-\frac{e^{d_1\tau }{\beta }_1{d}_{1}I}{{\lambda }_1}\right)-1=0.\hfill \end{array}$$

(9)

Let

$${\sigma }_0=min\left\{\frac{\lambda }{e^{d_1\tau }{d}_1},\frac{{\lambda }_1}{e^{d_1\tau }{\beta }_1{d}_1}\right\}.$$

If $I>{\sigma }_0$, then $min\{T,D\}<0$. Thus, we only need to consider whether $F_1\left(I\right)=0$ has a positive root on interval $[0,{\sigma }_0]$. Clearly, $F_1\left(I\right)$ is monotonically decreasing with respect to I on the interval $[0,{\sigma }_0]$, and

$$F_1\left(0\right)=R_0-1,F\left({\sigma }_0\right)=-1<0.$$

If $R_0\le 1$, then $F_1\left(I\right)=0$ has no positive root, and if $R_0>1$, then $F_1\left(I\right)=0$ has a unique positive root $I=I_1\in (0,{\sigma }_0)$, where

$$\begin{array}{cc}\hfill & I_1=\frac{M\beta ({\lambda }_1+\lambda {\beta }_1)-\sqrt{\Delta }}{2d_{1}M\beta {\beta }_1{e}^{d_1\tau }}=\frac{2cd\gamma (R_0-1)}{d_1[M\beta ({\lambda }_1+\lambda {\beta }_1)+\sqrt{\Delta }]},\hfill \\ \hfill & \Delta =M^2{\beta }^2{({\lambda }_1-\lambda {\beta }_1)}^2+4M{\beta }_1\beta cd\gamma e^{d_1\tau }.\hfill \end{array}$$

Thus, if $R_0>1$, model (2) has a unique immunity-inactivated equilibrium $E_1=(T_1,I_1,v_1,D_1,0)$, where

$$\begin{array}{cc}\hfill & T_1=\frac{1}{d}(\lambda -e^{d_1\tau }{d}_1{I}_1)=\frac{M\beta (\lambda {\beta }_1-{\lambda }_1)+\sqrt{\Delta }}{2dM\beta {\beta }_1}{D}_1=\frac{e^{d_1\tau }c}{M\beta T_1},v_1=\frac{d_{1}M}{c}{I}_1.\hfill \end{array}$$

Define the immune-activated reproduction number as follows:

$$R_1=\frac{2q{e}^{-d_1\tau }\beta \lambda {\lambda }_{1}M}{2qcd\gamma +b{d}_1[M\beta ({\lambda }_1+\lambda {\beta }_1)+\sqrt{▵}]}.$$

Then, it is not difficult to obtain the following lemma.

Lemma 1.

Assume that $R_0>1$, then the following statements are true:

(i) Clearly,

$$R_1=\frac{2qcd\gamma R_0}{2qcd\gamma +b{d}_1[M\beta ({\lambda }_1+\lambda {\beta }_1)+\sqrt{▵}]}<R_0.$$

(ii) Note that

$$I_1-\frac{b}{q}=\frac{2qcd\gamma (R_0-1)-b{d}_1[M\beta ({\lambda }_1+\lambda {\beta }_1)+\sqrt{▵}]}{q{d}_1[M\beta ({\lambda }_1+\lambda {\beta }_1)+\sqrt{▵}]}=\Theta (R_1-1),$$

where

$$\Theta =\frac{2qcd\gamma +b{d}_1[M\beta ({\lambda }_1+\lambda {\beta }_1)+\sqrt{▵}]}{q{d}_1[M\beta ({\lambda }_1+\lambda {\beta }_1)+\sqrt{▵}]}>0.$$

Thus, $R_1-1$ has the same sign as $I_1-\frac{b}{q}$.

Suppose $(\tilde{T},\tilde{I},\tilde{v},\tilde{D},\tilde{Z})$ ($\tilde{T}>0$, $\tilde{I}>0$, $\tilde{v}>0$, $\tilde{D}>0$, $\tilde{Z}>0$) is an immunity-activated equilibrium (positive equilibrium) of model (2). From model (2), it is not difficult to obtain the following relationships:

$$\begin{array}{cc}\hfill & \tilde{I}=\frac{b}{q},\tilde{v}=\frac{d_{1}Mb}{cq},\hfill \\ \hfill & \tilde{T}=\frac{1}{d}(\lambda -e^{d_1\tau }{d}_1\tilde{I}-e^{d_1\tau }p\tilde{I}\tilde{Z})=\frac{1}{d}\left(\lambda -\frac{e^{d_1\tau }{d}_{1}b}{q}-\frac{e^{d_1\tau }pb}{q}\tilde{Z}\right):={\mathrm{\Gamma }}_1\left(\tilde{Z}\right),\hfill \\ \hfill & \tilde{D}=\frac{1}{\gamma }({\lambda }_1-e^{d_1\tau }{\beta }_1{d}_1\tilde{I}-e^{d_1\tau }{\beta }_{1}p\tilde{I}\tilde{Z})=\frac{1}{\gamma }\left({\lambda }_1-\frac{e^{d_1\tau }{\beta }_1{d}_{1}b}{q}-\frac{e^{d_1\tau }{\beta }_{1}pb}{q}\tilde{Z}\right):={\mathrm{\Gamma }}_2\left(\tilde{Z}\right).\hfill \end{array}$$

(10)

Define the following condition:

$$\left(\mathrm{H}\right)\frac{b}{q}<{\sigma }_0.$$

Obviously, condition (H) is necessary to ensure that $\tilde{T}>0$ and $\tilde{D}>0$. Thus, if condition (H) does not hold, model (2) has no positive equilibrium.

If condition (H) holds, let

$$\begin{array}{cc}\hfill & {\sigma }_1=min\left\{\frac{q}{e^{d_1\tau }pb}\left(\lambda -\frac{e^{d_1\tau }{d}_{1}b}{q}\right),\frac{q}{e^{d_1\tau }{\beta }_{1}pb}\left({\lambda }_1-\frac{e^{d_1\tau }{\beta }_1{d}_{1}b}{q}\right)\right\},\hfill \\ \hfill & \widehat{R_1}=R_0\left(1-\frac{e^{d_1\tau }{d}_{1}b}{\lambda q}\right)\left(1-\frac{e^{d_1\tau }{\beta }_1{d}_{1}b}{q{\lambda }_1}\right),\hfill \end{array}$$

then ${\sigma }_1>0$ and $\widehat{R_1}>0$. From (8) and the second equation of model (2), we can ascertain that $\tilde{Z}$ satisfies the following equation:

$$\begin{array}{cc}\hfill F_2\left(\tilde{Z}\right)\equiv & \frac{e^{-d_1\tau }\beta M}{cd\gamma }\left(\lambda -\frac{e^{d_1\tau }{d}_{1}b}{q}-\frac{e^{d_1\tau }pb}{q}\tilde{Z}\right)\left({\lambda }_1-\frac{e^{d_1\tau }{\beta }_1{d}_{1}b}{q}-\frac{e^{d_1\tau }{\beta }_{1}pb}{q}\tilde{Z}\right)-1-\frac{p}{d_1}\tilde{Z}\hfill \\ \hfill =& 0.\hfill \end{array}$$

If $\tilde{Z}>{\sigma }_1$, then $min\{\tilde{T},\tilde{D}\}<0$. Thus, we only need to consider whether $F_2(\tilde{Z})=0$ has a positive root on interval $[0,{\sigma }_1]$. Clearly, if condition (H) holds, then $F_2(\tilde{Z})$ is monotonically decreasing with respect to $\tilde{Z}$ on the interval $[0,{\sigma }_1]$, and

$$F_2\left(0\right)=\widehat{R_1}-1,F_2\left({\sigma }_1\right)=-1-\frac{p}{d_1}{\sigma }_1<0.$$

If $\widehat{R_1}\le 1$, then $F_2(\tilde{Z})=0$ has no positive root, and if $\widehat{R_1}>1$, then $F_2(\tilde{Z})=0$ has a unique positive root $Z=Z^{*}\in (0,{\sigma }_1)$.

Therefore, if condition (H) holds and $\widehat{R_1}>1$, then model (2) has a unique immunity-activated equilibrium $E^{*}=(T^{*},I^{*},V^{*},D^{*},Z^{*})$, where

$$\begin{array}{cc}\hfill & T^{*}={\mathrm{\Gamma }}_1\left(Z^{*}\right),I^{*}=\frac{b}{q},v^{*}=\frac{d_{1}Mb}{cq},D^{*}={\mathrm{\Gamma }}_2\left(Z^{*}\right).\hfill \end{array}$$

(11)

Note that

$$M\beta ({\lambda }_1+\lambda {\beta }_1)+\sqrt{▵}\ge M\beta ({\lambda }_1+\lambda {\beta }_1)+M\beta |{\lambda }_1-\lambda {\beta }_1|=2M\beta max\{{\lambda }_1,\lambda {\beta }_1\}.$$

If $R_1>1$, then $2q{e}^{-d_1\tau }\beta \lambda {\lambda }_{1}M>2b{d}_{1}M\beta max\{{\lambda }_1,\lambda {\beta }_1\}$, which implies that condition (H) holds. From Lemma 1, if $R_1>1$, then

$$I_1>\frac{b}{q},\widehat{R_1}-1=F_1\left(\frac{b}{q}\right)>0.$$

Thus, if $R_1>1$, then condition (H) holds and $\widehat{R_1}>1$, then model (2) has a unique immunity-activated equilibrium $E^{*}=(T^{*},I^{*},v^{*},D^{*},Z^{*})$.

Based on the above discussion, we have the following results.

Theorem 2.

The following statements are true:

(i)

Model (2) always has an infection-free equilibrium $E_0=(T_0,0,0,D_0,0)$.

(ii)

If $R_0>1$, then model (2) has a unique immunity-inactivated equilibrium $E_1=(T_1,I_1,v_1,D_1,0)$.

(iii)

If $R_1>1$, then model (2) has a unique immunity-activated equilibrium $E^{*}=(T^{*},I^{*},v^{*},D^{*},Z^{*})$.

3. Local Stability

In this section, we will study the local stability of the infection-free equilibrium $E_0$, the immunity-inactivated equilibrium $E_1$, and the immunity-activated equilibrium $E^{*}$.

3.1. Local Stability of the Infection-Free Equilibrium $E_0$

Theorem 3.

The following statements are true:

(i)

If $R_0<1$, then the infection-free equilibrium $E_0$ of model (2) is locally asymptotically stable.

(ii)

If $R_0>1$, then the infection-free equilibrium $E_0$ of model (2) is unstable.

Proof. 

The characteristic equation of model (2) at the infection-free equilibrium $E_0$ is given by the following:

$$(\rho +b)(\rho +\gamma )(\rho +d)[(\rho +c)(\rho +d_1)-d_{1}M{e}^{-d_1\tau }\beta T_0{D}_0{e}^{-\rho \tau }]=0.$$

(12)

We can see that Equation (12) has 3 negative real roots ${\rho }_1=-b<0$, ${\rho }_2=-\gamma <0$, ${\rho }_3=-d<0$, and the other roots are determined by the following equation:

$$L_0\left(\rho \right)\equiv (\rho +c)(\rho +d_1)-d_{1}M{e}^{-d_1\tau }\beta T_0{D}_0{e}^{-\rho \tau }=0,$$

(13)

which is equivalent to

$$(\rho +c)(\rho +d_1)=d_{1}M{e}^{-d_1\tau }\beta T_0{D}_0{e}^{-\rho \tau }.$$

(14)

Assume that $\rho =x_0+i{y}_0$ is a root of Equation (14). We will show that $x_0<0$ if $R_0<1$. Otherwise, if $x_0\ge 0$, then we have the following:

$$|(\rho +c)(\rho +d_1)|=|\rho +c\left|\right|\rho +d_1|\ge c{d}_1.$$

Note that

$$c{d}_1\le \left|(\rho +c)(\rho +d_1)\right|=\left|d_{1}M{e}^{-d_1\tau }\beta T_0{D}_0{e}^{-\rho \tau }\right|\le d_{1}M{e}^{-d_1\tau }\beta T_0{D}_0=c{d}_1{R}_0<c{d}_1,$$

which is a contradiction. Thus, all roots of Equation (12) have negative real parts if $R_0<1$. This implies that the infection-free equilibrium $E_0$ is locally asymptotically stable.

Assume that $R_0>1$. We note that

$$L_0\left(0\right)=c{d}_1(1-R_0)<0,\underset{\rho \to +\infty }{lim}{L}_0\left(\rho \right)=+\infty .$$

Thus, there exists at least a positive real constant ${\rho }^{*}$ such that $L_1\left({\rho }^{*}\right)=0$. This implies that Equation (12) has a positive root, and then the infection-free equilibrium $E_0$ is unstable. □

3.2. Local Stability of the Immunity-Inactivated Equilibrium $E_1$

Theorem 4.

The following statements are true:

(i)

If $R_0>1>R_1$, then the immunity-inactivated equilibrium $E_1$ of model (2) is locally asymptotically stable.

(ii)

If $R_1>1$, then the immunity-inactivated equilibrium $E_1$ of model (2) is unstable.

Proof. 

For simplicity of presentation, we define the following matrices:

$$\begin{array}{cc}\hfill & {\mathbf{A}}_{\mathbf{1}}=\left(\begin{array}{cccc}\rho +\beta D_1{v}_1+d& 0& \beta D_1{T}_1& \beta v_1{T}_1\\ -\beta D_1{v}_1{e}^{-d_1\tau }{e}^{-\rho \tau }& \rho +d_1& -\beta D_1{T}_1{e}^{-d_1\tau }{e}^{-\rho \tau }& -\beta v_1{T}_1{e}^{-d_1\tau }{e}^{-\rho \tau }\\ 0& -d_{1}M& \rho +c& 0\\ \beta {\beta }_1{D}_1{v}_1& 0& \beta {\beta }_1{D}_1{T}_1& \rho +\beta {\beta }_1{v}_1{T}_1+\gamma \end{array}\right),\hfill \end{array}$$

$$\begin{array}{cc}\hfill & {\mathbf{A}}_{\mathbf{2}}=\left(\begin{array}{cccc}\rho +d& (\rho +d_1)e^{d_1\tau }{e}^{\rho \tau }& 0& 0\\ -\beta D_1{v}_1{e}^{-d_1\tau }{e}^{-\rho \tau }& \rho +d_1& -\beta D_1{T}_1{e}^{-d_1\tau }{e}^{-\rho \tau }& -\beta v_1{T}_1{e}^{-d_1\tau }{e}^{-\rho \tau }\\ 0& -d_{1}M& \rho +c& 0\\ 0& (\rho +d_1){\beta }_1{e}^{d_1\tau }{e}^{\rho \tau }& 0& \rho +\gamma \end{array}\right).\hfill \end{array}$$

The characteristic equation of model (2) at the immunity-inactivated equilibrium $E_1$ is given as follows:

$$(\rho +b-q{I}_1)\times \det \left({\mathbf{A}}_{\mathbf{1}}\right)=(\rho +b-q{I}_1)\times \det \left({\mathbf{A}}_{\mathbf{2}}\right)=0.$$

(15)

We can see that Equation (15) has one real root ${\tilde{\rho }}_1=q{I}_1-b$, and the other roots are determined by the following equation:

$$\begin{array}{cc}\hfill & (\rho +d)\left[(\rho +d_1)(\rho +c)(\rho +\gamma )-d_{1}M\beta D_1{T}_1{e}^{-d_1\tau }{e}^{-\rho \tau }(\rho +\gamma )+\beta {\beta }_1{v}_1{T}_1(\rho +d_1)(\rho +c)\right]\hfill \\ \hfill & +\beta D_1{v}_1(\rho +d_1)(\rho +c)(\rho +\gamma )=0.\hfill \end{array}$$

(16)

Assume that $R_0>1>R_1$. From Lemma 1, we have ${\tilde{\rho }}_1=q{I}_1-b<0$. Note that as $D_1{T}_1=\frac{e^{d_1\tau }c}{M\beta }$, then Equation (16) can be rewritten as follows:

$$\left(\rho +d_1\right)\left(\rho +c\right)\left(1+\frac{\beta {\beta }_1{v}_1{T}_1}{\rho +\gamma }+\frac{\beta D_1{v}_1}{\rho +d}\right)=d_{1}M\beta D_1{T}_1{e}^{-d_1\tau }{e}^{-\rho \tau }=d_{1}c{e}^{-\rho \tau }.$$

(17)

Clearly, $\rho =0$ is not a root of Equation (17). Assume that $\rho =x_1+i{y}_1$$(x_1^2+y_1^2>0)$ is a root of Equation (17). We will show that $x_1<0$. Otherwise, if $x_1\ge 0$, then we have $|\rho +d_1|>d_1$, $|\rho +c|>c$ and

$$\begin{array}{cc}\hfill {\mathrm{\Xi }}_1:=& \left|1+\frac{\beta {\beta }_1{v}_1{T}_1}{\rho +\gamma }+\frac{\beta D_1{v}_1}{\rho +d}\right|\hfill \\ \hfill =& \left|1+\frac{\beta {\beta }_1{v}_1{T}_1\left[(x_1+\gamma )-i{y}_1\right]}{{(x_1+\gamma )}^2+y_1^2}+\frac{\beta D_1{v}_1\left[(x_1+d)-i{y}_1\right]}{{(x_1+d)}^2+y_1^2}\right|\hfill \\ \hfill >& 1.\hfill \end{array}$$

Then, from Equation (17), we have

$$\begin{array}{cc}\hfill d_1{c}_1\ge & |d_{1}c{e}^{-\rho \tau }|\hfill \\ \hfill =& \left|\left(\rho +d_1\right)\left(\rho +c\right)\left(1+\frac{\beta {\beta }_1{v}_1{T}_1}{\rho +\gamma }+\frac{\beta D_1{v}_1}{\rho +d}\right)\right|\hfill \\ \hfill =& |\rho +d_1\left|\right|\rho +c|{\mathrm{\Xi }}_1\hfill \\ \hfill >& d_1{c}_1,\hfill \end{array}$$

which is a contradiction. Thus, all roots of Equation (15) have negative real parts if $R_0>1>R_1$. This implies that the immunity-inactivated equilibrium $E_1$ is locally asymptotically stable.

From Lemma 1, if $R_1>1$, then ${\tilde{\rho }}_1=q{I}_1-b>0$, which implies that the immunity-inactivated equilibrium $E_1$ is unstable. □

3.3. Local Stability of the Immunity-Activated Equilibrium $E^{*}$

We give the following lemma, which will be used in the proof of Theorem 5.

Lemma 2.

Assume $\rho =x+iy$ ($x\ge 0$, $x^2+y^2>0$). For any ${\theta }_1>0,{\theta }_2>0$, it the case that

$$|\rho (\rho +{\theta }_1)+{\theta }_2|\ge {\theta }_1\left|\rho \right|.$$

Proof. 

We note that

$$\rho (\rho +{\theta }_1)+{\theta }_2=(x+iy)(x+{\theta }_1+iy)+{\theta }_2=x(x+{\theta }_1)-y^2+{\theta }_2+iy(2x+{\theta }_1).$$

Thus, we have

$$\begin{array}{cc}\hfill |\rho (\rho +{\theta }_1)+{\theta }_2{|}^2=& {(x^2-y^2+{\theta }_2+x{\theta }_1)}^2+y^2{(2x+{\theta }_1)}^2\hfill \\ \hfill =& {(x^2-y^2+{\theta }_2)}^2+x^2{\theta }_1^2+2x{\theta }_1(x^2-y^2+{\theta }_2)+y^2{\theta }_1^2+y^2(4x^2+4x{\theta }_1)\hfill \\ \hfill \ge & (x^2+y^2){\theta }_1^2={\left|\rho \right|}^2{\theta }_1^2.\hfill \end{array}$$

□

Theorem 5.

If $R_1>1$, then the immunity-activated equilibrium $E^{*}$ of model (2) is locally asymptotically stable.

Proof. 

For simplicity of presentation, we define the following matrices:

$$\begin{array}{cc}\hfill & {\mathbf{A}}_{\mathbf{3}}=\left(\begin{array}{ccccc}\rho +\beta D^{*}{v}^{*}+d& 0& \beta D^{*}{T}^{*}& \beta v^{*}{T}^{*}& 0\\ -\beta D^{*}{v}^{*}{e}^{-d_1\tau }{e}^{-\rho \tau }& \rho +d_1+p{Z}^{*}& -\beta D^{*}{T}^{*}{e}^{-d_1\tau }{e}^{-\rho \tau }& -\beta v^{*}{T}^{*}{e}^{-d_1\tau }{e}^{-\rho \tau }& p{I}^{*}\\ 0& -d_{1}M& \rho +c& 0& 0\\ \beta {\beta }_1{D}^{*}{v}^{*}& 0& \beta {\beta }_1{D}^{*}{T}^{*}& \rho +\beta {\beta }_1{v}^{*}{T}^{*}+\gamma & 0\\ 0& -q{Z}^{*}& 0& 0& \rho \end{array}\right),\hfill \end{array}$$

$$\begin{array}{cc}\hfill & {\mathbf{A}}_{\mathbf{4}}=\left(\begin{array}{ccccc}\rho +d& (\rho +d_1+p{Z}^{*})e^{d_1\tau }{e}^{\rho \tau }& 0& 0& p{I}^{*}{e}^{d_1\tau }{e}^{\rho \tau }\\ -\beta D^{*}{v}^{*}{e}^{-d_1\tau }{e}^{-\rho \tau }& \rho +d_1+p{Z}^{*}& -\beta D^{*}{T}^{*}{e}^{-d_1\tau }{e}^{-\rho \tau }& -\beta v^{*}{T}^{*}{e}^{-d_1\tau }{e}^{-\rho \tau }& p{I}^{*}\\ 0& -d_{1}M& \rho +c& 0& 0\\ 0& {\beta }_1(\rho +d_1+p{Z}^{*})e^{d_1\tau }{e}^{\rho \tau }& 0& \rho +\gamma & {\beta }_{1}p{I}^{*}{e}^{d_1\tau }{e}^{\rho \tau }\\ 0& -q{Z}^{*}& 0& 0& \rho \end{array}\right),\hfill \end{array}$$

$$\begin{array}{cc}\hfill & {\mathbf{A}}_{\mathbf{5}}=\left(\begin{array}{cccc}\rho +d& (\rho +d_1+p{Z}^{*})e^{d_1\tau }{e}^{\rho \tau }& 0& 0\\ -\beta D^{*}{v}^{*}{e}^{-d_1\tau }{e}^{-\rho \tau }& \rho +d_1+p{Z}^{*}& -\beta D^{*}{T}^{*}{e}^{-d_1\tau }{e}^{-\rho \tau }& -\beta v^{*}{T}^{*}{e}^{-d_1\tau }{e}^{-\rho \tau }\\ 0& -d_{1}M& \rho +c& 0\\ 0& {\beta }_1(\rho +d_1+p{Z}^{*})e^{d_1\tau }{e}^{\rho \tau }& 0& \rho +\gamma \end{array}\right),\hfill \end{array}$$

$$\begin{array}{cc}\hfill & {\mathbf{A}}_{\mathbf{6}}=\left(\begin{array}{ccc}\rho +d& 0& e^{d_1\tau }{e}^{\rho \tau }\\ -\beta D^{*}{v}^{*}{e}^{-d_1\tau }{e}^{-\rho \tau }& -\beta v^{*}{T}^{*}{e}^{-d_1\tau }{e}^{-\rho \tau }& 1\\ 0& \rho +\gamma & {\beta }_1{e}^{d_1\tau }{e}^{\rho \tau }\end{array}\right).\hfill \end{array}$$

Note that if $I^{*}=\frac{b}{q}$, then the characteristic equation of model (2) at the immunity-activated equilibrium $E^{*}$ is given as follows:

$$\det \left({\mathbf{A}}_{\mathbf{3}}\right)=\det \left({\mathbf{A}}_{\mathbf{4}}\right)=\rho \times \det \left({\mathbf{A}}_{\mathbf{5}}\right)-q{Z}^{*}\left(\rho +c\right)p{I}^{*}\times \det \left({\mathbf{A}}_{\mathbf{6}}\right)=0.$$

(18)

Equation (18) can be rewritten as follows:

$$\begin{array}{cc}\hfill & \left[\rho (\rho +d_1+p{Z}^{*})+q{Z}^{*}p{I}^{*}\right](\rho +c)\left[(\rho +d)(\rho +\gamma )+(\rho +d)\beta {\beta }_1{v}^{*}{T}^{*}+\beta D^{*}{v}^{*}(\rho +\gamma )\right]\hfill \\ \hfill & -d_{1}M\beta D^{*}{T}^{*}(\rho +d)(\rho +\gamma )\rho e^{-d_1\tau }{e}^{-\rho \tau }=0,\hfill \end{array}$$

(19)

which is equivalent to

$$\left[\rho (\rho +d_1+p{Z}^{*})+q{Z}^{*}p{I}^{*}\right](\rho +c)\left\{1+\frac{\beta {\beta }_1{v}^{*}{T}^{*}}{\rho +\gamma }+\frac{\beta D^{*}{v}^{*}}{\rho +d}\right\}=d_{1}M\beta D^{*}{T}^{*}\rho e^{-d_1\tau }{e}^{-\rho \tau }.$$

(20)

From the second equation of model (2), we have $d_{1}M\beta D^{*}{T}^{*}{e}^{-d_1\tau }=c(d_1+p{Z}^{*})$. Then, Equation (20) can be rewritten as follows:

$$1+\frac{\beta {\beta }_1{v}^{*}{T}^{*}}{\rho +\gamma }+\frac{\beta D^{*}{v}^{*}}{\rho +d}=\frac{c(d_1+p{Z}^{*})\rho e^{-\rho \tau }}{[\rho (\rho +d_1+p{Z}^{*})+q{Z}^{*}p{I}^{*}](\rho +c)}.$$

(21)

Clearly, $\rho =0$ is not a root of Equation (21). Assume that $\rho =x_2+i{y}_2$$(x_2^2+y_2^2>0)$ is a root of Equation (21). We will show that $x_2<0$. Otherwise, if $x_2\ge 0$, then we have the following:

$$\begin{array}{cc}\hfill {\mathrm{\Xi }}_2:=& \left|1+\frac{\beta {\beta }_1{v}^{*}{T}^{*}}{\rho +\gamma }+\frac{\beta D^{*}{v}^{*}}{\rho +d}\right|\hfill \\ \hfill =& \left|1+\frac{\beta {\beta }_1{v}^{*}{T}^{*}\left[(x_2+\gamma )-i{y}_2\right]}{{(x_2+\gamma )}^2+y_2^2}+\frac{\beta D^{*}{v}^{*}\left[(x_2+d)-i{y}_2\right]}{{(x_2+d)}^2+y_2^2}\right|\hfill \\ \hfill >& 1.\hfill \end{array}$$

In addition, from Lemma 2, we have the following:

$$\begin{array}{cc}\hfill {\mathrm{\Xi }}_3:=& \left|\frac{c(d_1+p{Z}^{*})\rho e^{-\rho \tau }}{[\rho (\rho +d_1+p{Z}^{*})+q{Z}^{*}p{I}^{*}](\rho +c)}\right|\hfill \\ \hfill \le & \left|\frac{(d_1+p{Z}^{*})\rho }{[\rho (\rho +d_1+p{Z}^{*})+q{Z}^{*}p{I}^{*}]}\right|\hfill \\ \hfill \le & 1.\hfill \end{array}$$

Thus, ${\mathrm{\Xi }}_2>{\mathrm{\Xi }}_3$, which contradicts Equation (21). This implies that all roots of Equation (21) have negative real parts if $R_1>1$, and the immunity-activated equilibrium $E^{*}$ is locally asymptotically stable. □

4. Global Stability

In this section, we will study the global stability of the infection-free equilibrium $E_0$, the immunity-inactivated equilibrium $E_1$, and the immunity-activated equilibrium $E^{*}$ by constructing appropriate Lyapunov functionals. Some construction techniques of Lyapunov functionals (functions) can be found in [22,23,24,25,26,27,30,31,32,33] and the references therein.

For convenience, let

$$f\left(x\right)=x-1-lnx,x>0.$$

The function $f\left(x\right)\ge 0$ for any $x>0$ and $f\left(x\right)=0$ if and only if $x=1$.

4.1. Global Stability of the Infection-Free Equilibrium $E_0$

Theorem 6.

If $R_0\le 1$, then the infection-free equilibrium $E_0$ of model (2) is globally asymptotically stable in $X_1:=\{\varphi \in X^{+}|0<{\varphi }_1\left(0\right)\le T_0,0<{\varphi }_4\left(0\right)\le D_0\}$.

Proof. 

According to Theorem 3, we only need to prove that the infection-free equilibrium $E_0$ is globally attractive. It is easy to show that the set $X_1$ is positively invariant for model (2).

Define the Lyapunov functional as follows:

$$\begin{array}{cc}\hfill U_0=& T_{0}f\left(\frac{{\varphi }_1\left(0\right)}{T_0}\right)+e^{d_1\tau }(1+{\beta }_1){\varphi }_2\left(0\right)+\frac{e^{d_1\tau }(1+{\beta }_1)}{M}{\varphi }_3\left(0\right)+D_{0}f\left(\frac{{\varphi }_4\left(0\right)}{D_0}\right)\hfill \\ \hfill & +(1+{\beta }_1)\frac{e^{d_1\tau }p}{q}{\varphi }_5\left(0\right)+(1+{\beta }_1)\beta {\int }_{-\tau }^0{\varphi }_1\left(s\right){\varphi }_3\left(s\right){\varphi }_4\left(s\right)ds.\hfill \end{array}$$

It is clear that $U_0$ is continuous on $X_1$ and satisfies the condition (ii) of Lemma 3.1 in [34] on $\partial X_1=\overline{X_1}\setminus X_1$.

In calculating the derivative of $U_0$ along the solution ${(T\left(t\right),I\left(t\right),v\left(t\right),D\left(t\right),Z\left(t\right))}^T$ of model (2), it follows that for $t\ge 0$,

$$\begin{array}{cc}\hfill {\dot{U}}_0=& \left(1-\frac{T_0}{T\left(t\right)}\right)(\lambda -\beta D\left(t\right)v\left(t\right)T\left(t\right)-dT\left(t\right))+(1+{\beta }_1)\beta D(t-\tau )v(t-\tau )T(t-\tau )\hfill \\ \hfill & -e^{d_1\tau }(1+{\beta }_1)d_{1}I\left(t\right)-e^{d_1\tau }(1+{\beta }_1)pI\left(t\right)Z\left(t\right)+\frac{e^{d_1\tau }(1+{\beta }_1)}{M}(d_{1}MI\left(t\right)-cv\left(t\right))\hfill \\ \hfill & +\left(1-\frac{D_0}{D\left(t\right)}\right)({\lambda }_1-{\beta }_1\beta D\left(t\right)v\left(t\right)T\left(t\right)-\gamma D\left(t\right))\hfill \\ \hfill & +(1+{\beta }_1)\frac{e^{d_1\tau }p}{q}(qI\left(t\right)Z\left(t\right)-bZ\left(t\right))\hfill \\ \hfill & +(1+{\beta }_1)\beta [D\left(t\right)v\left(t\right)T\left(t\right)-D(t-\tau )v(t-\tau )T(t-\tau )]\hfill \\ \hfill =& \left(1-\frac{T_0}{T\left(t\right)}\right)(\lambda -dT\left(t\right))+\beta D\left(t\right)v\left(t\right)T_0-\frac{e^{d_1\tau }(1+{\beta }_1)c}{M}v\left(t\right)\hfill \\ \hfill & +\left(1-\frac{D_0}{D\left(t\right)}\right)({\lambda }_1-\gamma D\left(t\right))+{\beta }_1\beta D_{0}v\left(t\right)T\left(t\right)-\frac{e^{d_1\tau }(1+{\beta }_1)pb}{q}Z\left(t\right).\hfill \end{array}$$

Note that

$$T\left(t\right)\le \frac{\lambda }{d}=T_0,D\left(t\right)\le \frac{{\lambda }_1}{\gamma }=D_0,R_0=\frac{e^{-d_1\tau }\beta T_0{D}_{0}M}{c},$$

then we have

$$\begin{array}{cc}\hfill {\dot{U}}_0\le & -\frac{d}{T\left(t\right)}{(T\left(t\right)-T_0)}^2+\beta D_0{T}_{0}v\left(t\right)-\frac{e^{d_1\tau }(1+{\beta }_1)c}{M}v\left(t\right)\hfill \\ \hfill & -\frac{\gamma }{D\left(t\right)}{(D\left(t\right)-D_0)}^2+{\beta }_1\beta D_0{T}_{0}v\left(t\right)-\frac{e^{d_1\tau }(1+{\beta }_1)pb}{q}Z\left(t\right)\hfill \\ \hfill =& -\frac{d}{T\left(t\right)}{(T\left(t\right)-T_0)}^2-\frac{\gamma }{D\left(t\right)}{(D\left(t\right)-D_0)}^2\hfill \\ \hfill & +\frac{e^{d_1\tau }(1+{\beta }_1)c}{M}(R_0-1)v\left(t\right)-\frac{e^{d_1\tau }(1+{\beta }_1)pb}{q}Z\left(t\right).\hfill \end{array}$$

(22)

It follows from $R_0\le 1$ that ${\dot{U}}_0\le 0$ for $t\ge 0$. Thus, the infection-free equilibrium $E_0$ is stable. Moreover, ${\dot{U}}_0=0$ implies $T\left(t\right)=T_0$, $D\left(t\right)=D_0$, and $Z\left(t\right)=0$.

Let $M_0$ be the largest invariant set in the set ${\mathrm{\Omega }}_0:=\{\varphi \in \overline{X_1}:U_0<\infty \mathrm{and}{\dot{U}}_0=0\}$. From model (2) and the invariance of $M_0$, we can easily see that $M_0=\left\{E_0\right\}$. Thus, it follows from Lemma 3.1 in [34] that the infection-free equilibrium $E_0$ is globally attractive. □

4.2. Global Stability of the Immunity-Inactivated Equilibrium $E_1$

In order to prove Theorems 7 and 8, we need the following Lemma 3. In this subsection, we assume that $R_0>1$.

Let ${\theta }_1>1$ be a positive constant. Note that

$$R_0=\frac{e^{-d_1\tau }\beta \lambda {\lambda }_{1}M}{cd\gamma },T_1=\frac{\lambda }{d+\beta D_1{v}_1}<\frac{\lambda }{d},D_1=\frac{{\lambda }_1}{\gamma +{\beta }_1\beta T_1{v}_1}<\frac{{\lambda }_1}{\gamma },$$

then there exists a $\delta >0$ such that

$$T_1<\frac{\lambda }{d+\beta {\theta }_1{D}_0\delta }:=T_1^{\delta },D_1<\frac{{\lambda }_1}{\gamma +{\beta }_1\beta {\theta }_1{T}_0\delta }:=D_1^{\delta },I_1^{\delta }:=\frac{c}{2d_{1}M}\delta <\frac{b}{q},$$

(23)

and

$$e^{-d_1\tau }\beta D_1^{\delta }{T}_1^{\delta }\frac{M}{c}>1.$$

(24)

Clearly, $\delta <v_1\le v_{max}$. For convenience, we define the following:

$${\alpha }_0=-\frac{1}{c}ln\left(\frac{\delta }{4{\theta }_1{v}_{max}-2\delta }\right)>0,{\alpha }_1=max\{\widehat{{\alpha }_1},\widetilde{{\alpha }_1}\},$$

where

$$\begin{array}{cc}\hfill \widehat{{\alpha }_1}=& -\frac{1}{d+\frac{3}{4}\beta {\theta }_1{D}_0\delta }ln\left(1-\frac{d+\frac{3}{4}\beta {\theta }_1{D}_0\delta }{d+\beta {\theta }_1{D}_0\delta }\right)>0,\hfill \\ \hfill \widetilde{{\alpha }_1}=& -\frac{1}{\gamma +\frac{3}{4}{\beta }_1\beta {\theta }_1{T}_0\delta }ln\left(1-\frac{\gamma +\frac{3}{4}{\beta }_1\beta {\theta }_1{T}_0\delta }{\gamma +{\beta }_1\beta {\theta }_1{T}_0\delta }\right)>0.\hfill \end{array}$$

Furthermore, from (24), there exists a positive constant ${\alpha }_2>0$ such that

$$e^{-d_1\tau }\beta D_1^{\delta }{T}_1^{\delta }\frac{d_{1}M}{c}(1-e^{-c{\alpha }_2})>d_1+p{\theta }_1{Z}_{max}{e}^{-(b-q{I}_1^{\delta }){\alpha }_2}.$$

(25)

Using the methods and techniques in [35,36], we can obtain the following conclusion.

Lemma 3.

If $R_0>1$, then the solution ${(T\left(t\right),I\left(t\right),v\left(t\right),D\left(t\right),Z\left(t\right))}^T$ of model (2) with with any $\varphi \in X_2:=\{\varphi \in X^{+}|{\varphi }_i\left(0\right)>0,i=1,2,3,4\}$ satisfies

$$\begin{array}{cc}\hfill & \underset{t\to +\infty }{lim\; inf}T\left(t\right)\ge \frac{\lambda \gamma }{d\gamma +\beta {\lambda }_1{v}_{max}}:=T_{min},\hfill \\ \hfill & \underset{t\to +\infty }{lim\; inf}D\left(t\right)\ge \frac{{\lambda }_{1}d}{\gamma d+{\beta }_1\beta \lambda v_{max}}:=D_{min},\hfill \\ \hfill & \underset{t\to +\infty }{lim\; inf}I\left(t\right)\ge I_1^{\delta }{e}^{-(d_1+p{\theta }_1{Z}_{max})\alpha }:=I_{min},\hfill \\ \hfill & \underset{t\to +\infty }{lim\; inf}v\left(t\right)\ge \frac{d_{1}M}{c}{I}_{min}=\frac{\delta }{2}{e}^{-(d_1+p{\theta }_1{Z}_{max})\alpha }:=v_{min},\hfill \end{array}$$

where $\alpha =max\{{\alpha }_0+{\alpha }_1,{\alpha }_2\}+\tau$.

Proof. 

Let ${(T\left(t\right),I\left(t\right),v\left(t\right),D\left(t\right),Z\left(t\right))}^T$ be the solution of model (2) with any initial function $\varphi \in X_2$. Clearly, $T\left(t\right)>0$, $I\left(t\right)>0$, $v\left(t\right)>0$ and $D\left(t\right)>0$. By Theorem 1, we can easily obtain the following:

$$\underset{t\to +\infty }{lim\; inf}T\left(t\right)\ge \frac{\lambda }{d+\beta D_0{v}_{max}}=T_{min}>0,\underset{t\to +\infty }{lim\; inf}D\left(t\right)\ge \frac{{\lambda }_1}{\gamma +{\beta }_1\beta T_0{v}_{max}}=D_{min}>0.$$

From (4), it follows that there exists a $T_1>\tau$ such that for $t\ge T_1$,

$$T\left(t\right)<{\theta }_1{T}_0,D\left(t\right)<{\theta }_1{D}_0,v\left(t\right)<{\theta }_1{v}_{max}Z\left(t\right)<{\theta }_1{Z}_{max}.$$

(26)

Claim For any $t_0\ge T_1$, it is impossible to satisfy $I\left(t\right)\le I_1^{\delta }$ for $t\ge t_0$.

In the following, let us prove that the claim is true. If the claim is not true, then there exists a $t_0\ge T_1$ such that $I\left(t\right)\le I_1^{\delta }$ for $t\ge t_0$. From the third equation of model (2), we have

$$\underset{t\to +\infty }{lim\; sup}v\left(t\right)\le \frac{d_{1}M}{c}{I}_1^{\delta }=\frac{\delta }{2}.$$

Then, from (26) and the first and fourth equations of model (2), we have the following:

$$\underset{t\to +\infty }{lim\; inf}T\left(t\right)\ge \frac{\lambda }{d+\beta {\theta }_1{D}_0\frac{\delta }{2}}>T_1^{\delta },\underset{t\to +\infty }{lim\; inf}D\left(t\right)\ge \frac{{\lambda }_1}{\gamma +{\beta }_1\beta {\theta }_1{T}_0\frac{\delta }{2}}>D_1^{\delta }.$$

(27)

From (23), we note that $q{I}_1^{\delta }<b$. From the last equation of model (2), we have the following:

$$\dot{Z}\left(t\right)\le q{I}_1^{\delta }Z\left(t\right)-bZ\left(t\right)=-(b-q{I}_1^{\delta })Z\left(t\right),$$

which implies that

$$\underset{t\to +\infty }{lim}Z\left(t\right)=0.$$

(28)

By (24), there exists a sufficiently small positive constant $\epsilon >0$ such that

$${\mathrm{\Pi }}_1:=e^{-d_1\tau }\beta D_1^{\delta }{T}_1^{\delta }-\frac{c}{d_{1}M}(d_1+p\epsilon )>0.$$

By (27) and (28), there exists a $T_2>t_0$ such that for $t\ge T_2$,

$$T\left(t\right)>T_1^{\delta },D\left(t\right)>D_1^{\delta },Z\left(t\right)<\epsilon .$$

(29)

Define the following auxiliary function:

$$A\left(t\right)=I\left(t\right)+\frac{d_1+p\epsilon }{d_{1}M}v\left(t\right)+e^{-d_1\tau }\beta {\int }_{t-\tau }^{t}D\left(s\right)v\left(s\right)T\left(s\right)ds.$$

From (29), we have for $t\ge T_2$,

$$\begin{array}{cc}\hfill \dot{A}\left(t\right)=& e^{-d_1\tau }\beta D\left(t\right)v\left(t\right)T\left(t\right)-pI\left(t\right)Z\left(t\right)+pI\left(t\right)\epsilon -\frac{c}{d_{1}M}(d_1+p\epsilon )v\left(t\right)\hfill \\ \hfill \ge & e^{-d_1\tau }\beta D_1^{\delta }{T}_1^{\delta }v\left(t\right)-\frac{c}{d_{1}M}(d_1+p\epsilon )v\left(t\right)\hfill \\ \hfill =& {\mathrm{\Pi }}_{1}v\left(t\right)\ge 0.\hfill \end{array}$$

This implies that for $t\ge T_2$, $A\left(t\right)$ is monotonically increasing with respect to t. It follows from Theorem 1 that $A\left(t\right)$ is bounded for $t\ge T_2$. Thus, there exists a positive constant $A^{*}\ge I\left(T_2\right)>0$ such that ${lim}_{t\to +\infty }A\left(t\right)=A^{*}>0$. Moreover, according to Theorem 1, we have that for $t\ge T_2$, $\ddot{A}\left(t\right)$ is also bounded. This implies that $A\left(t\right)$ is uniformly continuous for $t>T_2$. Thus, it follows from Barbalat’s lemma [37] that ${lim}_{t\to +\infty }\dot{A}\left(t\right)=0$, which implies that ${lim}_{t\to +\infty }v\left(t\right)=0$. Then, from the second equation of model (2), it is not difficult to obtain ${lim}_{t\to +\infty }I\left(t\right)=0$. Thus, ${lim}_{t\to +\infty }A\left(t\right)=0$, which is a contradiction to ${lim}_{t\to +\infty }A\left(t\right)=A^{*}>0$. This proves the claim.

By the claim, there are two cases that need to be considered:

(i)

$I\left(t\right)\ge I_1^{\delta }$ for all sufficiently large t;

(ii)

$I\left(t\right)$ oscillates about $I_1^{\delta }$ for all sufficiently large t.

Clearly, we only need to consider case (ii). Let $t_1,t_2>T_1+\tau$ be sufficiently large such that

$$I\left(t_1\right)=I\left(t_2\right)=I_1^{\delta },I\left(t\right)<I_1^{\delta }(t_1<t<t_2).$$

If $t_2-t_1\le \alpha$, from (26) and the second equation of model (2), we have for $t_1\le t\le t_2$,

$$\dot{I}\left(t\right)\ge -(d_1+pZ\left(t\right))I\left(t\right)\ge -(d_1+p{\theta }_1{Z}_{max})I\left(t\right),$$

which implies that for $t_1\le t\le t_2$,

$$I\left(t\right)\ge I_1^{\delta }{e}^{-(d_1+p{\theta }_1{Z}_{max})(t-t_1)}\ge I_1^{\delta }{e}^{-(d_1+p{\theta }_1{Z}_{max})(t_2-t_1)}\ge I_1^{\delta }{e}^{-(d_1+p{\theta }_1{Z}_{max})\alpha }=I_{min}.$$

Assume that $t_2-t_1>\alpha$. It is easy to obtain $t_1\le t\le t_1+\alpha$, $I\left(t\right)\ge I_{min}$. Then, let us prove that, for $t_1+\alpha \le t\le t_2$, $I\left(t\right)\ge I_{min}$. Otherwise, there exists a ${\alpha }^{*}\ge 0$ such that for $t_1+\alpha \le t\le t_1+\alpha +{\alpha }^{*}:=\omega$, $I\left(t\right)\ge I_{min}$, $I\left(\omega \right)=I_{min}$ and $\dot{I}\left(\omega \right)\le 0$. From the third equation of model (2), we have for $t_1\le t\le t_2$,

$$\dot{v}\left(t\right)\le d_{1}M{I}_1^{\delta }-cv\left(t\right)=\frac{c\delta }{2}-cv\left(t\right),$$

which implies that for $t_1\le t\le t_2$,

$$v\left(t\right)\le \frac{\delta }{2}+\left(v\left(t_1\right)-\frac{\delta }{2}\right)e^{-c(t-t_1)}\le \frac{\delta }{2}+\left({\theta }_1{v}_{max}-\frac{\delta }{2}\right)e^{-c(t-t_1)}.$$

(30)

From (30), we have for $t_1+{\alpha }_0\le t\le t_2$,

$$v\left(t\right)\le \frac{\delta }{2}+\left({\theta }_1{v}_{max}-\frac{\delta }{2}\right)e^{-c{\alpha }_0}=\frac{3}{4}\delta .$$

(31)

From (26) and (31), we have for $t_1+{\alpha }_0\le t\le t_2$,

$$\dot{T}\left(t\right)\ge \lambda -\left(\beta {\theta }_1{D}_0\frac{3}{4}\delta +d\right)T\left(t\right),\dot{D}\left(t\right)\ge {\lambda }_1-\left({\beta }_1\beta {\theta }_1{T}_0\frac{3}{4}\delta +\gamma \right)T\left(t\right),$$

which implies that for $t_1+{\alpha }_0\le t\le t_2$,

$$\begin{array}{cc}\hfill T\left(t\right)\ge & \frac{\lambda }{\beta {\theta }_1{D}_0\frac{3}{4}\delta +d}+\left(T(t_1+{\alpha }_0)-\frac{\lambda }{\beta {\theta }_1{D}_0\frac{3}{4}\delta +d}\right)e^{-\left(\beta {\theta }_1{D}_0\frac{3}{4}\delta +d\right)(t-t_1-{\alpha }_0)}\hfill \\ \hfill \ge & \frac{\lambda }{\beta {\theta }_1{D}_0\frac{3}{4}\delta +d}\left(1-e^{-\left(\beta {\theta }_1{D}_0\frac{3}{4}\delta +d\right)(t-t_1-{\alpha }_0)}\right),\hfill \end{array}$$

(32)

$$\begin{array}{cc}\hfill D\left(t\right)\ge & \frac{{\lambda }_1}{{\beta }_1\beta {\theta }_1{T}_0\frac{3}{4}\delta +\gamma }+\left(D(t_1+{\alpha }_0)-\frac{{\lambda }_1}{{\beta }_1\beta {\theta }_1{T}_0\frac{3}{4}\delta +\gamma }\right)e^{-\left({\beta }_1\beta {\theta }_1{T}_0\frac{3}{4}\delta +\gamma \right)(t-t_1-{\alpha }_0)}\hfill \\ \hfill \ge & \frac{{\lambda }_1}{{\beta }_1\beta {\theta }_1{T}_0\frac{3}{4}\delta +\gamma }\left(1-e^{-\left({\beta }_1\beta {\theta }_1{T}_0\frac{3}{4}\delta +\gamma \right)(t-t_1-{\alpha }_0)}\right).\hfill \end{array}$$

(33)

From (32) and (33), we have for $t_1+{\alpha }_0+{\alpha }_1\le t\le t_2$,

$$\begin{array}{cc}\hfill T\left(t\right)\ge & \frac{\lambda }{\beta {\theta }_1{D}_0\frac{3}{4}\delta +d}\left(1-e^{-\left(\beta {\theta }_1{D}_0\frac{3}{4}\delta +d\right){\alpha }_1}\right)\hfill \\ \hfill \ge & \frac{\lambda }{\beta {\theta }_1{D}_0\frac{3}{4}\delta +d}\left(1-e^{-\left(\beta {\theta }_1{D}_0\frac{3}{4}\delta +d\right)\widehat{{\alpha }_1}}\right)=T_1^{\delta },\hfill \end{array}$$

(34)

$$\begin{array}{cc}\hfill D\left(t\right)\ge & \frac{{\lambda }_1}{{\beta }_1\beta {\theta }_1{T}_0\frac{3}{4}\delta +\gamma }\left(1-e^{-\left({\beta }_1\beta {\theta }_1{T}_0\frac{3}{4}\delta +\gamma \right){\alpha }_1}\right)\hfill \\ \hfill \ge & \frac{{\lambda }_1}{{\beta }_1\beta {\theta }_1{T}_0\frac{3}{4}\delta +\gamma }\left(1-e^{-\left({\beta }_1\beta {\theta }_1{T}_0\frac{3}{4}\delta +\gamma \right)\widetilde{{\alpha }_1}}\right)=D_1^{\delta }.\hfill \end{array}$$

(35)

From the last equation of model (2), it follows that for $t_1\le t\le t_2$,

$$\dot{Z}\left(t\right)\le q{I}_1^{\delta }Z\left(t\right)-bZ\left(t\right)=-(b-q{I}_1^{\delta })Z\left(t\right),$$

which implies that

$$Z\left(t\right)\le Z\left(t_1\right)e^{-(b-q{I}_1^{\delta })(t-t_1)}.$$

(36)

From (26) and (36), we have for $t_1+{\alpha }_2\le t\le t_2$,

$$Z\left(t\right)\le {\theta }_1{Z}_{max}{e}^{-(b-q{I}_1^{\delta }){\alpha }_2}.$$

(37)

Note that for $t_1\le t\le \omega$, $I\left(t\right)\ge I_{min}$. From the third equation of model (2), we have for $t_1\le t\le \omega$,

$$\dot{v}\left(t\right)\ge d_{1}M{I}_{min}-cv\left(t\right),$$

which implies that for $t_1\le t\le \omega$,

$$v\left(t\right)\ge \frac{d_{1}M{I}_{min}}{c}+\left(v\left(t_1\right)-\frac{d_{1}M{I}_{min}}{c}\right)e^{-c(t-t_1)}\ge \frac{d_{1}M{I}_{min}}{c}(1-e^{-c(t-t_1)}).$$

(38)

From (38), we have for $t_1+{\alpha }_2\le t\le \omega$,

$$v\left(t\right)\ge \frac{d_{1}M{I}_{min}}{c}(1-e^{-c{\alpha }_2}).$$

(39)

From (25), (34), (35), (37), (39), and the second equation of model (2), we have the following:

$$\begin{array}{cc}\hfill \dot{I}\left(\omega \right)=& e^{-d_1\tau }\beta D(\omega -\tau )v(\omega -\tau )T(\omega -\tau )-d_1{I}_{min}-p{I}_{min}Z\left(\omega \right)\hfill \\ \hfill \ge & e^{-d_1\tau }\beta D_1^{\delta }{T}_1^{\delta }\frac{d_{1}M{I}_{min}}{c}(1-e^{-c{\alpha }_2})-d_1{I}_{min}-p{I}_{min}{\theta }_1{Z}_{max}{e}^{-(b-q{I}_1^{\delta }){\alpha }_2}\hfill \\ \hfill =& \left(e^{-d_1\tau }\beta D_1^{\delta }{T}_1^{\delta }\frac{d_{1}M}{c}(1-e^{-c{\alpha }_2})-d_1-p{\theta }_1{Z}_{max}{e}^{-(b-q{I}_1^{\delta }){\alpha }_2}\right)I_{min}\hfill \\ \hfill >& 0,\hfill \end{array}$$

which is a contradiction. Thus, for $t_1\le t\le t_2$, $I\left(t\right)\ge I_{min}$.

Since the interval $t_1\le t\le t_2$ is arbitrary chosen, we can conclude that $I\left(t\right)\ge I_{min}$ for all sufficiently large t. Thus, ${lim\; inf}_{t\to +\infty }I\left(t\right)\ge I_{min}$. Then, from the third equation of model (2), the following can easily be obtained: ${lim\; inf}_{t\to +\infty }v\left(t\right)\ge \frac{d_{1}M{I}_{min}}{c}=v_{min}$. □

Remark 1.

According to Lemma 3, the viruses are persistent in the host if the immunity-inactivated reproduction number $R_0>1$. Lemma 3 gives an explicit estimate of the ultimate lower bound on the viral load.

For convenience, we define the following condition:

$$\left(\mathrm{H}1\right)4\left(\frac{d^2}{\lambda }+\frac{d{\beta }_1\beta v_{max}}{d+\gamma }\right)\left(\frac{{\gamma }^2}{{\lambda }_1}+\frac{\gamma \beta v_{max}}{d+\gamma }\right)>{\beta }_1{\beta }^2{v}_{max}^2.$$

Theorem 7.

If $R_0>1>R_1$ and condition (H1) holds, then the immunity-inactivated equilibrium $E_1$ of model (2) is globally asymptotically stable in $X_2$.

Proof. 

According to Theorem 4, we only need to prove that the immunity-inactivated equilibrium $E_1$ is globally attractive. Let ${(T\left(t\right),I\left(t\right),v\left(t\right),D\left(t\right),Z\left(t\right))}^T$ be the solution of model (2) with any initial function $\varphi \in X_2$. Clearly, $T\left(t\right)>0$, $I\left(t\right)>0$, $v\left(t\right)>0$, and $D\left(t\right)>0$ for $t\ge 0$.

Let $U_1\left(t\right)=W_{10}\left(t\right)+W_{11}\left(t\right)+W_{12}\left(t\right)$, where

$$\begin{array}{cc}\hfill W_{10}\left(t\right)=& {\beta }_1{e}^{-d_1\tau }{T}_{1}f\left(\frac{T\left(t\right)}{T_1}\right)+e^{-d_1\tau }{D}_{1}f\left(\frac{D\left(t\right)}{D_1}\right)\hfill \\ \hfill & +{\beta }_1{I}_1\left(\frac{I\left(t\right)}{I_1}\right)+\frac{{\beta }_1}{M}{v}_{1}f\left(\frac{v\left(t\right)}{v_1}\right)+{\beta }_1\frac{p}{q}Z\left(t\right),\hfill \\ \hfill W_{11}\left(t\right)=& {\beta }_1\beta D_1{v}_1{T}_1{e}^{-d_1\tau }{\int }_{t-\tau }^{t}f\left(\frac{D\left(s\right)v\left(s\right)T\left(s\right)}{D_1{v}_1{T}_1}\right)ds,\hfill \\ \hfill W_{12}\left(t\right)=& \frac{\beta v_{max}{e}^{-d_1\tau }}{2(d+\gamma )}{[{\beta }_1(T\left(t\right)-T_1)-(D\left(t\right)-D_1)]}^2.\hfill \end{array}$$

In calculating the derivative of $W_{10}\left(t\right)$ along the solution ${(T\left(t\right),I\left(t\right),v\left(t\right),D\left(t\right),Z\left(t\right))}^T$ of model (2), it follows that for $t\ge 0$,

$$\begin{array}{cc}\hfill {\dot{W}}_{10}\left(t\right)=& {\beta }_1{e}^{-d_1\tau }\left(1-\frac{T_1}{T\left(t\right)}\right)\left[\beta D_1{v}_1{T}_1-\beta D\left(t\right)v\left(t\right)T\left(t\right)-d(T\left(t\right)-T_1)\right]\hfill \\ \hfill & +e^{-d_1\tau }\left(1-\frac{D_1}{D\left(t\right)}\right)\left[{\beta }_1\beta D_1{v}_1{T}_1-{\beta }_1\beta D\left(t\right)v\left(t\right)T\left(t\right)-\gamma (D\left(t\right)-D_1)\right]\hfill \\ \hfill & +{\beta }_1\left(1-\frac{I_1}{I\left(t\right)}\right)\left[\beta e^{-d_1\tau }D(t-\tau )v(t-\tau )T(t-\tau )-d_{1}I\left(t\right)-pI\left(t\right)Z\left(t\right)\right]\hfill \\ \hfill & +\frac{{\beta }_1}{M}\left(1-\frac{v_1}{v\left(t\right)}\right)\left(d_{1}MI\left(t\right)-cv\left(t\right)\right)+{\beta }_{1}pI\left(t\right)Z\left(t\right)-{\beta }_1\frac{pb}{q}Z\left(t\right)\hfill \\ \hfill =& -\frac{d{\beta }_1{e}^{-d_1\tau }}{T\left(t\right)}{(T\left(t\right)-T_1)}^2-\frac{\gamma e^{-d_1\tau }}{D\left(t\right)}{(D\left(t\right)-D_1)}^2\hfill \\ \hfill & +2{\beta }_1\beta e^{-d_1\tau }(D_1{v}_1{T}_1-D\left(t\right)v\left(t\right)T\left(t\right))-{\beta }_1\beta e^{-d_1\tau }{D}_1{v}_1{T}_1\frac{T_1}{T\left(t\right)}\hfill \\ \hfill & +{\beta }_1\beta e^{-d_1\tau }D\left(t\right)v\left(t\right)T_1-\beta {\beta }_1{e}^{-d_1\tau }{D}_1{v}_1{T}_1\frac{D_1}{D\left(t\right)}+{\beta }_1\beta e^{-d_1\tau }{D}_{1}v\left(t\right)T\left(t\right)\hfill \\ \hfill & +{\beta }_1\beta e^{-d_1\tau }D(t-\tau )v(t-\tau )T(t-\tau )-{\beta }_1\beta e^{-d_1\tau }\frac{I_1}{I\left(t\right)}D(t-\tau )v(t-\tau )T(t-\tau )\hfill \\ \hfill & +{\beta }_1{d}_1{I}_1-{\beta }_1{d}_{1}I\left(t\right)\frac{v_1}{v\left(t\right)}+\frac{{\beta }_1}{M}c{v}_1-\frac{{\beta }_1}{M}cv\left(t\right)+{\beta }_{1}p\left(I_1-\frac{b}{q}\right)Z\left(t\right),\hfill \end{array}$$

(40)

where $\lambda =\beta D_1{v}_1{T}_1+d{T}_1$ and ${\lambda }_1={\beta }_1\beta D_1{v}_1{T}_1+\gamma D_1$ are used. Note that

$${\beta }_1\beta e^{-d_1\tau }{D}_1{v}_1{T}_1={\beta }_1{d}_1{I}_1=\frac{{\beta }_1}{M}c{v}_1,$$

thus, we have for $t\ge 0$,

$$\begin{array}{cc}\hfill {\dot{W}}_{10}\left(t\right)=& -\frac{d{\beta }_1{e}^{-d_1\tau }}{T\left(t\right)}{(T\left(t\right)-T_1)}^2-\frac{\gamma e^{-d_1\tau }}{D\left(t\right)}{(D\left(t\right)-D_1)}^2+{\beta }_{1}p\left(I_1-\frac{b}{q}\right)Z\left(t\right)\hfill \\ \hfill & +{\beta }_1\beta e^{-d_1\tau }{D}_1{v}_1{T}_1\left[4-\frac{T_1}{T\left(t\right)}-\frac{D_1}{D\left(t\right)}-\frac{D(t-\tau )v(t-\tau )T(t-\tau )I_1}{D_1{v}_1{T}_{1}I\left(t\right)}-\frac{I\left(t\right)v_1}{I_{1}v\left(t\right)}\right]\hfill \\ \hfill & -{\beta }_1\beta e^{-d_1\tau }[D\left(t\right)T\left(t\right)+D_1{T}_1-D_{1}T\left(t\right)-D\left(t\right)T_1]v\left(t\right)\hfill \\ \hfill & +{\beta }_1\beta e^{-d_1\tau }(D(t-\tau )v(t-\tau )T(t-\tau )-D\left(t\right)v\left(t\right)T\left(t\right)).\hfill \end{array}$$

(41)

In calculating the derivative of $W_{11}\left(t\right)$ and $W_{12}\left(t\right)$ along the solution $\left(T\right(t),I(t),v(t),$${D\left(t\right),Z\left(t\right))}^T$ of model (2), it follows that for $t\ge 0$,

$$\begin{array}{cc}\hfill {\dot{W}}_{11}\left(t\right)=& {\beta }_1\beta e^{-d_1\tau }[D\left(t\right)v\left(t\right)T\left(t\right)-D(t-\tau )v(t-\tau )T(t-\tau )]\hfill \\ \hfill & +{\beta }_1\beta e^{-d_1\tau }{D}_1{v}_1{T}_{1}ln\left(\frac{D(t-\tau )v(t-\tau )T(t-\tau )}{D\left(t\right)v\left(t\right)T\left(t\right)}\right)\hfill \\ \hfill =& {\beta }_1\beta e^{-d_1\tau }[D\left(t\right)v\left(t\right)T\left(t\right)-D(t-\tau )v(t-\tau )T(t-\tau )]\hfill \\ \hfill & +{\beta }_1\beta e^{-d_1\tau }{D}_1{v}_1{T}_1\{ln\left(\frac{T_1}{T\left(t\right)}\right)+ln\left(\frac{D_1}{D\left(t\right)}\right)\hfill \\ \hfill & +ln\left(\frac{D(t-\tau )v(t-\tau )T(t-\tau )I_1}{D_1{v}_1{T}_{1}I\left(t\right)}\right)+ln\left(\frac{I\left(t\right)v_1}{I_{1}v\left(t\right)}\right)\},\hfill \end{array}$$

(42)

$$\begin{array}{cc}\hfill {\dot{W}}_{12}\left(t\right)=& \frac{\beta v_{max}{e}^{-d_1\tau }}{d+\gamma }[{\beta }_1(T\left(t\right)-T_1)-(D\left(t\right)-D_1)]\{{\beta }_1(\beta D_1{v}_1{T}_1-\beta D\left(t\right)v\left(t\right)T\left(t\right)+d{T}_1-dT\left(t\right))\hfill \\ \hfill & -({\beta }_1\beta D_1{v}_1{T}_1-{\beta }_1\beta D\left(t\right)v\left(t\right)T\left(t\right)+\gamma D_1-\gamma D\left(t\right))\}\hfill \\ \hfill =& \frac{\beta v_{max}{e}^{-d_1\tau }}{d+\gamma }[{\beta }_1(T\left(t\right)-T_1)-(D\left(t\right)-D_1)][{\beta }_{1}d(T_1-T\left(t\right))-\gamma (D_1-D\left(t\right))]\hfill \\ \hfill =& -\frac{d{\beta }_1^2\beta v_{max}{e}^{-d_1\tau }}{d+\gamma }{(T\left(t\right)-T_1)}^2-\frac{\gamma \beta v_{max}{e}^{-d_1\tau }}{d+\gamma }{(D\left(t\right)-D_1)}^2\hfill \\ \hfill & +{\beta }_1\beta e^{-d_1\tau }{v}_{max}(T\left(t\right)-T_1)(D\left(t\right)-D_1).\hfill \end{array}$$

(43)

From (41), (42), and (43), we have for $t\ge 0$,

$$\begin{array}{cc}\hfill {\dot{U}}_1\left(t\right)=& -\left(\frac{d{\beta }_1{e}^{-d_1\tau }}{T\left(t\right)}+\frac{d{\beta }_1^2\beta v_{max}{e}^{-d_1\tau }}{d+\gamma }\right){(T\left(t\right)-T_1)}^2\hfill \\ \hfill & -\left(\frac{\gamma e^{-d_1\tau }}{D\left(t\right)}+\frac{\gamma \beta v_{max}{e}^{-d_1\tau }}{d+\gamma }\right){(D\left(t\right)-D_1)}^2\hfill \\ \hfill & +{\beta }_1\beta e^{-d_1\tau }(v_{max}-v\left(t\right))(T\left(t\right)-T_1)(D\left(t\right)-D_1)\hfill \\ \hfill & +{\beta }_{1}p\left(I_1-\frac{b}{q}\right)Z\left(t\right)-{\beta }_1\beta e^{-d_1\tau }{D}_1{v}_1{T}_1\mathrm{\Pi }\left(t\right),\hfill \end{array}$$

(44)

where

$$\mathrm{\Pi }\left(t\right)=f\left(\frac{T_1}{T\left(t\right)}\right)+f\left(\frac{D_1}{D\left(t\right)}\right)+f\left(\frac{D(t-\tau )v(t-\tau )T(t-\tau )I_1}{D_1{v}_1{T}_{1}I\left(t\right)}\right)+f\left(\frac{I\left(t\right)v_1}{I_{1}v\left(t\right)}\right)\ge 0.$$

If condition (H1) holds, then there exists a sufficiently small $\epsilon >0$ such that

$$4\left(\frac{d}{T_0+\epsilon }+\frac{d{\beta }_1\beta v_{max}}{d+\gamma }\right)\left(\frac{\gamma }{D_0+\epsilon }+\frac{\gamma \beta v_{max}}{d+\gamma }\right)>{\beta }_1{\beta }^2{v}_{max}^2,$$

which implies that the matrix

$$J:=\left(\begin{array}{cc}\frac{d{\beta }_1}{T_0+\epsilon }+\frac{d{\beta }_1^2\beta v_{max}}{d+\gamma }& -\frac{{\beta }_1\beta v_{max}}{2}\\ -\frac{{\beta }_1\beta v_{max}}{2}& \frac{\gamma }{D_0+\epsilon }+\frac{\gamma \beta v_{max}}{d+\gamma }\end{array}\right)$$

is positive definite. It follows from Theorem 1 that for the above $\epsilon >0$, there exists a $\widehat{T_1}\left(\epsilon \right)>\tau$ such that for $t>\widehat{T_1}\left(\epsilon \right)$,

$$T\left(t\right)<T_0+\epsilon ,D\left(t\right)<D_0+\epsilon ,|v_{max}-v\left(t\right)|<v_{max}.$$

(45)

From Lemma 1, (44), and (45), we have for $t>T_1\left(\epsilon \right)$,

$$\begin{array}{cc}\hfill {\dot{U}}_1\left(t\right)\le & -\left(\frac{d{\beta }_1{e}^{-d_1\tau }}{T_0+\epsilon }+\frac{d{\beta }_1^2\beta v_{max}{e}^{-d_1\tau }}{d+\gamma }\right){(T\left(t\right)-T_1)}^2\hfill \\ \hfill & -\left(\frac{\gamma e^{-d_1\tau }}{D_0+\epsilon }+\frac{\gamma \beta v_{max}{e}^{-d_1\tau }}{d+\gamma }\right){(D\left(t\right)-D_1)}^2\hfill \\ \hfill & +{\beta }_1\beta e^{-d_1\tau }{v}_{max}|(T\left(t\right)-T_1)\left|\right|(D\left(t\right)-D_1)|\hfill \\ \hfill & +{\beta }_{1}p\Theta (R_1-1)Z\left(t\right)-{\beta }_1\beta e^{-d_1\tau }{D}_1{v}_1{T}_1\mathrm{\Pi }\left(t\right)\hfill \\ \hfill =& -e^{-d_1\tau }\left(\right|T\left(t\right)-T_1|,|D\left(t\right)-D_1\left|\right)J\left(\begin{array}{c}|T\left(t\right)-T_1|\\ |D\left(t\right)-D_1|\end{array}\right)\hfill \\ \hfill & +{\beta }_{1}p\Theta (R_1-1)Z\left(t\right)-{\beta }_1\beta e^{-d_1\tau }{D}_1{v}_1{T}_1\mathrm{\Pi }\left(t\right).\hfill \end{array}$$

(46)

Note that as the matrix J is positive definite and $R_1<1$, we have for $t>\widehat{T_1}\left(\epsilon \right)$, $U_1\left(t\right)\le 0$, which implies that ${lim}_{t\to +\infty }{U}_1\left(t\right)$ exists. Moreover, according to Theorem 1 and Lemma 3, we have that for $t\ge \widehat{T_1}\left(\epsilon \right)$, $\ddot{A}\left(t\right)$ is bounded. This implies that $U_1\left(t\right)$ is uniformly continuous for $t>\widehat{T_1}\left(\epsilon \right)$. Then, it follows from Barbalat’s lemma [37] that

$$\underset{t\to +\infty }{lim}T\left(t\right)=T_1,\underset{t\to +\infty }{lim}D\left(t\right)=D_1,\underset{t\to +\infty }{lim}Z\left(t\right)=0,\underset{t\to +\infty }{lim}\frac{I\left(t\right)v_1}{I_{1}v\left(t\right)}=1.$$

Furthermore, from the first and third equations of model (2), we can obtain ${lim}_{t\to +\infty }v\left(t\right)=v_1$ and ${lim}_{t\to +\infty }I\left(t\right)=I_1$. Thus, the immunity-inactivated equilibrium $E_1$ is globally attractive. □

Remark 2.

Assume that $\tau =0$. Let $a\in (0,1)$, ${\mu }_3=min\{d,(1-a)d_1,c\}$, and ${\widehat{v}}_{max}=\frac{\lambda M}{a{\mu }_3}$. In [14], Tang et al. proved that the set

$$\mathrm{\Omega }:=\left\{(T,I,v,D)|0<T\le T_0,I\ge 0,v\ge 0,0<D\le D_0,T+I+\frac{a}{M}v\le \frac{\lambda }{{\mu }_3}\right\}$$

is attractive and positively invariant with respect to model (1), and the infected equilibrium $(T_1,I_1,v_1,D_1)$ of model (1) is globally asymptotically stable in Ω if $R_0>1$ and ${\beta }_1{\beta }^2{\lambda }_1\lambda {\left({\widehat{v}}_{max}\right)}^2\le 4d^2{\gamma }^2$. By using the Lyapuonv function

$$\begin{array}{cc}\hfill \widehat{U_1}\left(t\right)=& {\beta }_1{T}_{1}f\left(\frac{T\left(t\right)}{T_1}\right)+D_{1}f\left(\frac{D\left(t\right)}{D_1}\right)+{\beta }_1{I}_1\left(\frac{I\left(t\right)}{I_1}\right)+\frac{{\beta }_1}{M}{v}_{1}f\left(\frac{v\left(t\right)}{v_1}\right)\hfill \\ \hfill & +\frac{\beta {\widehat{v}}_{max}}{2(d+\gamma )}{[{\beta }_1(T\left(t\right)-T_1)-(D\left(t\right)-D_1)]}^2,\hfill \end{array}$$

the result for the global stability of the infected equilibrium $(T_1,I_1,v_1,D_1)$ of model (1) in [14] can be greatly improved (see Theorem 2 in [14]).

4.3. Global Stability of the Immunity-Activated Equilibrium $E^{*}$

Theorem 8.

If $R_1>1$ and condition (H1) holds, then the immunity-activated equilibrium $E^{*}$ of model (2) is globally asymptotically stable in $X_3:=\{\varphi \in X^{+}|{\varphi }_i\left(0\right)>0,i=1,2,3,4,5\}$.

Proof. 

By Theorem 5, we only need to prove that the immunity-activated equilibrium $E^{*}$ is globally attractive. Let ${(T\left(t\right),I\left(t\right),v\left(t\right),D\left(t\right),Z\left(t\right))}^T$ be the solution of model (2) with any initial function $\varphi \in X_3$. Clearly, $T\left(t\right)>0$, $I\left(t\right)>0$, $v\left(t\right)>0$, $D\left(t\right)>0$, and $Z\left(t\right)>0$ for $t\ge 0$.

Let $U_2\left(t\right)=W_{20}\left(t\right)+W_{21}\left(t\right)+W_{22}\left(t\right)$, where

$$\begin{array}{cc}\hfill W_{20}\left(t\right)=& {\beta }_1{e}^{-d_1\tau }{T}^{*}f\left(\frac{T\left(t\right)}{T^{*}}\right)+e^{-d_1\tau }{D}^{*}f\left(\frac{D\left(t\right)}{D^{*}}\right)+{\beta }_1{I}^{*}\left(\frac{I\left(t\right)}{I^{*}}\right)\hfill \\ \hfill & +\frac{{\beta }_1(d_1+p{Z}^{*})}{d_{1}M}{v}^{*}f\left(\frac{v\left(t\right)}{v^{*}}\right)+{\beta }_1\frac{p}{q}{Z}^{*}f\left(\frac{Z\left(t\right)}{Z^{*}}\right),\hfill \\ \hfill W_{21}\left(t\right)=& \beta {\beta }_1{D}^{*}{v}^{*}{T}^{*}{e}^{-d_1\tau }{\int }_{t-\tau }^{t}f\left(\frac{D\left(s\right)v\left(s\right)T\left(s\right)}{D^{*}{v}^{*}{T}^{*}}\right)ds,\hfill \\ \hfill W_{22}\left(t\right)=& \frac{\beta v_{max}{e}^{-d_1\tau }}{2(d+\gamma )}{[{\beta }_1(T\left(t\right)-T^{*})-(D\left(t\right)-D^{*})]}^2.\hfill \end{array}$$

Calculating the derivative of $W_{20}\left(t\right)$ along the solution ${(T\left(t\right),I\left(t\right),v\left(t\right),D\left(t\right),Z\left(t\right))}^T$ of model (2), it follows that for $t\ge 0$,

$$\begin{array}{cc}{\dot{W}}_{20}\left(t\right)=\hfill & {\beta }_1{e}^{-d_1\tau }\left(1-\frac{T^{*}}{T\left(t\right)}\right)\left[\lambda -\beta D\left(t\right)v\left(t\right)T\left(t\right)-dT\left(t\right)\right]\hfill \\ & +e^{-d_1\tau }\left(1-\frac{D^{*}}{D\left(t\right)}\right)\left[{\lambda }_1-\beta {\beta }_{1}D\left(t\right)v\left(t\right)T\left(t\right)-\gamma D\left(t\right)\right]\hfill \\ & +{\beta }_1\left(1-\frac{I^{*}}{I\left(t\right)}\right)\left[\beta e^{-d_1\tau }D(t-\tau )v(t-\tau )T(t-\tau )-(d_1+p{Z}^{*})I\left(t\right)+p(Z^{*}-Z\left(t\right))I\left(t\right)\right]\hfill \\ & +\frac{{\beta }_1(d_1+p{Z}^{*})}{d_{1}M}\left(1-\frac{v^{*}}{v\left(t\right)}\right)\left(d_{1}MI\left(t\right)-cv\left(t\right)\right)+{\beta }_1\frac{p}{q}\left(1-\frac{Z^{*}}{Z\left(t\right)}\right)(qI\left(t\right)-b)Z\left(t\right).\hfill \end{array}$$

(47)

From $\lambda =\beta D^{*}{v}^{*}{T}^{*}+d{T}^{*}$, ${\lambda }_1={\beta }_1\beta D^{*}{v}^{*}{T}^{*}+\gamma D^{*}$, $I^{*}=\frac{b}{q}$ and (47), we have for $t\ge 0$,

$$\begin{array}{cc}\hfill {\dot{W}}_{20}\left(t\right)=& {\beta }_1{e}^{-d_1\tau }\left(1-\frac{T^{*}}{T\left(t\right)}\right)\left[\beta D^{*}{v}^{*}{T}^{*}-\beta D\left(t\right)v\left(t\right)T\left(t\right)-d(T\left(t\right)-T^{*})\right]\hfill \\ \hfill & +e^{-d_1\tau }\left(1-\frac{D^{*}}{D\left(t\right)}\right)\left[\beta {\beta }_1{D}^{*}{v}^{*}{T}^{*}-\beta {\beta }_{1}D\left(t\right)v\left(t\right)T\left(t\right)-\gamma (D\left(t\right)-D^{*})\right]\hfill \\ \hfill & +{\beta }_1\left(1-\frac{I^{*}}{I\left(t\right)}\right)\left[\beta e^{-d_1\tau }D(t-\tau )v(t-\tau )T(t-\tau )-(d_1+p{Z}^{*})I\left(t\right)\right]\hfill \\ \hfill & +{\beta }_{1}p(I\left(t\right)-I^{*})(Z^{*}-Z\left(t\right))+\frac{{\beta }_1(d_1+p{Z}^{*})}{d_{1}M}\left(1-\frac{v^{*}}{v\left(t\right)}\right)\left(d_{1}MI\left(t\right)-cv\left(t\right)\right)\hfill \\ \hfill & +{\beta }_1\frac{p}{q}\left(Z\left(t\right)-Z^{*}\right)(qI\left(t\right)-q{I}^{*})\hfill \\ \hfill =& -\frac{d{\beta }_1{e}^{-d_1\tau }}{T\left(t\right)}{(T\left(t\right)-T^{*})}^2-\frac{\gamma e^{-d_1\tau }}{D\left(t\right)}{(D\left(t\right)-D^{*})}^2\hfill \\ \hfill & +2{\beta }_1\beta e^{-d_1\tau }(D^{*}{v}^{*}{T}^{*}-D\left(t\right)v\left(t\right)T\left(t\right))\hfill \\ \hfill & -{\beta }_1\beta e^{-d_1\tau }{D}^{*}{v}^{*}{T}^{*}\frac{T^{*}}{T\left(t\right)}+{\beta }_1\beta e^{-d_1\tau }D\left(t\right)v\left(t\right)T^{*}-\beta {\beta }_1{e}^{-d_1\tau }{D}^{*}{v}^{*}{T}^{*}\frac{D^{*}}{D\left(t\right)}\hfill \\ \hfill & +{\beta }_1\beta e^{-d_1\tau }{D}^{*}v\left(t\right)T\left(t\right)+{\beta }_1\beta e^{-d_1\tau }D(t-\tau )v(t-\tau )T(t-\tau )\hfill \\ \hfill & -{\beta }_1\beta e^{-d_1\tau }\frac{I^{*}}{I\left(t\right)}D(t-\tau )v(t-\tau )T(t-\tau )+{\beta }_1(d_1+p{Z}^{*})I^{*}\hfill \\ \hfill & -{\beta }_1(d_1+p{Z}^{*})I\left(t\right)\frac{v^{*}}{v\left(t\right)}+\frac{{\beta }_1(d_1+p{Z}^{*})}{d_{1}M}c{v}^{*}-\frac{{\beta }_1(d_1+p{Z}^{*})}{d_{1}M}cv\left(t\right).\hfill \end{array}$$

(48)

Note that

$${\beta }_1\beta e^{-d_1\tau }{D}^{*}{v}^{*}{T}^{*}={\beta }_1(d_1+p{Z}^{*})I^{*}=\frac{{\beta }_1(d_1+p{Z}^{*})}{d_{1}M}c{v}^{*},$$

then we have for $t\ge 0$,

$$\begin{array}{cc}\hfill {\dot{W}}_{20}\left(t\right)=& -\frac{d{\beta }_1{e}^{-d_1\tau }}{T\left(t\right)}{(T\left(t\right)-T^{*})}^2-\frac{\gamma e^{-d_1\tau }}{D\left(t\right)}{(D\left(t\right)-D^{*})}^2\hfill \\ \hfill & +{\beta }_1\beta e^{-d_1\tau }{D}^{*}{v}^{*}{T}^{*}\left[4-\frac{T^{*}}{T\left(t\right)}-\frac{D^{*}}{D\left(t\right)}-\frac{D(t-\tau )v(t-\tau )T(t-\tau )I^{*}}{D^{*}{v}^{*}{T}^{*}I\left(t\right)}-\frac{I\left(t\right)v^{*}}{I^{*}v\left(t\right)}\right]\hfill \\ \hfill & -{\beta }_1\beta e^{-d_1\tau }[D\left(t\right)T\left(t\right)+D^{*}{T}^{*}-D^{*}T\left(t\right)-D\left(t\right)T^{*}]v\left(t\right)\hfill \\ \hfill & +{\beta }_1\beta e^{-d_1\tau }(D(t-\tau )v(t-\tau )T(t-\tau )-D\left(t\right)v\left(t\right)T\left(t\right)).\hfill \end{array}$$

(49)

In calculating the derivative of $W_{21}\left(t\right)$ and $W_{22}\left(t\right)$ along the solution ${(T\left(t\right),I\left(t\right),v\left(t\right),D\left(t\right),Z\left(t\right))}^T$ of model (2), it follows that for $t\ge 0$,

$$\begin{array}{cc}\hfill {\dot{W}}_{11}\left(t\right)=& {\beta }_1\beta e^{-d_1\tau }[D\left(t\right)v\left(t\right)T\left(t\right)-D(t-\tau )v(t-\tau )T(t-\tau )]\hfill \\ \hfill & +{\beta }_1\beta e^{-d_1\tau }{D}^{*}{v}^{*}{T}^{*}ln\left(\frac{D(t-\tau )v(t-\tau )T(t-\tau )}{D\left(t\right)v\left(t\right)T\left(t\right)}\right)\hfill \\ \hfill =& \beta {\beta }_1{e}^{-d_1\tau }[D\left(t\right)v\left(t\right)T\left(t\right)-D(t-\tau )v(t-\tau )T(t-\tau )]\hfill \\ \hfill & +\beta {\beta }_1{e}^{-d_1\tau }{D}^{*}{v}^{*}{T}^{*}\{ln\left(\frac{T^{*}}{T\left(t\right)}\right)+ln\left(\frac{D^{*}}{D\left(t\right)}\right)\hfill \\ \hfill & +ln\left(\frac{D(t-\tau )v(t-\tau )T(t-\tau )I^{*}}{D^{*}{v}^{*}{T}^{*}I\left(t\right)}\right)+ln\left(\frac{I\left(t\right)v^{*}}{I^{*}v\left(t\right)}\right)\},\hfill \end{array}$$

(50)

$$\begin{array}{cc}\hfill {\dot{W}}_{12}\left(t\right)=& \frac{\beta v_{max}{e}^{-d_1\tau }}{d+\gamma }[{\beta }_1(T\left(t\right)-T^{*})-(D\left(t\right)-D^{*})]\{{\beta }_1(\beta D^{*}{v}^{*}{T}^{*}-\beta D\left(t\right)v\left(t\right)T\left(t\right)+d{T}^{*}-dT\left(t\right))\hfill \\ \hfill & -({\beta }_1\beta D^{*}{v}^{*}{T}^{*}-{\beta }_1\beta D\left(t\right)v\left(t\right)T\left(t\right)+\gamma D^{*}-\gamma D\left(t\right))\}\hfill \\ \hfill =& \frac{\beta v_{max}{e}^{-d_1\tau }}{d+\gamma }[{\beta }_1(T\left(t\right)-T^{*})-(D\left(t\right)-D^{*})][{\beta }_{1}d(T^{*}-T\left(t\right))-\gamma (D^{*}-D\left(t\right))]\hfill \\ \hfill =& -\frac{d{\beta }_1^2\beta v_{max}{e}^{-d_1\tau }}{d+\gamma }{(T\left(t\right)-T^{*})}^2-\frac{\gamma \beta v_{max}{e}^{-d_1\tau }}{d+\gamma }{(D\left(t\right)-D^{*})}^2\hfill \\ \hfill & +{\beta }_1\beta e^{-d_1\tau }{v}_{max}(T\left(t\right)-T^{*})(D\left(t\right)-D^{*}).\hfill \end{array}$$

(51)

From (49), (50), and (51), we have for $t\ge 0$,

$$\begin{array}{cc}\hfill {\dot{U}}_2\left(t\right)=& -\left(\frac{d{\beta }_1{e}^{-d_1\tau }}{T\left(t\right)}+\frac{d{\beta }_1^2\beta v_{max}{e}^{-d_1\tau }}{d+\gamma }\right){(T\left(t\right)-T^{*})}^2\hfill \\ \hfill & -\left(\frac{\gamma e^{-d_1\tau }}{D\left(t\right)}+\frac{\gamma \beta v_{max}{e}^{-d_1\tau }}{d+\gamma }\right){(D\left(t\right)-D^{*})}^2\hfill \\ \hfill & +{\beta }_1\beta e^{-d_1\tau }(v_{max}-v\left(t\right))(T\left(t\right)-T^{*})(D\left(t\right)-D^{*})-{\beta }_1\beta e^{-d_1\tau }{D}^{*}{v}^{*}{T}^{*}\widehat{\mathrm{\Pi }}\left(t\right),\hfill \end{array}$$

(52)

where

$$\widehat{\mathrm{\Pi }}\left(t\right)=f\left(\frac{T^{*}}{T\left(t\right)}\right)+f\left(\frac{D^{*}}{D\left(t\right)}\right)+f\left(\frac{D(t-\tau )v(t-\tau )T(t-\tau )I^{*}}{D^{*}{v}^{*}{T}^{*}I\left(t\right)}\right)+f\left(\frac{I\left(t\right)v^{*}}{I^{*}v\left(t\right)}\right)\ge 0.$$

We claim that the immunity-inactivated equilibrium $E_1$ is globally attractive if $R_1>1$ and condition (H1) holds. The rest of the proof is very similar to the proof of Theorem 7. We omit the details here to avoid repetition. □

5. Numerical Simulations and Conclusions

We present some numerical simulations to illustrate our main theoretical results. Here, we fix $\lambda =4$, $\beta =0.0001$, $d=0.015$, $d_1=0.1$, $M=100$, $p=0.0025$, $c=6.5$, ${\lambda }_1=1$, ${\beta }_1=0.01$, $b=0.5$, and $\tau =1$ and change the values of $\gamma$ and q.

If we choose $q=0.01$ and $\gamma =0.42$, then we have $R_0\approx 0.8838460738<1$, and model (2) has an infection-free equilibrium $E_0\approx (266.6666667,0,0,2.380952381,0)$. According to Theorem 6, the infection-free equilibrium $E_0\approx (266.6666667,0,0,2.380952381,0)$ is globally asymptotically stable (see Figure 2), and the viral load eventually converges to 0.

If we choose $q=0.01$ and $\gamma =0.16$, then we have $R_0\approx 2.320095944>1$ and $R_1\approx 0.5444124848<1$, and the model (2) has two equilibria: the infection-free equilibrium $E_0\approx (266.6666667,0,0,6.25,0)$ and the immunity-inactivated equilibrium $E_1\approx (117.5671563$, $20.2366224,31.13326522,6.110219209,0)$. The calculation shows that condition (H1) is satisfied. Thus, it follows from Theorem 3 and Theorem 7 that the infection-free equilibrium $E_0\approx (266.6666667,0,0,6.25,0)$ is unstable, the immunity-inactivated equilibrium $E_1\approx (117.5671563,20.2366224,31.13326522,6.110219209,0)$ is globally asymptotically stable (see Figure 3), and the viral load eventually converges to $v_1\approx 31.13326522$.

If we choose $q=0.2$ and $\gamma =0.16$, then we have $R_1\approx 1.994781846>1$, and model (2) has three equilibria: the infection-free equilibrium $E_0\approx (266.6666667,0,0,6.25,0)$, the immunity-inactivated equilibrium $E_1\approx (117.5671563,20.2366224,31.13326522,6.110219209$, $0)$ and the immunity-activated equilibrium $E^{*}\approx (230.0089421,2.5,3.846153846,6.215633383$, $39.60627404)$. The calculation shows that condition (H1) is satisfied. Thus, it follows from Theorems 3, 4, and 8 that the infection-free equilibrium $E_0\approx (266.6666667,0,0,6.25,0)$ and the immunity-inactivated equilibrium $E_1\approx (117.5671563,20.2366224,31.13326522$, $6.110219209,0)$ are unstable, and the immunity-activated equilibrium $E^{*}\approx (230.0089421$, $2.5,3.846153846,6.215633383,39.60627404)$ is globally asymptotically stable (see Figure 4), and the viral load eventually converges to $v^{*}\approx 3.846153846<v_1\approx 31.13326522$.

In this paper, we propose a MERS-CoV infection model with CTL immune response and intracellular delay based on model (1). By analyzing the characteristic equations of the infection-free equilibrium $E_0$, the immunity-inactivated equilibrium $E_1$, and the immunity-activated equilibrium $E^{*}$ of model (2), we establish the complete results of local stability for three types of equilibria (see Theorems 3–5). The results also show that the intracellular delay $\tau$ does not change the local stability of the equilibria of model (2) (i.e., intracellular delay $\tau$ is harmless) and do not cause Hopf bifurcation.

Moreover, we investigated the global properties of model (2). The main results of this paper improve and extend the main results in [14]. Our main result shows that the infection-free equilibrium $E_0$ is globally asymptotically stable if the immunity-inactivated reproduction number

$$R_0=\frac{e^{-d_1\tau }\beta \lambda {\lambda }_{1}M}{cd\gamma }\le 1.$$

The infection-free equilibrium $E_0$ is globally asymptotically stable, meaning that the viruses in the host are eventually cleared. Note that $R_0$ is monotonically decreasing with respect to the delay $\tau$. According to the expression for $R_0$ and the global stability results for the infection-free equilibrium $E_0$, increasing the intracellular delay $\tau$, reducing the expression rate ${\lambda }_1$ of DPP4, and increasing the hydrolysis rate $\gamma$ of DPP4 are beneficial for controlling MERS-CoV infection. By constructing suitable Lyapunov functionals and using Barbalat’s lemma, we ascertain that if condition (H1) holds, the immunity-inactivated equilibrium $E_1$ is globally asymptotically stable when $R_0>1>R_1$, and the immunity-activated equilibrium $E^{*}$ is globally asymptotically stable when $R_1>1$. The immunity-inactivated equilibrium $E_1$ and immunity-activated equilibrium $E^{*}$ are globally asymptotically stable indicating that the viral load in the host eventually converges to the positive constant values $v_1$ and $v^{*}$, respectively. Condition (H1) in Theorems 7 and 8 is a technical assumption that may be weakened or even eliminated if more suitable Lyapunov functionals can be constructed.

Note that the immune-activated reproduction number

$$R_1=\frac{2q{e}^{-d_1\tau }\beta \lambda {\lambda }_{1}M}{2qcd\gamma +b{d}_1[M\beta ({\lambda }_1+\lambda {\beta }_1)+\sqrt{M^2{\beta }^2{({\lambda }_1-\lambda {\beta }_1)}^2+4M{\beta }_1\beta cd\gamma e^{d_1\tau }}]}$$

is positively correlated to the parameter q and

$$v^{*}=\frac{d_{1}Mb}{cq}\left(v^{*}<v_1=\frac{2Md\gamma (R_0-1)}{M\beta ({\lambda }_1+\lambda {\beta }_1)+\sqrt{M^2{\beta }^2{({\lambda }_1-\lambda {\beta }_1)}^2+4M{\beta }_1\beta cd\gamma e^{d_1\tau }}}\right)$$

is negatively correlated to the parameter q. If the immunity-inactivated reproduction number $R_0>1$, then it is advantageous to reduce the viral load by increasing the parameter q such that the immune-activated reproduction number $R_1>1$. Therefore, both intracellular delay $\tau$ and CTL immune response play critical roles in controlling MERS-CoV infection. Our Lemma 3 shows that if the immunity-inactivated reproduction number $R_0>1$, the viruses are persistent in the host. Moreoever, Lemma 3 gives a specific estimate of the ultimate lower bound on viral load. The result also suggests that the CTL immune response, while reducing the viral load in the host, does not ultimately clear the virus.