# Global Properties of Virus Dynamics Models with Multitarget Cells and Discrete-Time Delays

Authors: A. M. Elaiw; M. A. Alghamdi
Journal: Discrete Dynamics in Nature and Society
 (2011)

---

## Abstract

We propose a class of virus dynamics models with multitarget cells and multiple
intracellular delays and study their global properties. The first model is a 5-dimensional system of nonlinear delay differential equations (DDEs) that describes the interaction of the virus with two classes of target cells. The second model is a (2n+1)-dimensional system of nonlinear DDEs that describes the dynamics of the virus, n classes of uninfected target cells, and n classes of infected target cells. The third model generalizes the second one by assuming that the incidence rate of infection is given by saturation functional response. Two types of discrete time delays are incorporated into these models to describe (i) the latent period between the time the target cell is contacted by the virus particle and the time the virus enters the cell, (ii) the latent period between the time the virus has penetrated into a cell and the time of the emission of infectious (mature) virus particles. Lyapunov functionals are constructed to establish the global asymptotic stability of the uninfected and infected steady states of these models. We have proven that if the basic reproduction number R0 is less than unity, then the uninfected steady state is globally asymptotically stable, and if R0>1 (or if the infected steady state exists), then the infected steady state is globally asymptotically stable.


---

## Body

## 1. Introduction

Nowadays, various types of viruses infect the human body and cause serious and dangerous diseases. Mathematical modeling and model analysis of virus dynamics have attracted the interests of mathematicians during the recent years, due to their importance in understanding the associated characteristics of the virus dynamics and guiding in developing efficient antiviral drug therapies. Several mathematical models have been proposed in the literature to describe the interaction of the virus with the target cells [1]. Some of these models are given by a system of nonlinear ordinary differential equations (ODEs). Others are given by a system of nonlinear delay differential equations (DDEs) to account the intracellular time delays. The basic virus dynamics model with intracellular discrete time delay has been proposed in [2] and given by(1.1)ẋ(t)=λ-dx(t)-βx(t)v(t),(1.2)ẏ(t)=e-mτβx(t-τ)v(t-τ)-ay(t),(1.3)v̇(t)=py(t)-cv(t),
where x(t),y(t), and v(t) represent the populations of uninfected target cells, infected cells, and free virus particles at time t, respectively. Here, λ represents the rate of which new target cells are generated from sources within the body, d is the death rate constant, and β is the infection rate constant. Equation (1.2) describes the population dynamics of the infected cells and shows that they die with rate constant a. The virus particles are produced by the infected cells with rate constant p, and are removed from the system with rate constant c. The parameter τ accounts for the time between viral entry into the target cell and the production of new virus particles. The recruitment of virus-producing cells at time t  is given by the number of cells that were newly infected cells at time t-τ and are still alive at time t. The probability of surviving the time period from t-τ to t is e-mτ, where m is the constant death rate of infected cells but not yet virus-producing cells.A great effort has been made in developing various mathematical models of viral infections with discrete or distributed delays and studying their basic and global properties, such as positive invariance properties, boundedness of the model solutions and stability analysis [3–19]. In [20–24], multiple inracellular delays have been incroporated into the virus dynamics model. Most of the existing models are based on the assumption that the virus attacks one class of target cells (e.g., CD4+ T cells in case of HIV or hepatic cells in case of HCV and HBV). Since the interactions of some types of viruses inside the human body is not very clear and complicated, therefore, the virus may attack more than one class of target cells. Hence, virus dynamics models describing the interaction of the virus with more than one class of target cells are needed. In case of HIV infection, Perelson et al. [25] observed that the HIV attack two classes of target cells, CD4+ T cells and macrophages. In [26, 27], an HIV model with two target cells has been proposed. In very recent works [28–30], we have proposed several HIV models with two target cells and investigated the global asymptotic stability of their steady states. In [31], we have proposed a class of virus dynamics models with multitarget cells. However, the intracellular time delay has been neglected in [26–31].The purpose of this paper is to propose a class of virus dynamics models with multitarget cells and establish the global stability of their steady states. The first model considers the interaction of the virus with two classes of target cells. In the second model, we assume that the virus attacksn classes of target cells. The third model generalizes the second one by assuming that the infection rate is given by saturation functional response. We incorporate two types of discrete time delays into these models describing (i) the time between the target cell is contacted by the virus particle and the contacting virus enters the cell, (ii) the time between the virus has penetrated into a cell and the emission of infectious (mature) virus particles. The global stability of these models is established using Lyapunov functionals, which are similar in nature to those used in [11, 20]. We prove that the global dynamics of these models are determined by the basic reproduction number R0. If R0≤1, then the uninfected steady state is globally asymptotically stable (GAS). If R0>1 (or if the infected steady state exists), then the infected steady state is GAS for all time delays.


## 2. Virus Dynamics Model with Two Target Cells and Delays

In this section, we introduce a mathematical model of virus infection with two classes of target cells. This model can describe the HIV dynamics with two classes of target cells, CD4+ T cells and macrophages [26, 27]. This model can be considered as an extension of the models given in [11, 26, 27]: (2.1)ẋ1(t)=λ1-d1x1(t)-β1x1(t)v(t),(2.2)ẏ1(t)=e-m1τ1β1x1(t-τ1)v(t-τ1)-a1y1(t),(2.3)ẋ2(t)=λ2-d2x2(t)-β2x2(t)v(t),(2.4)ẏ2(t)=e-m2τ2β2x2(t-τ2)v(t-τ2)-a2y2(t),(2.5)v̇(t)=e-n1ω1p1y1(t-ω1)+e-n2ω2p2y2(t-ω2)-cv(t),
where x1and x2 represent the populations of the two classes of uninfected target cells; y1and y2 are the populations of the infected cells. The population of the target cells are described by (2.1) and (2.3), where λ1 and λ2 represent the rates of which new target cells are generated, d1 and d2 are the death rate constants, and β1 and β2 are the infection rate constants. Equations (2.2) and (2.4) describe the population dynamics of the two classes of infected cells and show that they die with rate constants a1 and a2. The virus particles are produced by the two classes of infected cells with rate constants p1 and p2 and are cleared with rate constant c. Here the parameter τi accounts for the time between the target cells of class i are contacted by the virus particle and the contacting virus enters the cells. The recruitment of virus-producing cells at time t  is given by the number of cells that were newly infected cells at times t-τi and are still alive at time t. Also, mi is assumed to be a constant death rate for infected target cells, but not yet virus-producing cells. Thus, the probability of surviving the time period from t-τi to t is e-miτi, i=1,2. The time between the virus has penetrated into a target cell of class i and the emission of infectious (matures) virus particles is represented by ωi. The probability of survival of an immature virus is given by e-niωi, where ni is constant.


### 2.1. Initial Conditions

The initial conditions for system (2.1)–(2.5) take the form(2.6)x1(θ)=φ1(θ),y1(θ)=φ2(θ),x2(θ)=φ3(θ),y2(θ)=φ4(θ),v(θ)=φ5(θ),φi(θ)≥0,θ∈[-max⁡{τ1,τ2,ω1,ω2},0],i=1,…,5,
where (φ1(θ),…,φ5(θ))∈C([-max{τ1,τ2,ω1,ω2},0],ℝ+5), the Banach space of continuous functions mapping the interval [-max{τ1,τ2,ω1,ω2},0]into ℝ+5, where ℝ+5={(z1,z2,…,z5):zi≥0}.By the fundamental theory of functional differential equations [32], system (2.1)–(2.5) has a unique solution(x1(t),y1(t),x2(t),y2(t),v(t)) satisfying the initial conditions (2.6).


### 2.2. Nonnegativity and Boundedness of Solutions

In the following, we establish the nonnegativity and boundedness of solutions of (2.1)–(2.5) with initial conditions (2.6).Proposition 2.1.
Let(x1(t),y1(t),x2(t),y2(t),v(t)) be any solution of (2.1)–(2.5) satisfying the initial conditions (2.6), then x1(t), y1(t), x2(t), y2(t),and v(t) are all nonnegative for t≥0 and ultimately bounded.Proof.
From (2.1) and (2.3), we have
(2.7)xi(t)=xi(0)e-∫0t(di+βiv(ξ))dξ+λi∫0te-∫ηt(di+βiv(ξ))dξdη,i=1,2,
which indicates that x1(t)≥0, x2(t)≥0 for all t≥0. Now from (2.2), (2.4), and (2.5), we have
(2.8)yi(t)=yi(0)e-ait+e-miτi∫0tβixi(η-τi)v(η-τi)e-ai(t-η)dη,i=1,2,v(t)=v(0)e-ct+∫0t[e-n1ω1p1y1(η-ω1)+e-n2ω2p2y2(η-ω2)]e-c(t-η)dη,
confiming that y1(t)≥0, y2(t)≥0, v(t)≥0 for all t∈[0,max{τ1,τ2,ω1,ω2}]. By a recursive argument, we obtain y1(t)≥0, y2(t)≥0, v(t)≥0 for all t≥0.
To show the boundedness of the solutions, we letX1(t)=e-m1τ1x1(t-τ1-ω1)+y1(t-ω1) and X2(t)=e-m2τ2x2(t-τ2-ω2)+y2(t-ω2), then
(2.9)Ẋ1(t)≤λ1e-m1τ1-σ1X1(t),Ẋ2(t)≤λ2e-m2τ2-σ2X2(t),
where σ1=min{d1,a1} and σ2=min{d2,a2}. Hence, limsupt→∞X1(t)≤L1, and limsupt→∞X2(t)≤L2, where L1=λ1e-m1τ1/σ1 and L2=λ2e-m2τ2/σ2. On the other hand,
(2.10)v̇(t)≤e-n1ω1p1L1+e-n2ω2p2L2-cv,
then limsupt→∞v(t)≤L3, where L3=(e-n1ω1p1L1+e-n2ω2p2L2)/c. It follows that the solution (x1(t),y1(t),x2(t),y2(t),v(t)) is ultimately bounded.


