Global Properties of Latent Virus Dynamics Models with Immune Impairment and Two Routes of Infection

Abstract

This paper studies the global stability of viral infection models with CTL immune impairment. We incorporate both productively and latently infected cells. The models integrate two routes of transmission, cell-to-cell and virus-to-cell. In the second model, saturated virus–cell and cell–cell incidence rates are considered. The basic reproduction number is derived and two steady states are calculated. We first establish the nonnegativity and boundedness of the solutions of the system, then we investigate the global stability of the steady states. We utilize the Lyapunov method to prove the global stability of the two steady states. We support our theorems by numerical simulations.

1. Introduction

In the literature, several mathematical models of within-host virus dynamics have been constructed and analyzed [1,2,3,4,5,6,7,8,9]. The cytotoxic T Lymphocyte (CTL) is one of the central components of the immune system against viral infections. CTLs lyse the viral-infected cells which participate in reducing or clearing the viruses from the body. Several mathematical models have been presented which integrate the effect of the CTL immune response on viral dynamics (see e.g., [10,11,12]). Nowak and Bangham [10] have presented a mathematical model to characterize the dynamics of the virus (J) with uninfected cells (G), infected cells (I) and CTLs (K) as:

$$\begin{array}{cc}\hfill \dot{G}\left(t\right)& =\theta -\mu G\left(t\right)-\xi G\left(t\right)J\left(t\right),\hfill \end{array}$$

(1)

$$\begin{array}{cc}\hfill \dot{I}\left(t\right)& =\xi G\left(t\right)J\left(t\right)-\varrho I\left(t\right)-\beta I\left(t\right)K\left(t\right),\hfill \end{array}$$

(2)

$$\begin{array}{cc}\hfill \dot{J}\left(t\right)& =\vartheta I\left(t\right)-cJ\left(t\right),\hfill \end{array}$$

(3)

$$\begin{array}{cc}\hfill \dot{K}\left(t\right)& =\rho I\left(t\right)K\left(t\right)-ϵK\left(t\right).\hfill \end{array}$$

(4)

The uninfected cells are replenished at rate $\theta$, die at rate $\mu G$ and become infected at rate $\xi GJ$, where $\xi$ is the virus–cell incidence rate constant. $\beta IK$ is the killer rate of infected cells by CTL and $\varrho I$ is the death rate of the infected cells, where $\beta$ and $\varrho$ are constants. The CTLs are proliferated and die at rates $\rho IK$ and $ϵK$, respectively, where $\rho$ and $ϵ$ are constants.

Models (1)–(4) assume that the presence of antigen can activate the CTL immune response, however, the CTL immune impairment is negelcted. To model the immune impairment, Regoes et al. [13] have modified models (1)–(4) as:

$$\begin{array}{cc}\hfill \dot{G}\left(t\right)& =\theta -\mu G\left(t\right)-\xi G\left(t\right)J\left(t\right),\hfill \end{array}$$

(5)

$$\begin{array}{cc}\hfill \dot{I}\left(t\right)& =\xi G\left(t\right)J\left(t\right)-\varrho I\left(t\right)-\beta I\left(t\right)K\left(t\right),\hfill \end{array}$$

(6)

$$\begin{array}{cc}\hfill \dot{J}\left(t\right)& =\vartheta I\left(t\right)-cJ\left(t\right),\hfill \end{array}$$

(7)

$$\begin{array}{cc}\hfill \dot{K}\left(t\right)& =\rho I\left(t\right)-ϵK\left(t\right)-hI\left(t\right)K\left(t\right),\hfill \end{array}$$

(8)

where the terms $\rho I$ and $hIK$ represents the proliferation rate and the immune impairment, respectively, and h is a constant. Mathematical models of virus dynamics with impairment of CTL functions have been constructed in seveal papers (see e.g., [13,14,15]). The works presented in [13,14,15] assume that the virus infects the uninfected cells by virus-to-cell transmission.

The uninfected target cells can be infected via two ways of transmissions, namely, the diffusion-limited virus-to-cell transmission and the direct cell-to-cell transfer using virological synapses [16]. The cell-to-cell transmission has been recognized in several works (see e.g., [17,18,19,20]). Recent studies have revealed that over 50% of viral infection is due to cell-to-cell transmission [21] and even with an antiretroviral therapy, the cell-to-cell spread of the virus can still permit ongoing replication [22]. Thus, for some viruses, cell-to-cell transmission seems to be a more powerful means of virus propagation than the virus-to-cell transmission [23,24]. Several mathematical models of virus dynamics with two ways of infection have been developed by many researchers (see [25,26,27,28,29,30]). However, in these papers, the impairment of CTL functions is not included. In a very recent work, Elaiw et al. [31] have studied the dynamic behavior of virus infection with impairment of CTL functions and two routes of infection, but with one class of infected cells, productively infected cells.

In case of human immunodeficiency virus (HIV) infection, current treatment consisting of several antiretroviral drugs can suppress viral replication to a low level but cannot completely eradicate the HIV [29]. An important reason is that HIV provirus can reside in latently infected cells [32,33]. Latently infected cells live long, are not affected by antiretroviral drugs or immune responses, but can be activated to produce HIV by relevant antigens.

The aim of the present paper is to propose and analyze viral infection models which include (i) both productively infected cells and latently infected cells, (ii) both virus-to-cell and cell-to-cell transmissions, and (iii) impairment of CTL functions. We first show that the solutions of the models are nonnegative and bounded, then we derive the basic reproduction number which determines the existence and global stability of the steady states. We utilize the Lyapunov method to prove the global stability of the two steady states. We support our theorems by numerical simulations.

2. The Model

We study the following model:

$$\begin{array}{cc}\hfill \dot{G}\left(t\right)& =\theta -\mu G\left(t\right)-{\xi }_{1}G\left(t\right)J\left(t\right)-{\xi }_{2}G\left(t\right)I\left(t\right),\hfill \end{array}$$

(9)

$$\begin{array}{cc}\hfill \dot{L}\left(t\right)& =(1-\nu )({\xi }_{1}G\left(t\right)J\left(t\right)+{\xi }_{2}G\left(t\right)I\left(t\right))-(b+d)L\left(t\right),\hfill \end{array}$$

(10)

$$\begin{array}{cc}\hfill \dot{I}\left(t\right)& =\nu ({\xi }_{1}G\left(t\right)J\left(t\right)+{\xi }_{2}G\left(t\right)I\left(t\right))-\varrho I\left(t\right)-\beta I\left(t\right)K\left(t\right)+bL\left(t\right),\hfill \end{array}$$

(11)

$$\begin{array}{cc}\hfill \dot{J}\left(t\right)& =\vartheta I\left(t\right)-cJ\left(t\right),\hfill \end{array}$$

(12)

$$\begin{array}{cc}\hfill \dot{K}\left(t\right)& =\rho I\left(t\right)-ϵK\left(t\right)-hI\left(t\right)K\left(t\right),\hfill \end{array}$$

(13)

where, L is the concentration of the latently infected cells. The uninfected cells become infected at rates ${\xi }_{1}GJ$ and ${\xi }_{2}GI$ due to virus-to-cell and cell-to-cell infections, respectively, where ${\xi }_1$ and ${\xi }_2$ are the incidence rate constants. The fractions $1-\nu$ and $\nu$ with $0<\nu \le 1$ are the probabilities that upon infection, an uninfected cell will becomes either latently infected or productively infected, respectively. Parameter b denotes the average number of latently infected cells cells that become productively infected cells, and d denotes the death rate constant of the latently infected cells.

2.1. Nonnegativity and Boundedness

Let us define

$$\mathsf{\Omega }=\left\{\left(G,L,I,J,K\right)\in {\mathbb{R}}_{\ge 0}^5:0\le G,L,I\le N_1,0\le J\le N_2,0\le K\le N_3\right\}.$$

(14)

Lemma 1. 

The compact set Ω is positively invariant for system (9)–(13).

Proof. 

We observe that

$$\begin{array}{cc}\hfill {\left.\dot{G}\right|}_{\left(G=0\right)}& =\theta >0,\hfill \\ \hfill {\left.\dot{L}\right|}_{(L=0)}& =(1-\nu )({\xi }_{1}GJ+{\xi }_{2}GI)\ge 0,\forall G,J,I\ge 0,\hfill \\ \hfill {\left.\dot{I}\right|}_{\left(I=0\right)}& =\nu {\xi }_{1}GJ+bL\ge 0,\forall G,J,L\ge 0,\hfill \\ \hfill {\left.\dot{J}\right|}_{\left(J=0\right)}& =\vartheta I\ge 0,\forall I\ge 0,\hfill \\ \hfill {\left.\dot{K}\right|}_{\left(K=0\right)}& =\rho I\ge 0,\forall I\ge 0.\hfill \end{array}$$

