Global Stability of Within-Host Virus Dynamics Models with Multitarget Cells

Abstract

In this paper, we study the stability analysis of two within-host virus dynamics models with antibody immune response. We assume that the virus infects n classes of target cells. The second model considers two types of infected cells: (i) latently infected cells; and (ii) actively infected cells that produce the virus particles. For each model, we derive a biological threshold number ${\mathcal{R}}_0$. Using the method of Lyapunov function, we establish the global stability of the steady states of the models. The theoretical results are confirmed by numerical simulations.

1. Introduction

Recently, many mathematicians have proposed several mathematical models to describe the interaction between viruses (such as HIV, HCV, HBV, HTLV and CHIKV) and human target cells (see, e.g., [1,2,3,4,5,6,7,8,9,10,11,12,13,14,15,16,17,18,19,20,21,22,23,24,25,26,27,28,29,30,31,32,33,34,35,36]). Mathematical models of human viruses can lead to the development of efficient antiviral drugs and to understand the interaction of the viruses with target cells [2]. Studying the stability analysis of the models is also important to understand the behavior of the virus. Immune response plays an important role in controlling the infection of several viruses. Cytotoxic T Lymphocyte (CTL) and antibodies are the two effector responses of the immune system. CTL cells attack and kill the infected cells. The B cell produces antibodies to neutralize the viruses. The antibody immune response is more effective than CTL immune response in some infection processes [37]. The basic virus dynamics model with antibody immune response has been presented in [29,30] as:

$$\begin{array}{cc}\hfill \dot{S}& =\mu -dS-bSV,\hfill \end{array}$$

(1)

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

(2)

$$\begin{array}{cc}\hfill \dot{V}& =mI-rV-qBV,\hfill \end{array}$$

(3)

$$\begin{array}{cc}\hfill \dot{B}& =cBV-\delta B,\hfill \end{array}$$

(4)

where S, I, V, and B are the concentrations of uninfected target cells, infected cells, virus particles and B cells, respectively. Parameters d and $\mu$ represent the death rate and birth rate constants of the uninfected cells, respectively. The uninfected cells become infected at rate $bSV$, where b is rate constant of the virus-target incidence. The infected cells and free virus particles die at rates $ϵI$ and $rV$, respectively. An actively infected cell produces an average number m of virus particles per unit time. The virus particles are attacked by the B cells at rate $qVB$. The term $cBV$ represents the growth rate of B cells after encountering the virus. The B cells die at rate $\delta B$. In 2017, Wang and Liu [36] presented a mathematical model for in host virus infection by considering a constant production rate of the B cells in addition to their proliferation rate. Equation (4) has been modified as:

$$\dot{B}=\eta +cBV-\delta B,$$

(5)

where $\eta$ is the production rate of the B cells.

The model in Equations (1)–(4) assumes that the virus infects only one category of target cells. There are several models for viral infections that have included two categories of target cells (see, e.g., [38,39,40,41,42,43,44,45,46,47]). It has been reported in [48] that HIV infects vital cells in the human immune system such as helper T cells (specifically CD4${}^{+}$ T cells), macrophages, and dendritic cells. In the case of CHIKV infection, the CHIKV inoculates into the body via bites from infected mosquitoes and replicates in a variety of cells, such as skeletal muscle satellite cells, fibroblasts, macrophages, monocytes and other skin cells [49,50,51,52]. To model the virus dynamics with multiplecategories of target cells, Elaiw [53] proposed the following viral infection model:

$$\begin{array}{cc}\hfill {\dot{S}}_i& ={\mu }_i-d_i{S}_i-b_i{S}_{i}V,i=1,\dots ,n,\hfill \end{array}$$

(6)

$$\begin{array}{cc}\hfill {\dot{I}}_i& =b_i{S}_{i}V-{ϵ}_i{I}_i,i=1,\dots ,n,\hfill \end{array}$$

(7)

$$\begin{array}{rr}\dot{V}& =\sum _{i=1}^n{m}_i{I}_i-rV,\end{array}$$

(8)

where n is the number categories of target cells. This model was generalized by Xia et al. [54] by considering general nonlinear rates for viral infection and cell death. In [53,54], the antibody immune response has been neglected. Several mathematical models have been proposed which take the antibody immune response into account (see, e.g., [29,30,31,32,33,34,35]). However, these models have included one target cell population. Therefore, our aim in this paper is to introduce two virus infection models to describe the dynamics of the virus with n classes of target cells. The antibody immune response is considered where the population dynamics of the B cells is described by Equation (5). In the second model, we incorporate both latently and actively infected cells. We investigate the nonnegativity and boundedness of the solutions of the models. We establish the existence of the steady states and analyze their global stability. We construct Lyapunov function using the method of Korobeinikov [55].

2. Virus Dynamics Model

We consider a within-host virus dynamics model with n classes of uninfected target cells.

$$\begin{array}{cc}\hfill {\dot{S}}_i& ={\mu }_i-d_i{S}_i-b_i{S}_{i}V,i=1,\dots ,n,\hfill \end{array}$$

(9)

$$\begin{array}{cc}\hfill {\dot{I}}_i& =b_i{S}_{i}V-{ϵ}_i{I}_i,i=1,\dots ,n,\hfill \end{array}$$

(10)

$$\begin{array}{cc}\hfill \dot{V}& =\stackrel{n}{\sum _{i=1}}{m}_i{I}_i-rV-qBV,\hfill \end{array}$$