### 2.3. Steady States

It will be explained in the following that the global behavior of model (2.1)–(2.5) crucially depends on the basic reproduction number given by(2.11)R0=e-(m1τ1+n1ω1)p1β1a2x10+e-(m2τ2+n2ω2)p2β2a1x20a1a2c,
where x10=λ1/d1 and x20=λ2/d2. We observe that R0 can be written as(2.12)R0=R1+R2,
where(2.13)R1=e-(m1τ1+n1ω1)p1β1λ1a1d1c,R2=e-(m2τ2+n2ω2)p2β2λ2a2d2c
are the basic reproduction numbers of each class of target cell dynamics separately (see [28]).Following the same line as in [28], we can show that if R0≤1, then system (2.1)–(2.5) has only one steady state E0=(x10,0,x20,0,0) which is called uninfected steady state, and if R0>1, then system (2.1)–(2.5) has two steady states E0 and infected steady state E1=(x1*,y1*,x2*,y2*,v*). The coordinates of the infected steady state are given by (2.14)x1*={α5cα1α5+α2α3,ifα4=0,-(α1α5+α2α3-α4c)+(α1α5+α2α3-α4c)2+4α1α4α5c2α1α4,ifα4≠0,x2*={α3cα1α5+α2α3,ifα4=0,cα2+(α1α5+α2α3-α4c)-(α1α5+α2α3-α4c)2+4α1α4α5c2α2α4,ifα4≠0,y1*=d1a1em1τ1(x10x1*-1)x1*,y2*=d2a2em2τ2(x20x2*-1)x2*,v*=d1β1(x10x1*-1),
where (2.15)α1=e-(n1ω1+m1τ1)p1β1a1,α2=e-(n2ω2+m2τ2)p2β2a2,α3=λ2β1,α4=β1d2-β2d1,α5=λ1β2.


### 2.4. Global Stability

In this section, we prove the global stability of the uninfected and infected steady states of system (2.1)–(2.5). The strategy of the proof is to use suitable Lyapunov functionals which are similar in nature to those used in [11, 20]. Next we will use the following notation: z=z(t), for any z∈{x1,y1,x2,y2,v}. We also define a function H:ℝ>0→ℝ≥0 as (2.16)H(z)=z-1-ln⁡z.
It is clear that H(z)≥0 for any z>0 and H has the global minimum H(1)=0.Theorem 2.2.
(i) IfR0≤1, then E0 is GAS for any τ1,τ2,ω1,ω2≥0.
(ii) IfR0>1, then E1 is GAS for any τ1,τ2,ω1,ω2≥0.Proof.
(i) We consider a Lyapunov functional(2.17)W1=W11+W12+W13,
where
(2.18)W11=e-m1τ1x10H(x1x10)+y1+γ(e-m2τ2x20H(x2x20)+y2)+a1p1en1ω1v,W12=e-m1τ1∫0τ1β1x1(t-θ)v(t-θ)dθ+γe-m2τ2∫0τ2β2x2(t-θ)v(t-θ)dθ,W13=a1∫0ω1y1(t-θ)dθ+γa2∫0ω2y2(t-θ)dθ,
where γ=(p2a1/p1a2)en1ω1-n2ω2. We note that W1 is defined, continuous, and positive definite for all (x1,y1,x2,y2,v)>0. Also, the global minimum W1=0 occurs at the uninfected steady state E0.
The time derivatives ofW11 and W12 are given by
(2.19)dW11dt=e-m1τ1(1-x10x1)ẋ1+ẏ1+γ(e-m2τ2(1-x20x2)ẋ2+ẏ2)+a1p1en1ω1v̇,dW12dt=e-m1τ1∫0τ1ddtβ1x1(t-θ)v(t-θ)dθ+γe-m2τ2∫0τ2ddtβ2x2(t-θ)v(t-θ)dθ=-e-m1τ1∫0τ1ddθβ1x1(t-θ)v(t-θ)dθ-γe-m2τ2∫0τ2ddθβ2x2(t-θ)v(t-θ)dθ=e-m1τ1(β1x1v-β1x1(t-τ1)v(t-τ1))+γe-m2τ2[β2x2v-β2x2(t-τ2)v(t-τ2)].
Similarly, dW13/dt is given by
(2.20)dW13dt=a1(y1-y1(t-ω1))+γa2(y2-y2(t-ω2)).
It follows that
(2.21)dW1dt=e-m1τ1(1-x10x1)(λ1-d1x1-β1x1v)+e-m1τ1β1x1(t-τ1)v(t-τ1)-a1y1+γ[e-m2τ2(1-x20x2)(λ2-d2x2-β2x2v)+e-m2τ2β2x2(t-τ2)v(t-τ2)-a2y2]+a1p1en1ω1[e-n1ω1p1y1(t-ω1)+e-n2ω2p2y2(t-ω2)-cv]+e-m1τ1[β1x1v-β1x1(t-τ1)v(t-τ1)]+γe-m2τ2[β2x2v-β2x2(t-τ2)v(t-τ2)]+a1[y1-y1(t-ω1)]+γa2[y2-y2(t-ω2)]=e-m1τ1λ1[2-x10x1-x1x10]+e-m2τ2γλ2[2-x20x2-x2x20]+e-m1τ1β1x10v+e-m2τ2γβ2x20v-a1cp1en1ω1v=e-m1τ1λ1[2-x10x1-x1x10]+e-m2τ2γλ2[2-x20x2-x2x20]+a1cen1ω1p1(R0-1)v.
Since the arithmetical mean is greater than or equal to the geometrical mean, then the first two terms of (2.21) are less than or equal to zero. Therefore, if R0≤1, then dW1/dt≤0 for all x1,x2,v>0. By Theorem  5.3.1 in [32], the solutions of system (2.1)–(2.5) limit to M, the largest invariant subset of {dW1/dt=0}. Clearly, it follows from (2.21) that dW1/dt=0 if and only if x1=x10, x2=x20, v=0. Noting that M is invariant, for each element of M we have v=0, v̇=0. From (2.5) we drive that
(2.22)0=v̇=e-n1ω1p1y1(t-ω1)+e-n2ω2p2y2(t-ω2).
Since y1(t-θ)≥0 and y2(t-θ)≥0 for all θ∈[0,max{τ1,τ2,ω1,ω2}], then e-n1ω1p1y1(t-ω1)+e-n2ω2p2y2(t-ω2)=0 if and only if y1(t-ω1)=y2(t-ω2)=0. Hence, dW1/dt=0 if and only if x1=x10, x2=x20, y1=y2=v=0. From LaSalle’s invariance principle, E0 is GAS for any τ1,τ2,ω1,ω2≥0.
(ii) Define a Lyapunov functional as(2.23)W2=e-m1τ1x1*H(x1x1*)+y1*H(y1y1*)+γ[e-m2τ2x2*H(x2x2*)+y2*H(y2y2*)]+a1p1en1ω1v*H(vv*)+e-m1τ1β1x1*v*∫0τ1H(x1(t-θ)v(t-θ)x1*v*)dθ+e-m2τ2γβ2x2*v*∫0τ2H(x2(t-θ)v(t-θ)x2*v*)dθ+a1y1*∫0ω1H(y1(t-θ)y1*)dθ+γa2y2*∫0ω2H(y2(t-θ)y2*)dθ.
Differentiating with respect to time yields
(2.24)dW2dt=e-m1τ1(1-x1*x1)(λ1-d1x1-β1x1v)+(1-y1*y1)(e-m1τ1β1x1(t-τ1)v(t-τ1)-a1y1)+γ[e-m2τ2(1-x2*x2)(λ2-d2x2-β2x2v)+(1-y2*y2)(e-m2τ2β2x2(t-τ2)v(t-τ2)-a2y2)]+a1p1en1ω1(1-v*v)(e-n1ω1p1y1(t-ω1)+e-n2ω2p2y2(t-ω2)-cv)+e-m1τ1[β1x1v-β1x1(t-τ1)v(t-τ1)+β1x1*v*ln⁡(x1(t-τ1)v(t-τ1)x1v)]+γe-m2τ2[β2x2v-β2x2(t-τ2)v(t-τ2)+β2x2*v*ln⁡(x2(t-τ2)v(t-τ2)x2v)]+a1[y1-y1(t-ω1)+y1*ln⁡(y1(t-ω1)y1)]+a2γ[y2-y2(t-ω2)+y2*ln⁡(y2(t-ω2)y2)]=e-m1τ1(1-x1*x1)(λ1-d1x1)+e-m1τ1β1x1*v-e-m1τ1y1*β1x1(t-τ1)v(t-τ1)y1+a1y1*+γ[e-m2τ2(1-x2*x2)(λ2-d2x2)+e-m2τ2β2x2*v-e-m2τ2y2*β2x2(t-τ2)v(t-τ2)y2+a2y2*y2*β2x2(t-τ2)v(t-τ2)y2]-a1v*y1(t-ω1)v-γa2v*y2(t-ω2)v-a1cp1en1ω1v+a1cp1en1ω1v*+e-m1τ1β1x1*v*ln⁡(x1(t-τ1)v(t-τ1)x1v)+γe-m2τ2β2x2*v*ln⁡(x2(t-τ2)v(t-τ2)x2v)+a1y1*ln⁡(y1(t-ω1)y1)+γa2y2*ln⁡(y2(t-ω2)y2).
Using the infected steady state E1 conditions
(2.25)λ1=d1x1*+β1x1*v*,λ2=d2x2*+β2x2*v*,a1y1*em1τ1=β1x1*v*,a2y2*em2τ2=β2x2*v*,cv*=p1e-n1ω1y1*+p2e-n2ω2y2*,
we obtain
(2.26)dW2dt=e-m1τ1[d1x1*+a1y1*em1τ1-d1x1-x1*x1(d1x1*+a1y1*em1τ1-d1x1)]+e-m1τ1β1x1*v*(vv*)-a1y1*y1*x1(t-τ1)v(t-τ1)y1x1*v*+2a1y1*+γ[y2*x2(t-τ2)v(t-τ2)y2x2*v*e-m2τ2(d2x2*+a2y2*em2τ2-d2x2)-e-m2τ2x2*x2(d2x2*+a2y2*em2τ2-d2x2)+e-m2τ2β2x2*v*(vv*)-a2y2*y2*x2(t-τ2)v(t-τ2)y2x2*v*+2a2y2*]-a1y1*v*y1(t-ω1)vy1*-γa2y2*v*y2(t-ω2)vy2*-a1cp1en1ω1v*(vv*)+a1y1*ln⁡(x1(t-τ1)v(t-τ1)x1v)+γa2y2*ln⁡(x2(t-τ2)v(t-τ2)x2v)+a1y1*ln⁡(y1(t-ω1)y1)+γa2y2*ln⁡(y2(t-ω2)y2)=e-m1τ1d1x1*(2-x1x1*-x1*x1)-a1y1*x1*x1-a1y1*y1*x1(t-τ1)v(t-τ1)y1x1*v*-a1y1*v*y1(t-ω1)vy1*+3a1y1*+a1y1*ln⁡(x1(t-τ1)v(t-τ1)y1(t-ω1)x1vy1)+γ[e-m2τ2d2x2*(2-x2x2*-x2*x2)-a2y2*x2*x2-a2y2*y2*x2(t-τ2)v(t-τ2)y2x2*v*-a2y2*v*y2(t-ω2)vy2*+3a2y2*+a2y2*ln⁡(x2(t-τ2)v(t-τ2)y2(t-ω2)x2vy2)]+a1p1en1ω1(e-n1ω1p1y1*+e-n2ω2p2y2*-cv*)vv*.
From (2.25), we see that the last term in (2.26) vanishes. Then, using the following equalities:
(2.27)ln⁡(x1(t-τ1)v(t-τ1)y1(t-ω1)x1vy1)=ln⁡(x1*x1)+ln⁡(y1(t-ω1)v*y1*v)+ln⁡(y1*x1(t-τ1)v(t-τ1)y1x1*v*),ln⁡(x2(t-τ2)v(t-τ2)y2(t-ω2)x2vy2)=ln⁡(x2*x2)+ln⁡(y2(t-ω2)v*y2*v)+ln⁡(y2*x2(t-τ2)v(t-τ2)y2x2*v*),
we can rewrite (2.26) as
(2.28)dW2dt=e-m1τ1d1x1*(2-x1x1*-x1*x1)+e-m2τ2γd2x2*(2-x2x2*-x2*x2)-a1y1*[H(x1*x1)+H(y1(t-ω1)v*y1*v)+H(y1*x1(t-τ1)v(t-τ1)y1x1*v*)]-γa2y2*[H(x2*x2)+H(y2(t-ω2)v*y2*v)+H(y2*x2(t-τ2)v(t-τ2)y2x2*v*)].
Since the arithmetical mean is greater than or equal to the geometrical mean, then the first two terms of (2.28) are less than or equal to zero. It is easy to see that if x1*,y1*,x2*,y2*,v*>0, then dW2/dt≤0. By [32,Theorem  5.3.1], the solutions of system (2.1)–(2.5) limit to M, the largest invariant subset of {dW2/dt=0}. It can be seen that dW2/dt=0 if and only if x1=x1*,x2=x2*,v=v*, and H=0, that is,
(2.29)y1(t-ω1)v*y1*v=y2(t-ω2)v*y2*v=y1*x1(t-τ1)v(t-τ1)y1x1*v*=y2*x2(t-τ2)v(t-τ2)y2x2*v*=1.
If v=v*, then from (2.29), we have y1=y1* and y2=y2*, and hence dW2/dt equal to zero at E1. LaSalle’s invariance principle implies global stability of E1.


## 2.1. Initial Conditions

The initial conditions for system (2.1)–(2.5) take the form(2.6)x1(θ)=φ1(θ),y1(θ)=φ2(θ),x2(θ)=φ3(θ),y2(θ)=φ4(θ),v(θ)=φ5(θ),φi(θ)≥0,θ∈[-max⁡{τ1,τ2,ω1,ω2},0],i=1,…,5,
where (φ1(θ),…,φ5(θ))∈C([-max{τ1,τ2,ω1,ω2},0],ℝ+5), the Banach space of continuous functions mapping the interval [-max{τ1,τ2,ω1,ω2},0]into ℝ+5, where ℝ+5={(z1,z2,…,z5):zi≥0}.By the fundamental theory of functional differential equations [32], system (2.1)–(2.5) has a unique solution(x1(t),y1(t),x2(t),y2(t),v(t)) satisfying the initial conditions (2.6).


## 2.2. Nonnegativity and Boundedness of Solutions

In the following, we establish the nonnegativity and boundedness of solutions of (2.1)–(2.5) with initial conditions (2.6).Proposition 2.1.
Let(x1(t),y1(t),x2(t),y2(t),v(t)) be any solution of (2.1)–(2.5) satisfying the initial conditions (2.6), then x1(t), y1(t), x2(t), y2(t),and v(t) are all nonnegative for t≥0 and ultimately bounded.Proof.
From (2.1) and (2.3), we have
(2.7)xi(t)=xi(0)e-∫0t(di+βiv(ξ))dξ+λi∫0te-∫ηt(di+βiv(ξ))dξdη,i=1,2,
which indicates that x1(t)≥0, x2(t)≥0 for all t≥0. Now from (2.2), (2.4), and (2.5), we have
(2.8)yi(t)=yi(0)e-ait+e-miτi∫0tβixi(η-τi)v(η-τi)e-ai(t-η)dη,i=1,2,v(t)=v(0)e-ct+∫0t[e-n1ω1p1y1(η-ω1)+e-n2ω2p2y2(η-ω2)]e-c(t-η)dη,
confiming that y1(t)≥0, y2(t)≥0, v(t)≥0 for all t∈[0,max{τ1,τ2,ω1,ω2}]. By a recursive argument, we obtain y1(t)≥0, y2(t)≥0, v(t)≥0 for all t≥0.
To show the boundedness of the solutions, we letX1(t)=e-m1τ1x1(t-τ1-ω1)+y1(t-ω1) and X2(t)=e-m2τ2x2(t-τ2-ω2)+y2(t-ω2), then
(2.9)Ẋ1(t)≤λ1e-m1τ1-σ1X1(t),Ẋ2(t)≤λ2e-m2τ2-σ2X2(t),
where σ1=min{d1,a1} and σ2=min{d2,a2}. Hence, limsupt→∞X1(t)≤L1, and limsupt→∞X2(t)≤L2, where L1=λ1e-m1τ1/σ1 and L2=λ2e-m2τ2/σ2. On the other hand,
(2.10)v̇(t)≤e-n1ω1p1L1+e-n2ω2p2L2-cv,
then limsupt→∞v(t)≤L3, where L3=(e-n1ω1p1L1+e-n2ω2p2L2)/c. It follows that the solution (x1(t),y1(t),x2(t),y2(t),v(t)) is ultimately bounded.


## 2.3. Steady States

It will be explained in the following that the global behavior of model (2.1)–(2.5) crucially depends on the basic reproduction number given by(2.11)R0=e-(m1τ1+n1ω1)p1β1a2x10+e-(m2τ2+n2ω2)p2β2a1x20a1a2c,
where x10=λ1/d1 and x20=λ2/d2. We observe that R0 can be written as(2.12)R0=R1+R2,
where(2.13)R1=e-(m1τ1+n1ω1)p1β1λ1a1d1c,R2=e-(m2τ2+n2ω2)p2β2λ2a2d2c
are the basic reproduction numbers of each class of target cell dynamics separately (see [28]).Following the same line as in [28], we can show that if R0≤1, then system (2.1)–(2.5) has only one steady state E0=(x10,0,x20,0,0) which is called uninfected steady state, and if R0>1, then system (2.1)–(2.5) has two steady states E0 and infected steady state E1=(x1*,y1*,x2*,y2*,v*). The coordinates of the infected steady state are given by (2.14)x1*={α5cα1α5+α2α3,ifα4=0,-(α1α5+α2α3-α4c)+(α1α5+α2α3-α4c)2+4α1α4α5c2α1α4,ifα4≠0,x2*={α3cα1α5+α2α3,ifα4=0,cα2+(α1α5+α2α3-α4c)-(α1α5+α2α3-α4c)2+4α1α4α5c2α2α4,ifα4≠0,y1*=d1a1em1τ1(x10x1*-1)x1*,y2*=d2a2em2τ2(x20x2*-1)x2*,v*=d1β1(x10x1*-1),
where (2.15)α1=e-(n1ω1+m1τ1)p1β1a1,α2=e-(n2ω2+m2τ2)p2β2a2,α3=λ2β1,α4=β1d2-β2d1,α5=λ1β2.


## 2.4. Global Stability

In this section, we prove the global stability of the uninfected and infected steady states of system (2.1)–(2.5). The strategy of the proof is to use suitable Lyapunov functionals which are similar in nature to those used in [11, 20]. Next we will use the following notation: z=z(t), for any z∈{x1,y1,x2,y2,v}. We also define a function H:ℝ>0→ℝ≥0 as (2.16)H(z)=z-1-ln⁡z.
It is clear that H(z)≥0 for any z>0 and H has the global minimum H(1)=0.Theorem 2.2.
(i) IfR0≤1, then E0 is GAS for any τ1,τ2,ω1,ω2≥0.
(ii) IfR0>1, then E1 is GAS for any τ1,τ2,ω1,ω2≥0.Proof.
(i) We consider a Lyapunov functional(2.17)W1=W11+W12+W13,
where
(2.18)W11=e-m1τ1x10H(x1x10)+y1+γ(e-m2τ2x20H(x2x20)+y2)+a1p1en1ω1v,W12=e-m1τ1∫0τ1β1x1(t-θ)v(t-θ)dθ+γe-m2τ2∫0τ2β2x2(t-θ)v(t-θ)dθ,W13=a1∫0ω1y1(t-θ)dθ+γa2∫0ω2y2(t-θ)dθ,
where γ=(p2a1/p1a2)en1ω1-n2ω2. We note that W1 is defined, continuous, and positive definite for all (x1,y1,x2,y2,v)>0. Also, the global minimum W1=0 occurs at the uninfected steady state E0.
The time derivatives ofW11 and W12 are given by
(2.19)dW11dt=e-m1τ1(1-x10x1)ẋ1+ẏ1+γ(e-m2τ2(1-x20x2)ẋ2+ẏ2)+a1p1en1ω1v̇,dW12dt=e-m1τ1∫0τ1ddtβ1x1(t-θ)v(t-θ)dθ+γe-m2τ2∫0τ2ddtβ2x2(t-θ)v(t-θ)dθ=-e-m1τ1∫0τ1ddθβ1x1(t-θ)v(t-θ)dθ-γe-m2τ2∫0τ2ddθβ2x2(t-θ)v(t-θ)dθ=e-m1τ1(β1x1v-β1x1(t-τ1)v(t-τ1))+γe-m2τ2[β2x2v-β2x2(t-τ2)v(t-τ2)].
Similarly, dW13/dt is given by
(2.20)dW13dt=a1(y1-y1(t-ω1))+γa2(y2-y2(t-ω2)).
It follows that
(2.21)dW1dt=e-m1τ1(1-x10x1)(λ1-d1x1-β1x1v)+e-m1τ1β1x1(t-τ1)v(t-τ1)-a1y1+γ[e-m2τ2(1-x20x2)(λ2-d2x2-β2x2v)+e-m2τ2β2x2(t-τ2)v(t-τ2)-a2y2]+a1p1en1ω1[e-n1ω1p1y1(t-ω1)+e-n2ω2p2y2(t-ω2)-cv]+e-m1τ1[β1x1v-β1x1(t-τ1)v(t-τ1)]+γe-m2τ2[β2x2v-β2x2(t-τ2)v(t-τ2)]+a1[y1-y1(t-ω1)]+γa2[y2-y2(t-ω2)]=e-m1τ1λ1[2-x10x1-x1x10]+e-m2τ2γλ2[2-x20x2-x2x20]+e-m1τ1β1x10v+e-m2τ2γβ2x20v-a1cp1en1ω1v=e-m1τ1λ1[2-x10x1-x1x10]+e-m2τ2γλ2[2-x20x2-x2x20]+a1cen1ω1p1(R0-1)v.
Since the arithmetical mean is greater than or equal to the geometrical mean, then the first two terms of (2.21) are less than or equal to zero. Therefore, if R0≤1, then dW1/dt≤0 for all x1,x2,v>0. By Theorem  5.3.1 in [32], the solutions of system (2.1)–(2.5) limit to M, the largest invariant subset of {dW1/dt=0}. Clearly, it follows from (2.21) that dW1/dt=0 if and only if x1=x10, x2=x20, v=0. Noting that M is invariant, for each element of M we have v=0, v̇=0. From (2.5) we drive that
(2.22)0=v̇=e-n1ω1p1y1(t-ω1)+e-n2ω2p2y2(t-ω2).
Since y1(t-θ)≥0 and y2(t-θ)≥0 for all θ∈[0,max{τ1,τ2,ω1,ω2}], then e-n1ω1p1y1(t-ω1)+e-n2ω2p2y2(t-ω2)=0 if and only if y1(t-ω1)=y2(t-ω2)=0. Hence, dW1/dt=0 if and only if x1=x10, x2=x20, y1=y2=v=0. From LaSalle’s invariance principle, E0 is GAS for any τ1,τ2,ω1,ω2≥0.
(ii) Define a Lyapunov functional as(2.23)W2=e-m1τ1x1*H(x1x1*)+y1*H(y1y1*)+γ[e-m2τ2x2*H(x2x2*)+y2*H(y2y2*)]+a1p1en1ω1v*H(vv*)+e-m1τ1β1x1*v*∫0τ1H(x1(t-θ)v(t-θ)x1*v*)dθ+e-m2τ2γβ2x2*v*∫0τ2H(x2(t-θ)v(t-θ)x2*v*)dθ+a1y1*∫0ω1H(y1(t-θ)y1*)dθ+γa2y2*∫0ω2H(y2(t-θ)y2*)dθ.
Differentiating with respect to time yields
(2.24)dW2dt=e-m1τ1(1-x1*x1)(λ1-d1x1-β1x1v)+(1-y1*y1)(e-m1τ1β1x1(t-τ1)v(t-τ1)-a1y1)+γ[e-m2τ2(1-x2*x2)(λ2-d2x2-β2x2v)+(1-y2*y2)(e-m2τ2β2x2(t-τ2)v(t-τ2)-a2y2)]+a1p1en1ω1(1-v*v)(e-n1ω1p1y1(t-ω1)+e-n2ω2p2y2(t-ω2)-cv)+e-m1τ1[β1x1v-β1x1(t-τ1)v(t-τ1)+β1x1*v*ln⁡(x1(t-τ1)v(t-τ1)x1v)]+γe-m2τ2[β2x2v-β2x2(t-τ2)v(t-τ2)+β2x2*v*ln⁡(x2(t-τ2)v(t-τ2)x2v)]+a1[y1-y1(t-ω1)+y1*ln⁡(y1(t-ω1)y1)]+a2γ[y2-y2(t-ω2)+y2*ln⁡(y2(t-ω2)y2)]=e-m1τ1(1-x1*x1)(λ1-d1x1)+e-m1τ1β1x1*v-e-m1τ1y1*β1x1(t-τ1)v(t-τ1)y1+a1y1*+γ[e-m2τ2(1-x2*x2)(λ2-d2x2)+e-m2τ2β2x2*v-e-m2τ2y2*β2x2(t-τ2)v(t-τ2)y2+a2y2*y2*β2x2(t-τ2)v(t-τ2)y2]-a1v*y1(t-ω1)v-γa2v*y2(t-ω2)v-a1cp1en1ω1v+a1cp1en1ω1v*+e-m1τ1β1x1*v*ln⁡(x1(t-τ1)v(t-τ1)x1v)+γe-m2τ2β2x2*v*ln⁡(x2(t-τ2)v(t-τ2)x2v)+a1y1*ln⁡(y1(t-ω1)y1)+γa2y2*ln⁡(y2(t-ω2)y2).
Using the infected steady state E1 conditions
(2.25)λ1=d1x1*+β1x1*v*,λ2=d2x2*+β2x2*v*,a1y1*em1τ1=β1x1*v*,a2y2*em2τ2=β2x2*v*,cv*=p1e-n1ω1y1*+p2e-n2ω2y2*,
we obtain
(2.26)dW2dt=e-m1τ1[d1x1*+a1y1*em1τ1-d1x1-x1*x1(d1x1*+a1y1*em1τ1-d1x1)]+e-m1τ1β1x1*v*(vv*)-a1y1*y1*x1(t-τ1)v(t-τ1)y1x1*v*+2a1y1*+γ[y2*x2(t-τ2)v(t-τ2)y2x2*v*e-m2τ2(d2x2*+a2y2*em2τ2-d2x2)-e-m2τ2x2*x2(d2x2*+a2y2*em2τ2-d2x2)+e-m2τ2β2x2*v*(vv*)-a2y2*y2*x2(t-τ2)v(t-τ2)y2x2*v*+2a2y2*]-a1y1*v*y1(t-ω1)vy1*-γa2y2*v*y2(t-ω2)vy2*-a1cp1en1ω1v*(vv*)+a1y1*ln⁡(x1(t-τ1)v(t-τ1)x1v)+γa2y2*ln⁡(x2(t-τ2)v(t-τ2)x2v)+a1y1*ln⁡(y1(t-ω1)y1)+γa2y2*ln⁡(y2(t-ω2)y2)=e-m1τ1d1x1*(2-x1x1*-x1*x1)-a1y1*x1*x1-a1y1*y1*x1(t-τ1)v(t-τ1)y1x1*v*-a1y1*v*y1(t-ω1)vy1*+3a1y1*+a1y1*ln⁡(x1(t-τ1)v(t-τ1)y1(t-ω1)x1vy1)+γ[e-m2τ2d2x2*(2-x2x2*-x2*x2)-a2y2*x2*x2-a2y2*y2*x2(t-τ2)v(t-τ2)y2x2*v*-a2y2*v*y2(t-ω2)vy2*+3a2y2*+a2y2*ln⁡(x2(t-τ2)v(t-τ2)y2(t-ω2)x2vy2)]+a1p1en1ω1(e-n1ω1p1y1*+e-n2ω2p2y2*-cv*)vv*.
From (2.25), we see that the last term in (2.26) vanishes. Then, using the following equalities:
(2.27)ln⁡(x1(t-τ1)v(t-τ1)y1(t-ω1)x1vy1)=ln⁡(x1*x1)+ln⁡(y1(t-ω1)v*y1*v)+ln⁡(y1*x1(t-τ1)v(t-τ1)y1x1*v*),ln⁡(x2(t-τ2)v(t-τ2)y2(t-ω2)x2vy2)=ln⁡(x2*x2)+ln⁡(y2(t-ω2)v*y2*v)+ln⁡(y2*x2(t-τ2)v(t-τ2)y2x2*v*),
we can rewrite (2.26) as
(2.28)dW2dt=e-m1τ1d1x1*(2-x1x1*-x1*x1)+e-m2τ2γd2x2*(2-x2x2*-x2*x2)-a1y1*[H(x1*x1)+H(y1(t-ω1)v*y1*v)+H(y1*x1(t-τ1)v(t-τ1)y1x1*v*)]-γa2y2*[H(x2*x2)+H(y2(t-ω2)v*y2*v)+H(y2*x2(t-τ2)v(t-τ2)y2x2*v*)].
Since the arithmetical mean is greater than or equal to the geometrical mean, then the first two terms of (2.28) are less than or equal to zero. It is easy to see that if x1*,y1*,x2*,y2*,v*>0, then dW2/dt≤0. By [32,Theorem  5.3.1], the solutions of system (2.1)–(2.5) limit to M, the largest invariant subset of {dW2/dt=0}. It can be seen that dW2/dt=0 if and only if x1=x1*,x2=x2*,v=v*, and H=0, that is,
(2.29)y1(t-ω1)v*y1*v=y2(t-ω2)v*y2*v=y1*x1(t-τ1)v(t-τ1)y1x1*v*=y2*x2(t-τ2)v(t-τ2)y2x2*v*=1.
If v=v*, then from (2.29), we have y1=y1* and y2=y2*, and hence dW2/dt equal to zero at E1. LaSalle’s invariance principle implies global stability of E1.


## 3. Basic Virus Dynamics Model with Multitarget Cells and Delays

In this section, we propose a virus dynamics model which describes the interaction of the virus withn classes of target cells. Two types of discrete-time delays (τi,ωi,i=1,…,n) are incorporated into the model. The model is a generalization of those of one class of target cells and two classes of target cells models presented, respectively, in [26, 33]. Moreover, it can be seen that when n=1 and ω1=0, then the following model leads to the model presented in [11]. (3.1)ẋi=λi-dixi-βixiv,i=1,…,n,ẏi=e-miτiβixi(t-τi)v(t-τi)-aiyi,i=1,…,n,v̇=∑i=1ne-niωipiyi(t-ωi)-cv,
where xi and yi represent the populations of the uninfected target cells and infected cells of class i, respectively, v is the population of the virus particles. All the parameters of the model have the same biological meaning as given in the previous section.The initial conditions for system (3.1) take the form(3.2)xj(θ)=φj(θ),j=1,…,n,yj(θ)=ϕj+n(θ),j=1,…,n,v(θ)=φ2n+1(θ),φj(θ)≥0,j=1,…,2n+1,θ∈[-max⁡{τ1,…,τn,ω1,…,ωn},0],
where (φ1(θ),φ2(θ),…,φ2n+1(θ))∈C and C=C([-max{τ1,…,τn,ω1,…,ωn},0],ℝ+2n+1) is the Banach space of continuous functions mapping the interval [-max{τ1,…,τn,ω1,…,ωn},0] into ℝ+2n+1.Similar to the previous section, the nonnegativity and the boundedness of the solutions of system (3.1) can be shown.


### 3.1. Steady States

It is clear that system (3.1) has an uninfected steady state E0=(x10,…,xn0,y10,…,yn0,v0), where xi0=λi/di, yi0=0, i=1,…,n, and v0=0. The system can also have a positive infected steady state E1(x1*,…,xn*,y1*,…,yn*,v*). The coordinates of the infected steady state, if they exist, satisfy the equalities:(3.3)λi=dixi*+βixi*v*,i=1,…,n,(3.4)aiyi*=e-miτiβixi*v*,i=1,…,n,(3.5)cv*=∑i=1ne-niωipiyi*.
The basic reproduction number of system (3.1) is given by(3.6)R0=∑i=1nRi=∑i=1ne-(miτi+niωi)βipiλiaidic,
where Ri is the basic reproduction number for the dynamics of the interaction of the virus only with the target cells of class i.


### 3.2. Global Stability

In the following theorem, the global stability of the uninfected and infected steady states of system (3.1) will be established.Theorem 3.1.
(i) IfR0≤1, then E0 is GAS for any τi,ωi≥0, i=1,…,n.
(ii) IfE1 exists, then it is GAS for any τi,ωi≥0, i=1,…,n.Proof.
(i) Define a Lyapunov functionalW1 as follows:
(3.7)W1=∑i=1nγi[e-miτixi0H(xixi0)+yi+e-miτiβi∫0τixi(t-θ)v(t-θ)dθ+ai∫0ωiyi(t-θ)dθ]+a1p1en1ω1v,
where γi=(a1pi/aip1)en1ω1-niωi. The time derivative of W1 along the solution of system (3.1) satisfies
(3.8)dW1dt=∑i=1nγi[e-miτi(1-xi0xi)(λi-dixi-βixiv)+e-miτiβixi(t-τi)v(t-τi)-aiyi+e-miτi(βixiv-βixi(t-τi)v(t-τi))+ai(yi-yi(t-ωi))(1-xi0xi)]+a1p1en1ω1(∑i=1ne-niωipiyi(t-ωi)-cv)=∑i=1ne-miτiγiλi[2-xixi0-xi0xi]+a1cp1en1ω1(∑i=1ne-(miτi+niωi)piβixi0aic-1)v=∑i=1ne-miτiγiλi[2-xixi0-xi0xi]+a1cp1en1ω1(R0-1)v.
Since the arithmetical mean is greater than or equal to the geometrical mean, then the first term of (3.8) is less than or equal to zero. Therefore, if R0≤1, then dW1/dt≤0 for all xi,yi,v>0. Similar to the previous section, one can show that the maximal compact invariant set in {dW1/dt=0} is the singleton {E0} when R0≤1. The global stability of E0 follows from LaSalle’s invariance principle.
To prove (ii), we consider the Lyapunov functional(3.9)W2=∑i=1nγi[e-miτixi*H(xixi*)+yi*H(yiyi*)+e-miτiβixi*v*∫0τiH(xi(t-θ)v(t-θ)xi*v*)dθ+aiyi*∫0ωiH(yi(t-θ)yi*)dθ]+a1p1en1ω1v*H(vv*).
Differentiating with respect to time yields
(3.10)dW2dt=∑i=1nγi[e-miτi(1-xi*xi)(λi-dixi-βixiv)+(1-yi*yi)(e-miτiβixi(t-τi)v(t-τi)-aiyi)+e-miτi(βixiv-βixi(t-τi)v(t-τi))+e-miτiβixi*v*ln⁡(xi(t-τi)v(t-τi)xiv)+ai(yi-yi(t-ωi))+aiyi*ln⁡(yi(t-ωi)yi)]+a1en1ω1p1(1-v*v)(∑i=1ne-niωipiyi(t-ωi)-cv)=∑i=1nγi[e-miτi(λi-dixi-λixi*xi+dixi*+βixi*v)-e-miτiβixi(t-τi)v(t-τi)yi*yi+aiyi*+e-miτiβixi*v*ln⁡(xi(t-τi)v(t-τi)xiv)+aiyi*ln⁡(yi(t-ωi)yi)]-a1cen1ω1p1v-v*v∑i=1nγiaiyi(t-ωi)+a1cen1ω1p1v*.
Using the infected steady state conditions (3.3)–(3.5), and the following equality:
(3.11)a1cen1ω1p1v*=a1en1ω1p1∑i=1ne-niωipiyi*=∑i=1nγiaiyi*,
we obtain
(3.12)dW2dt=∑i=1nγi[e-miτi(dixi*+emiτiaiyi*-dixi-xi*xi(dixi*+emiτiaiyi*)+dixi*)-aiyi*yi*xi(t-τi)v(t-τi)yixi*v*+2aiyi*-aiyi*v*yi(t-ωi)vyi*+aiyi*ln⁡(xi(t-τi)v(t-τi)xiv)+aiyi*ln⁡(yi(t-ωi)yi)]+(∑i=1ne-miτiγiβixi*-a1cen1ω1p1)v=∑i=1nγi[e-miτidixi*(2-xi*xi-xixi*)-aiyi*xi*xi-aiyi*v*yi(t-ωi)vyi*-aiyi*yi*xi(t-τi)v(t-τi)yixi*v*+3aiyi*+aiyi*ln⁡(xi(t-τi)v(t-τi)yi(t-ωi)xivyi)]+a1en1ω1p1(∑i=1ne-niωipiyi*-cv*)vv*.
From (3.5), we can see that the last term in (3.12) vanishes. Then, by using the following equality:
(3.13)ln⁡(xi(t-τi)v(t-τi)yi(t-ωi)xivyi)=ln⁡(xi*xi)+ln⁡(v*yi(t-ωi)vyi*)+ln⁡(yi*xi(t-τi)v(t-τi)yixi*v*),
we can rewrite (3.12) as
(3.14)dW2dt=∑i=1nγi[e-miτidixi*(2-xi*xi-xixi*)-aiyi*(H(xi*xi)+H(v*yi(t-ωi)vyi*)+H(yi*xi(t-τi)v(t-τi)yixi*v*))].
It is easy to see that if xi*,yi*,v*>0, i=1,…,n, then dW2/dt≤0 for all (xi,yi,v)>0 (the arithmetical mean is greater than or equal to the geometrical mean and H≥0). Clearly, the singleton {E1} is the only invariant set in {dW2/dt=0}. LaSalle's invariance principle implies global stability of E1.


## 3.1. Steady States

It is clear that system (3.1) has an uninfected steady state E0=(x10,…,xn0,y10,…,yn0,v0), where xi0=λi/di, yi0=0, i=1,…,n, and v0=0. The system can also have a positive infected steady state E1(x1*,…,xn*,y1*,…,yn*,v*). The coordinates of the infected steady state, if they exist, satisfy the equalities:(3.3)λi=dixi*+βixi*v*,i=1,…,n,(3.4)aiyi*=e-miτiβixi*v*,i=1,…,n,(3.5)cv*=∑i=1ne-niωipiyi*.
The basic reproduction number of system (3.1) is given by(3.6)R0=∑i=1nRi=∑i=1ne-(miτi+niωi)βipiλiaidic,
where Ri is the basic reproduction number for the dynamics of the interaction of the virus only with the target cells of class i.


## 3.2. Global Stability

In the following theorem, the global stability of the uninfected and infected steady states of system (3.1) will be established.Theorem 3.1.
(i) IfR0≤1, then E0 is GAS for any τi,ωi≥0, i=1,…,n.
(ii) IfE1 exists, then it is GAS for any τi,ωi≥0, i=1,…,n.Proof.
(i) Define a Lyapunov functionalW1 as follows:
(3.7)W1=∑i=1nγi[e-miτixi0H(xixi0)+yi+e-miτiβi∫0τixi(t-θ)v(t-θ)dθ+ai∫0ωiyi(t-θ)dθ]+a1p1en1ω1v,
where γi=(a1pi/aip1)en1ω1-niωi. The time derivative of W1 along the solution of system (3.1) satisfies
(3.8)dW1dt=∑i=1nγi[e-miτi(1-xi0xi)(λi-dixi-βixiv)+e-miτiβixi(t-τi)v(t-τi)-aiyi+e-miτi(βixiv-βixi(t-τi)v(t-τi))+ai(yi-yi(t-ωi))(1-xi0xi)]+a1p1en1ω1(∑i=1ne-niωipiyi(t-ωi)-cv)=∑i=1ne-miτiγiλi[2-xixi0-xi0xi]+a1cp1en1ω1(∑i=1ne-(miτi+niωi)piβixi0aic-1)v=∑i=1ne-miτiγiλi[2-xixi0-xi0xi]+a1cp1en1ω1(R0-1)v.
Since the arithmetical mean is greater than or equal to the geometrical mean, then the first term of (3.8) is less than or equal to zero. Therefore, if R0≤1, then dW1/dt≤0 for all xi,yi,v>0. Similar to the previous section, one can show that the maximal compact invariant set in {dW1/dt=0} is the singleton {E0} when R0≤1. The global stability of E0 follows from LaSalle’s invariance principle.
To prove (ii), we consider the Lyapunov functional(3.9)W2=∑i=1nγi[e-miτixi*H(xixi*)+yi*H(yiyi*)+e-miτiβixi*v*∫0τiH(xi(t-θ)v(t-θ)xi*v*)dθ+aiyi*∫0ωiH(yi(t-θ)yi*)dθ]+a1p1en1ω1v*H(vv*).
Differentiating with respect to time yields
(3.10)dW2dt=∑i=1nγi[e-miτi(1-xi*xi)(λi-dixi-βixiv)+(1-yi*yi)(e-miτiβixi(t-τi)v(t-τi)-aiyi)+e-miτi(βixiv-βixi(t-τi)v(t-τi))+e-miτiβixi*v*ln⁡(xi(t-τi)v(t-τi)xiv)+ai(yi-yi(t-ωi))+aiyi*ln⁡(yi(t-ωi)yi)]+a1en1ω1p1(1-v*v)(∑i=1ne-niωipiyi(t-ωi)-cv)=∑i=1nγi[e-miτi(λi-dixi-λixi*xi+dixi*+βixi*v)-e-miτiβixi(t-τi)v(t-τi)yi*yi+aiyi*+e-miτiβixi*v*ln⁡(xi(t-τi)v(t-τi)xiv)+aiyi*ln⁡(yi(t-ωi)yi)]-a1cen1ω1p1v-v*v∑i=1nγiaiyi(t-ωi)+a1cen1ω1p1v*.
Using the infected steady state conditions (3.3)–(3.5), and the following equality:
(3.11)a1cen1ω1p1v*=a1en1ω1p1∑i=1ne-niωipiyi*=∑i=1nγiaiyi*,
we obtain
(3.12)dW2dt=∑i=1nγi[e-miτi(dixi*+emiτiaiyi*-dixi-xi*xi(dixi*+emiτiaiyi*)+dixi*)-aiyi*yi*xi(t-τi)v(t-τi)yixi*v*+2aiyi*-aiyi*v*yi(t-ωi)vyi*+aiyi*ln⁡(xi(t-τi)v(t-τi)xiv)+aiyi*ln⁡(yi(t-ωi)yi)]+(∑i=1ne-miτiγiβixi*-a1cen1ω1p1)v=∑i=1nγi[e-miτidixi*(2-xi*xi-xixi*)-aiyi*xi*xi-aiyi*v*yi(t-ωi)vyi*-aiyi*yi*xi(t-τi)v(t-τi)yixi*v*+3aiyi*+aiyi*ln⁡(xi(t-τi)v(t-τi)yi(t-ωi)xivyi)]+a1en1ω1p1(∑i=1ne-niωipiyi*-cv*)vv*.
From (3.5), we can see that the last term in (3.12) vanishes. Then, by using the following equality:
(3.13)ln⁡(xi(t-τi)v(t-τi)yi(t-ωi)xivyi)=ln⁡(xi*xi)+ln⁡(v*yi(t-ωi)vyi*)+ln⁡(yi*xi(t-τi)v(t-τi)yixi*v*),
we can rewrite (3.12) as
(3.14)dW2dt=∑i=1nγi[e-miτidixi*(2-xi*xi-xixi*)-aiyi*(H(xi*xi)+H(v*yi(t-ωi)vyi*)+H(yi*xi(t-τi)v(t-τi)yixi*v*))].
It is easy to see that if xi*,yi*,v*>0, i=1,…,n, then dW2/dt≤0 for all (xi,yi,v)>0 (the arithmetical mean is greater than or equal to the geometrical mean and H≥0). Clearly, the singleton {E1} is the only invariant set in {dW2/dt=0}. LaSalle's invariance principle implies global stability of E1.


## 4. Virus Dynamics Model with Saturation Infection Rate

In this section, we proposed a virus dynamics model which describes the interaction of the virus withn classes of target cells taking into account the saturation infection rate and multiple intracellular delays: (4.1)ẋi(t)=λi-dixi(t)-βixi(t)v(t)1+αiv(t),i=1,…,n,ẏi(t)=e-mtτiβixi(t-τi)v(t-τi)1+αiv(t)-aiyi(t),i=1,…,n,v̇(t)=∑i=1ne-niωipiyi(t-ωi)-cv(t),
where αi,i=1,…,n are positive constants. The variables and parameters of the model have the same definitions as given in Section 2. We mention that if n=1 and ω1=0, then model (4.1) leads to the model presented in [9], and if n=2 and ω1=ω2=0, α1=α2=1, then model (4.1) leads to the model presented in [34].


### 4.1. Steady States

It is clear that system (4.1) has an uninfected steady state E0=(x10,…,xn0,y10,…,yn0,v0), where xi0=λi/di, yi0=0, and v0=0. The system can also have a positive infected steady state E1(x1*,…,xn*,y1*,…,yn*,v*). The coordinates of the infected steady state, if they exist, satisfy the equalities:(4.2)λi=dixi*+βixi*v*1+αiv*,i=1,...,n,aiyi*=e-miτiβixi*v*1+αiv*,i=1,...,n,cv*=∑i=1ne-niωipiyi*.
The basic reproduction number R0 for system (4.1) is the same as given by (3.6).


### 4.2. Global Stability

In this section, we study the global stability of the uninfected and infected steady states of system (4.1).Theorem 4.1.
(i) IfR0≤1, then E0 is GAS for any τi,ωi≥0, i=1,…,n.
(ii) IfE1 exists, then it is GAS for any τi,ωi≥0, i=1,…,n.Proof.
(i) Define a Lyapunov functionalW1 as follows:
(4.3)W1=∑i=1nγi[e-miτixi0H(xixi0)+yi+e-miτiβi∫0τixi(t-θ)v(t-θ)1+αiv(t-θ)dθ+ai∫0ωiyi(t-θ)dθ]+a1en1ω1p1v,
The time derivative of W1 along the trajectories of (4.1) satisfies
(4.4)dW1dt=∑i=1nγi[e-miτi(1-xi0xi)(λi-dixi-βixiv1+αiv)+e-miτiβixi(t-τi)v(t-τi)1+αiv(t-τi)-aiyi+e-miτi(βixiv1+αiv-βixi(t-τi)v(t-τi)1+αiv(t-τi))+ai(yi-yi(t-ωi))(1-xi0xi)]+a1en1ω1p1(∑i=1ne-niωipiyi(t-ωi)-cv)=∑i=1nγie-miτi[λi-dixi-λixi0xi+dixi0+βixi0v1+αiv]-a1cen1ω1p1v=∑i=1nγie-miτiλi[2-xixi0-xi0xi]-a1cen1ω1p1v+a1cen1ω1p1∑i=1ne-(miτi+niωi)piβixi0vaic(1+αiv)=∑i=1nγie-miτiλi[2-xixi0-xi0xi]-a1cen1ω1p1v+a1cen1ω1p1∑i=1nRiv1+αiv=∑i=1nγie-miτiλi[2-xixi0-xi0xi]+a1cen1ω1p1(R0-1)v-a1cen1ω1p1∑i=1nRiαiv21+αiv.
It is clear that if R0≤1, then dW1/dt≤0 for all xi,yi,v>0, where equality occurs at E0. The global stability of E0 follows from LaSalle’s invariance principle.
To prove (ii), we consider the Lyapunov functional:(4.5)W2=∑i=1nγi[e-miτixi*H(xixi*)+yi*H(yiyi*)+e-miτiβixi*v*1+αiv*∫0τiH(xi(t-θ)v(t-θ)(1+αiv*)xi*v*(1+αiv(t-θ)))dθ+aiyi*∫0ωiH(yi(t-θ)yi*)dθ]+a1en1ω1p1v*H(vv*).
Differentiating with respect to time yields
(4.6)dW2dt=∑i=1nγi[e-miτi(1-xi*xi)(λi-dixi-βixiv1+αiv)+(1-yi*yi)(e-miτiβixi(t-τi)v(t-τi)1+αiv(t-τi)-aiyi)+e-miτi(+βixi*v*1+αiv*ln⁡(xi(t-τi)v(t-τi)(1+αiv)xiv(1+αiv(t-τi))))βixiv1+αiv-βixi(t-τi)v(t-τi)1+αiv(t-τi)+βixi*v*1+αiv*ln⁡(xi(t-τi)v(t-τi)(1+αiv)xiv(1+αiv(t-τi))))+ai(yi-yi(t-ωi))+aiyi*ln⁡(yi(t-ωi)yi)]+a1en1ω1p1(1-v*v)(∑i=1ne-niωipiyi(t-ωi)-cv)=∑i=1nγi[e-miτi(λi-dixi-λixi*xi+dixi*+βixi*v1+αiv)-e-miτiβixi(t-τi)v(t-τi)1+αiv(t-τi)yi*yi+aiyi*+e-miτiβixi*v*1+αiv*ln⁡(xi(t-τi)v(t-τi)(1+αiv)xiv(1+αiv(t-τi)))+aiyi*ln⁡(yi(t-ωi)yi)]-a1cen1ω1p1v-a1en1ω1p1v*v∑i=1ne-niωipiyi(t-ωi)+a1cen1ω1p1v*.
Using the infected steady state conditions (4.2), and the following equality:
(4.7)a1cen1ω1p1v=a1cen1ω1p1v*vv*=a1en1ω1p1vv*∑i=1ne-niωipiyi*=vv*∑i=1nγiaiyi*,
we obtain
(4.8)dW2dt=∑i=1nγi[e-miτi(dixi*+emiτiaiyi*-dixi-xi*xi(dixi*+emiτiaiyi*)+dixi*)+aiyi*v(1+αiv*)v*(1+αiv)-aiyi*yi*xi(t-τi)v(t-τi)(1+αiv*)yixi*v*(1+αiv(t-τi))+2aiyi*-aiyi*v*yi(t-ωi)vyi*+aiyi*ln⁡(xi(t-τi)v(t-τi)yi(t-ωi)(1+αiv)xivyi(1+αiv(t-τi)))-aiyi*vv*]=∑i=1nγi[e-miτidixi*(2-xi*xi-xixi*)-aiyi*xi*xi+3aiyi*+aiyi*(v(1+αiv*)v*(1+αiv)-vv*)-aiyi*yi*xi(t-τi)v(t-τi)(1+αiv*)yixi*v*(1+αiv(t-τi))-aiyi*v*yi(t-ωi)vyi*+aiyi*ln⁡(xi(t-τi)v(t-τi)yi(t-ωi)(1+αiv)xivyi(1+αiv(t-τi)))]=∑i=1nγi[e-miτidixi*(2-xi*xi-xixi*)-aiyi*xi*xi+4aiyi*-aiyi*yi*xi(t-τi)v(t-τi)(1+αiv*)yixi*v*(1+αiv(t-τi))+aiyi*(-1+v(1+αiv*)v*(1+αiv)-vv*+1+αiv1+αiv*-1+αiv1+αiv*)-aiyi*v*yi(t-ωi)vyi*+aiyi*ln⁡(xi(t-τi)v(t-τi)yi(t-ωi)(1+αiv)xivyi(1+αiv(t-τi)))(2-xi*xi-xixi*)].
Then using the following equalities:
(4.9)ln⁡(xi(t-τi)v(t-τi)yi(t-ωi)(1+αiv)xivyi(1+αiv(t-τi)))=ln⁡(xi*xi)+ln⁡(v*yi(t-ωi)vyi*)+ln⁡(1+αiv1+αiv*)+ln⁡(yi*xi(t-τi)v(t-τi)(1+αiv*)yixi*v*(1+αiv(t-τi))),-1+v(1+αiv*)v*(1+αiv)-vv*+1+αiv1+αiv*=-αi(v-v*)2v*(1+αiv*)(1+αiv),
we obtain
(4.10)dW2dt=∑i=1nγi[e-miτidixi*(2-xi*xi-xixi*)-aiyi*αi(v-v*)2v*(1+αiv*)(1+αiv)-aiyi*((yi*xi(t-τi)v(t-τi)(1+αiv*)yixi*v*(1+αiv(t-τi)))+H(yi*xi(t-τi)v(t-τi)(1+αiv*)yixi*v*(1+αiv(t-τi))))H(xi*xi)+H(v*yi(t-ωi)vyi*)+H(1+αiv1+αiv*)+H(yi*xi(t-τi)v(t-τi)(1+αiv*)yixi*v*(1+αiv(t-τi))))].
It is easy to see that if xi*,yi*,v*>0,i=1,…,n, then dW2/dt≤0 for all (xi,yi,v)>0. Clearly, the singleton {E1} is the only invariant set in {dW2/dt=0}. LaSalle’s invariance principle implies global stability of E1.


## 4.1. Steady States

It is clear that system (4.1) has an uninfected steady state E0=(x10,…,xn0,y10,…,yn0,v0), where xi0=λi/di, yi0=0, and v0=0. The system can also have a positive infected steady state E1(x1*,…,xn*,y1*,…,yn*,v*). The coordinates of the infected steady state, if they exist, satisfy the equalities:(4.2)λi=dixi*+βixi*v*1+αiv*,i=1,...,n,aiyi*=e-miτiβixi*v*1+αiv*,i=1,...,n,cv*=∑i=1ne-niωipiyi*.
The basic reproduction number R0 for system (4.1) is the same as given by (3.6).


## 4.2. Global Stability

In this section, we study the global stability of the uninfected and infected steady states of system (4.1).Theorem 4.1.
(i) IfR0≤1, then E0 is GAS for any τi,ωi≥0, i=1,…,n.
(ii) IfE1 exists, then it is GAS for any τi,ωi≥0, i=1,…,n.Proof.
(i) Define a Lyapunov functionalW1 as follows:
(4.3)W1=∑i=1nγi[e-miτixi0H(xixi0)+yi+e-miτiβi∫0τixi(t-θ)v(t-θ)1+αiv(t-θ)dθ+ai∫0ωiyi(t-θ)dθ]+a1en1ω1p1v,
The time derivative of W1 along the trajectories of (4.1) satisfies
(4.4)dW1dt=∑i=1nγi[e-miτi(1-xi0xi)(λi-dixi-βixiv1+αiv)+e-miτiβixi(t-τi)v(t-τi)1+αiv(t-τi)-aiyi+e-miτi(βixiv1+αiv-βixi(t-τi)v(t-τi)1+αiv(t-τi))+ai(yi-yi(t-ωi))(1-xi0xi)]+a1en1ω1p1(∑i=1ne-niωipiyi(t-ωi)-cv)=∑i=1nγie-miτi[λi-dixi-λixi0xi+dixi0+βixi0v1+αiv]-a1cen1ω1p1v=∑i=1nγie-miτiλi[2-xixi0-xi0xi]-a1cen1ω1p1v+a1cen1ω1p1∑i=1ne-(miτi+niωi)piβixi0vaic(1+αiv)=∑i=1nγie-miτiλi[2-xixi0-xi0xi]-a1cen1ω1p1v+a1cen1ω1p1∑i=1nRiv1+αiv=∑i=1nγie-miτiλi[2-xixi0-xi0xi]+a1cen1ω1p1(R0-1)v-a1cen1ω1p1∑i=1nRiαiv21+αiv.
It is clear that if R0≤1, then dW1/dt≤0 for all xi,yi,v>0, where equality occurs at E0. The global stability of E0 follows from LaSalle’s invariance principle.
To prove (ii), we consider the Lyapunov functional:(4.5)W2=∑i=1nγi[e-miτixi*H(xixi*)+yi*H(yiyi*)+e-miτiβixi*v*1+αiv*∫0τiH(xi(t-θ)v(t-θ)(1+αiv*)xi*v*(1+αiv(t-θ)))dθ+aiyi*∫0ωiH(yi(t-θ)yi*)dθ]+a1en1ω1p1v*H(vv*).
Differentiating with respect to time yields
(4.6)dW2dt=∑i=1nγi[e-miτi(1-xi*xi)(λi-dixi-βixiv1+αiv)+(1-yi*yi)(e-miτiβixi(t-τi)v(t-τi)1+αiv(t-τi)-aiyi)+e-miτi(+βixi*v*1+αiv*ln⁡(xi(t-τi)v(t-τi)(1+αiv)xiv(1+αiv(t-τi))))βixiv1+αiv-βixi(t-τi)v(t-τi)1+αiv(t-τi)+βixi*v*1+αiv*ln⁡(xi(t-τi)v(t-τi)(1+αiv)xiv(1+αiv(t-τi))))+ai(yi-yi(t-ωi))+aiyi*ln⁡(yi(t-ωi)yi)]+a1en1ω1p1(1-v*v)(∑i=1ne-niωipiyi(t-ωi)-cv)=∑i=1nγi[e-miτi(λi-dixi-λixi*xi+dixi*+βixi*v1+αiv)-e-miτiβixi(t-τi)v(t-τi)1+αiv(t-τi)yi*yi+aiyi*+e-miτiβixi*v*1+αiv*ln⁡(xi(t-τi)v(t-τi)(1+αiv)xiv(1+αiv(t-τi)))+aiyi*ln⁡(yi(t-ωi)yi)]-a1cen1ω1p1v-a1en1ω1p1v*v∑i=1ne-niωipiyi(t-ωi)+a1cen1ω1p1v*.
Using the infected steady state conditions (4.2), and the following equality:
(4.7)a1cen1ω1p1v=a1cen1ω1p1v*vv*=a1en1ω1p1vv*∑i=1ne-niωipiyi*=vv*∑i=1nγiaiyi*,
we obtain
(4.8)dW2dt=∑i=1nγi[e-miτi(dixi*+emiτiaiyi*-dixi-xi*xi(dixi*+emiτiaiyi*)+dixi*)+aiyi*v(1+αiv*)v*(1+αiv)-aiyi*yi*xi(t-τi)v(t-τi)(1+αiv*)yixi*v*(1+αiv(t-τi))+2aiyi*-aiyi*v*yi(t-ωi)vyi*+aiyi*ln⁡(xi(t-τi)v(t-τi)yi(t-ωi)(1+αiv)xivyi(1+αiv(t-τi)))-aiyi*vv*]=∑i=1nγi[e-miτidixi*(2-xi*xi-xixi*)-aiyi*xi*xi+3aiyi*+aiyi*(v(1+αiv*)v*(1+αiv)-vv*)-aiyi*yi*xi(t-τi)v(t-τi)(1+αiv*)yixi*v*(1+αiv(t-τi))-aiyi*v*yi(t-ωi)vyi*+aiyi*ln⁡(xi(t-τi)v(t-τi)yi(t-ωi)(1+αiv)xivyi(1+αiv(t-τi)))]=∑i=1nγi[e-miτidixi*(2-xi*xi-xixi*)-aiyi*xi*xi+4aiyi*-aiyi*yi*xi(t-τi)v(t-τi)(1+αiv*)yixi*v*(1+αiv(t-τi))+aiyi*(-1+v(1+αiv*)v*(1+αiv)-vv*+1+αiv1+αiv*-1+αiv1+αiv*)-aiyi*v*yi(t-ωi)vyi*+aiyi*ln⁡(xi(t-τi)v(t-τi)yi(t-ωi)(1+αiv)xivyi(1+αiv(t-τi)))(2-xi*xi-xixi*)].
Then using the following equalities:
(4.9)ln⁡(xi(t-τi)v(t-τi)yi(t-ωi)(1+αiv)xivyi(1+αiv(t-τi)))=ln⁡(xi*xi)+ln⁡(v*yi(t-ωi)vyi*)+ln⁡(1+αiv1+αiv*)+ln⁡(yi*xi(t-τi)v(t-τi)(1+αiv*)yixi*v*(1+αiv(t-τi))),-1+v(1+αiv*)v*(1+αiv)-vv*+1+αiv1+αiv*=-αi(v-v*)2v*(1+αiv*)(1+αiv),
we obtain
(4.10)dW2dt=∑i=1nγi[e-miτidixi*(2-xi*xi-xixi*)-aiyi*αi(v-v*)2v*(1+αiv*)(1+αiv)-aiyi*((yi*xi(t-τi)v(t-τi)(1+αiv*)yixi*v*(1+αiv(t-τi)))+H(yi*xi(t-τi)v(t-τi)(1+αiv*)yixi*v*(1+αiv(t-τi))))H(xi*xi)+H(v*yi(t-ωi)vyi*)+H(1+αiv1+αiv*)+H(yi*xi(t-τi)v(t-τi)(1+αiv*)yixi*v*(1+αiv(t-τi))))].
It is easy to see that if xi*,yi*,v*>0,i=1,…,n, then dW2/dt≤0 for all (xi,yi,v)>0. Clearly, the singleton {E1} is the only invariant set in {dW2/dt=0}. LaSalle’s invariance principle implies global stability of E1.


## 5. Conclusion

In this paper, we have studied the global properties of a class of virus dynamics models with multitarget cells and multiple delays. First, we have introduced a model with two classes of target cells (CD4+ T and macrophages in case of HIV). Then, we have proposed a model describing the interaction of the virus with n classes of target cells. A model with multitarget cells taking into account the saturation infection rate is also studied. Two types of discrete time delays have been incorporated into these models to take into account (i) the latent period between the time the target cell is contacted by the virus particle and the time virus enters the cell, (ii) the latent period between the time the virus has penetrated into a cell and the time of the emission of infectious (mature) virus particles. The global stability of the uninfected and infected steady states has been established by using suitable Lyapunov functionals and LaSalle invariant principle. We have proven that, if the basic reproduction number R0 is less than unity, then the uninfected steady state is GAS and if R0>1 (or the infected steady state exists) then the infected steady state is GAS for all time delays.



---

*Source: 201274-2011-12-20.xml*