Uniform Persistence and Global Attractivity in a Delayed Virus Dynamic Model with Apoptosis and Both Virus-to-Cell and Cell-to-Cell Infections

Abstract

In this paper, we study the global dynamics of a delayed virus dynamics model with apoptosis and both virus-to-cell and cell-to-cell infections. When the basic reproduction number $R_0>1$, we obtain the uniform persistence of the model, and give some explicit expressions of the ultimate upper and lower bounds of any positive solution of the model. In addition, by constructing the appropriate Lyapunov functionals, we obtain some sufficient conditions for the global attractivity of the disease-free equilibrium and the chronic infection equilibrium of the model. Our results extend existing related works.

1. Introduction

It is well known that human health and safety have been seriously threatened by known or emerging new viral infections, such as human immunodeficiency virus (HIV), severe acute respiratory syndrome coronavirus 2 (SARS-CoV-2), etc. The mechanisms of transmission of viral infections have become increasingly complex, due to the mutation and evolution of viruses caused by changes in the physical environment, the use of drugs, and other factors. Since the 1980s, differential equations have been widely used in the study of important issues, such as the transmission mechanisms and control strategies of virus infection, and have gradually developed and established the important interdisciplinary research branch of virus dynamics [1,2,3,4,5]. In particular, in [1,2], the authors proposed the following classical viral infection dynamics model, describing the interactions among uninfected cells, infected cells and free viruses:

$$\left\{\begin{array}{cc}\hfill \dot{x}(t)=& s-dx(t)-\beta x(t)v(t),\hfill \\ \hfill \dot{y}(t)=& \beta x(t)v(t)-py(t),\hfill \\ \hfill \dot{v}(t)=& ky(t)-uv(t),\hfill \end{array}\right.$$

(1)

where $x(t)$, $y(t)$ and $v(t)$ denote the concentrations of uninfected cells, infected cells and free viruses at time t, respectively. The constant $s>0$ is the rate at which new uninfected cells are generated. The constant $d>0$ is the death rate of uninfected cells. The constant $\beta \ge 0$ is the characterizing infection of the cells. The constant $p>0$ is the death rate of infected cells. The infected cells produce virus particles at the constant rate $k\ge 0$, and the constant $u>0$ is the rate at which the virus is cleared. The term $\beta x(t)v(t)$ denotes the rate at which uninfected cells become infected cells through their contact with free viruses.

Based on model (1), many scholars have extended the linear growth rate $s-dx(t)$ of uninfected cells to classical logistic growth or more general nonlinear functions, and extended the bilinear functional response function $\beta x(t)v(t)$ to the following classical forms, $\beta x(t)v(t)/(1+v(t))$, $\beta x(t)v(t)/(1+ax(t)+bv(t))$, $\beta x(t)v(t)/(1+av(t))$, and more general nonlinear functions (see, for example, [6,7,8,9,10,11,12,13,14,15,16,17] and the references therein). Additionally, important factors, such as delay and immune response, were considered, and a series of important results on global stability and existence of Hopf bifurcations were obtained [16,18,19,20,21,22,23].

In addition, recent studies have also shown that a large number of viral particles can also be transferred from infected cells to uninfected cells through the formation of virally induced structures termed virological synapses [24,25]. The direct fusion between infected cells and uninfected cells can also lead to cell infection, which is also known as cell-to-cell infection [25,26,27]. Based on this important fact, many scholars have proposed several important classes of viral dynamics models, considering the more general nonlinear growth rate of uninfected cells (which can include linear and logistic growth), while introducing important factors such as virus-to-cell infection, cell-to-cell infection, and immune response and delay, and have thoroughly investigated the local and global dynamics of equilibria and the existence of Hopf bifurcations. For details, see, for example, [18,19,20,21,25,26,27,28,29,30,31] and the references cited therein. In particular, based on the studies in [7,25,26], the authors [28] proposed and studied the following virus infection dynamic model:

$$\left\{\begin{array}{cc}\hfill \dot{x}(t)=& rx(t)\left(1-\frac{x(t)+\alpha y(t)}{K}\right)-{\beta }_{1}x(t)y(t)-{\beta }_{2}x(t)v(t),\hfill \\ \hfill \dot{y}(t)=& {\beta }_{1}x(t)y(t)+{\beta }_{2}x(t)v(t)-py(t),\hfill \\ \hfill \dot{v}(t)=& ky(t)-uv(t).\hfill \end{array}\right.$$

(2)

In model (2), the constant $K>0$ denotes the effective carrying capacity of the environment of uninfected cells and infected cells. The term $rx(1-(x+\alpha y)/K)$ indicates that the growth of uninfected cells conforms to the logistic growth function and takes into account the effect of infected cells on the maximum carrying capacity of the environment, the constant $r>0$ is the growth rate, and $\alpha \ge 0$ is a constant; the constant ${\beta }_1\ge 0$ is the cell-to-cell infection rate, and the constant ${\beta }_2\ge 0$ is the virus-to-cell infection rate. All other parameters in model (2) have the same biological meaning as that in model (1). In [28], for model (2), the authors obtained the local stability of the equilibria, uniform persistence and the existence of Hopf bifurcations caused by the cell-to-cell infection rate ${\beta }_1$ or the virus-to-cell infection rate ${\beta }_2$.

HIV gene expression products can produce toxicity, which directly or indirectly induces apoptosis in uninfected cells [32]. Studies have shown that viral proteins interact with uninfected cells and induce an apoptotic signal, which induces the death of uninfected cells [33]. In [34], the authors considered the following virus infection dynamic model with delay:

$$\left\{\begin{array}{cc}\hfill \dot{x}(t)=& s-dx(t)-cx(t)y(t)-\beta x(t)v(t),\hfill \\ \hfill \dot{y}(t)=& \delta x(t-\tau )v(t-\tau )-py(t),\hfill \\ \hfill \dot{v}(t)=& ky(t)-uv(t),\hfill \end{array}\right.$$

(3)

where the constant $\delta =\beta e^{-m_0\tau }>0$ denotes the surviving rate of infected cells before it becomes productively infected, $m_0\ge 0$ is a constant, and $\tau \ge 0$ is a delay. The constant $c\ge 0$ is the rate of apoptosis at which infected cells induce uninfected cells [32,33]. All other parameters in model (3) have the same biological meaning as that in model (1). Based on models (2) and (3), in [22,23,35], the authors further considered virus dynamics models with the logistic growth of uninfected cells, nonlinear infection rate, cell-to-cell infection, virus-to-cell infection and delay, and investigated the permanence, the global stability of the disease-free equilibrium, the local stability of the chronic infection equilibrium and the existence of Hopf bifurcations.

In this paper, based on [18,19,20,21,22,23,25,28,30,34,35], etc., we continue to consider the following delayed virus dynamic model with apoptosis and both virus-to-cell and cell-to-cell infections:

$$\left\{\begin{array}{cc}\hfill \dot{x}(t)=& s-dx(t)+rx(t)\left(1-\frac{x(t)+\alpha y(t)}{K}\right)-cx(t)y(t)-{\beta }_{1}x(t)y(t)-{\beta }_{2}x(t)v(t),\hfill \\ \hfill \dot{y}(t)=& {\beta }_1{e}^{-m_1{\tau }_1}x(t-{\tau }_1)y(t-{\tau }_1)+{\beta }_2{e}^{-m_2{\tau }_2}x(t-{\tau }_2)v(t-{\tau }_2)-py(t),\hfill \\ \hfill \dot{v}(t)=& ky(t)-uv(t).\hfill \end{array}\right.$$

(4)

In model (4), the delay ${\tau }_1\ge 0$ represents the time between infected cells spreading viruses into uninfected cells and the production of new free viruses; the delay ${\tau }_2\ge 0$ represents the time between viral entry into an uninfected cell and the production of new free viruses. $m_1\ge 0$ and $m_2\ge 0$ are constants, and ${\delta }_1:={\beta }_1{e}^{-m_1{\tau }_1}$ and ${\delta }_2:={\beta }_2{e}^{-m_2{\tau }_2}$ denote the survival rates of uninfected cells during successful infection with infected cells and free viruses, respectively. All other parameters have the same biological meaning as that in models (1)–(3).

Let $\tau =max\{{\tau }_1,{\tau }_2\}$. The initial condition of model (4) is given as follows, $x(\theta )={\varphi }_1(\theta )$, $y(\theta )={\varphi }_2(\theta )$, $v(\theta )={\varphi }_3(\theta )$ $(\theta \in [-\tau ,0])$, where $\varphi =({\varphi }_1,{\varphi }_2,{\varphi }_3)\in C^{+}$ $:=\{\varphi \in C|{\varphi }_i\ge 0,i=1,2,3\}$, $C=C([-\tau ,0],{\mathbb{R}}^3)$ is the Banach space of continuous functions from $[-\tau ,0]$ to ${\mathbb{R}}^3$ equipped with the supremum norm.

By using the standard theory of functional differential equations (see [36,37,38,39]), it is easy to show that the solution $(x(t),y(t),v(t))$ of model (4) with the above initial condition is existent, unique, non-negative on $[0,+\infty )$, and satisfies

$$\underset{t\to \infty }{lim\; sup}x(t)\le x_0,\underset{t\to \infty }{lim\; sup}y(t)\le M_2,\underset{t\to \infty }{lim\; sup}v(t)\le M_3,$$

(5)

where

$$x_0=\frac{K}{2r}\left(r-d+\sqrt{{(d-r)}^2+\frac{4rs}{K}}\right),M_2=\frac{2n_0(e^{-m_1{\tau }_1}+e^{-m_2{\tau }_2})}{\overline{p}},M_3=\frac{k{M}_2}{u},$$

$$n_0=\underset{x\in [0,x_0]}{max}\left(s-dx+rx\left(1-\frac{x}{K}\right)\right),\overline{p}=min\left\{p,\frac{n_0}{x_0}\right\}.$$

Obviously, model (4) always has a disease-free equilibrium (boundary equilibrium) $E_0=(x_0,0,0)$. We can easily derive the expression of the basic reproduction number of model (4) as

$$R_0=\frac{x_0(u{\beta }_1{e}^{-m_1{\tau }_1}+k{\beta }_2{e}^{-m_2{\tau }_2})}{pu}=\frac{x_0(u{\delta }_1+k{\delta }_2)}{pu}$$

by the method of the next generation matrix [40,41]. The basic reproduction number $R_0$ is positively correlated with respect to the cell-to-cell infection rate ${\beta }_1$ and the virus-to-cell infection rate ${\beta }_2$. Hence, when only one route of infection is considered, the evolution of the disease infection may be underestimated.

The function is defined as

$$\mathsf{\Gamma }(z)=s-dz+rz\left(1-\frac{z}{K}\right)=\frac{r}{K}(z+x_1)(x_0-z)(z\ge 0),x_1=-\frac{K}{2r}\left(r-d-\sqrt{{(d-r)}^2+\frac{4rs}{K}}\right)>0.$$

Note that if $R_0=\frac{x_0}{x^{*}}>1$, then model (4) has a unique chronic infection equilibrium (positive equilibrium) $E^{*}=(x^{*},y^{*},v^{*})$, where

$$x^{*}=\frac{pu}{u{\delta }_1+k{\delta }_2},y^{*}=\frac{\mathsf{\Gamma }(x^{*})}{x^{*}\xi }=\frac{r(x^{*}+x_1)(x_0-x^{*})}{K{x}^{*}\xi }>0,v^{*}=\frac{k{y}^{*}}{u},\xi =\frac{r\alpha }{K}+c+{\beta }_1+\frac{k{\beta }_2}{u}.$$

(6)

It is noted that the apoptosis rate c has effects in reducing the loads of both infected cells and free viruses. In addition, it is easy to show that the set $G:=\{\varphi =({\varphi }_1,{\varphi }_2,{\varphi }_3)\in C^{+}|0\le {\varphi }_1\le x_0\}$ is attractive and positively invariant with respect to model (4).

For the global asymptotic stability (global attractivity) of the disease-free equilibrium $E_0$ of model (4), using the method similar to that in [18,23,34,35], the following conclusion can be obtained (the proof is omitted): if $R_0<1$ ($R_0=1$), then the disease-free equilibrium $E_0$ is globally asymptotically stable (globally attractive) in G.

As far as we know, the global attractivity of the chronic infection equilibrium $E^{*}$ of model (4) is still a difficult mathematical question and worthy of further study. This paper has the following two main purposes. First of all, we study the uniform persistence of model (4) in Section 2, and give explicit expressions of the ultimate upper and lower bounds of any positive solution of model (4). Second, by constructing some appropriate Lyapunov functionals and combining inequality analysis, some sufficient conditions for the global attractivity of the chronic infection equilibrium $E^{*}$ of model (4) are given in Section 3. The brief summary of the conclusions of this paper is given in Section 4.

2. Uniform Persistence

In this section, we assume that $R_0>1$. It is not difficult to find that the function

$$f_1(z):=\frac{K}{2r}\left(z+\sqrt{z^2+\frac{4rs}{K}}\right)$$

is strictly monotonically increasing with respect to z on $\mathbb{R}$. According to the first equation of model (4), $x^{*}$ can be rewritten as $x^{*}=f_1(l_0)$, where $l_0=r-d-\left(\frac{r\alpha }{K}+c+{\beta }_1\right)y^{*}-{\beta }_2{v}^{*}$. For convenience, let us define the following parameters:

$$\begin{array}{c}{\nu }_1=f_1(l_1),x_1^{*}=f_1(l_2),\widehat{x_1^{*}}=\frac{K}{2r}\left(l_2-\sqrt{l_2^2+\frac{4rs}{K}}\right),\hfill \\ l_1=r-d-\left(\frac{r\alpha }{K}+c+{\beta }_1\right)M_2-{\beta }_2{M}_3,l_2=r-d-\frac{1}{2}\left(\frac{r\alpha }{K}+c+{\beta }_1\right)y^{*}-{\beta }_2{v}^{*}.\hfill \end{array}$$

Note that $x_0=f_1(r-d)$ and $r-d>l_2>l_0>l_1$, we can obtain $x_0>x_1^{*}>x^{*}>{\nu }_1>0$. For the uniform persistence of model (4), we have the following main results.

Theorem 1.

If $R_0>1$, then model (4) is uniformly persistent in $X^{+}:=\{\varphi =({\varphi }_1,{\varphi }_2,{\varphi }_3)\in C^{+}|{\varphi }_2(0)>0,{\varphi }_3(0)>0\}$, and the solution $(x(t),y(t),v(t))$ of model (4) with any $\varphi \in X^{+}$ satisfies

$$\begin{array}{cc}\hfill & \underset{t\to \infty }{lim\; inf}x(t)\ge {\nu }_1,\underset{t\to \infty }{lim\; inf}y(t)\ge \frac{y^{*}}{2}{e}^{-p\varpi }\equiv {\nu }_2,\underset{t\to \infty }{lim\; inf}v(t)\ge \frac{k{y}^{*}}{2u}{e}^{-p\varpi }=\frac{k}{u}{\nu }_2\equiv {\nu }_3,\hfill \end{array}$$

(7)

where $\varpi =T_0+T_1+T_2+{\tau }_1+{\tau }_2$,

$$\begin{array}{c}{T}_0=-\frac{1}{u}ln\left(\frac{y^{*}}{2M_2}\right)>0,T_1=\frac{-K}{r(x_1^{*}-\widehat{x_1^{*}})}ln\left[\left(\frac{x_1^{*}-x^0}{x^0-\widehat{x_1^{*}}}\right)\left(\frac{\gamma {\nu }_1-\widehat{x_1^{*}}}{x_1^{*}-\gamma {\nu }_1}\right)\right]>0,\hfill \\ T_2=\frac{q}{u(1-q)}>0,\gamma \in (0,1),x^0\in (x^{*},x_1^{*}),q=\frac{x^{*}}{x^0}<1.\hfill \end{array}$$

Proof. 

Let $(x(t),y(t),v(t))$ be any solution of model (4) with any $\varphi \in X^{+}$. By (5), for any $\epsilon >0$, there exists a sufficiently large $\widehat{t}>\tau$ such that, for $t>\widehat{t}$, $y(t)\le M_2+\epsilon$ and $v(t)\le M_3+\epsilon$. From the first equation of model (4), we have, for $t>\widehat{t}$,

$$\begin{array}{cc}\hfill \dot{x}(t)\ge & s-dx(t)+rx(t)\left(1-\frac{x(t)+\alpha (M_2+\epsilon )}{K}\right)-(c+{\beta }_1)\left(M_2+\epsilon \right)x(t)-{\beta }_2\left(M_3+\epsilon \right)x(t)\hfill \\ \hfill =& -\frac{r}{K}(x(t)-{\nu }_1(\epsilon ))(x(t)-\widehat{{\nu }_1}(\epsilon )),\hfill \end{array}$$

where

$$\begin{array}{c}{\nu }_1(\epsilon )=\frac{K}{2r}\left(l_1(\epsilon )+\sqrt{l_1^2(\epsilon )+\frac{4rs}{K}}\right)>0,\widehat{{\nu }_1}(\epsilon )=\frac{K}{2r}\left(l_1(\epsilon )-\sqrt{l_1^2(\epsilon )+\frac{4rs}{K}}\right)<0,\hfill \\ l_1(\epsilon )=r-d-\left(\frac{r\alpha }{K}+c+{\beta }_1\right)(M_2+\epsilon )-{\beta }_2(M_3+\epsilon ).\hfill \end{array}$$

Using the arbitrariness of $\epsilon$, we have ${lim\; inf}_{t\to \infty }x(t)\ge {\nu }_1(0)={\nu }_1$.

Next, using a method similar to that in [34,35], let us prove that $lim{inf}_{t\to \infty }y(t)\ge {\nu }_2$.

Note that $\gamma \in (0,1)$, then there exists a sufficiently large $T>0$ such that for $t\ge T$,

$$x(t)>\gamma {\nu }_1,v(t)\le \frac{k{M}_2}{u}+\frac{k{y}^{*}}{2u}.$$

Let us first claim that, for any $t_0\ge T$, when $t\ge t_0$, the inequality $y(t)\le \frac{y^{*}}{2}$ cannot always hold.

If this claim is not true, then there exists a $t_0\ge T$ such that $y(t)\le \frac{y^{*}}{2}$ for all $t\ge t_0$. Then, from the third equation of model (4), we have for $t\ge t_0$, $\dot{v}(t)\le \frac{1}{2}k{y}^{*}-uv(t),$ which implies that, for $t\ge t_0$,

$$v(t)\le \frac{k{y}^{*}}{2u}+\left(v(t_0)-\frac{k{y}^{*}}{2u}\right)e^{-u(t-t_0)}\le \frac{k{y}^{*}}{2u}+\frac{k{M}_2}{u}{e}^{-u(t-t_0)}.$$

Hence, we have, for $t\ge t_0+T_0$, $v(t)\le \frac{k}{u}{y}^{*}=v^{*}$. Then, from the first equation of model (4), we have, for $t\ge t_0+T_0$,

$$\begin{array}{cc}\hfill \dot{x}(t)\ge & s-dx(t)+rx(t)\left(1-\frac{x(t)}{K}\right)-\frac{r\alpha y^{*}}{2K}x(t)-\frac{c{y}^{*}}{2}x(t)-\frac{{\beta }_1{y}^{*}}{2}x(t)-{\beta }_2{v}^{*}x(t)\hfill \\ \hfill =& -\frac{r}{K}(x(t)-x_1^{*})(x(t)-\widehat{x_1^{*}}),\hfill \end{array}$$

which implies that, for $t\ge t_0+T_0$,

$$\begin{array}{c}\hfill x(t)\ge \frac{x_1^{*}-\widehat{x_1^{*}}\left(\frac{x(t_0+T_0)-x_1^{*}}{x(t_0+T_0)-\widehat{x_1^{*}}}\right)e^{-\frac{r}{K}(x_1^{*}-\widehat{x_1^{*}})(t-t_0-T_0)}}{1-\left(\frac{x(t_0+T_0)-x_1^{*}}{x(t_0+T_0)-\widehat{x_1^{*}}}\right)e^{-\frac{r}{K}(x_1^{*}-\widehat{x_1^{*}})(t-t_0-T_0)}}>\frac{x_1^{*}+\widehat{x_1^{*}}\left(\frac{x_1^{*}-\gamma {\nu }_1}{\gamma {\nu }_1-\widehat{x_1^{*}}}\right)e^{-\frac{r}{K}(x_1^{*}-\widehat{x_1^{*}})(t-t_0-T_0)}}{1+\left(\frac{x_1^{*}-\gamma {\nu }_1}{\gamma {\nu }_1-\widehat{x_1^{*}}}\right)e^{-\frac{r}{K}(x_1^{*}-\widehat{x_1^{*}})(t-t_0-T_0)}}.\end{array}$$

Hence, we have, for $t\ge t_0+T_0+T_1$,

$$x(t)>\frac{x_1^{*}+\widehat{x_1^{*}}\left(\frac{x_1^{*}-\gamma {\nu }_1}{\gamma {\nu }_1-\widehat{x_1^{*}}}\right)e^{-\frac{r}{K}(x_1^{*}-\widehat{x_1^{*}})T_1}}{1+\left(\frac{x_1^{*}-\gamma {\nu }_1}{\gamma {\nu }_1-\widehat{x_1^{*}}}\right)e^{-\frac{r}{K}(x_1^{*}-\widehat{x_1^{*}})T_1}}=x^0>x^{*}.$$

(8)

Define $m=min\{\overline{y},\frac{u\overline{v}}{k}\}>0$, where

$$\begin{array}{c}\hfill \overline{y}=\underset{\theta \in [-({\tau }_1+{\tau }_2),0]}{min}y(T_{*}+\theta )>0,\overline{v}=\underset{\theta \in [-({\tau }_1+{\tau }_2),0]}{min}v(T_{*}+\theta )>0,T_{*}=t_0+T_0+T_1+{\tau }_1+{\tau }_2.\end{array}$$

Next, we show that $y(t)\ge m$ for $t\ge t_0+T_0+T_1$. In fact, otherwise, there exists a $\widehat{T_1}\ge 0$ such that $y(t)\ge m$ for $t_0+T_0+T_1\le t\le T_{*}+\widehat{T_1}$, $y(T_{*}+\widehat{T_1})=m$ and $\dot{y}(T_{*}+\widehat{T_1})\le 0$. Then, from the third equation of model (4), we have, for $t_0+T_0+T_1\le t\le T_{*}+\widehat{T_1}$, $\dot{v}(t)\ge km-uv(t)$, which implies that, for $t_0+T_0+T_1\le t\le T_{*}+\widehat{T_1}$,

$$\begin{array}{cc}\hfill v(t)& \ge \left(v(t_0+T_0+T_1)-\frac{km}{u}\right)e^{-u(t-t_0-T_0-T_1)}+\frac{km}{u}\hfill \\ \hfill & \ge (v(t_0+T_0+T_1)-\overline{v})e^{-u(t-t_0-T_0-T_1)}+\frac{km}{u}\hfill \\ \hfill & \ge \frac{km}{u}.\hfill \end{array}$$

(9)

Therefore, from (8) and (9), we have, for $t=T_{*}+\widehat{T_1}$,

$$\begin{array}{cc}\hfill \dot{y}(t)& ={\delta }_{1}x(t-{\tau }_1)y(t-{\tau }_1)+{\delta }_{2}x(t-{\tau }_2)v(t-{\tau }_2)-pm\hfill \\ \hfill & \ge {\delta }_1{x}^{0}m+{\delta }_2{x}^0\frac{km}{u}-pm\hfill \\ \hfill & =pm\left(\frac{x^0}{x^{*}}-1\right)>0.\hfill \end{array}$$

This is a contradiction to $\dot{y}(T_{*}+\widehat{T_1})\le 0$. This shows that for $t\ge t_0+T_0+T_1$, $y(t)\ge m$.

Using the derivation completely similar to (9), we have, for $t\ge t_0+T_0+T_1$, $v(t)\ge \frac{km}{u}$.

Consider the following auxiliary function:

$$V(t)=y(t)+\frac{{\delta }_2}{u}{x}^{*}v(t)+{\delta }_1{\int }_{t-{\tau }_1}^{t}x(\theta )y(\theta )d\theta +{\delta }_2{\int }_{t-{\tau }_2}^{t}x(\theta )v(\theta )d\theta .$$

Then, we have, for $t\ge t_0+T_0+T_1$,

$$\dot{V}(t)={\delta }_1(x(t)-x^{*})y(t)+{\delta }_2(x(t)-x^{*})v(t)\ge m\left({\delta }_1+\frac{k}{u}{\delta }_2\right)(x^0-x^{*})>0,$$

which leads to, for $t\ge t_0+T_0+T_1$,

$$V(t)\ge V(t_0+T_0+T_1)+m\left({\delta }_1+\frac{k}{u}{\delta }_2\right)(x^0-x^{*})(t-t_0-T_0-T_1),$$

which implies that $V(t)\to +\infty (t\to +\infty )$. This is a contradiction with the boundedness of $V(t)$. Therefore, the claim is proved.

Below, there are two remaining cases that need to be discussed.

(i) $y(t)\ge \frac{y^{*}}{2}$ for sufficiently large t. (ii) $y(t)$ oscillates about $\frac{y^{*}}{2}$ for sufficiently large t.

For the case (i), it clearly has $lim{inf}_{t\to +\infty }y(t)\ge {\nu }_2$.

For the case (ii), let $t_1$ and $t_2$ be sufficiently large such that $y(t_1)=y(t_2)=\frac{y^{*}}{2},y(t)<\frac{y^{*}}{2}(t_1<t<t_2).$

If $t_2-t_1\le \varpi$, from the second equation of model (4), we have, for $t_1\le t\le t_2$, $\dot{y}(t)\ge -py(t)$, which implies that for $t_1\le t\le t_2$,

$$y(t)\ge y(t_1)e^{-p(t-t_1)}\ge \frac{y^{*}}{2}{e}^{-p\varpi }={\nu }_2.$$

If $t_2-t_1>\varpi$, it is easily obtained that, for $t_1\le t\le t_1+\varpi$, $y(t)\ge {\nu }_2$. Further, we prove that for $t_1+\varpi \le t\le t_2$, $y(t)\ge {\nu }_2$.

If not, there exists a $\widehat{T_2}\ge 0$ such that for $t_1\le t\le t_1+\varpi +\widehat{T_2}(<t_2)$, $y(t)\ge {\nu }_2$, $y(t_1+\varpi +\widehat{T_2})={\nu }_2$ and $\dot{y}(t_1+\varpi +\widehat{T_2})\le 0$. Using the derivation method similar to (8), treating $t_1$ as $t_0$, we have, for $t_1+T_0+T_1\le t\le t_1+\varpi +\widehat{T_2}$, $x(t)>x^0>x^{*}$. Then, let us prove that there exists $\overline{t}\in [t_1,t_1+T_0+T_1+T_2]$ such that $v(\overline{t})\ge \frac{qk{\nu }_2}{u}$. If not, then we have, for $t_1\le t\le t_1+T_0+T_1+T_2$, $v(t)<\frac{qk{\nu }_2}{u}$. From the third equation of model (4), we have, for $t_1\le t\le t_1+T_0+T_1+T_2$,

$$\dot{v}(t)\ge k{\nu }_2-uv(t)\ge k{\nu }_2(1-q),$$

which implies that, for $t=t_1+T_0+T_1+T_2$,

$$v(t)\ge v(t_1)+k{\nu }_2(1-q)(t-t_1)>v(t_1)+k{\nu }_2(1-q)T_2>\frac{qk{\nu }_2}{u},$$

which is a contradiction. Hence, we conclude that there exists $\overline{t}\in [t_1,t_1+T_0+T_1+T_2]$ such that $v(\overline{t})\ge \frac{qk{\nu }_2}{u}$.

Note that, for $\overline{t}\le t\le t_1+\varpi +\widehat{T_2}$, $\dot{v}(t)\ge qk{\nu }_2-uv(t)$, which implies that, for $\overline{t}\le t\le t_1+\varpi +\widehat{T_2}$,

$$v(t)\ge \left(v(\overline{t})-\frac{qk{\nu }_2}{u}\right)e^{-u(t-t_1)}+\frac{qk{\nu }_2}{u}\ge \frac{qk{\nu }_2}{u}.$$

Hence, we have from the second equation of model (4) that for $t=t_1+\varpi +\widehat{T_2}$,

$$\begin{array}{cc}\hfill \dot{y}(t)=& {\delta }_{1}x(t-{\tau }_1)y(t-{\tau }_1)+{\delta }_{2}x(t-{\tau }_2)v(t-{\tau }_2)-p{\nu }_2\hfill \\ \hfill \ge & p{\nu }_2\left(\frac{{\delta }_1}{p}x(t-{\tau }_1)+\frac{qk{\delta }_2}{up}x(t-{\tau }_2)-1\right)\hfill \\ \hfill >& p{\nu }_2\left(\frac{q{x}^0}{x^{*}}-1\right)=0.\hfill \end{array}$$

(10)

This is a contradiction to $\dot{y}(t_1+\varpi +\widehat{T_2})\le 0$. Based on the above analysis, we have $y(t)\ge {\nu }_2$ for $t\in [t_1,t_2]$. Since this kind of interval $[t_1,t_2]$ is chosen in an arbitrary way, we conclude that $y(t)\ge {\nu }_2$ for any sufficiently large t. Thus, $lim{inf}_{t\to \infty }y(t)\ge {\nu }_2$.

Finally, according to the third equation of model (4), we have $lim{inf}_{t\to \infty }v(t)\ge \frac{k{\nu }_2}{u}={\nu }_3$. □

3. Global Attractivity of the Chronic Infection Equilibrium

Now, we continue to discuss global attractivity of the chronic infection equilibrium $E^{*}$. The following lemma is used.

Lemma 1.

(Barbalat’s lemma [42]) Let $q(t)$ be a real valued differentiable function defined on some half line $[\vartheta ,+\infty )$, $\vartheta \in (-\infty ,+\infty )$. If $li{m}_{t\to +\infty }q(t)={\vartheta }_1$ ($|{\vartheta }_1|<+\infty$) and $\dot{q}(t)$ is uniformly continuous for $t>\vartheta$, then $li{m}_{t\to +\infty }\dot{q}(t)=0$.

Due to the technical requirements of the proof, we assume that $m_1{\tau }_1=m_2{\tau }_2$.

For any sufficient small $0<\epsilon <{\nu }_1$, we define

$$\begin{array}{cc}\hfill {\nu }_1(\epsilon )=& {\nu }_1-\epsilon ,x_0(\epsilon )=x_0+\epsilon ,M_2(\epsilon )=M_2+\epsilon ,M_3(\epsilon )=M_3+\epsilon ,{\mathsf{\Lambda }}_1=\frac{r\alpha }{K}+c+{\beta }_1,\hfill \\ \hfill {\mathsf{{\rm Y}}}_1(\epsilon )=& d-r+\frac{r}{K}(x_0(\epsilon )+x^{*})+{\mathsf{\Lambda }}_1{y}^{*}+{\beta }_2{v}^{*},{\mathsf{{\rm Y}}}_2(\epsilon )={\beta }_1{x}_0(\epsilon )+{\beta }_1{y}^{*}+{\beta }_2{x}_0(\epsilon )+{\beta }_2{v}^{*}+p{e}^{m_1{\tau }_1},\hfill \\ \hfill {\mathsf{\Psi }}_1(\epsilon )=& \frac{1}{2}{\beta }_1[x_0(\epsilon )e^{-m_1{\tau }_1}{\mathsf{{\rm Y}}}_2(\epsilon )+M_2(\epsilon )({\mathsf{{\rm Y}}}_1(\epsilon )+{\mathsf{\Lambda }}_1{x}_0(\epsilon )+{\beta }_2{x}_0(\epsilon ))],\hfill \\ \hfill {\mathsf{\Psi }}_2(\epsilon )=& \frac{1}{2}{\beta }_2[M_3(\epsilon )({\mathsf{{\rm Y}}}_1(\epsilon )+{\mathsf{\Lambda }}_1{x}_0(\epsilon )+{\beta }_2{x}_0(\epsilon ))+(k+u)x_0(\epsilon )],\hfill \\ \hfill {\mathsf{\Psi }}_3(\epsilon )=& \frac{1}{2}{\beta }_1^2{x}_0(\epsilon )y^{*}(1+e^{-m_1{\tau }_1}),{\mathsf{\Psi }}_4(\epsilon )=\frac{1}{2}{\beta }_1{\beta }_2{x}_0(\epsilon )v^{*}(1+e^{-m_1{\tau }_1}),\hfill \\ \hfill {\mathsf{\Psi }}_5(\epsilon )=& \frac{1}{2}{\mathsf{{\rm Y}}}_1(\epsilon ){\beta }_1{M}_2(\epsilon )(1+e^{m_1{\tau }_1}),{\mathsf{\Psi }}_6(\epsilon )=\frac{1}{2}{\beta }_1^2{x}_0^2(\epsilon )(1+e^{-m_1{\tau }_1}),\hfill \\ \hfill {\mathsf{\Psi }}_7(\epsilon )=& \frac{1}{2}p{\beta }_1{x}_0(\epsilon )(1+e^{m_1{\tau }_1})+\frac{1}{2}{\mathsf{\Lambda }}_1{x}_0(\epsilon ){\beta }_1{M}_2(\epsilon )(1+e^{m_1{\tau }_1}),{\mathsf{\Psi }}_8(\epsilon )=\frac{1}{2}{\beta }_1{\beta }_2{x}_0^2(\epsilon )(1+e^{-m_1{\tau }_1}),\hfill \\ \hfill {\mathsf{\Psi }}_9(\epsilon )=& \frac{1}{2}{\beta }_2{x}_0(\epsilon ){\beta }_1{M}_2(\epsilon )(1+e^{m_1{\tau }_1}),{\mathsf{\Psi }}_{10}(\epsilon )=\frac{1}{2}{\mathsf{{\rm Y}}}_1(\epsilon ){\beta }_2{M}_3(\epsilon )(1+e^{m_1{\tau }_1}),\hfill \\ \hfill {\mathsf{\Psi }}_{11}(\epsilon )=& \frac{1}{2}k{\beta }_2{x}_0(\epsilon )(1+e^{m_1{\tau }_1})+\frac{1}{2}{\mathsf{\Lambda }}_1{x}_0(\epsilon ){\beta }_2{M}_3(\epsilon )(1+e^{m_1{\tau }_1}),\hfill \\ \hfill {\mathsf{\Psi }}_{12}(\epsilon )=& \frac{1}{2}u{\beta }_2{x}_0(\epsilon )(1+e^{m_1{\tau }_1})+\frac{1}{2}{\beta }_2{x}_0(\epsilon ){\beta }_2{M}_3(\epsilon )(1+e^{m_1{\tau }_1}).\hfill \end{array}$$

Let us define the real symmetric matrices as follows,

$$J(\epsilon )=\left(\begin{array}{ccc}{A}_{11}(\epsilon )& -A_{12}(\epsilon )& 0\\ -A_{12}(\epsilon )& A_{22}(\epsilon )& -A_{23}(\epsilon )\\ 0& -A_{23}(\epsilon )& A_{33}(\epsilon )\end{array}\right),$$

where

$$\begin{array}{cc}\hfill A_{11}(\epsilon )=& \frac{1}{x_0(\epsilon )}\left[d-r+\frac{r{x}^{*}}{K}+\left(\frac{r\alpha }{K}+c\right)y^{*}\right]+\frac{r}{K}+{\theta }_1\{d-r+\frac{r}{K}({\nu }_1(\epsilon )+x^{*})+\left(c+\frac{r\alpha }{K}\right)y^{*}\hfill \\ \hfill & -[({\mathsf{\Psi }}_1(\epsilon )+{\mathsf{\Psi }}_3(\epsilon )+{\mathsf{\Psi }}_4(\epsilon )+{\mathsf{\Psi }}_5(\epsilon )){\tau }_1+({\mathsf{\Psi }}_2(\epsilon )+{\mathsf{\Psi }}_{10}(\epsilon )){\tau }_2]\},\hfill \\ \hfill A_{12}(\epsilon )=& \frac{1}{2}\left(c+\frac{r\alpha }{K}\right)+\frac{1}{2}{\theta }_1\left\{e^{m_1{\tau }_1}\left[d-r+\frac{r}{K}(x_0(\epsilon )+x^{*})+\left(c+\frac{r\alpha }{K}\right)y^{*}+p\right]+\left(c+\frac{r\alpha }{K}\right)x_0(\epsilon )\right\},\hfill \\ \hfill A_{22}(\epsilon )=& {\theta }_1\{e^{m_1{\tau }_1}\left[e^{m_1{\tau }_1}p+\left(c+\frac{r\alpha }{K}\right){\nu }_1(\epsilon )\right]\hfill \\ \hfill & -[({\mathsf{\Psi }}_1(\epsilon )e^{m_1{\tau }_1}+{\mathsf{\Psi }}_6(\epsilon )+{\mathsf{\Psi }}_7(\epsilon )){\tau }_1+({\mathsf{\Psi }}_2(\epsilon )e^{m_1{\tau }_1}+{\mathsf{\Psi }}_{11}(\epsilon )){\tau }_2]\},\hfill \\ \hfill A_{23}(\epsilon )=& \frac{1}{2}{\theta }_{2}k,A_{33}(\epsilon )={\theta }_{2}u-{\theta }_1[({\mathsf{\Psi }}_8(\epsilon )+{\mathsf{\Psi }}_9(\epsilon )){\tau }_1+{\mathsf{\Psi }}_{12}(\epsilon ){\tau }_2],\hfill \end{array}$$

where ${\theta }_1$ and ${\theta }_2$ are arbitrary positive constants.

Theorem 2.

If $R_0>1$, $d-r+\frac{r{x}^{*}}{K}\ge 0$, $m_1{\tau }_1=m_2{\tau }_2$ and matrix $J(0)$ is positive definite, then the chronic infection equilibrium $E^{*}$ is globally attractive in $X^{+}$.

Proof. 

Let $(x(t),y(t),v(t))\in X^{+}$$(t\ge 0)$ be any solution of model (4). If $J(0)$ is positive definite, then $J(\epsilon )$ is also positive definite for any sufficiently small $0<\epsilon <{\nu }_1$. For the above $\epsilon$, there exists a sufficiently large $T(\epsilon )>{\tau }_1+{\tau }_2$ such that, for $t>T(\epsilon )$,

$$\begin{array}{c}\hfill 0<{\nu }_1(\epsilon )<x(t)<x_0(\epsilon ),y(t)<M_2(\epsilon ),v(t)<M_3(\epsilon ).\end{array}$$

Let $g(z)=z-1-lnz$$(z>0)$. Clearly, $g(z)\ge 0$$(z>0)$, and $g(z)=0$ if and only if $z=1$. Define

$$\begin{array}{cc}\hfill U_1=& x^{*}g\left(\frac{x(t)}{x^{*}}\right)+e^{m_1{\tau }_1}{y}^{*}g\left(\frac{y(t)}{y^{*}}\right)+\frac{{\beta }_2{x}^{*}{v}^{*}}{k{y}^{*}}{v}^{*}g\left(\frac{v(t)}{v^{*}}\right)\hfill \\ \hfill & +{\beta }_1{x}^{*}{y}^{*}{\int }_{t-{\tau }_1}^{t}g\left(\frac{x(s)y(s)}{x^{*}{y}^{*}}\right)ds+{\beta }_2{x}^{*}{v}^{*}{\int }_{t-{\tau }_2}^{t}g\left(\frac{x(s)v(s)}{x^{*}{v}^{*}}\right)ds.\hfill \end{array}$$

Hence, $U_1$ is positive definite with respect to the chronic infection equilibrium $E^{*}=(x^{*},y^{*},v^{*})$. Similar to the calculation in [18,19], for $t\ge 0$, the derivative along the solution of model (4) satisfies

$$\begin{array}{cc}\hfill \frac{d{U}_1}{dt}=& -\frac{1}{x(t)}\left(d-r+\frac{r{x}^{*}}{K}\right){(x(t)-x^{*})}^2-\frac{r}{K}{(x(t)-x^{*})}^2\hfill \\ \hfill & +\left(\frac{r\alpha }{K}+c\right)\left(1-\frac{x^{*}}{x(t)}\right)(x^{*}{y}^{*}-x(t)y(t))\hfill \\ \hfill & +{\beta }_1{x}^{*}{y}^{*}\left\{-g\left(\frac{x(t-{\tau }_1)y(t-{\tau }_1)}{x^{*}y(t)}\right)-g\left(\frac{x^{*}}{x(t)}\right)\right\}\hfill \\ \hfill & +{\beta }_2{x}^{*}{v}^{*}\left\{-g\left(\frac{x(t-{\tau }_2)v(t-{\tau }_2)y^{*}}{x^{*}{v}^{*}y(t)}\right)-g\left(\frac{x^{*}}{x(t)}\right)-g\left(\frac{y(t)v^{*}}{y^{*}v(t)}\right)\right\}\hfill \\ \hfill =& -\frac{1}{x(t)}\left[d-r+\frac{r{x}^{*}}{K}+\left(\frac{r\alpha }{K}+c\right)y^{*}\right]{(x(t)-x^{*})}^2-\frac{r}{K}{(x(t)-x^{*})}^2\hfill \\ \hfill & -\left(\frac{r\alpha }{K}+c\right)(x(t)-x^{*})(y(t)-y^{*})\hfill \\ \hfill & +{\beta }_1{x}^{*}{y}^{*}\left\{-g\left(\frac{x(t-{\tau }_1)y(t-{\tau }_1)}{x^{*}y(t)}\right)-g\left(\frac{x^{*}}{x(t)}\right)\right\}\hfill \\ \hfill & +{\beta }_2{x}^{*}{v}^{*}\left\{-g\left(\frac{x(t-{\tau }_2)v(t-{\tau }_2)y^{*}}{x^{*}{v}^{*}y(t)}\right)-g\left(\frac{x^{*}}{x(t)}\right)-g\left(\frac{y(t)v^{*}}{y^{*}v(t)}\right)\right\}.\hfill \end{array}$$

(11)

It is worth mentioning that if $\alpha =c=0$ and $d-r+\frac{r{x}^{*}}{K}\ge 0$, then $\frac{d{U}_1}{dt}\le 0$, which leads to $E^{*}$ is stable (see [37,38]). Then, it follows from [19] that $E^{*}$ is globally attractive. Thus, $E^{*}$ is globally asymptotically stable.

If $\alpha$ and c are not 0 at the same time, inspired by [43,44], we define

$$U_2=\frac{1}{2}{[(x(t)-x^{*})+e^{m_1{\tau }_1}(y(t)-y^{*})]}^2.$$

From model (4), we have, for $t>T(\epsilon )$,

$$\begin{array}{cc}\hfill \dot{x}(t)+e^{m_1{\tau }_1}\dot{y}(t)=& s-dx(t)+rx(t)\left(1-\frac{x(t)}{K}\right)-\left(c+\frac{r\alpha }{K}\right)x(t)y(t)-e^{m_1{\tau }_1}py(t)\hfill \\ \hfill & +{\beta }_{1}x(t-{\tau }_1)y(t-{\tau }_1)-{\beta }_{1}x(t)y(t)+{\beta }_{2}x(t-{\tau }_2)v(t-{\tau }_2)-{\beta }_{2}x(t)v(t)\hfill \\ \hfill =& -\left[d-r+\frac{r}{K}(x(t)+x^{*})+\left(c+\frac{r\alpha }{K}\right)y^{*}\right](x(t)-x^{*})\hfill \\ \hfill & -\left[e^{m_1{\tau }_1}p+\left(c+\frac{r\alpha }{K}\right)x(t)\right](y(t)-y^{*})\hfill \\ \hfill & +{\beta }_{1}x(t-{\tau }_1)(y(t-{\tau }_1)-y(t))+{\beta }_{1}y(t)(x(t-{\tau }_1)-x(t))\hfill \\ \hfill & +{\beta }_{2}x(t-{\tau }_2)(v(t-{\tau }_2)-v(t))+{\beta }_{2}v(t)(x(t-{\tau }_2)-x(t)).\hfill \end{array}$$

Further, we have, for $t>T(\epsilon )$,

$$\begin{array}{cc}\hfill \frac{d{U}_2}{dt}=& [(x(t)-x^{*})+e^{m_1{\tau }_1}(y(t)-y^{*})](\dot{x}(t)+e^{m_1{\tau }_1}\dot{y}(t))\hfill \\ \hfill =& -\left[d-r+\frac{r}{K}(x(t)+x^{*})+\left(c+\frac{r\alpha }{K}\right)y^{*}\right]{(x(t)-x^{*})}^2\hfill \\ \hfill & -e^{m_1{\tau }_1}\left[e^{m_1{\tau }_1}p+\left(c+\frac{r\alpha }{K}\right)x(t)\right]{(y(t)-y^{*})}^2\hfill \\ \hfill & -\left\{e^{m_1{\tau }_1}\left[d-r+\frac{r}{K}(x(t)+x^{*})+\left(c+\frac{r\alpha }{K}\right)y^{*}+p\right]+\left(c+\frac{r\alpha }{K}\right)x(t)\right\}\hfill \\ \hfill & \times (x(t)-x^{*})(y(t)-y^{*})+{\mathsf{\Gamma }}_1(t)+{\mathsf{\Gamma }}_2(t)+{\mathsf{\Gamma }}_3(t)+{\mathsf{\Gamma }}_4(t),\hfill \end{array}$$

(12)

where

$$\begin{array}{cc}\hfill {\mathsf{\Gamma }}_1(t)=& -{\beta }_{1}x(t-{\tau }_1)[(x(t)-x^{*})+e^{m_1{\tau }_1}(y(t)-y^{*})]{\int }_{t-{\tau }_1}^t\dot{y}(s)ds,\hfill \\ \hfill {\mathsf{\Gamma }}_2(t)=& -{\beta }_{1}y(t)[(x(t)-x^{*})+e^{m_1{\tau }_1}(y(t)-y^{*})]{\int }_{t-{\tau }_1}^t\dot{x}(s)ds,\hfill \\ \hfill {\mathsf{\Gamma }}_3(t)=& -{\beta }_{2}x(t-{\tau }_2)[(x(t)-x^{*})+e^{m_1{\tau }_1}(y(t)-y^{*})]{\int }_{t-{\tau }_2}^t\dot{v}(s)ds,\hfill \\ \hfill {\mathsf{\Gamma }}_4(t)=& -{\beta }_{2}v(t)[(x(t)-x^{*})+e^{m_1{\tau }_1}(y(t)-y^{*})]{\int }_{t-{\tau }_2}^t\dot{x}(s)ds.\hfill \end{array}$$

Since $E^{*}$ is a positive equilibrium of model (4), $\dot{x}(t)$, $\dot{y}(t)$ and $\dot{v}(t)$ can be rewritten as

$$\begin{array}{cc}\hfill \dot{x}(t)=& -\left[d-r+\frac{r}{K}(x(t)+x^{*})\right](x(t)-x^{*})+{\mathsf{\Lambda }}_1{x}^{*}{y}^{*}-{\mathsf{\Lambda }}_{1}x(t)y(t)+{\beta }_2{x}^{*}{v}^{*}-{\beta }_{2}x(t)v(t)\hfill \\ \hfill =& -\left[d-r+\frac{r}{K}(x(t)+x^{*})\right](x(t)-x^{*})+{\mathsf{\Lambda }}_1(x^{*}-x(t))y^{*}+{\mathsf{\Lambda }}_{1}x(t)(y^{*}-y(t))\hfill \\ \hfill & +{\beta }_2(x^{*}-x(t))v^{*}+{\beta }_{2}x(t)(v^{*}-v(t))\hfill \\ \hfill =& -\tilde{{\mathsf{{\rm Y}}}_1}(t)(x(t)-x^{*})+{\mathsf{\Lambda }}_{1}x(t)(y^{*}-y(t))+{\beta }_{2}x(t)(v^{*}-v(t)),\hfill \end{array}$$

(13)

where $\tilde{{\mathsf{{\rm Y}}}_1}(t)=d-r+\frac{r}{K}(x(t)+x^{*})+{\mathsf{\Lambda }}_1{y}^{*}+{\beta }_2{v}^{*}$,

$$\begin{array}{cc}\hfill \dot{y}(t)=& {\beta }_1{e}^{-m_1{\tau }_1}x(t-{\tau }_1)(y(t-{\tau }_1)-y^{*})+{\beta }_1{e}^{-m_1{\tau }_1}{y}^{*}(x(t-{\tau }_1)-x^{*})\hfill \\ \hfill & +{\beta }_2{e}^{-m_2{\tau }_2}x(t-{\tau }_2)(v(t-{\tau }_2)-v^{*})+{\beta }_2{e}^{-m_2{\tau }_2}{v}^{*}(x(t-{\tau }_2)-x^{*})+p(y^{*}-y(t)),\hfill \end{array}$$

(14)

$$\begin{array}{c}\hfill \dot{v}(t)=k(y(t)-y^{*})+u(v^{*}-v(t)).\end{array}$$

(15)

From (14), we have, for $t>T(\epsilon )+2({\tau }_1+{\tau }_2)$,

$$\begin{array}{cc}\hfill |{\mathsf{\Gamma }}_1(t)|=& {\beta }_{1}x(t-{\tau }_1)\left|[(x(t)-x^{*})+e^{m_1{\tau }_1}(y(t)-y^{*})]\right|\hfill \\ \hfill & \times \left|{\int }_{t-{\tau }_1}^t\right\{{\beta }_1{e}^{-m_1{\tau }_1}x(s-{\tau }_1)(y(s-{\tau }_1)-y^{*})+{\beta }_1{e}^{-m_1{\tau }_1}{y}^{*}(x(s-{\tau }_1)-x^{*})\hfill \\ \hfill & +{\beta }_2{e}^{-m_2{\tau }_2}x(s-{\tau }_2)(v(s-{\tau }_2)-v^{*})+{\beta }_2{e}^{-m_2{\tau }_2}{v}^{*}(x(s-{\tau }_2)-x^{*})+p(y^{*}-y(s))\left\}ds\right|\hfill \\ \hfill \le & {\beta }_1{x}_0(\epsilon ){\int }_{t-{\tau }_1}^t\{\frac{{\beta }_1{e}^{-m_1{\tau }_1}{x}_0(\epsilon )}{2}[{(x(t)-x^{*})}^2+{(y(s-{\tau }_1)-y^{*})}^2]\hfill \\ \hfill & +\frac{{\beta }_1{x}_0(\epsilon )}{2}[{(y(t)-y^{*})}^2+{(y(s-{\tau }_1)-y^{*})}^2]\hfill \\ \hfill & +\frac{{\beta }_1{e}^{-m_1{\tau }_1}{y}^{*}}{2}[{(x(t)-x^{*})}^2+{(x(s-{\tau }_1)-x^{*})}^2]+\frac{{\beta }_1{y}^{*}}{2}[{(y(t)-y^{*})}^2+{(x(s-{\tau }_1)-x^{*})}^2]\hfill \\ \hfill & +\frac{{\beta }_2{e}^{-m_2{\tau }_2}{x}_0(\epsilon )}{2}[{(x(t)-x^{*})}^2+{(v(s-{\tau }_2)-v^{*})}^2]+\frac{{\beta }_2{x}_0(\epsilon )}{2}[{(y(t)-y^{*})}^2+{(v(s-{\tau }_2)-v^{*})}^2]\hfill \\ \hfill & +\frac{{\beta }_2{e}^{-m_2{\tau }_2}{v}^{*}}{2}[{(x(t)-x^{*})}^2+{(x(s-{\tau }_2)-x^{*})}^2]+\frac{{\beta }_2{v}^{*}}{2}[{(y(t)-y^{*})}^2+{(x(s-{\tau }_2)-x^{*})}^2]\hfill \\ \hfill & +\frac{p}{2}[{(x(t)-x^{*})}^2+{(y(s)-y^{*})}^2]+\frac{p{e}^{m_1{\tau }_1}}{2}[{(y(t)-y^{*})}^2+{(y(s)-y^{*})}^2]\}ds\hfill \\ \hfill =& \frac{1}{2}{\beta }_1{x}_0(\epsilon )e^{-m_1{\tau }_1}{\mathsf{{\rm Y}}}_2(\epsilon ){(x(t)-x^{*})}^2{\tau }_1+\frac{1}{2}{\beta }_1{x}_0(\epsilon ){\mathsf{{\rm Y}}}_2(\epsilon ){(y(t)-y^{*})}^2{\tau }_1\hfill \\ \hfill & +\frac{1}{2}{\beta }_1^2{x}_0^2(\epsilon )(1+e^{-m_1{\tau }_1}){\int }_{t-{\tau }_1}^t{(y(s-{\tau }_1)-y^{*})}^{2}ds\hfill \\ \hfill & +\frac{1}{2}{\beta }_1^2{x}_0(\epsilon )y^{*}(1+e^{-m_1{\tau }_1}){\int }_{t-{\tau }_1}^t{(x(s-{\tau }_1)-x^{*})}^{2}ds\hfill \\ \hfill & +\frac{1}{2}{\beta }_1{\beta }_2{x}_0^2(\epsilon )(1+e^{-m_1{\tau }_1}){\int }_{t-{\tau }_1}^t{(v(s-{\tau }_2)-v^{*})}^{2}ds\hfill \\ \hfill & +\frac{1}{2}{\beta }_1{\beta }_2{x}_0(\epsilon )v^{*}(1+e^{-m_1{\tau }_1}){\int }_{t-{\tau }_1}^t{(x(s-{\tau }_2)-x^{*})}^{2}ds\hfill \\ \hfill & +\frac{1}{2}p{\beta }_1{x}_0(\epsilon )(1+e^{m_1{\tau }_1}){\int }_{t-{\tau }_1}^t{(y(s)-y^{*})}^{2}ds.\hfill \end{array}$$

Similarly, from (13) and (15), we have, for $t>T(\epsilon )+2({\tau }_1+{\tau }_2)$,

$$\begin{array}{cc}\hfill |{\mathsf{\Gamma }}_2(t)|\le & {\beta }_1{M}_2(\epsilon )\left|[(x(t)-x^{*})+e^{m_1{\tau }_1}(y(t)-y^{*})]\right|\hfill \\ \hfill & \times {\int }_{t-{\tau }_1}^t\left\{{\mathsf{{\rm Y}}}_1(\epsilon )|(x(s)-x^{*})|+{\mathsf{\Lambda }}_1{x}_0(\epsilon )|(y(s)-y^{*})|+{\beta }_2{x}_0(\epsilon )\left|(v(s)-v^{*})\right|\right\}ds\hfill \\ \hfill \le & \frac{1}{2}{\beta }_1{M}_2(\epsilon )({\mathsf{{\rm Y}}}_1(\epsilon )+{\mathsf{\Lambda }}_1{x}_0(\epsilon )+{\beta }_2{x}_0(\epsilon )){(x(t)-x^{*})}^2{\tau }_1\hfill \\ \hfill & +\frac{1}{2}{\beta }_1{M}_2(\epsilon )e^{m_1{\tau }_1}({\mathsf{{\rm Y}}}_1(\epsilon )+{\mathsf{\Lambda }}_1{x}_0(\epsilon )+{\beta }_2{x}_0(\epsilon )){(y(t)-y^{*})}^2{\tau }_1\hfill \\ \hfill & +\frac{1}{2}{\mathsf{{\rm Y}}}_1(\epsilon ){\beta }_1{M}_2(\epsilon )(1+e^{m_1{\tau }_1}){\int }_{t-{\tau }_1}^t{(x(s)-x^{*})}^{2}ds\hfill \\ \hfill & +\frac{1}{2}{\mathsf{\Lambda }}_1{x}_0(\epsilon ){\beta }_1{M}_2(\epsilon )(1+e^{m_1{\tau }_1}){\int }_{t-{\tau }_1}^t{(y(s)-y^{*})}^{2}ds\hfill \\ \hfill & +\frac{1}{2}{\beta }_2{x}_0(\epsilon ){\beta }_1{M}_2(\epsilon )(1+e^{m_1{\tau }_1}){\int }_{t-{\tau }_1}^t{(v(s)-v^{*})}^{2}ds,\hfill \end{array}$$

$$\begin{array}{cc}\hfill |{\mathsf{\Gamma }}_3(t)|\le & {\beta }_2{x}_0(\epsilon )\left[\right|(x(t)-x^{*})|+e^{m_1{\tau }_1}|(y(t)-y^{*})\left|\right]{\int }_{t-{\tau }_2}^t\left[k\right|(y(s)-y^{*})|+u|(v^{*}-v(s))\left|\right]ds,\hfill \\ \hfill \le & \frac{1}{2}{\beta }_2{x}_0(\epsilon )(k+u){(x(t)-x^{*})}^2{\tau }_2+\frac{1}{2}{\beta }_2{x}_0(\epsilon )(k+u)e^{m_1{\tau }_1}{(y(t)-y^{*})}^2{\tau }_2\hfill \\ \hfill & +\frac{1}{2}k{\beta }_2{x}_0(\epsilon )(1+e^{m_1{\tau }_1}){\int }_{t-{\tau }_2}^t{(y(s)-y^{*})}^{2}ds+\frac{1}{2}u{\beta }_2{x}_0(\epsilon )(1+e^{m_1{\tau }_1}){\int }_{t-{\tau }_2}^t{(v(s)-v^{*})}^{2}ds,\hfill \end{array}$$

$$\begin{array}{cc}\hfill |{\mathsf{\Gamma }}_4(t)|\le & {\beta }_2{M}_3(\epsilon )\left[\right|(x(t)-x^{*})|+e^{m_1{\tau }_1}|(y(t)-y^{*})\left|\right]\hfill \\ \hfill & \times {\int }_{t-{\tau }_2}^t\left\{{\mathsf{{\rm Y}}}_1(\epsilon )|(x(s)-x^{*})|+{\mathsf{\Lambda }}_1{x}_0(\epsilon )|(y(s)-y^{*})|+{\beta }_2{x}_0(\epsilon )\left|(v(s)-v^{*})\right|\right\}ds\hfill \\ \hfill \le & \frac{1}{2}{\beta }_2{M}_3(\epsilon )({\mathsf{{\rm Y}}}_1(\epsilon )+{\mathsf{\Lambda }}_1{x}_0(\epsilon )+{\beta }_2{x}_0(\epsilon )){(x(t)-x^{*})}^2{\tau }_2\hfill \\ \hfill & +\frac{1}{2}{\beta }_2{M}_3(\epsilon )e^{m_1{\tau }_1}({\mathsf{{\rm Y}}}_1(\epsilon )+{\mathsf{\Lambda }}_1{x}_0(\epsilon )+{\beta }_2{x}_0(\epsilon )){(y(t)-y^{*})}^2{\tau }_2\hfill \\ \hfill & +\frac{1}{2}{\mathsf{{\rm Y}}}_1(\epsilon ){\beta }_2{M}_3(\epsilon )(1+e^{m_1{\tau }_1}){\int }_{t-{\tau }_2}^t{(x(s)-x^{*})}^{2}ds\hfill \\ \hfill & +\frac{1}{2}{\mathsf{\Lambda }}_1{x}_0(\epsilon ){\beta }_2{M}_3(\epsilon )(1+e^{m_1{\tau }_1}){\int }_{t-{\tau }_2}^t{(y(s)-y^{*})}^{2}ds\hfill \\ \hfill & +\frac{1}{2}{\beta }_2{x}_0(\epsilon ){\beta }_2{M}_3(\epsilon )(1+e^{m_1{\tau }_1}){\int }_{t-{\tau }_2}^t{(v(s)-v^{*})}^{2}ds.\hfill \end{array}$$

Hence, we have, for $t>T(\epsilon )+2({\tau }_1+{\tau }_2)$,

$$\begin{array}{cc}\hfill \mathsf{\Gamma }(t):=& |{\mathsf{\Gamma }}_1(t)|+|{\mathsf{\Gamma }}_2(t)|+|{\mathsf{\Gamma }}_3(t)|+|{\mathsf{\Gamma }}_4(t)|\hfill \\ \hfill \le & \frac{1}{2}{\beta }_1[x_0(\epsilon )e^{-m_1{\tau }_1}{\mathsf{{\rm Y}}}_2(\epsilon )+M_2(\epsilon )({\mathsf{{\rm Y}}}_1(\epsilon )+{\mathsf{\Lambda }}_1{x}_0(\epsilon )+{\beta }_2{x}_0(\epsilon ))]{\tau }_1{(x(t)-x^{*})}^2\hfill \\ \hfill & +\frac{1}{2}{\beta }_1[x_0(\epsilon ){\mathsf{{\rm Y}}}_2(\epsilon )+e^{m_1{\tau }_1}{M}_2(\epsilon )({\mathsf{{\rm Y}}}_1(\epsilon )+{\mathsf{\Lambda }}_1{x}_0(\epsilon )+{\beta }_2{x}_0(\epsilon ))]{\tau }_1{(y(t)-y^{*})}^2\hfill \\ \hfill & +\frac{1}{2}{\beta }_2[M_3(\epsilon )({\mathsf{{\rm Y}}}_1(\epsilon )+{\mathsf{\Lambda }}_1{x}_0(\epsilon )+{\beta }_2{x}_0(\epsilon ))+x_0(\epsilon )(k+u)]{\tau }_2{(x(t)-x^{*})}^2\hfill \\ \hfill & +\frac{1}{2}{\beta }_2{e}^{m_1{\tau }_1}[M_3(\epsilon )({\mathsf{{\rm Y}}}_1(\epsilon )+{\mathsf{\Lambda }}_1{x}_0(\epsilon )+{\beta }_2{x}_0(\epsilon ))+x_0(\epsilon )(k+u)]{\tau }_2{(y(t)-y^{*})}^2\hfill \\ \hfill & +\frac{1}{2}{\beta }_1^2{x}_0(\epsilon )y^{*}(1+e^{-m_1{\tau }_1}){\int }_{t-{\tau }_1}^t{(x(s-{\tau }_1)-x^{*})}^{2}ds\hfill \\ \hfill & +\frac{1}{2}{\beta }_1{\beta }_2{x}_0(\epsilon )v^{*}(1+e^{-m_1{\tau }_1}){\int }_{t-{\tau }_1}^t{(x(s-{\tau }_2)-x^{*})}^{2}ds\hfill \\ \hfill & +\frac{1}{2}{\mathsf{{\rm Y}}}_1(\epsilon ){\beta }_1{M}_2(\epsilon )(1+e^{m_1{\tau }_1}){\int }_{t-{\tau }_1}^t{(x(s)-x^{*})}^{2}ds\hfill \\ \hfill & +\frac{1}{2}{\beta }_1^2{x}_0^2(\epsilon )(1+e^{-m_1{\tau }_1}){\int }_{t-{\tau }_1}^t{(y(s-{\tau }_1)-y^{*})}^{2}ds\hfill \\ \hfill & +\left[\frac{1}{2}p{\beta }_1{x}_0(\epsilon )(1+e^{m_1{\tau }_1})+\frac{1}{2}{\mathsf{\Lambda }}_1{x}_0(\epsilon ){\beta }_1{M}_2(\epsilon )(1+e^{m_1{\tau }_1})\right]{\int }_{t-{\tau }_1}^t{(y(s)-y^{*})}^{2}ds\hfill \\ \hfill & +\frac{1}{2}{\beta }_1{\beta }_2{x}_0^2(\epsilon )(1+e^{-m_1{\tau }_1}){\int }_{t-{\tau }_1}^t{(v(s-{\tau }_2)-v^{*})}^{2}ds\hfill \\ \hfill & +\frac{1}{2}{\beta }_2{x}_0(\epsilon ){\beta }_1{M}_2(\epsilon )(1+e^{m_1{\tau }_1}){\int }_{t-{\tau }_1}^t{(v(s)-v^{*})}^{2}ds\hfill \\ \hfill & +\frac{1}{2}{\mathsf{{\rm Y}}}_1(\epsilon ){\beta }_2{M}_3(\epsilon )(1+e^{m_1{\tau }_1}){\int }_{t-{\tau }_2}^t{(x(s)-x^{*})}^{2}ds\hfill \\ \hfill & +\left[\frac{1}{2}k{\beta }_2{x}_0(\epsilon )(1+e^{m_1{\tau }_1})+\frac{1}{2}{\mathsf{\Lambda }}_1{x}_0(\epsilon ){\beta }_2{M}_3(\epsilon )(1+e^{m_1{\tau }_1})\right]{\int }_{t-{\tau }_2}^t{(y(s)-y^{*})}^{2}ds\hfill \\ \hfill & +\left[\frac{1}{2}u{\beta }_2{x}_0(\epsilon )(1+e^{m_1{\tau }_1})+\frac{1}{2}{\beta }_2{x}_0(\epsilon ){\beta }_2{M}_3(\epsilon )(1+e^{m_1{\tau }_1})\right]{\int }_{t-{\tau }_2}^t{(v(s)-v^{*})}^{2}ds\hfill \\ \hfill =& {\mathsf{\Psi }}_1(\epsilon ){\tau }_1{(x(t)-x^{*})}^2+{\mathsf{\Psi }}_1(\epsilon )e^{m_1{\tau }_1}{\tau }_1{(y(t)-y^{*})}^2+{\mathsf{\Psi }}_2(\epsilon ){\tau }_2{(x(t)-x^{*})}^2+{\mathsf{\Psi }}_2(\epsilon )e^{m_1{\tau }_1}{\tau }_2{(y(t)-y^{*})}^2\hfill \\ \hfill & +{\mathsf{\Psi }}_3(\epsilon ){\int }_{t-{\tau }_1}^t{(x(s-{\tau }_1)-x^{*})}^{2}ds+{\mathsf{\Psi }}_4(\epsilon ){\int }_{t-{\tau }_1}^t{(x(s-{\tau }_2)-x^{*})}^{2}ds+{\mathsf{\Psi }}_5(\epsilon ){\int }_{t-{\tau }_1}^t{(x(s)-x^{*})}^{2}ds\hfill \\ \hfill & +{\mathsf{\Psi }}_6(\epsilon ){\int }_{t-{\tau }_1}^t{(y(s-{\tau }_1)-y^{*})}^{2}ds+{\mathsf{\Psi }}_7(\epsilon ){\int }_{t-{\tau }_1}^t{(y(s)-y^{*})}^{2}ds\hfill \\ \hfill & +{\mathsf{\Psi }}_8(\epsilon ){\int }_{t-{\tau }_1}^t{(v(s-{\tau }_2)-v^{*})}^{2}ds+{\mathsf{\Psi }}_9(\epsilon ){\int }_{t-{\tau }_1}^t{(v(s)-v^{*})}^{2}ds\hfill \\ \hfill & +{\mathsf{\Psi }}_{10}(\epsilon ){\int }_{t-{\tau }_2}^t{(x(s)-x^{*})}^{2}ds+{\mathsf{\Psi }}_{11}(\epsilon ){\int }_{t-{\tau }_2}^t{(y(s)-y^{*})}^{2}ds+{\mathsf{\Psi }}_{12}(\epsilon ){\int }_{t-{\tau }_2}^t{(v(s)-v^{*})}^{2}ds.\hfill \end{array}$$

For $t>T(\epsilon )+2({\tau }_1+{\tau }_2)$, we define

$$\begin{array}{cc}\hfill U_3=& {\mathsf{\Psi }}_3(\epsilon )\left[{\int }_{t-{\tau }_1}^t{\int }_{\theta }^t{(x(s-{\tau }_1)-x^{*})}^{2}dsd\theta +{\tau }_1{\int }_{t-{\tau }_1}^t{(x(s)-x^{*})}^{2}ds\right]\hfill \\ \hfill & +{\mathsf{\Psi }}_4(\epsilon )\left[{\int }_{t-{\tau }_1}^t{\int }_{\theta }^t{(x(s-{\tau }_2)-x^{*})}^{2}dsd\theta +{\tau }_1{\int }_{t-{\tau }_2}^t{(x(s)-x^{*})}^{2}ds\right]+{\mathsf{\Psi }}_5(\epsilon ){\int }_{t-{\tau }_1}^t{\int }_{\theta }^t{(x(s)-x^{*})}^{2}dsd\theta \hfill \\ \hfill & +{\mathsf{\Psi }}_6(\epsilon )\left[{\int }_{t-{\tau }_1}^t{\int }_{\theta }^t{(y(s-{\tau }_1)-y^{*})}^{2}dsd\theta +{\tau }_1{\int }_{t-{\tau }_1}^t{(y(s)-y^{*})}^{2}ds\right]+{\mathsf{\Psi }}_7(\epsilon ){\int }_{t-{\tau }_1}^t{\int }_{\theta }^t{(y(s)-y^{*})}^{2}dsd\theta \hfill \\ \hfill & +{\mathsf{\Psi }}_8(\epsilon )\left[{\int }_{t-{\tau }_1}^t{\int }_{\theta }^t{(v(s-{\tau }_2)-v^{*})}^{2}dsd\theta +{\tau }_1{\int }_{t-{\tau }_2}^t{(v(s)-v^{*})}^{2}ds\right]+{\mathsf{\Psi }}_9(\epsilon ){\int }_{t-{\tau }_1}^t{\int }_{\theta }^t{(v(s)-v^{*})}^{2}dsd\theta \hfill \\ \hfill & +{\mathsf{\Psi }}_{10}(\epsilon ){\int }_{t-{\tau }_2}^t{\int }_{\theta }^t{(x(s)-x^{*})}^{2}dsd\theta +{\mathsf{\Psi }}_{11}(\epsilon ){\int }_{t-{\tau }_2}^t{\int }_{\theta }^t{(y(s)-y^{*})}^{2}dsd\theta \hfill \\ \hfill & +{\mathsf{\Psi }}_{12}(\epsilon ){\int }_{t-{\tau }_2}^t{\int }_{\theta }^t{(v(s)-v^{*})}^{2}dsd\theta .\hfill \end{array}$$

Computing the derivative of $U_3$, we have, for $t>T(\epsilon )+2({\tau }_1+{\tau }_2)$,

$$\begin{array}{cc}\hfill \frac{d{U}_3}{dt}=& {\mathsf{\Psi }}_3(\epsilon )\left[-{\int }_{t-{\tau }_1}^t{(x(s-{\tau }_1)-x^{*})}^{2}ds+{\tau }_1{(x(t)-x^{*})}^2\right]\hfill \\ \hfill & +{\mathsf{\Psi }}_4(\epsilon )\left[-{\int }_{t-{\tau }_1}^t{(x(s-{\tau }_2)-x^{*})}^{2}ds+{\tau }_1{(x(t)-x^{*})}^2\right]+{\mathsf{\Psi }}_5(\epsilon )\left[-{\int }_{t-{\tau }_1}^t{(x(s)-x^{*})}^{2}ds+{\tau }_1{(x(t)-x^{*})}^2\right]\hfill \\ \hfill & +{\mathsf{\Psi }}_6(\epsilon )\left[-{\int }_{t-{\tau }_1}^t{(y(s-{\tau }_1)-y^{*})}^{2}ds+{\tau }_1{(y(t)-y^{*})}^2\right]+{\mathsf{\Psi }}_7(\epsilon )\left[-{\int }_{t-{\tau }_1}^t{(y(s)-y^{*})}^{2}ds+{\tau }_1{(y(t)-y^{*})}^2\right]\hfill \\ \hfill & +{\mathsf{\Psi }}_8(\epsilon )\left[-{\int }_{t-{\tau }_1}^t{(v(s-{\tau }_2)-v^{*})}^{2}ds+{\tau }_1{(v(t)-v^{*})}^2\right]+{\mathsf{\Psi }}_9(\epsilon )\left[-{\int }_{t-{\tau }_1}^t{(v(s)-v^{*})}^{2}ds+{\tau }_1{(v(t)-v^{*})}^2\right]\hfill \\ \hfill & +{\mathsf{\Psi }}_{10}(\epsilon )\left[-{\int }_{t-{\tau }_2}^t{(x(s)-x^{*})}^{2}ds+{\tau }_2{(x(t)-x^{*})}^2\right]+{\mathsf{\Psi }}_{11}(\epsilon )\left[-{\int }_{t-{\tau }_2}^t{(y(s)-y^{*})}^{2}ds+{\tau }_2{(y(t)-y^{*})}^2\right]\hfill \\ \hfill & +{\mathsf{\Psi }}_{12}(\epsilon )\left[-{\int }_{t-{\tau }_2}^t{(v(s)-v^{*})}^{2}ds+{\tau }_2{(v(t)-v^{*})}^2\right].\hfill \end{array}$$

Hence, we have, for $t>T(\epsilon )+2({\tau }_1+{\tau }_2)$,

$$\begin{array}{cc}\hfill \frac{d{U}_3}{dt}+\mathsf{\Gamma }(t)\le & [({\mathsf{\Psi }}_1(\epsilon )+{\mathsf{\Psi }}_3(\epsilon )+{\mathsf{\Psi }}_4(\epsilon )+{\mathsf{\Psi }}_5(\epsilon )){\tau }_1+({\mathsf{\Psi }}_2(\epsilon )+{\mathsf{\Psi }}_{10}(\epsilon )){\tau }_2]{(x(t)-x^{*})}^2\hfill \\ \hfill +& [({\mathsf{\Psi }}_1(\epsilon )e^{m_1{\tau }_1}+{\mathsf{\Psi }}_6(\epsilon )+{\mathsf{\Psi }}_7(\epsilon )){\tau }_1+({\mathsf{\Psi }}_2(\epsilon )e^{m_1{\tau }_1}+{\mathsf{\Psi }}_{11}(\epsilon )){\tau }_2]{(y(t)-y^{*})}^2\hfill \\ \hfill +& [({\mathsf{\Psi }}_8(\epsilon )+{\mathsf{\Psi }}_9(\epsilon )){\tau }_1+{\mathsf{\Psi }}_{12}(\epsilon ){\tau }_2]{(v(t)-v^{*})}^2.\hfill \end{array}$$

Further, we have, for $t>T(\epsilon )+2({\tau }_1+{\tau }_2)$,

$$\begin{array}{cc}\hfill & \frac{d{U}_2}{dt}+\frac{d{U}_3}{dt}\hfill \\ \hfill \le & -\{d-r+\frac{r}{K}({\nu }_1(\epsilon )+x^{*})+\left(c+\frac{r\alpha }{K}\right)y^{*}\hfill \\ \hfill & -[({\mathsf{\Psi }}_1(\epsilon )+{\mathsf{\Psi }}_3(\epsilon )+{\mathsf{\Psi }}_4(\epsilon )+{\mathsf{\Psi }}_5(\epsilon )){\tau }_1+({\mathsf{\Psi }}_2(\epsilon )+{\mathsf{\Psi }}_{10}(\epsilon )){\tau }_2]\}{(x(t)-x^{*})}^2\hfill \\ \hfill & -\{e^{m_1{\tau }_1}\left[e^{m_1{\tau }_1}p+\left(c+\frac{r\alpha }{K}\right){\nu }_1(\epsilon )\right]\hfill \\ \hfill & -[({\mathsf{\Psi }}_1(\epsilon )e^{m_1{\tau }_1}+{\mathsf{\Psi }}_6(\epsilon )+{\mathsf{\Psi }}_7(\epsilon )){\tau }_1+({\mathsf{\Psi }}_2(\epsilon )e^{m_1{\tau }_1}+{\mathsf{\Psi }}_{11}(\epsilon )){\tau }_2]\}{(y(t)-y^{*})}^2\hfill \\ \hfill & +\left\{e^{m_1{\tau }_1}\left[d-r+\frac{r}{K}(x_0(\epsilon )+x^{*})+\left(c+\frac{r\alpha }{K}\right)y^{*}+p\right]+\left(c+\frac{r\alpha }{K}\right)x_0(\epsilon )\right\}\left|(x(t)-x^{*})(y(t)-y^{*})\right|\hfill \\ \hfill & +[({\mathsf{\Psi }}_8(\epsilon )+{\mathsf{\Psi }}_9(\epsilon )){\tau }_1+{\mathsf{\Psi }}_{12}(\epsilon ){\tau }_2]{(v(t)-v^{*})}^2.\hfill \end{array}$$

(16)

Define

$$U_4=\frac{1}{2}{(v(t)-v^{*})}^2,$$

then we have, for $t>T(\epsilon )+2({\tau }_1+{\tau }_2)$,

$$\frac{d{U}_4}{dt}=(v(t)-v^{*})[k(y(t)-y^{*})+u(v^{*}-v(t))]=k(y(t)-y^{*})(v(t)-v^{*})-u{(v(t)-v^{*})}^2.$$

(17)

Finally, we define

$$U=U_1+{\theta }_1(U_2+U_3)+{\theta }_2{U}_4.$$

From (11), (16) and (17), we have, for $t>T(\epsilon )+2({\tau }_1+{\tau }_2)$,

$$\begin{array}{cc}\hfill \frac{dU}{dt}=& \frac{d{U}_1}{dt}+{\theta }_1\left(\frac{d{U}_2}{dt}+\frac{d{U}_3}{dt}\right)+{\theta }_2\frac{d{U}_4}{dt}\hfill \\ \hfill \le & -A_{11}(\epsilon ){(x(t)-x^{*})}^2-A_{22}(\epsilon ){(y(t)-y^{*})}^2-A_{33}(\epsilon ){(v(t)-v^{*})}^2\hfill \\ \hfill & +2A_{12}(\epsilon )|x(t)-x^{*}\left|\right|y(t)-y^{*}|+2A_{23}(\epsilon )|y(t)-y^{*}\left|\right|v(t)-v^{*}|\hfill \\ \hfill =& -\left(|x(t)-x^{*}|,|y(t)-y^{*}|,|v(t)-v^{*}|\right)J(\epsilon ){\left(|x(t)-x^{*}|,|y(t)-y^{*}|,|v(t)-v^{*}|\right)}^T.\hfill \end{array}$$

(18)

Since $J(\epsilon )$ is positive definite, then using the classic Barbalat’s lemma [42], we have

$$\underset{t\to +\infty }{lim}|x(t)-x^{*}|=\underset{t\to +\infty }{lim}|y(t)-y^{*}|=\underset{t\to +\infty }{lim}|v(t)-v^{*}|=0.$$

Thus, the chronic infection equilibrium $E^{*}$ is globally attractive. □

4. Conclusions

In this paper, we mainly study the uniform persistence and global attractivity of chronic infection equilibrium $E^{*}$ of model (4). For the uniform persistence of model (4), Theorem 1 and (5) give explicit expressions of the ultimate upper and lower bounds of any positive solution of model (4). In fact, the global attractivity of the chronic infection equilibrium $E^{*}$ of model (4) is still a question worthy of further discussion. Using standard analytical methods, it is not difficult to find that when delay ${\tau }_1$ or ${\tau }_2$ changes, Hopf bifurcations can appear near the chronic infection equilibrium $E^{*}$. Under the condition $m_1{\tau }_1=m_2{\tau }_2$, Theorem 2 gives a class of sufficient conditions to determine the global attractivity of the chronic infection equilibrium $E^{*}$. Of course, it is also easy to see that the verification of these sufficient conditions are complex and highly conservative, relying strongly on the construction of the Lyapunov functional $U=U_1+{\theta }_1(U_2+U_3)+{\theta }_2{U}_4$.

To illustrate the feasibility of the application of Theorem 2, we specifically choose the following parameter values: $s=1,d=1,r=0.01,K=1000,c=0.01,p=2,k=1,u=1,\alpha =2,{\beta }_1=1,{\beta }_2=1.06,{\tau }_1=0.015,{\tau }_2=0.01,m_1=0.02,m_2=0.03$. Then, we have $R_0\approx 1.040081$ $>1$, $m_1{\tau }_1=m_2{\tau }_2$, and model (4) has a unique chronic infection equilibrium $E^{*}\approx (0.971165,0.0191695,0.0191695)$. We further choose ${\theta }_1=1$ and ${\theta }_2=0.5$, then we can obtain $A_{11}(0)\approx 1.660181>0$, $A_{12}(0)\approx 1.505625$, $A_{22}(0)\approx 1.637131>0$, $A_{23}(0)=0.25$, and $A_{33}(0)\approx 0.3623425>0$. It is easy to verify that $J(0)$ is positive definite and $d-r+\frac{r{x}^{*}}{K}\approx 0.99001>0$. Therefore, the conditions of Theorem 2 are satisfied, and chronic infection equilibrium $E^{*}$ is globally attractive.