(11)

$$\begin{array}{cc}\hfill \dot{B}& =\eta +cBV-\delta B,\hfill \end{array}$$

(12)

where $S_i,I_i$ represent the concentrations of the uninfected target cells and infected cells of class $i,$ respectively.

2.1. Properties of Solutions

To show that the model in Equations (9)–(12) is biologically acceptable in the sense that no population goes negative or infinity, we establish the nonnegativity and boundedness of solutions of the model. Let ${S}=(S_1,S_2,\dots ,S_n)$ and ${I}=(I_1,I_2,\dots ,I_n).$

Proposition 1.

For the system in Equations (9)–(12), there exists a positively invariant compact set

$$\mathsf{\Omega }=\left\{\left({S},{I},V,B\right)\in {\mathbb{R}}_{\ge 0}^{2n+2}:0\le S_i,I_i\le N_i,i=1,\dots ,n.0\le V\le M_1,0\le B\le M_2\right\}.$$

Proof. 

Since

$$\begin{array}{cc}\hfill {\left.{\dot{S}}_i\right|}_{S_i=0}& ={\mu }_i>0,i=1,\dots ,n,{\left.{\dot{I}}_i\right|}_{I_i=0}=b_i{S}_{i}V\ge 0\mathrm{for}\mathrm{all}{S}_i,V\ge 0,i=1,\dots ,n,\hfill \\ \hfill {\left.\dot{V}\right|}_{V=0}& =\sum _{i=1}^n{m}_i{I}_i\ge 0\mathrm{for}\mathrm{all}{I}_i\ge 0,{\left.\dot{B}\right|}_{B=0}=\eta >0,\hfill \end{array}$$

then, ${\mathbb{R}}_{\ge 0}^{2n+2}$ is positively invariant for the system in Equations (9)–(12).

Next, let $T_i\left(t\right)=S_i+I_i,$$i=1,\dots ,n,$ then

$${\dot{T}}_i\le {\mu }_i-{\sigma }_i{T}_i,$$

where ${\sigma }_i=min\{d_i,{ϵ}_i\}.$ Hence, $0\le$$T_i\left(t\right)\le \frac{{\mu }_i}{{\sigma }_i}$ for all $t\ge 0$ if $T_i\left(0\right)\le \frac{{\mu }_i}{{\sigma }_i}.$ It follows that $0\le S_i\left(t\right),I_i\left(t\right)\le N_i$ for all $t\ge 0,$ if $0\le S_i\left(0\right)+I_i\left(0\right)\le N_i,$ where, $N_i=\frac{{\mu }_i}{{\sigma }_i}.$ Moreover, let $X\left(t\right)=V\left(t\right)+\frac{q}{c}B\left(t\right),$ then

$$\dot{X}\left(t\right)=\sum _{i=1}^n{m}_i{I}_i-rV+\frac{q}{c}\eta -\frac{\delta q}{c}B\le \sum _{i=1}^n{m}_i{N}_i+\frac{q}{c}\eta -{\rho }_1(V+\frac{q}{c}B)=\sum _{i=1}^n{m}_i{N}_i+\frac{q}{c}\eta -{\rho }_{1}X\left(t\right),$$

where ${\rho }_1=min\{r,\delta \}.$ Hence, $X\left(t\right)\le M_1,$ where $M_1=\sum _{i=1}^n\frac{m_i{N}_i+\frac{q}{c}\eta }{{\rho }_1}.$ Since $V\left(t\right)\ge 0$ and $B\left(t\right)\ge 0,$ then $0\le V\left(t\right)\le M_1$ and $0\le B\left(t\right)\le M_2$ if $0\le V\left(0\right)+\frac{q}{c}B\left(0\right)\le M_1,$ where $M_2=\frac{c{M}_1}{q}.☐$

2.2. Steady States

In this subsection, we show the existence of the steady states of the model in Equations (9)–(12). The basic reproduction number of the system in Equations (9)–(12) is defined as:

$${\mathcal{R}}_0=\stackrel{n}{\sum _{i=1}}{\mathcal{R}}_{0i}=\stackrel{n}{\sum _{i=1}}\frac{m_i{b}_i{\mu }_i\delta }{d_i{ϵ}_i(r\delta +q\eta )}.$$

Lemma 1.

(i) If ${\mathcal{R}}_0\le 1,$ then the virus-free steady state $Q_0$ is the only steady state for the system; and (ii) if ${\mathcal{R}}_0>1,$ then the system has a unique endemic steady state $Q_1$ and $Q_1\in \mathring{\mathsf{\Omega }},$ where $\mathring{\mathsf{\Omega }}$ is the interior of $\mathsf{\Omega }.$

Proof. 

Let the R.H.S of the system in Equations (9)–(12) equal zero; then, we get:

$$\begin{array}{cc}\hfill S_i& =\frac{{\mu }_i}{d_i+b_{i}V},i=1,\dots ,n,\hfill \end{array}$$

(13)

$$\begin{array}{cc}\hfill I_i& =\frac{b_i{S}_{i}V}{{ϵ}_i},i=1,\dots ,n,\hfill \end{array}$$

(14)

$$\begin{array}{cc}\hfill \sum _{i=1}^n{m}_i{I}_i& =\left(r+qB\right)V,\hfill \end{array}$$

(15)

$$\begin{array}{cc}\hfill B& =\frac{\eta }{\delta -cV}.\hfill \end{array}$$

(16)