This confirms that $\left(G\left(t\right),L\left(t\right),I\left(t\right),J\left(t\right),K\left(t\right)\right)\in {\mathbb{R}}_{\ge 0}^5$ with $\left(G\left(0\right),L\left(0\right),I\left(0\right),J\left(0\right),K\left(0\right)\right)\in {\mathbb{R}}_{\ge 0}^5$. Let $F=G+L+I+\frac{\varrho }{2\vartheta }J+\frac{\varrho }{4\rho }K.$ Then

$$\begin{array}{cc}\hfill \dot{F}& =\theta -\mu G-{\xi }_{1}GJ-{\xi }_{2}GI+(1-\nu )({\xi }_{1}GJ+{\xi }_{2}GI)-(b+d)L+\nu ({\xi }_{1}GJ+{\xi }_{2}GI)-\varrho I\hfill \\ & -\beta IK+bL+\frac{\varrho }{2\vartheta }\left(\vartheta I-cJ\right)+\frac{\varrho }{4\rho }\left(\rho I-ϵK-hIK\right)\hfill \\ & =\theta -\mu G-dL-\frac{\varrho }{4}I-\left(\beta +\frac{\varrho h}{4\rho }\right)IK-\frac{\varrho c}{2\vartheta }J-\frac{\varrho ϵ}{4\rho }K\hfill \\ & \le \theta -\mu G-dL-\frac{\varrho }{4}I-\frac{\varrho c}{2\vartheta }J-\frac{\varrho ϵ}{4\rho }K\hfill \\ & \le \theta -\sigma \left(G+L+I+\frac{\varrho }{2\vartheta }J+\frac{\varrho }{4\rho }K\right)=\theta -\sigma F,\hfill \end{array}$$

where, $\sigma =min\{\mu ,d,\frac{\varrho }{4},c,ϵ\}$. Hence, $0\le F\left(t\right)\le N_1$ for all $t\ge 0$ if $F\left(0\right)\le N_1,$ where $N_1=\frac{\theta }{\sigma }$. Consequently, $0\le G\left(t\right),L\left(t\right),I\left(t\right)\le N_1,0\le J\left(t\right)\le N_2$ and $0\le K\left(t\right)\le N_3$ for all $t\ge 0$ if $G\left(0\right)+L\left(0\right)+I\left(0\right)+\frac{\varrho }{2\vartheta }J\left(0\right)+\frac{\varrho }{4\vartheta }$ $K\left(0\right)\le N_1$, where $N_2=\frac{2\vartheta \theta }{\varrho \sigma }$ and $N_3=\frac{4\rho \theta }{\varrho \sigma }$. This establishes the bondedness of $G\left(t\right),L\left(t\right),I\left(t\right),J\left(t\right)$ and $K\left(t\right)$. □

Let us define the basic reproduction number of system (9)–(13) as:

$${\mathcal{R}}_0=\frac{\theta \left(d\nu +b\right)\left(\vartheta {\xi }_1+c{\xi }_2\right)}{\varrho c\mu (b+d)}.$$

(15)

Lemma2. 

For system (9)–(13),

(i)

if ${\mathcal{R}}_0\le 1$ then there exists a disease-free steady state ${\mathsf{\Delta }}_0$,

(ii)

if ${\mathcal{R}}_0>1$, then there exist two steady states ${\mathsf{\Delta }}_0$ and endemic steady state ${\mathsf{\Delta }}_1$.

Proof. 

The steady states of the system satisfy

$$\begin{array}{cc}\hfill 0& =\theta -\mu G-{\xi }_{1}GJ-{\xi }_{2}GI,\hfill \end{array}$$

(16)

$$\begin{array}{cc}\hfill 0& =(1-\nu )({\xi }_{1}GJ+{\xi }_{2}GI)-(b+d)L,\hfill \end{array}$$

(17)

$$\begin{array}{cc}\hfill 0& =\nu ({\xi }_{1}GJ+{\xi }_{2}GI)-\varrho I-\beta IK+bL,\hfill \end{array}$$

(18)

$$\begin{array}{cc}\hfill 0& =\vartheta I-cJ,\hfill \end{array}$$

(19)

$$\begin{array}{cc}\hfill 0& =\rho I-ϵK-hIK.\hfill \end{array}$$

(20)

By solving Equations (16)–(20) we get two steady states, disease-free steady state ${\mathsf{\Delta }}_0=(G_0,0,0,0,0)$ where $G_0={\displaystyle \frac{\theta }{\mu }}$. In addition, we have

$$A_1{I}^2+B_{1}I+C_1=0,$$

where

$$\begin{array}{cc}\hfill A_1& =(h\varrho +\beta \rho )(b+d)(\vartheta {\xi }_1+c{\xi }_2),\hfill \\ \hfill B_1& =((\vartheta {\xi }_1+c{\xi }_2)ϵ+c\mu h)(b+d)\varrho -(\vartheta {\xi }_1+c{\xi }_2)\theta dh\nu -(\vartheta {\xi }_1+c{\xi }_2)b\theta h+(b+d)\beta \rho c\mu ,\hfill \\ \hfill C_1& =ϵ(bc\mu \varrho +cd\mu \varrho )\left(1-{\mathcal{R}}_0\right).\hfill \end{array}$$

Define a function ${\psi }_1$ by

$${\psi }_1\left(I\right)=A_1{I}^2+B_{1}I+C_1=0.$$

Then, ${\psi }_1\left(0\right)=ϵ(bc\mu \varrho +cd\mu \varrho )\left(1-{\mathcal{R}}_0\right)<0$ when ${\mathcal{R}}_0>1$ and ${\displaystyle \underset{I\to \infty }{lim}{\psi }_1\left(I\right)=\infty }$. Hence, there exists $I_1\in (0,\infty )$ such that ${\psi }_1\left(I_1\right)=0$. Hence, when ${\mathcal{R}}_0>1$, then

$$\begin{array}{cc}\hfill G_1& =\frac{\theta c}{{\xi }_1\vartheta I_1+{\xi }_{2}c{I}_1+c\mu }>0,J_1={\displaystyle \frac{\vartheta I_1}{c}}>0,\hfill \\ \hfill L_1& =\frac{(1-\nu )\theta I_1({\xi }_1\vartheta +{\xi }_{2}c)}{({\xi }_1\vartheta I_1+{\xi }_{2}c{I}_1+c\mu )(d+b)}>0,K_1={\displaystyle \frac{\rho I_1}{h{I}_1+ϵ}}>0.\hfill \end{array}$$

It follows that, an endemic steady state ${\mathsf{\Delta }}_1(G_1,L_{1,}{I}_1,J_1,K_1),$ exists if ${\mathcal{R}}_0>1$. □

2.2. Global Stability

We define $\mathsf{\Gamma }\left(\ell \right)=\ell -1-ln\ell$. We note that $\mathsf{\Gamma }\left(\ell \right)\ge 0$ for any $\ell >0$ and $\mathsf{\Gamma }\left(1\right)=0$. To investigate the global stability of the steady states, we construct Lyapunov functions using the method presented [4] and followed by [5,6,7].

Theorem 1. 

Let ${\mathcal{R}}_0<1$, then ${\mathsf{\Delta }}_0$ of models (9)–(13), is globally asymptotically stable and it is unstable if ${\mathcal{R}}_0>1$.

Proof. 

Constructing a function ${\mathsf{\Lambda }}_0(G,L,I,J,K)$ as:

$${\mathsf{\Lambda }}_0(G,L,I,J,K)=G_0\mathsf{\Gamma }\left(\frac{G}{G_0}\right)+\left(\frac{b}{\nu d+b}\right)L+\left(\frac{b+d}{\nu d+b}\right)I+\frac{{\xi }_1{G}_0}{c}J+\frac{\varrho (1-{\mathcal{R}}_0)}{\rho }\left(\frac{b+d}{\nu d+b}\right)K.$$

Clearly, ${\mathsf{\Lambda }}_0(G,L,I,J,K)$ for all $G,L,I,J,K>0$, while ${\mathsf{\Lambda }}_0(G,L,I,J,K)$ reaches its global minimum at ${\mathsf{\Delta }}_0$. Calculating $\frac{d{\mathsf{\Lambda }}_0}{dt}$ along the trajectories of (9)–(13) we get