Substituting Equations (13), (14) and (16) into Equation (15), we get

$$\begin{array}{cc}\hfill \sum _{i=1}^n\frac{m_i{b}_i{\mu }_{i}V}{{ϵ}_i\left(d_i+b_{i}V\right)}-\left(r+qB\right)V& =0,\hfill \\ \hfill \left(\stackrel{n}{\sum _{i=1}}\frac{m_i{b}_i{\mu }_i\delta }{d_i{ϵ}_i\left(1+{\alpha }_{i}V\right)(r\delta +q\eta )}-\frac{\left(r+qB\right)\delta }{(r\delta +q\eta )}\right)V& =0,\hfill \end{array}$$

$$\begin{array}{cc}\hfill \left(\stackrel{n}{\sum _{i=1}}\frac{{\mathcal{R}}_{0i}}{1+{\alpha }_{i}V}-\frac{\delta \left(\gamma -rcV\right)}{\gamma \left(\delta -cV\right)}\right)V& =0,\hfill \end{array}$$

(17)

where ${\alpha }_i=\frac{b_i}{d_i}$ and $\gamma =$$r\delta +q\eta .$ Equation (17) admits $V=0$ as a solution. Substituting $V=0$ in Equations (13),(14) and (16), we get the virus-free steady state $Q_0=({{S}}^0,{{I}}^0,V^0,B^0),$ where $S_i^0=\frac{{\mu }_i}{d_i},$$I_i^0=0,$$V^0=0,$ and $B^0=\frac{\eta }{\delta }.$ The other possibility of Equation (17) is

$$\sum _{i=1}^n\frac{{\mathcal{R}}_{0i}\gamma \left(\delta -cV\right)}{1+{\alpha }_{i}V}-\delta \left(\gamma -rcV\right)=0,$$

(18)

where $V\ne \frac{\delta }{c}$.

Let us define a function $F\left(V\right)$ as:

$$F\left(V\right)=\sum _{i=1}^n\frac{{\mathcal{R}}_{0i}\gamma \left(\delta -cV\right)}{1+{\alpha }_{i}V}-\delta \left(\gamma -rcV\right)$$

Then, we get

$$\begin{array}{cc}\hfill F\left(0\right)& =\gamma \delta \left({\mathcal{R}}_0-1\right),\hfill \\ \hfill F\left(\frac{\delta }{c}\right)& =-q\delta \eta .\hfill \end{array}$$

Thus, if ${\mathcal{R}}_0>1,$ then $F\left(0\right)>0,$$F\left(\frac{\delta }{c}\right)<0$ and there exists $V^{*}\in \left(0,\frac{\delta }{c}\right)$ such that $F\left(V^{*}\right)=0.$ Moreover, from Equations (13), (14) and (16), we obtain $S_i^{*}>0,I_i^{*}>0$ and $B^{*}>0$. Then, $Q_1=({{S}}^{*},{{L}}^{*},{{I}}^{*},V^{*},B^{*})$ exists when ${\mathcal{R}}_0>1$.

Clearly, $Q_0\in \mathsf{\Omega }.$ Now, we show that if ${\mathcal{R}}_0>1,$ then $Q_1\in \mathring{\mathsf{\Omega }}.$ From Equations (13) and (14), we have $d_i{S}_i^{*}+{ϵ}_i{I}_i^{*}={\mu }_i,$$i=1,\dots ,n.$ Since $S_i^{*}>0$ and $I_i^{*}>0,$ then

$$\begin{array}{cc}\hfill d_i{S}_i^{*}& <{\mu }_i\Rightarrow S_i^{*}<\frac{{\mu }_i}{d_i}\le N_i,\hfill \\ \hfill {ϵ}_i{I}_i^{*}& <{\mu }_i\Rightarrow I_i^{*}<\frac{{\mu }_i}{{ϵ}_i}\le N_i.\hfill \end{array}$$

Moreover, from Equations (15)–(16), we have

$$\sum _{i=1}^n{m}_i{I}_i^{*}-r{V}^{*}-q{V}^{*}{B}^{*}+\frac{q}{c}(\eta +c{B}^{*}{V}^{*}-\delta B^{*})=0,$$

$$\begin{array}{c}\Rightarrow r{V}^{*}+\frac{\delta q}{c}{B}^{*}=\sum _{i=1}^n{m}_i{I}_i^{*}+\frac{q}{c}\eta <\sum _{i=1}^n{m}_i{N}_i+\frac{q}{c}\eta \hfill \\ \Rightarrow V^{*}<\sum _{i=1}^n\frac{m_i{N}_i+\frac{q}{c}\eta }{r}\le M_1,B_1<\frac{c}{q}\sum _{i=1}^n\frac{m_i{N}_i+\frac{q}{c}\eta }{\delta }\le \frac{c{M}_1}{q}=M_2.\hfill \end{array}$$

It follows that $Q_1\in \mathring{\mathsf{\Omega }}.☐$

2.3. Global Stability

In the following theorems, we establish the global stability of the two steady states of the system in Equations (9)–(12) by constructing suitable Lyapunov functions. Let us define

$$H\left(x\right)=x-lnx-1.$$

Clearly, $H\left(x\right)\ge 0$ for $x>0$ and $H\left(1\right)=0$.

Theorem 1.

For the system in Equations (9)–(12), suppose that ${\mathcal{R}}_0\le 1$, then $Q_0$ is globally asymptotically stable in $\mathsf{\Omega }.$