$$\begin{array}{cc}\hfill {\displaystyle \frac{d{\mathsf{\Lambda }}_0}{dt}}& =\left(1-\frac{G_0}{G}\right)\left(\theta -\mu G-{\xi }_{1}GJ-{\xi }_{2}GI\right)+\left(\frac{b}{\nu d+b}\right)\left((1-\nu )({\xi }_{1}GJ+{\xi }_{2}GI)-(d+b)L\right)\hfill \\ & +\left(\frac{b+d}{\nu d+b}\right)(\nu ({\xi }_{1}GJ+{\xi }_{2}GI)-\varrho I+bL-\beta IK)+\frac{{\xi }_1{G}_0}{c}(\vartheta I-cJ)+\frac{\varrho (1-{\mathcal{R}}_0)}{\rho }\left(\frac{b+d}{\nu d+b}\right)(\rho I-ϵK-hIK)\hfill \\ & =\left(1-\frac{G_0}{G}\right)(\theta -\mu G)+\varrho \left(\frac{b+d}{\nu d+b}\right)\left(\frac{{\xi }_2{G}_0(\nu d+b)}{\varrho (b+d)}-1+\frac{{\xi }_1{G}_0\vartheta (\nu d+b)}{\varrho c(b+d)}+(1-{\mathcal{R}}_0)\right)I\hfill \\ & -\left(\frac{b+d}{\nu d+b}\right)\left(\beta +\frac{\varrho h(1-{\mathcal{R}}_0)}{\rho }\right)IK-\frac{\varrho ϵ(1-{\mathcal{R}}_0)}{\rho }\left(\frac{b+d}{\nu d+b}\right)K\hfill \\ & =-\mu \frac{{(G-G_0)}^2}{G}-\left(\frac{b+d}{\nu d+b}\right)\left(\beta +\frac{\varrho h(1-{\mathcal{R}}_0)}{\rho }\right)IK-\left(\frac{b+d}{\nu d+b}\right)\frac{\varrho ϵ(1-{\mathcal{R}}_0)}{\rho }K.\hfill \end{array}$$

Since ${\mathcal{R}}_0<1,$ then for all $G,L,I,J,K>0$ we have ${\displaystyle \frac{d{\mathsf{\Lambda }}_0}{dt}}\le 0$. The solutions of the system tend to the largest invariant subset of $\{(G,L,I,J,K):\frac{d{\mathsf{\Lambda }}_0}{dt}=0\}$ [34]. It can be easily show that $\frac{d{\mathsf{\Lambda }}_0}{dt}=0$ at ${\mathsf{\Delta }}_0$. Applying LaSalle’s invariance principle (LIP), we get that ${\mathsf{\Delta }}_0$ is globally asymptotically stable.

We calculate the characteristic equation at the steady state ${\mathsf{\Delta }}_0$ as:

$$(\lambda +\mu )(\lambda +ϵ)({\lambda }^3+a_1{\lambda }^2+a_2\lambda +a_3)=0,$$

(21)

where

$$\begin{array}{cc}\hfill a_1& =-{\xi }_2\nu G_0+\varrho +b+c+d,\hfill \end{array}$$

(22)

$$\begin{array}{cc}\hfill a_2& =((-{\xi }_2\nu G_0+\varrho +b+d)c-({\xi }_{2}b+\nu ({\xi }_1\vartheta +{\xi }_{2}d))G_0+\varrho (d+b),\hfill \end{array}$$

(23)

$$\begin{array}{cc}\hfill a_3& =(-{\xi }_2(\nu d+b)c{G}_0-{\xi }_1\vartheta G_0\left((\nu d+b)\right))+\varrho c(d+b)=\varrho c(d+b)(1-{\mathcal{R}}_0).\hfill \end{array}$$

(24)

Define

$${\psi }_2\left(\lambda \right)={\lambda }^3+a_1{\lambda }^2+a_2\lambda +a_3.$$

We have ${\psi }_2\left(0\right)=\varrho c(d+b)(1-{\mathcal{R}}_0)$. Hence, ${\psi }_2\left(0\right)<0$ when ${\mathcal{R}}_0$$>1$. We have also ${lim}_{\lambda \to \infty }{\psi }_2\left(\lambda \right)=\infty ,$ which shows that ${\psi }_2$ has a positive real root and then, ${\mathsf{\Delta }}_0$ is unstable for ${\mathcal{R}}_0$$>1$. □

Theorem 2. 

For system (9)–(13), if${\mathcal{R}}_0>1$, then Δ${}_1$ isglobally asymptotically stable.

Proof. 

Let a function ${\mathsf{\Lambda }}_1(G,L,I,J,K)$ be defined as:

$$\begin{array}{cc}\hfill {\mathsf{\Lambda }}_1(G,L,I,J,K)& =G_1\mathsf{\Gamma }\left(\frac{G}{G_1}\right)+\left(\frac{b}{\nu d+b}\right)L_1\mathsf{\Gamma }\left(\frac{L}{L_1}\right)+\left(\frac{b+d}{\nu d+b}\right)I_1\mathsf{\Gamma }\left(\frac{I}{I_1}\right)+\frac{{\xi }_1{G}_1}{c}{J}_1\mathsf{\Gamma }\left(\frac{J}{J_1}\right)\hfill \\ & +\frac{\beta }{2(\rho -h{K}_1)}\left(\frac{b+d}{\nu d+b}\right){(K-K_1)}^2.\hfill \end{array}$$

Clearly, ${\mathsf{\Lambda }}_1(G,L,I,J,K)>0$ for all $G,L,I,J,K>0$, and ${\mathsf{\Lambda }}_1(G_1,L_1,I_1,J_1,K_1)=0$. Calculating $\frac{d{\mathsf{\Lambda }}_1}{dt}$ along the trajectories of (9)–(13) we get

$$\begin{array}{cc}\hfill \frac{d{\mathsf{\Lambda }}_1}{dt}& =\left(1-\frac{G_1}{G}\right)\left(\theta -\mu G-{\xi }_{1}GJ-{\xi }_{2}GI\right)+\left(\frac{b}{\nu d+b}\right)\left(1-\frac{L_1}{L}\right)\left((1-\nu )({\xi }_{1}GJ+{\xi }_{2}GI)-(d+b)L\right)\hfill \\ & +\left(\frac{b+d}{\nu d+b}\right)\left(1-\frac{I_1}{I}\right)\left(\nu ({\xi }_{1}GJ+{\xi }_{2}GI)-\varrho I+bL-\beta IK\right)\hfill \\ & +\frac{{\xi }_1{G}_1}{c}\left(1-\frac{J_1}{J}\right)\left(\vartheta I-cJ\right)+\frac{\beta }{\rho -h{K}_1}\left(\frac{b+d}{\nu d+b}\right)(K-K_1)\left(\rho I-ϵK-hIK\right)\hfill \\ & =\left(1-\frac{G_1}{G}\right)(\theta -\mu G)+{\xi }_2{G}_{1}I-\left(\frac{b}{\nu d+b}\right)\frac{L_1}{L}(1-\nu )({\xi }_{1}GJ+{\xi }_{2}GI)+\left(\frac{b(d+b)}{\nu d+b}\right)L_1\hfill \\ & -\nu \left(\frac{b+d}{\nu d+b}\right)\frac{I_1}{I}({\xi }_{1}GJ+{\xi }_{2}GI)-\varrho \left(\frac{b+d}{\nu d+b}\right)(I-I_1)-b\left(\frac{b+d}{\nu d+b}\right)\frac{I_1}{I}L-\beta \left(\frac{b+d}{\nu d+b}\right)(I-I_1)K+\frac{\vartheta {\xi }_1{G}_1}{c}I\hfill \\ & -\frac{\vartheta {\xi }_1{G}_1}{c}\frac{J_1}{J}I+{\xi }_1{G}_1{J}_1+\frac{\beta }{\rho -h{K}_1}\left(\frac{b+d}{\nu d+b}\right)(K-K_1)\left(\rho I-ϵK-hIK\right).\hfill \end{array}$$

(25)

Simplifying Equation (25) and applying the following conditions for ${\mathsf{\Delta }}_1$:

$$\begin{array}{cc}\hfill \theta -\mu G_1& ={\xi }_1{G}_1{J}_1+{\xi }_2{G}_1{I}_1,(1-\nu )\left({\xi }_1{G}_1{J}_1+{\xi }_2{G}_1{I}_1\right)=(d+b)L_1,\hfill \\ \hfill \nu \left({\xi }_1{G}_1{J}_1+{\xi }_2{G}_1{I}_1\right)+b{L}_1& =\varrho I_1+\beta I_1{K}_1,\vartheta I_1=c{J}_1,\rho I_1=ϵK_1+h{I}_1{K}_1,\hfill \\ \hfill \left(\frac{b+d}{\nu d+b}\right)(\varrho I_1+\beta I_1{K}_1)& ={\xi }_1{G}_1{J}_1+{\xi }_2{G}_1{I}_1,\hfill \end{array}$$

we get

$$\begin{array}{cc}\hfill \frac{d{\mathsf{\Lambda }}_1}{dt}& =-\left(\mu +{\xi }_2{I}_1\frac{(b+d)\nu }{\nu d+b}\right)\frac{{\left(G-G_1\right)}^2}{G}-\beta \left(\frac{ϵ+hI}{\rho -h{K}_1}\right)\left(\frac{b+d}{\nu d+b}\right){(K-K_1)}^2\hfill \\ & +{\xi }_1{G}_1{J}_1\left(\frac{b(1-\nu )}{\nu d+b}\right)\left(4-\frac{G_1}{G}-\frac{L_{1}GJ}{L{G}_1{J}_1}-\frac{I_{1}L}{I{L}_1}-\frac{J_{1}I}{J{I}_1}\right)\hfill \\ & +{\xi }_1{G}_1{J}_1\left(\frac{(b+d)\nu }{\nu d+b}\right)\left(3-\frac{G_1}{G}-\frac{I_{1}GJ}{I{G}_1{J}_1}-\frac{J_{1}I}{J{I}_1}\right)\hfill \\ & +{\xi }_2{G}_1{I}_1\left(\frac{b(1-\nu )}{\nu d+b}\right)\left(3-\frac{G_1}{G}-\frac{L_{1}GI}{L{G}_1{I}_1}-\frac{I_{1}L}{I_1{L}_1}\right).\hfill \end{array}$$

(26)

We have if ${\mathcal{R}}_0>1,$ then $G_1,L_{1,}{I}_1,J_1,K_1>0$. The geometrical and arithmetical means relationship implies that

$$\begin{array}{cc}\hfill 4& \le \frac{G_1}{G}+\frac{L_{1}GJ}{L{G}_1{J}_1}+\frac{I_{1}L}{I{L}_1}+\frac{J_{1}I}{J{I}_1},\hfill \\ \hfill 3& \le \frac{G_1}{G}+\frac{I_{1}GJ}{I{G}_1{J}_1}+\frac{J_{1}I}{J{I}_1},\hfill \\ \hfill 3& \le \frac{G_1}{G}+\frac{L_{1}GI}{L{G}_1{I}_1}+\frac{I_{1}L}{I_1{L}_1}.\hfill \end{array}$$

Hence for all $G,L,I,J,K>0$ we have $\frac{d{\mathsf{\Lambda }}_1}{dt}\le 0$ and $\frac{d{\mathsf{\Lambda }}_1}{dt}=0$ when $G=G_1$, $L=L_1,$$I=I_1$, $J=J_1$ and $K=K_1$. Utilizing LIP we obtain that if ${\mathcal{R}}_0>1$, then ${\mathsf{\Delta }}_1$ is globally asymptotically stable. □

3. Model with Saturated Incidence Rate

The rate of infection in model (9)–(13) is bilinear in the virus and the uninfected cell. Actual incidence rates are probably not strictly linear. A less than linear response in viruses and infected cells could occur due to saturation at high virus or infected cell concentrations [35]. Therefore, it is reasonable for us to assume that the infection rate of modeling viral infection is given by saturated mass action. In this section, we study a vial infection model with saturation:

$$\begin{array}{cc}\hfill \dot{G}& =\theta -\mu G-\frac{{\xi }_{1}GJ}{1+{\alpha }_{1}J}-\frac{{\xi }_{2}GI}{1+{\alpha }_{2}I},\hfill \end{array}$$

(27)

$$\begin{array}{cc}\hfill \dot{L}& =(1-\nu )\left(\frac{{\xi }_{1}GJ}{1+{\alpha }_{1}J}+\frac{{\xi }_{2}GI}{1+{\alpha }_{2}I}\right)-(d+b)L,\hfill \end{array}$$

(28)

$$\begin{array}{cc}\hfill \dot{I}& =\nu \left(\frac{{\xi }_{1}GJ}{1+{\alpha }_{1}J}+\frac{{\xi }_{2}GI}{1+{\alpha }_{2}I)}\right)-\varrho I+bL-\beta IK,\hfill \end{array}$$