Theorem 2.

For the system in Equations (9)–(12), suppose that ${\mathcal{R}}_0>1$, then $Q_1$ is globally asymptotically stable in $\mathring{\mathsf{\Omega }}.$

The proofs of these theorems are given in the Appendix A.

Biologically, when ${\mathcal{R}}_0<1$, then each infected cell will produce less than one infected cell during its life at the beginning of the infection. The virus will be decreased and eliminated from the body. When ${\mathcal{R}}_0>1$, then, at the beginning of the infection, each infected cell will produce more than one infected cell during its life. The viruses will be increased and the infection becomes chronic.

3. Virus Model with Latency

In this section, we study the mathematical model of virus infection with n classes of uninfected target cells, taking into account the latently infected cells (such cells contain the viruses but are not producing it) and the actively infected cells (such cells are producing the viruses).

$$\begin{array}{cc}\hfill {\dot{S}}_i& ={\mu }_i-d_i{S}_i-b_i{S}_{i}V,i=1,\dots ,n,\hfill \end{array}$$

(19)

$$\begin{array}{cc}\hfill {\dot{L}}_i& =(1-p_i)b_i{S}_{i}V-({\theta }_i+{\lambda }_i)L_i,i=1,\dots ,n,\hfill \end{array}$$

(20)

$$\begin{array}{cc}\hfill {\dot{I}}_i& =p_i{b}_i{S}_{i}V+{\lambda }_i{L}_i-{ϵ}_i{I}_i,i=1,\dots ,n,\hfill \end{array}$$

(21)

$$\begin{array}{cc}\hfill \dot{V}& =\stackrel{n}{\sum _{i=1}}{m}_i{I}_i-rV-qBV,\hfill \end{array}$$

(22)

$$\begin{array}{cc}\hfill \dot{B}& =\eta +cBV-\delta B,\hfill \end{array}$$

(23)

where $L_i$ and $I_i$ are the concentrations of latently infected and actively infected target cells of class $i,$ respectively. A fraction $(1-p_i)$ of infected target cells is assumed to be latently infected cells and the remaining $p_i$ becomes actively infected cells, where $0<p_i<1,$ $i=1,\dots ,n.$ The latently infected cells are transmitted to actively infected cells at rate ${\lambda }_i{L}_i$ and die at rate ${\theta }_i{L}_i.$

3.1. Properties of Solutions

In the following, we establish the nonnegativity and boundedness of solutions of the model in Equations (19)–(23). Let ${L}=(L_1,L_2,\dots ,L_n).$

Proposition 2.

For the system in Equations (19)–(23), there exists a positively invariant compact set

$${\mathsf{\Omega }}^L=\left\{\left({S},{L},{I},V,B\right)\in {\mathbb{R}}_{\ge 0}^{3n+2}:0\le S_i,L_i,I_i\le N_i^L,i=1,\dots ,n.0\le V\le M_1^L,0\le B\le M_2^L\right\}$$

Proof. 

Since

$$\begin{array}{cc}\hfill {\left.{\dot{S}}_i\right|}_{S_i=0}& ={\mu }_i>0,i=1,\dots ,n,\hfill \\ \hfill {\left.{\dot{L}}_i\right|}_{L_i=0}& =(1-p_i)b_i{S}_{i}V\ge 0\mathrm{for}\mathrm{all}{S}_i,V\ge 0,i=1,\dots ,n,\hfill \\ \hfill {\left.{\dot{I}}_i\right|}_{I_i=0}& =p_i{b}_i{S}_{i}V+{\lambda }_i{L}_i\ge 0\mathrm{for}\mathrm{all}{S}_i,L_i,V\ge 0,i=1,\dots ,n,\hfill \\ \hfill {\left.\dot{V}\right|}_{V=0}& =\sum _{i=1}^n{m}_i{I}_i\ge 0\mathrm{for}\mathrm{all}{I}_i\ge 0,\hfill \\ \hfill {\left.\dot{B}\right|}_{B=0}& =\eta >0,\hfill \end{array}$$

then, ${\mathbb{R}}_{\ge 0}^{3n+2}$ is positively invariant for the system in Equations (19)–(23).

We consider $T_i\left(t\right)=S_i+L_i+I_i,$$i=1,\dots ,n,$ then

$${\dot{T}}_i\left(t\right)\le {\mu }_i-{\sigma }_i^L{T}_i\left(t\right),$$

where ${\sigma }_i^L=min\{d_i,{\theta }_i,{ϵ}_i\}.$ Hence, $0\le$ $T_i\left(t\right)\le \frac{{\mu }_i}{{\sigma }_i^L}$ for all $t\ge 0$ if $T_i\left(0\right)\le \frac{{\mu }_i}{{\sigma }_i^L}.$ It follows that $0\le S_i\left(t\right),L_i\left(t\right),I_i\left(t\right)\le N_i^L$ for all $t\ge 0,$ if $0\le S_i\left(0\right)+L_i\left(0\right)+I_i\left(0\right)\le N_i^L,$ where $N_i^L=\frac{{\mu }_i}{{\sigma }_i^L}.$ On the other hand, let $X\left(t\right)=V\left(t\right)+\frac{q}{c}B\left(t\right),$ then

$$\dot{X}\left(t\right)=\sum _{i=1}^n{m}_i{I}_i-rV+\frac{q}{c}\eta -\frac{\delta q}{c}B\le \sum _{i=1}^n{m}_i{N}_i+\frac{q}{c}\eta -{\rho }_2(V+\frac{q}{c}B)=\sum _{i=1}^n{m}_i{N}_i+\frac{q}{c}\eta -{\rho }_{1}X\left(t\right).$$

Hence, $X\left(t\right)\le M_1^L,$ where $M_1^L=\sum _{i=1}^n\frac{m_i{N}_i+\frac{q}{c}\eta }{{\rho }_1}.$ Since $V\left(t\right)\ge 0$ and $B\left(t\right)\ge 0,$ then $0\le V\left(t\right)\le M_1^L$ and $0\le B\left(t\right)\le M_2^L$ if $0\le V\left(0\right)+\frac{q}{c}B\left(0\right)\le M_1^L,$ where $M_2^L=\frac{c{M}_1^L}{q}.☐$

3.2. Steady States

In this subsection, we establish the existence of the steady states of the model in Equations (19)–(23). The basic reproduction number of the system in Equations (19)–(23) is defined as:

$${\mathcal{R}}_0^L=\stackrel{n}{\sum _{i=1}}{\mathcal{R}}_{0i}^L=\stackrel{n}{\sum _{i=1}}\frac{m_i{b}_i{\mu }_i\delta ({\theta }_i{p}_i+{\lambda }_i)}{d_i{ϵ}_i({\theta }_i+{\lambda }_i)(r\delta +q\eta )}.$$

Lemma 2.

(i) If ${\mathcal{R}}_0^L\le 1,$ then the virus-free steady state $Q_0^L$ is the only steady state for the system; and (ii) if ${\mathcal{R}}_0^L>1,$ then the system has a unique endemic steady state $Q_1^L$ and $Q_1^L\in {\mathring{\mathsf{\Omega }}}^L.$

Proof. 

Let the R.H.S of the system in Equations (19)–(23) equal zero; then, we get:

$$\begin{array}{cc}\hfill S_i& =\frac{{\mu }_i}{d_i+b_{i}V},i=1,\dots ,n,\hfill \end{array}$$

(24)

$$\begin{array}{cc}\hfill L_i& =\frac{(1-p_i)b_i{S}_{i}V}{{\theta }_i+{\lambda }_i},i=1,\dots ,n,\hfill \end{array}$$

(25)

$$\begin{array}{cc}\hfill I_i& =\frac{b_i{S}_{i}V({\theta }_i{p}_i+{\lambda }_i)}{{ϵ}_i({\theta }_i+{\lambda }_i)},i=1,\dots ,n,\hfill \end{array}$$

(26)

$$\begin{array}{cc}\hfill \sum _{i=1}^n{m}_i{I}_i& =\left(r+qB\right)V,\hfill \end{array}$$

(27)

$$\begin{array}{cc}\hfill B& =\frac{\eta }{\delta -cV}.\hfill \end{array}$$

(28)

Substituting Equations (24)–(26) and (28) into Equation (27), we get

$$\begin{array}{cc}\hfill \sum _{i=1}^n\frac{m_i{b}_i{\mu }_i({\theta }_i{p}_i+{\lambda }_i)V}{{ϵ}_i\left(d_i+b_{i}V\right)({\theta }_i+{\lambda }_i)}-\left(r+qB\right)V& =0,\hfill \\ \hfill \left(\stackrel{n}{\sum _{i=1}}\frac{m_i{b}_i{\mu }_i\delta ({\theta }_i{p}_i+{\lambda }_i)}{d_i{ϵ}_i\left(1+{\alpha }_{i}V\right)({\theta }_i+{\lambda }_i)(r\delta +q\eta )}-\frac{\left(r+qB\right)\delta }{(r\delta +q\eta )}\right)V& =0,\hfill \\ \hfill \left(\stackrel{n}{\sum _{i=1}}\frac{{\mathcal{R}}_{0i}^L}{1+{\alpha }_{i}V}-\frac{\delta \left(\gamma -rcV\right)}{\gamma \left(\delta -cV\right)}\right)V& =0,\hfill \end{array}$$

(29)

where ${\alpha }_i=\frac{b_i}{d_i}$ and $\gamma =r\delta +q\eta .$ Equation (29) gives two possible solutions. If $V=0,$ then we get virus-free steady state $Q_0^L=({{S}}^0,{{L}}^0,{{I}}^0,V^0,B^0)=(\frac{{\mu }_i}{d_i},0,0,0,\frac{\eta }{\delta }).$ If $V\ne 0,$ then Equation (29) becomes

$$\sum _{i=1}^n\frac{{\mathcal{R}}_{0i}^L\gamma \left(\delta -cV\right)}{1+{\alpha }_{i}V}-\delta \left(\gamma -rcV\right)=0.$$

(30)

We define a function $G\left(V\right)$ as:

$$G\left(V\right)=\sum _{i=1}^n\frac{{\mathcal{R}}_{0i}^L\gamma \left(\delta -cV\right)}{1+{\alpha }_{i}V}-\delta \left(\gamma -rcV\right).$$

If ${\mathcal{R}}_0^L>1,$then

$$\begin{array}{cc}\hfill G\left(0\right)& =\gamma \delta \left({\mathcal{R}}_0^L-1\right)>0,\hfill \\ \hfill G\left(\frac{\delta }{c}\right)& =-q\delta \eta <0.\hfill \end{array}$$