(29)

$$\begin{array}{cc}\hfill \dot{J}& =\vartheta I-cJ,\hfill \end{array}$$

(30)

$$\begin{array}{cc}\hfill \dot{K}& =\rho I-ϵK-hIK,\hfill \end{array}$$

(31)

where ${\alpha }_1,{\alpha }_2$ are saturation constants. All parameters and variables have the same meaning as (9)–(13).

3.1. Basic Properties

The next lemma shows the nonnegativity and boundedness of the solutions of system (27)–(31)

Lemma 3. 

The compact set Ω is positively invariant for system (27)–(31).

The proof is similar to that of Lemma 1.

The basic reproduction number of system (27)–(31) is the same as given by Equation (15).

Lemma 4. 

Consider models (27)–(31), then

(i) 

A disease-free steady state ${\mathsf{\Delta }}_0$ exists when ${\mathcal{R}}_0\le 1,$

(ii) 

An endemic steady state ${\mathsf{\Delta }}_1$ exists when ${\mathcal{R}}_0>1$.

Proof. 

Let

$$\begin{array}{cc}\hfill 0& =\theta -\mu G-\frac{{\xi }_{1}GJ}{1+{\alpha }_{1}J}-\frac{{\xi }_{2}GI}{1+{\alpha }_{2}I},\hfill \end{array}$$

(32)

$$\begin{array}{cc}\hfill 0& =(1-\nu )\left(\frac{{\xi }_{1}GJ}{1+{\alpha }_{1}J}+\frac{{\xi }_{2}GI}{1+{\alpha }_{2}I}\right)-(d+b)L,\hfill \end{array}$$

(33)

$$\begin{array}{cc}\hfill 0& =\nu \left(\frac{{\xi }_{1}GJ}{1+{\alpha }_{1}J}+\frac{{\xi }_{2}GI}{1+{\alpha }_{2}I}\right)+bL-\varrho I-\beta IK,\hfill \end{array}$$

(34)

$$\begin{array}{cc}\hfill 0& =\vartheta I-cJ,\hfill \end{array}$$

(35)

$$\begin{array}{cc}\hfill 0& =\rho I-ϵK-hIK.\hfill \end{array}$$

(36)

By solving the algebraic Equations (32)–(36) we obtain a disease-free steady state ${\mathsf{\Delta }}_0=(G_0,0,0,0,0)$. Moreover we have

$$A_2{I}^3+B_2{I}^2+C_{2}I+D_2=0,$$

$$\begin{array}{cc}\hfill A_2& =\vartheta (h\varrho +\beta \rho )(\mu {\alpha }_1{\alpha }_2+{\xi }_1{\alpha }_2+{\xi }_2{\alpha }_1)(b+d),\hfill \\ \hfill B_2& =(\vartheta (\mu {\alpha }_1+{\xi }_1)(ϵ{\alpha }_2\varrho +h\varrho +\beta \rho )+(\vartheta ϵ{\alpha }_1\varrho +ch\varrho +c\beta \rho ){\xi }_2+c\mu {\alpha }_2(h\varrho +\beta \rho ))(b+d)\hfill \\ & -({\xi }_1{\alpha }_2+{\xi }_2{\alpha }_1)\theta h\vartheta (d\nu +b),\hfill \\ \hfill C_2& =((\vartheta \mu {\alpha }_1+c\mu {\alpha }_2+\vartheta {\xi }_1+c{\xi }_2)\varrho ϵ+(h\varrho +\rho \beta )c\mu )(b+d)-(({\xi }_1{\alpha }_2+{\xi }_2{\alpha }_1)ϵ\vartheta +(\vartheta {\xi }_1+c{\xi }_2)h)\theta (d\nu +b),\hfill \\ \hfill D_2& =\frac{ϵ}{c\mu \varrho }(1-{\mathcal{R}}_0),\hfill \end{array}$$

where ${\mathcal{R}}_0$ is defined by Equation (15). Define

$${\psi }_3\left(I\right)=A_2{I}^3+B_2{I}^2+C_{2}I+D_2=0.$$