and there exists $V^{*}\in \left(0,\frac{\delta }{c}\right)$ such that $G\left(V^{*}\right)=0.$ Moreover, from Equations (24)–(26) and (28), we obtain $S_i^{*}>0,L_i^{*}>0,I_i^{*}>0$ and $B^{*}>0.$

Similar to the proof of Lemma 1, one can easily show that $Q_0^L\in {\mathsf{\Omega }}^L$and $Q_1^L\in {\mathring{\mathsf{\Omega }}}^L.☐$

3.3. Global Stability

In this section, we use Lyapunov method to prove the global stability of the two steady states of the system in Equations (19)–(23).

Theorem 3.

For the system in Equations (19)–(23), suppose that ${\mathcal{R}}_0^L\le 1$, then $Q_0^L$ is globally asymptotically stable in ${\mathsf{\Omega }}^L.$

Theorem 4.

For the system in Equations (19)–(23), suppose that ${\mathcal{R}}_0^L>1$, then $Q_1^L$ is globally asymptotically stable in ${\mathring{\mathsf{\Omega }}}^L.$

The proofs of these theorems are given in the Appendix A.

4. Numerical Simulations

To illustrate our theoretical results, we perform numerical simulations for the systems in Equations (9)–(12) and Equations (19)–(23). We consider the case $n=2$.

4.1. Simulations for Virus Dynamics Model

Using the values of the parameters given in

Table 1. The values of the parameters of the model in Equations (9)–(12).

Parameter	Value	Parameter	Value
${\mu }_1$	$1.826$	$d_1$	$0.7979$
${\mu }_2$	$3.198$	$d_2$	$0.5$
$b_1$	varied	q	$0.5964$
c	$0.5$	r	$0.4418$
$m_1=m_2$	$2.02$	$\eta$	$1.402$
$b_2$	varied	$\delta$	$1.251$
${ϵ}_1={ϵ}_2$	$0.4441$

, we show the dynamical behavior of the system states ${S}$, ${I}$, V and B, to confirm the theoretical results given in Theorems 1–2.

• Effect of $b_1$ and $b_2$ on the stability of steady states: To show the global stability results, we consider three different initial values as:

IV1: $S_1\left(0\right)=2.0,S_2\left(0\right)=6.0,I_1\left(0\right)=0.3,I_2\left(0\right)=0.4,V\left(0\right)=0.4$ and $B\left(0\right)=2.0,$

IV2: $S_1\left(0\right)=1.7,S_2\left(0\right)=5.5,I_1\left(0\right)=0.7,I_2\left(0\right)=0.7,V\left(0\right)=0.6$ and $B\left(0\right)=4.0,$

IV3: $S_1\left(0\right)=1.4,S_2\left(0\right)=5.0,I_1\left(0\right)=1.1,I_2\left(0\right)=1.0,V\left(0\right)=0.8$ and $B\left(0\right)=6.0.$

We consider two sets of the parameters $b_1$ and $b_2$ as follows:

Set (I): We choose $b_1=0.005$ and $b_2=0.02.$ Using these data, we compute ${\mathcal{R}}_0=0.5710<1$, then the system has one steady state $Q_0$. In Figure 1, we can see that the concentrations of the uninfected target cells and B cells return to their normal values $S_1^0=\frac{{\mu }_1}{d_1}=2.2885 ,$$S_2^0=\frac{{\mu }_2}{d_2}=6.3960,$ and $B^0=\frac{\eta }{\delta }=1.1207$. On the other hand, the concentrations of infected target cells and virus particles are decaying and approaching zero for all the three initial values IV1–IV3. It means that $Q_0$ is globally asymptotically stable and the virus will be cleared. This result confirms the result of Theorem 1.

Set (II): We take $b_1=0.5$ and $b_2=0.1.$ Then, we calculate ${\mathcal{R}}_0=7.3086>1.$ The system has two steady states $Q_0$ and $Q_1$. It is clear in Figure 1 that both the numerical results and the theoretical results of Theorem 2 are consistent. It is seen that the solutions of the system converge to the steady state $Q_1=(0.99595,4.52267,2.32229,2.10913,2.07104,6.50639)$ for all the three initial values IV1–IV3.

4.2. Simulations for Virus Model with Latency

In this subsection, we show the numerical results for the system in Equations (19)–(23) with parameters values given in

Table 2. The values of the parameters of the model in Equations (19)–(23).

Parameter	Value	Parameter	Value
${\mu }_1$	$1.826$	$d_1$	$0.7979$
${\mu }_2$	$3.198$	$d_2$	$0.5$
le $b_1=b_2$	$0.04$	q	$0.5964$
c	$0.5$	r	$0.4418$
$m_1$=$m_2$	$2.02$	$\eta$	$1.402$
${\lambda }_1={\lambda }_2$	$0.1$	$\delta$	$1.251$
${ϵ}_1$=${ϵ}_2$	$0.4441$	${\theta }_1={\theta }_2$	$0.5$
$p_1$=$p_2$	varied

. The effect of parameters $p_i$ and $p_2$ on the qualitative behavior of the system is discussed below. We take $p=p_1=p_2.$ The initial values are chosen as: $S_1\left(0\right)=2.25$,$S_2\left(0\right)=6.0$,$L_1\left(0\right)=0.1,L_2\left(0\right)=0.05,I_1\left(0\right)=0.1,I_2\left(0\right)=0.5,V\left(0\right)=0.4$ and $B\left(0\right)=1.4.$ In

Table 3. The values of steady states, ${\mathcal{R}}_0^L$ for the model in Equations (19)–(23) with different values of p.

$\mathit{p}={\mathit{p}}_1={\mathit{p}}_2$	Steady States	${\mathit{R}}_0^{\mathit{L}}$
$0.1$	$Q_0=(2.26262,6.28133,0,0,0,0,0,1.23318)$	$0.3558$
$0.2$	$Q_0=(2.26262,6.28133,0,0,0,0,0,1.23318)$	$0.4744$
$0.3$	$Q_0=(2.26262,6.28133,0,0,0,0,0,1.23318)$	$0.5930$
$0.4$	$Q_0=(2.26262,6.28133,0,0,0,0,0,1.23318)$	$0.7116$
$0.5$	$Q_0=(2.26262,6.28133,0,0,0,0,0,1.23318)$	$0.8302$
$0.6$	$Q_0=(2.26262,6.28133,0,0,0,0,0,1.23318)$	$0.9488$
$0.643144$	$Q_0=(2.26262,6.28133,0,0,0,0,0,1.23318)$	$1.0000$
$0.7$	$Q_1=(2.23995,6.18216,0.04480,0.12364,0.07115,0.19637,0.43238,1.35484)$	$1.0674$
$0.8$	$Q_1=(2.22562.\mathrm{},6.12005,0.02968,0.08160,0.09707,0.26693,0.56363,1.44657)$	$1.1860$
$0.9$	$Q_1=(2.20775.\mathrm{},6.04324,0.01472,0.04029,0.13390,0.36652,0.72967,1.58210)$	$1.3046$

, we have calculated the values of the steady states and ${\mathcal{R}}_0^L$ for different values of p. It is clearly seen that, as p is increased, ${\mathcal{R}}_0^L$ is also increased. Let $p^{cr}$ be the value of p, such that

$${\mathcal{R}}_0^L=\stackrel{2}{\sum _{i=1}}\frac{m_i{b}_i{\mu }_i\delta ({\theta }_i{p}^{cr}+{\lambda }_i)}{d_i{ϵ}_i({\theta }_i+{\lambda }_i)(r\delta +q\eta )}=1.$$

Using the data given in Table 2, we obtain $p^{cr}=0.643144.$ In Figure , we can see that, for $p\le 0.643144,$ the trajectory of the system will converge to $Q_0^L$ and, for $p>0.643144,$ the trajectory will converge to $Q_1^L.$ This shows that, the factor $1-p$ plays the role of a controller which can used to stabilize the system around $Q_0^L.$ Biologically, the factor $1-p$ plays the role of an antiviral treatment which can be applied to eradicate the virus from the body.

5. Conclusions and Discussion

Most of the existing mathematical models of viral infection study the viral infection and production in one or two classes of target cells.However, HIV and CHIKV can infect three and five types of target cells, respectively. In this paper, we have studied two within-host virus dynamics models with antibody immune response and with n classes of target cells. In the second model, we have considered two types of infected cells, latently infected cells (such cells contain the viruses but are not producing it) and the actively infected cells (such cells are producing the viruses). We have shown that the solutions of each model are nonnegative and bounded, which ensure the well-posedness of the models. For each model, we have derived a biological threshold number ${\mathcal{R}}_0$ (the basic reproduction number) which fully determines the existence and stability of the two steady states of the model. We have investigated the global stability of the steady states of the model by using Lyapunov method and LaSalle’s invariance principle. We have proven that: (i) if ${\mathcal{R}}_0\le 1$ (${\mathcal{R}}_0^L\le 1$), then the virus-free steady state $Q_0$ ($Q_0^L$) is globally asymptotically stable and the virus is predicted to be completely cleared from infected individuals; and (ii) if ${\mathcal{R}}_0>1$ (${\mathcal{R}}_0^L>1)$, then the endemic steady state $Q_1$ ($Q_1^L$) is globally asymptotically stable and a chronic virus infection is attained. We have conducted numerical simulations and have shown that both the theoretical and numerical results are consistent. Our analysis extends some existing results in the literature. For example, the global stability was analyzed for a model with one target cell population [36].

The model in Equations (1)–(4) has three steady states virus-free steady state $Q_0^H,$ endemic steady state without antibody immune response $Q_1^H$ and endemic steady state with antibody immune response $Q_2^H$. Moreover, the existence and stability of the steady states are determined by two threshold parameters, the basic reproduction number ${\mathcal{R}}_0^H$ (which determines whether or not the disease will progress) and the antibody immune response activation number ${\mathcal{R}}_1^H$ (which determines whether a persistent antibody immune response can be established ) [30], where

$${\mathcal{R}}_0^H=\frac{mb\mu }{dϵr},{\mathcal{R}}_1^H=\frac{{\mathcal{R}}_0^H}{1+\frac{b\delta }{dc}}.$$

We note that the values of the parameters q, c and $\delta$ have no impact on the values of ${\mathcal{R}}_0^H$. Thus, the model in Equations (1)–(4) implies that the antibody immune response do not play a role in clearing the viruses but can play a significant role in reducing the infection progress.

Our models have two steady states and their existence and stability are determined by one threshold parameter ${\mathcal{R}}_0$. This is because of considering the production rate of the B cells $\eta$. The basic reproduction number of the model in Equations (9)–(12) in the case of $n=1$ is given by:

$${\mathcal{R}}_0=\frac{m_1{b}_1{\mu }_1\delta }{d_1{ϵ}_1(r\delta +q\eta )}.$$