We have

$$\begin{array}{cc}\hfill {\psi }_3\left(0\right)& =\frac{ϵ}{c\mu \varrho }(1-{\mathcal{R}}_0),\hfill \\ \hfill \underset{I\to \infty }{lim}{\psi }_3\left(I\right)& =\infty .\hfill \end{array}$$

Since ${\mathcal{R}}_0>$1, then ${\psi }_3\left(0\right)<0$ and there exists $I_1\in (0,\infty )$ such that ${\psi }_3\left(I_1\right)=0$. Hence

$$\begin{array}{cc}\hfill G_1& =\frac{\theta ({\alpha }_1\vartheta I_1+c)({\alpha }_2{I}_1+1)}{({\alpha }_1{\alpha }_2\mu +{\xi }_1{\alpha }_2+{\xi }_2{\alpha }_1)\vartheta I_1^2+({\alpha }_1\mu \vartheta +{\alpha }_{2}c\mu +{\xi }_1\vartheta +{\xi }_{2}c)I_1+c\mu }>0,\hfill \end{array}$$

(37)

$$\begin{array}{cc}\hfill L_1& =\frac{(1-\nu )(\vartheta ({\xi }_1{\alpha }_2+{\xi }_2{\alpha }_1)I_1+{\xi }_{2}c+{\xi }_1\vartheta ))\theta I_1}{(({\alpha }_1{\alpha }_2\mu +{\xi }_1{\alpha }_2+{\xi }_2{\alpha }_1)\vartheta I_1^2+({\alpha }_1\mu \vartheta +{\alpha }_{2}c\mu +{\xi }_1\vartheta +{\xi }_{2}c)I_1+c\mu )(d+b)}>0,\hfill \end{array}$$

(38)

$$\begin{array}{cc}\hfill J_1& ={\displaystyle \frac{\vartheta I_1}{c}}>0,K_1={\displaystyle \frac{\rho I_1}{h{I}_1+ϵ}}>0.\hfill \end{array}$$

(39)

Hence, the endemic steady state ${\mathsf{\Delta }}_1(G_1,L_1,I_1,J_1,K_1)$ exists when ${\mathcal{R}}_0>1$. □

3.2. Global Properties

Theorem 3. 

Let ${\mathcal{R}}_0<1$, then ${\mathsf{\Delta }}_0$ of system (27)–(31) is globally asymptotically stable and it is unstable if ${\mathcal{R}}_0>1$.

Proof. 

Define ${\mathsf{\Lambda }}_0^G(G,L,I,J,K)$ as the following

$${\mathsf{\Lambda }}_0^G(G,L,I,J,K)=G_0\mathsf{\Gamma }\left(\frac{G}{G_0}\right)+\left(\frac{b}{\nu d+b}\right)L+\left(\frac{b+d}{\nu d+b}\right)I+\frac{{\xi }_1{G}_0}{c}J+\frac{\varrho (1-{\mathcal{R}}_0)}{\rho }\left(\frac{b+d}{\nu d+J}\right)K.$$

It is seen that ${\mathsf{\Lambda }}_0^G(G,L,I,J,K)>0$ for all $G,L,I,J,K>0$ while ${\mathsf{\Lambda }}_0^G(G,L,I,J,K)$ reaches its global minimum at ${\mathsf{\Delta }}_0$. We calculate ${\displaystyle \frac{d{\mathsf{\Lambda }}_0^G}{dt}}$ as:

$$\begin{array}{cc}\hfill {\displaystyle \frac{d{\mathsf{\Lambda }}_0^G}{dt}}& =\left(1-\frac{G_0}{G}\right)\left(\theta -\mu G-\frac{{\xi }_{1}GJ}{1+{\alpha }_{1}J}-\frac{{\xi }_{2}GI}{1+{\alpha }_{2}I}\right)+\left(\frac{b}{\nu d+b}\right)\left((1-\nu )\left(\frac{{\xi }_{1}GJ}{1+{\alpha }_{1}J}+\frac{{\xi }_{2}GI}{1+{\alpha }_{2}I}\right)-(d+b)L\right)\hfill \\ & +\left(\frac{b+d}{\nu d+b}\right)\left(\nu \left(\frac{{\xi }_{1}GJ}{1+{\alpha }_{1}J}+\frac{{\xi }_{2}GI}{1+{\alpha }_{2}I}\right)+bL-\varrho I-\beta IK\right)+\frac{{\xi }_1{G}_0}{c}\left(\vartheta I-cJ\right)\hfill \\ & +\frac{\varrho (1-{\mathcal{R}}_0)}{\rho }\left(\frac{b+d}{\nu d+b}\right)\left(\rho I-ϵK-hIK\right)\hfill \\ & =\mu \left(1-\frac{G_0}{G}\right)\left(G_0-G\right)+\frac{{\xi }_1{G}_{0}J}{1+{\alpha }_{1}J}+\frac{{\xi }_2{G}_{0}I}{1+{\alpha }_{2}I}+\varrho \left(\frac{b+d}{\nu d+b}\right)(1-1-{\mathcal{R}}_0)I+\frac{{\xi }_1{G}_0}{c}\left(\vartheta I-cJ\right)\hfill \\ & -\left(\frac{b+d}{\nu d+b}\right)\left(\left(\beta +\frac{\varrho h(1-{\mathcal{R}}_0)}{\rho }\right)IK+\frac{\varrho (1-{\mathcal{R}}_0)ϵ}{\rho }K\right)\hfill \\ & =-\left(\mu \frac{{(G-G_0)}^2}{G}+\frac{{\alpha }_1{\xi }_1{G}_0{J}^2}{1+{\alpha }_{1}J}+\frac{{\alpha }_2{\xi }_2{G}_0{I}^2}{1+{\alpha }_{2}I}\right)+\varrho \left(\frac{b+d}{\nu d+b}\right)\left(\frac{{\xi }_1{G}_0\vartheta (\nu d+b)}{\varrho c(b+d)}+\frac{{\xi }_2{G}_0(\nu d+b)}{\varrho (b+d)}-{\mathcal{R}}_0\right)I\hfill \\ & -\left(\frac{b+d}{\nu d+b}\right)\left(\left(\beta +\frac{\varrho h(1-{\mathcal{R}}_0)}{\rho }\right)IK+\frac{\varrho (1-{\mathcal{R}}_0)ϵ}{\rho }K\right)\hfill \\ & =-\mu \frac{{(G-G_0)}^2}{G}-\frac{{\alpha }_1{\xi }_1{G}_0{J}^2}{1+{\alpha }_{1}J}-\frac{{\alpha }_2{\xi }_2{G}_0{J}^2}{1+{\alpha }_{2}J}-\left(\frac{b+d}{\nu d+b}\right)\left((\beta +\frac{\varrho h(1-{\mathcal{R}}_0)}{\rho })IK+\frac{\varrho (1-{\mathcal{R}}_0)ϵ}{\rho }K\right).\hfill \end{array}$$

(40)

Clearly if ${\mathcal{R}}_0<1,$ then for all $G,L,I,J,K>0,$ we have ${\displaystyle \frac{d{\mathsf{\Lambda }}_0^G}{dt}}\le 0$, and ${\displaystyle \frac{d{\mathsf{\Lambda }}_0^G}{dt}}=0$ when $G=G_0,L=0,I=0,J=0$ and $K=0$. Applying LIP implies we get that if ${\mathcal{R}}_0<1,$ then ${\mathsf{\Delta }}_0$ is globally asymptotically stable. Similar to the previous section we can easily show that if ${\mathcal{R}}_0>1$, then ${\mathsf{\Delta }}_0$ is unstable. □

Theorem 4. 

Let ${\mathcal{R}}_0>1$ then Δ${}_1$ ofsystem (27)–(31) is globally asymptotically stable.

Proof. 

Define a function ${\mathsf{\Lambda }}_1^G(G,L,I,J,K)$ as:

$$\begin{array}{cc}\hfill {\mathsf{\Lambda }}_1^G(G,L,I,J,K)& =G_1\mathsf{\Gamma }\left(\frac{G}{G_1}\right)+\left(\frac{b}{\nu d+b}\right)L_1\mathsf{\Gamma }\left(\frac{L}{L_1}\right)+\left(\frac{b+d}{\nu d+b}\right)I_1\mathsf{\Gamma }\left(\frac{I}{I_1}\right)+\frac{{\xi }_1{G}_1}{c(1+{\alpha }_1{J}_1)}{J}_1\mathsf{\Gamma }\left(\frac{J}{J_1}\right)\hfill \\ & +\frac{\beta }{2(\rho -h{K}_1)}\left(\frac{b+d}{\nu d+b}\right){(K-K_1)}^2.\hfill \end{array}$$