We can see that ${\mathcal{R}}_0$ depends on the parameter $\eta$. Therefore, if the production rate of the B cells $\eta$ is increased such that ${\mathcal{R}}_0<1$, then $Q_0$ is globally asymptotically stable. Thus, the model in Equations (9)–(12) implies that the antibody immune response can clear the virus from body.

In our proposed models, we have only considered one arm of the immune system which is based on the antibodies. CTL cells play a prominent role in achieving the best representation of the dynamics of some types of viruses. The virus dynamics model with n categories of target cells and CTL immune response can be given as:

$$\begin{array}{cc}\hfill {\dot{S}}_i& ={\mu }_i-d_i{S}_i-b_i{S}_{i}V,i=1,\dots ,n,\hfill \\ \hfill {\dot{I}}_i& =b_i{S}_{i}V-{ϵ}_i{I}_i-{\overline{q}}_i{I}_{i}Z,i=1,\dots ,n,\hfill \\ \hfill \dot{V}& =\stackrel{n}{\sum _{i=1}}{m}_i{I}_i-rV,\hfill \\ \hfill \dot{Z}& =\overline{\eta }+\stackrel{n}{\sum _{i=1}}{\overline{c}}_i{I}_{i}Z-\overline{\delta }Z.\hfill \end{array}$$

In addition to this model, one can formulate a virus dynamics model with n categories of target cells and both antibodies and CTL immune response.

5.1. Effects of Latency on the Virus Dynamics

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

$$\begin{array}{cc}\hfill {\dot{S}}_i& ={\mu }_i-d_i{S}_i-(1-u)b_i{S}_{i}V,i=1,\dots ,n,\hfill \end{array}$$

(31)

$$\begin{array}{cc}\hfill {\dot{I}}_i& =(1-u)b_i{S}_{i}V-{ϵ}_i{I}_i,i=1,\dots ,n,\hfill \end{array}$$

(32)

$$\begin{array}{cc}\hfill \dot{V}& =\stackrel{n}{\sum _{i=1}}{m}_i{I}_i-rV-qBV,\hfill \end{array}$$

(33)

$$\begin{array}{cc}\hfill \dot{B}& =\eta +cBV-\delta B.\hfill \end{array}$$

(34)

Consequently, the parameter ${\mathcal{R}}_0$ for the system in Equations (31)–(34) is given by

$${\mathcal{R}}_0\left(u\right)=(1-u)\stackrel{n}{\sum _{i=1}}\frac{m_i{b}_i{\mu }_i\delta }{d_i{ϵ}_i(r\delta +q\eta )}.$$

The model in Equations (19)–(23) under the effect of treatment is given by:

$$\begin{array}{cc}\hfill {\dot{S}}_i& ={\mu }_i-d_i{S}_i-(1-u)b_i{S}_{i}V,i=1,\dots ,n,\hfill \end{array}$$

(35)

$$\begin{array}{cc}\hfill {\dot{L}}_i& =(1-u)(1-p_i)b_i{S}_{i}V-({\theta }_i+{\lambda }_i)L_i,i=1,\dots ,n,\hfill \end{array}$$

(36)

$$\begin{array}{cc}\hfill {\dot{I}}_i& =(1-u)p_i{b}_i{S}_{i}V+{\lambda }_i{L}_i-{ϵ}_i{I}_i,i=1,\dots ,n,\hfill \end{array}$$

(37)

$$\begin{array}{cc}\hfill \dot{V}& =\stackrel{n}{\sum _{i=1}}{m}_i{I}_i-rV-qBV,\hfill \end{array}$$

(38)

$$\begin{array}{cc}\hfill \dot{B}& =\eta +cBV-\delta B.\hfill \end{array}$$

(39)

The parameter ${\mathcal{R}}_0^L$ for the system in Equations (35)–(39) is given by

$${\mathcal{R}}_0^L\left(u\right)=(1-u)\stackrel{n}{\sum _{i=1}}\frac{m_i{b}_i{\mu }_i\delta ({\theta }_i{p}_i+{\lambda }_i)}{d_i{ϵ}_i({\theta }_i+{\lambda }_i)(r\delta +q\eta )}.$$

Since $0<p_i<1$, then

$${\mathcal{R}}_0^L\left(u\right)=(1-u)\stackrel{n}{\sum _{i=1}}\frac{m_i{b}_i{\mu }_i\delta ({\theta }_i{p}_i+{\lambda }_i)}{d_i{ϵ}_i({\theta }_i+{\lambda }_i)(r\delta +q\eta )}<(1-u)\stackrel{n}{\sum _{i=1}}\frac{m_i{b}_i{\mu }_i\delta }{d_i{ϵ}_i(r\delta +q\eta )}={\mathcal{R}}_0\left(u\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 able to clear the viruses from the body. We determine $u_{crit}$ and $u_{crit}^L$ that make

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

to stabilize the system in Equations (31)–(39) around $Q_0$ and $Q_0^L$, respectively. Now, we calculate $u_{crit}$ and $u_{crit}^L$ as:

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

Clearly, ${\mathcal{R}}_0^L\left(0\right)<{\mathcal{R}}_0\left(0\right)$ and thus $u_{crit}^L<u_{crit}$. Therefore, the drug efficacy necessary to drive the system to the virus-free steady state is actually less for the system in Equations (35)–(39) than that for the system in Equations (31)–(34).