It is seen that ${\mathsf{\Lambda }}_1^G(G,L,I,J,K)>0$ for all $G,L,I,J,K>0$ while ${\mathsf{\Lambda }}_1^G(G,L,I,J,K)$ reaches its global minimum at ${\mathsf{\Delta }}_1$. Calculating $\frac{d{\mathsf{\Lambda }}_1^G}{dt}$ as:

$$\begin{array}{cc}\hfill \frac{d{\mathsf{\Lambda }}_1^G}{dt}& =\left(1-\frac{G_1}{G}\right)\left(\theta -\mu G-\frac{{\xi }_{1}GJ}{1+{\alpha }_{1}J}-\frac{{\xi }_{2}GI}{1+{\alpha }_{2}I}\right)\hfill \\ & +\left(\frac{b}{\nu d+b}\right)\left(1-\frac{L_1}{L}\right)\left((1-\nu )\left(\frac{{\xi }_{1}GJ}{1+{\alpha }_{1}J}+\frac{{\xi }_{2}GI}{1+{\alpha }_{2}I}\right)-(d+b)L\right)\hfill \\ & +\left(\frac{b+d}{\nu d+b}\right)\left(1-\frac{I_1}{I}\right)\left(\nu \left(\frac{{\xi }_{1}GJ}{1+{\alpha }_{1}J}+\frac{{\xi }_{2}GI}{1+{\alpha }_{2}I}\right)+bL-\varrho I-\beta IK\right)\hfill \\ & +\frac{{\xi }_1{G}_1}{c(1+{\alpha }_1{J}_1)}\left(1-\frac{J_1}{J}\right)\left(\vartheta I-cJ\right)+\frac{\beta }{\rho -h{K}_1}\left(\frac{b+d}{\nu d+b}\right)(K-K_1)\left(\rho I-ϵK-hIK\right)\hfill \\ & =\left(1-\frac{G_1}{G}\right)\left(\theta -\mu G\right)+\frac{{\xi }_1{G}_{1}J}{1+{\alpha }_{1}J}+\frac{{\xi }_2{G}_{1}I}{1+{\alpha }_{2}I}-\left(\frac{(1-\nu )b}{\nu d+b}\right)\left(\frac{{\xi }_{1}GJ}{1+{\alpha }_{1}J}+\frac{{\xi }_{2}GI}{1+{\alpha }_{2}I}\right)\frac{L_1}{L_{}}\hfill \\ & +\frac{b(d+b)}{\nu d+b}{L}_1-\nu \left(\frac{b+d}{\nu d+b}\right)\left(\frac{{\xi }_{1}GJ}{1+{\alpha }_{1}J}+\frac{{\xi }_{2}GI}{1+{\alpha }_{2}I}\right)\frac{I_1}{I}-\varrho \left(\frac{b+d}{\nu d+b}\right)\left(I-I_1\right)\hfill \\ & -\beta \left(\frac{b+d}{\nu d+b}\right)\left(I-I_1\right)K-\left(\frac{b(b+d)}{\nu d+b}\right)\frac{I_1}{I}L+\frac{\vartheta {\xi }_1{G}_1}{c(1+{\alpha }_1{J}_1)}I-\frac{{\xi }_1{G}_{1}J}{1+{\alpha }_1{J}_1}\hfill \\ & -\frac{\vartheta {\xi }_1{G}_1}{c(1+{\alpha }_1{J}_1)}\frac{J_1}{J}I+\frac{{\xi }_1{G}_1{J}_1}{1+{\alpha }_1{J}_1}+\frac{\beta }{\rho -h{K}_1}\left(\frac{b+d}{\nu d+b}\right)(K-K_1)\left(\rho I-ϵK-hIK\right).\hfill \end{array}$$

(41)

The steady state conditions of ${\mathsf{\Delta }}_1$ implies that:

$$\begin{array}{cc}\hfill \theta -\mu G_1& =\frac{{\xi }_1{G}_1{J}_1}{1+{\alpha }_1{J}_1}+\frac{{\xi }_2{G}_1{I}_1}{1+{\alpha }_2{I}_1},(b+d)L_1=(1-\nu )\left(\frac{{\xi }_1{G}_1{J}_1}{1+{\alpha }_1{J}_1}+\frac{{\xi }_2{G}_1{I}_1}{1+{\alpha }_2{I}_1}\right),\hfill \\ \hfill \varrho I_1+\beta I_1{K}_1& =\nu \left(\frac{{\xi }_1{G}_1{J}_1}{1+{\alpha }_1{J}_1}+\frac{{\xi }_2{G}_1{I}_1}{1+{\alpha }_2{I}_1}\right)+b{L}_1,\vartheta I_1=c{J}_1,\rho I_1=ϵK_1+h{I}_1{K}_1,\hfill \\ \hfill \left(\frac{b+d}{\nu d+b}\right)(\varrho I_1+\beta I_1{K}_1)& =\frac{{\xi }_1{G}_1{J}_1}{1+{\alpha }_1{J}_{{}_1}}+\frac{{\xi }_2{G}_1{I}_1}{1+{\alpha }_2{I}_1},\hfill \end{array}$$

we get

$$\begin{array}{cc}\hfill \frac{d{\mathsf{\Lambda }}_1^G}{dt}& =-\mu \frac{{\left(G-G_1\right)}^2}{G}-\frac{{\xi }_1{G}_1{J}_1}{1+{\alpha }_1{J}_1}\left(\frac{{\alpha }_1{(J-J_1)}^2}{J_1(1+{\alpha }_{1}J)(1+{\alpha }_1{J}_1)}\right)-\frac{{\xi }_2{G}_1{I}_1}{1+{\alpha }_2{I}_1}\left(\frac{{\alpha }_2{(I-I_1)}^2}{I_1(1+{\alpha }_{2}I)(1+{\alpha }_2{I}_1)}\right)\hfill \\ & +\frac{{\xi }_1{G}_1{J}_1}{1+{\alpha }_1{J}_1}\left(\frac{b(1-\nu )}{\nu d+b}\right)\left(5-\frac{G_1}{G}-\frac{L_{1}GJ(1+{\alpha }_1{J}_1)}{L{G}_1{J}_1(1+{\alpha }_{1}J)}-\frac{I_{1}L}{I{L}_1}-\frac{J_{1}I}{J{I}_1}-\frac{1+{\alpha }_{1}J}{1+{\alpha }_1{J}_1}\right)\hfill \\ & +\frac{{\xi }_1{G}_1{J}_1}{1+{\alpha }_1{J}_1}\left(\frac{(b+d)\nu }{\nu d+b}\right)\left(4-\frac{G_1}{G}-\frac{I_{1}GJ(1+{\alpha }_1{J}_1)}{I{G}_1{J}_1(1+{\alpha }_{1}J)}-\frac{I{J}_1}{I_{1}J}-\frac{1+{\alpha }_{1}J}{1+{\alpha }_1{J}_1}\right)\hfill \\ & +\frac{{\xi }_2{G}_1{I}_1}{1+{\alpha }_2{I}_1}\left(\frac{b(1-\nu )}{\nu d+b}\right)\left(4-\frac{G_1}{G}-\frac{L_{1}GI(1+{\alpha }_2{I}_1)}{L{G}_1{I}_1(1+{\alpha }_{2}I)}-\frac{I_{1}L}{I{L}_1}-\frac{1+{\alpha }_{2}I}{1+{\alpha }_2{I}_1}\right)\hfill \\ & +\frac{{\xi }_2{G}_1{I}_1}{1+{\alpha }_2{I}_1}\left(\frac{(b+d)\nu }{\nu d+b}\right)\left(3-\frac{G_1}{G}-\frac{G(1+{\alpha }_2{I}_1)}{G_1(1+{\alpha }_{2}I)}-\frac{1+{\alpha }_{2}I}{1+{\alpha }_2{I}_1}\right)-\beta \left(\frac{ϵ+hI}{\rho -h{K}_1}\right)\left(\frac{b+d}{\nu d+b}\right){(K-K_1)}^2.\hfill \end{array}$$

The geometrical and arithmetical means relationship implies that

$$\begin{array}{cc}\hfill 5& \le \frac{G_1}{G}+\frac{L_{1}GJ(1+{\alpha }_1{J}_1)}{L{G}_1{J}_1(1+{\alpha }_{1}J)}+\frac{I_{1}L}{I{L}_1}+\frac{J_{1}I}{J{I}_1}+\frac{1+{\alpha }_{1}J}{1+{\alpha }_1{J}_1},\hfill \\ \hfill 4& \le \frac{G_1}{G}+\frac{I_{1}GJ(1+{\alpha }_1{J}_1)}{I{G}_1{J}_1(1+{\alpha }_{1}J)}+\frac{I{J}_1}{I_{1}J}+\frac{1+{\alpha }_{1}J}{1+{\alpha }_1{J}_1},\hfill \\ \hfill 4& \le \frac{G_1}{G}+\frac{L_{1}GI(1+{\alpha }_2{I}_1)}{L{G}_1{I}_1(1+{\alpha }_{2}I)}+\frac{I_{1}L}{I{L}_1}+\frac{1+{\alpha }_{2}I}{1+{\alpha }_2{I}_1},\hfill \\ \hfill 3& \le \frac{G_1}{G}+\frac{G(1+{\alpha }_2{I}_1)}{G_1(1+{\alpha }_{2}I)}+\frac{1+{\alpha }_{2}I}{1+{\alpha }_2{I}_1}.\hfill \end{array}$$

Thus $\frac{d{\mathsf{\Lambda }}_1^G}{dt}\le 0$ for all $G,L,I,J,K>0$ and $\frac{d{\mathsf{\Lambda }}_1^G}{dt}=0$ at ${\mathsf{\Delta }}_1$. Using LIP one can easily show that ${\mathsf{\Delta }}_1$ is globally asymptotically stable. □

4. Numerical Simulations

In this section, we solve system (27)–(31) numerically with values of the parameters given as: $\theta =270$, $\mu =0.2$, ${\xi }_2=0.005$, $b=0.1$, $d=0.2$, $\varrho =\vartheta =5.5$, $c=3$, $\rho =0.5$, $ϵ=0.1$ and $\nu =0.5$. The parameters ${\xi }_1$, ${\alpha }_1,{\alpha }_2$, $\beta$ and h will be varied. We take $\alpha ={\alpha }_1={\alpha }_2$ and choose different initial conditions as:

IC1: $G\left(0\right)=900,L\left(0\right)=200,I\left(0\right)=15,J\left(0\right)=30,K\left(0\right)=4$,

IC2: $G\left(0\right)=600,L\left(0\right)=150,I\left(0\right)=10,J\left(0\right)=20,K\left(0\right)=3$,

IC3: $G\left(0\right)=400,L\left(0\right)=75,I\left(0\right)=5,J\left(0\right)=10,K\left(0\right)=2,$

IC4: $G\left(0\right)=900,L\left(0\right)=200,I\left(0\right)=140,J\left(0\right)=100,K\left(0\right)=4.2,$

IC5: $G\left(0\right)=900,L\left(0\right)=140,I\left(0\right)=15,J\left(0\right)=100,K\left(0\right)=4.$

Case(1) Stability of steady states:

We take $\alpha =0,h=0.1,\beta =0.04$ and ${\xi }_1$ is varied as:

(i)${\xi }_1=0.0005,$ then ${\mathcal{R}}_0=0.9682<1$. Figure 1 shows that, the solution of the system with different initial conditions IC1–IC3 tends to ${\mathsf{\Delta }}_0$. This result implies that ${\mathsf{\Delta }}_0$ is globally asymptotically stable which confirms Theorem 3.

(ii)${\xi }_1=0.005$ then, ${\mathcal{R}}_0=2.3182>1$. The numerical results show that the solutions of the system tends to the steady state ${\mathsf{\Delta }}_1=(602.3861,249.2046,17.5212,32.1223,4.7300)$ for all IC1–IC3. This supports the global stability result of Theorem 4.

Case(2) Virus dynamics with variation of$\alpha$:

In this case, we fix ${\xi }_1=0.005,h=0.1,\beta =0.4$ and $\alpha$ is changed. We solve the system numerically with the initial condition IC4. In Figure 2, we show the effect of saturated incidence parameter $\alpha$. We can see that the concentration of the uninfected cells is increased as $\alpha$ is increased. Moreover, the concentration of latently infected cells, productively infected, viruses and CTLs are decreased as $\alpha$ is increased.

Case(3) Effect ofhon the virus dynamics:

Here, we fix ${\xi }_1=0.005,\alpha =0.05,\beta =0.4$ and h is changed. The system is solved with initial condition IC5, Figure 3 shows that the increasing of h will increase both $G\left(t\right)$ and $K\left(t\right)$ and decrease all of $L\left(t\right)$, $I\left(t\right)$ and $J\left(t\right)$.

5. Discussion and Conclusions

In this paper, we have proposed two virus dynamics models with impairment of CTL functions. We consider that the healthy cells are infected by two ways, viral and cellular infections. We have considered both latently and productively infected cells. The incidence rate is represented by bilinear and saturation in the first and second models, respectively. We have established the well-posedness of the model. We have derived the basic reproduction numbers ${\mathcal{R}}_0$ which determine the existence and stability of the disease-free steady state ${\mathsf{\Delta }}_0$ and endemic steady state ${\mathsf{\Delta }}_1$ of the model. We have investigated the global stability of the steady states of the model by using the Lyapunov method and LaSalle’s invariance principle. We have proven that (i) if ${\mathcal{R}}_0<1$, then ${\mathsf{\Delta }}_0$ is globally asymptotically stable and the viruses is cleared (ii) if ${\mathcal{R}}_0>1,$ then ${\mathsf{\Delta }}_1$ exists then it is globally asymptotically stable. This case corresponds to the persistence of the viruses. The effects of saturation and CTL impairment have been studied. We have supported the theoretical results by numerical simulations.

Models (1)–(4) have three steady states; disease-free steady state ${\mathsf{\Delta }}_0^C,$ endemic steady state without a CTL immune response ${\mathsf{\Delta }}_1^C$ and endemic steady state with a CTL immune response ${\mathsf{\Delta }}_2^C$. Moreover, the existence and stability of the steady states are determined by two threshold parameters, the basic reproduction number $R_0^C$ (which determines whether or not the disease will progress) and the CTL immune response activation number $R_1^C$ (which determines whether or not a persistent CTL immune response can be established), where

$${\mathcal{R}}_0^C=\frac{\theta \vartheta \xi }{\varrho c\mu },{\mathcal{R}}_1^C=\frac{{\mathcal{R}}_0^C}{1+\frac{ϵ\vartheta \xi }{c\mu \rho }}.$$

In contrast, models (5)–(8) as well as our proposed models (9)–(13) and (27)–(31) have two steady states (${\mathsf{\Delta }}_0$ and ${\mathsf{\Delta }}_1$) and their existence and stability are determined by only the basic reproduction number $R_0$.

It has been reported in several works (see e.g., [10,13,36]) that viruses mutate fast and there is a generation of quasi species that may vary in infectivity. In fact, mutations are one of the ways of immune evasion whereby viruses can evade CTL activity. The high mutation rate of viruses naturally leads to the study of the interplay between immune response and virus diversity for a number of different strains [36]. A viral infection model with CTL immune response and mutations has been proposed in [10] as:

$$\begin{array}{cc}\hfill \dot{G}\left(t\right)& =\theta -\mu G\left(t\right)-\sum _{i=1}^n{\xi }_{i}G\left(t\right)J_i\left(t\right),\hfill \end{array}$$

(42)

$$\begin{array}{cc}\hfill {\dot{I}}_i\left(t\right)& =\sum _{i=1}^n{\xi }_{i}G\left(t\right)J_i\left(t\right)-{\varrho }_i{I}_i\left(t\right)-{\beta }_i{I}_i\left(t\right)K_i\left(t\right),i=1,2,...,n\hfill \end{array}$$

(43)

$$\begin{array}{cc}\hfill {\dot{J}}_i\left(t\right)& ={\vartheta }_i{I}_i\left(t\right)-c_i{J}_i\left(t\right),i=1,2,...,n\hfill \end{array}$$

(44)

$$\begin{array}{cc}\hfill {\dot{K}}_i\left(t\right)& ={\rho }_i{I}_i\left(t\right)K_i\left(t\right)-{ϵ}_i{K}_i\left(t\right),i=1,2,...,n\hfill \end{array}$$

(45)

where, $I_i$ is the concentration of actively infected cells with virus strain i, $J_i$ denotes the concentration of different strains of virus particles and $K_i$ denotes the concentration of strain specific immune responses. It has been assumed that there are n diffierent strains of virus. Models (42)–(45) can be extended to take into account (i) cell-to-cell transmision, (ii) latently infected cells, (iii) immune impairment, and (iv) time delay as:

$$\begin{array}{cc}\hfill \dot{G}\left(t\right)& =\theta -\mu G\left(t\right)-\sum _{i=1}^{n}G\left(t\right)\left[{\xi }_{1,i}{J}_i\left(t\right)+{\xi }_{2,i}{I}_i\left(t\right)\right],\hfill \end{array}$$

(46)

$$\begin{array}{cc}\hfill {\dot{L}}_i\left(t\right)& =(1-{\nu }_i)\sum _{i=1}^n{e}^{-{\gamma }_i{\tau }_i}G(t-{\tau }_i)\left[{\xi }_{1,i}{J}_i(t-{\tau }_i)+{\xi }_{2,i}{I}_i(t-{\tau }_i)\right]-(b_i+d_i)L_i\left(t\right),i=1,2,...,n\hfill \end{array}$$

(47)

$$\begin{array}{cc}\hfill {\dot{I}}_i\left(t\right)& ={\nu }_i\sum _{i=1}^n{e}^{-{\kappa }_i{\omega }_i}G(t-{\omega }_i)\left[{\xi }_{1,i}{J}_i(t-{\omega }_i)+{\xi }_{2,i}{I}_i(t-{\omega }_i)\right]-{\rho }_i{I}_i\left(t\right)-{\beta }_i{I}_i\left(t\right)K_i\left(t\right)+b_i{L}_i\left(t\right),i=1,2,...,n\hfill \end{array}$$

(48)

$$\begin{array}{cc}\hfill {\dot{J}}_i\left(t\right)& ={\theta }_i{e}^{-{\varphi }_i{\kappa }_i}{I}_i(t-{\kappa }_i)-c_i{J}_i\left(t\right),i=1,2,...,n\hfill \end{array}$$

(49)

$$\begin{array}{cc}\hfill {\dot{K}}_i\left(t\right)& ={\rho }_i{I}_i\left(t\right)-{ϵ}_i{K}_i\left(t\right)-h_i{I}_i\left(t\right)K_i\left(t\right),i=1,2,...,n\hfill \end{array}$$

(50)

where $L_i$ is the concentration of latently infected cells with virus strain i. Here, ${\tau }_i$ is the time between a virus strain i entering an uninfected cell to become latently infected cell with virus strain i, and ${\omega }_i$ is the time between a virus strain i entering an uninfected cell and the production of immature viruses of type i. The immature viruses of type i need time ${\kappa }_i$ to be mature. The factors $e^{-{\gamma }_i{\tau }_i}$, $e^{-{\kappa }_i{\omega }_i}$ and $e^{-{\varphi }_i{\kappa }_i}$ represent the probability of surviving to the age of ${\tau }_i$, ${\omega }_i$ and ${\kappa }_i$, respectively, where ${\gamma }_i$, ${\kappa }_i$ and, ${\varphi }_i$ are positive constants. It is worth stressing that the role of the delay term does not only take into account the delay in the dynamical response of the interacting entities, but also their heterogeneity. This can be accounted for by modeling interactions as shown in [37].

Effects of Latent Infection on the Virus Dynamics

In this subsection, we show the effect of the presence of latently infected cells on virus dynamics. Let us incorporate an antiviral drug with efficacy $\eta$ where $\eta \in [0,1)$. The virus dynamics model (9)–(13) under the effect of treatment is given by:

$$\begin{array}{cc}\hfill \dot{G}\left(t\right)& =\theta -\mu G\left(t\right)-(1-\eta )\left[{\xi }_{1}J\left(t\right)+{\xi }_{2}I\left(t\right)\right]G\left(t\right),\hfill \end{array}$$

(51)

$$\begin{array}{cc}\hfill \dot{L}\left(t\right)& =(1-\nu )(1-\eta )\left[{\xi }_{1}J\left(t\right)+{\xi }_{2}I\left(t\right)\right]G\left(t\right)-(b+d)L\left(t\right),\hfill \end{array}$$

(52)

$$\begin{array}{cc}\hfill \dot{I}\left(t\right)& =\nu (1-\eta )\left[{\xi }_{1}J\left(t\right)+{\xi }_{2}I\left(t\right)\right]G\left(t\right)-\varrho I\left(t\right)-\beta I\left(t\right)K\left(t\right)+bL\left(t\right),\hfill \end{array}$$

(53)

$$\begin{array}{cc}\hfill \dot{J}\left(t\right)& =\vartheta I\left(t\right)-cJ\left(t\right),\hfill \end{array}$$

(54)

$$\begin{array}{cc}\hfill \dot{K}\left(t\right)& =\rho I\left(t\right)-ϵK\left(t\right)-hI\left(t\right)K\left(t\right).\hfill \end{array}$$

(55)

The basic reproduction number $R_0^L$ for system (51)–(55) is given by

$${\mathcal{R}}_0^L\left(\eta \right)=(1-\eta )\frac{\theta \left(d\nu +b\right)\left(\vartheta {\xi }_1+c{\xi }_2\right)}{\varrho c\mu (b+d)}.$$

When the population of the latently infected cells are not modeled then models (51)–(55) will become:

$$\begin{array}{cc}\hfill {\mathcal{R}}_0^L\left(\eta \right)& <1,\mathrm{for}\mathrm{all}{\eta }_{crit}^L<\eta <1,\hfill \\ \hfill {\mathcal{R}}_0^W\left(\eta \right)& <1,\mathrm{for}\mathrm{all}{\eta }_{crit}^W<\eta <1,\hfill \end{array}$$

and stabilize the disease-free steady state for systems (51)–(55) and (56)–(59). Now, we calculate ${\eta }_{crit}^W$ and ${\eta }_{crit}^L$ as:

$$\begin{array}{cc}\hfill \dot{G}\left(t\right)& =\theta -\mu G\left(t\right)-(1-\eta )\left[{\xi }_{1}J\left(t\right)+{\xi }_{2}I\left(t\right)\right]G\left(t\right),\hfill \end{array}$$

(56)

$$\begin{array}{cc}\hfill \dot{I}\left(t\right)& =(1-\eta )\left[{\xi }_{1}J\left(t\right)+{\xi }_{2}I\left(t\right)\right]G\left(t\right)-\varrho I\left(t\right)-\beta I\left(t\right)K\left(t\right),\hfill \end{array}$$

(57)

$$\begin{array}{cc}\hfill \dot{J}\left(t\right)& =\vartheta I\left(t\right)-cJ\left(t\right),\hfill \end{array}$$

(58)

$$\begin{array}{cc}\hfill \dot{K}\left(t\right)& =\rho I\left(t\right)-ϵK\left(t\right)-hI\left(t\right)K\left(t\right).\hfill \end{array}$$

(59)

The basic reproduction number $R_0^W$ for system (56)–(59) is given by

$${\mathcal{R}}_0^W\left(\eta \right)=(1-\eta )\frac{\theta \left(\vartheta {\xi }_1+c{\xi }_2\right)}{\varrho c\mu }.$$

Since $0<\nu <1$, then

$${\mathcal{R}}_0^L\left(\eta \right)=(1-\eta )\frac{\theta \left(d\nu +b\right)\left(\vartheta {\xi }_1+c{\xi }_2\right)}{\varrho c\mu (b+d)}<(1-\eta )\frac{\theta \left(\vartheta {\xi }_1+c{\xi }_2\right)}{\varrho c\mu }={\mathcal{R}}_0^W\left(\eta \right).$$

Clearly, the presence of latently infected cells deceases the basic reproduction number of the system. Now, we aim to determine the minimum drug efficacy that can clear the viruses from the body. We determine ${\eta }_{crit}^L$ and ${\eta }_{crit}^W$ that make

$$\begin{array}{cc}\hfill {\eta }_{crit}^L& =max\left\{0,\frac{{\mathcal{R}}_0^L\left(0\right)-1}{{\mathcal{R}}_0^L\left(0\right)}\right\},\hfill \\ \hfill {\eta }_{crit}^W& =max\left\{0,\frac{{\mathcal{R}}_0^W\left(0\right)-1}{{\mathcal{R}}_0^W\left(0\right)}\right\}.\hfill \end{array}$$

Clearly, $R_0^L\left(0\right)<R_0^W\left(0\right)$ and thus ${\eta }_{crit}^L<{\eta }_{crit}^W$. Therefore, the drug efficacy necessary to steer the states of the system to the disease-free steady state is actually less for system (51)–(55) than that for system (56)–(59).