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

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 (2𝑛+1)-dimensional system of nonlinear DDEs that describes the dynamics of the virus, 𝑛 classes of uninfected target cells, and 𝑛 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 𝑅₀ is less than unity, then the uninfected steady state is globally asymptotically stable, and if 𝑅₀>1 (or if the infected steady state exists), then the infected steady state is globally asymptotically stable.

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)𝑣̇(𝑡𝛽)𝑥=(𝑡𝑝−𝑦(𝜏𝑡))𝑣−(𝑡𝑐𝑣−(𝑡𝜏)),−𝑎𝑦(𝑡),((11..23)) where 𝑥(𝑡),𝑦(𝑡), and 𝑣(𝑡) represent the populations of uninfected target cells, infected cells, and free virus particles at time 𝑡, respectively. Here, 𝜆 represents the rate of which new target cells are generated from sources within the body, 𝑑 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 𝑎. The virus particles are produced by the infected cells with rate constant 𝑝, and are removed from the system with rate constant 𝑐. 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 𝑡 is given by the number of cells that were newly infected cells at time 𝑡−𝜏 and are still alive at time 𝑡. The probability of surviving the time period from 𝑡−𝜏 to 𝑡 is 𝑒^−𝑚𝜏, where 𝑚 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 attacks 𝑛 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 𝑅₀. If 𝑅₀≤1, then the uninfected steady state is globally asymptotically stable (GAS). If 𝑅₀>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)(𝑡)=𝑥̇₂𝑒^{−𝑚1𝜏1}𝛽₁𝑥₁𝑡−𝜏₁𝑣𝑡−𝜏₁−𝑎₁𝑦₁(𝑡),(2.2)𝑦̇₂(𝑡)=𝜆₂−𝑑₂𝑥₂(𝑡)−𝛽₂𝑥₂(𝑡)𝑣(𝑡),(2.3)𝑣̇(𝑡)=(𝑒𝑡)=_{−𝑛1𝜔1}𝑒^{−𝑚2𝜏2}𝛽₂𝑥₂𝑡−𝜏₂𝑣𝑡−𝜏₂−𝑎₂𝑦₂(𝑡),(2.4)𝑝₁𝑦₁𝑡−𝜔₁+𝑒^{−𝑛2𝜔2}𝑝₂𝑦₂𝑡−𝜔₂−𝑐𝑣(𝑡),(2.5) where 𝑥₁and 𝑥₂ represent the populations of the two classes of uninfected target cells; 𝑦₁and 𝑦₂ are the populations of the infected cells. The population of the target cells are described by (2.1) and (2.3), where 𝜆₁ and 𝜆₂ represent the rates of which new target cells are generated, 𝑑₁ and 𝑑₂ are the death rate constants, and 𝛽₁ and 𝛽₂ 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 𝑎₁ and 𝑎₂. The virus particles are produced by the two classes of infected cells with rate constants 𝑝₁ and 𝑝₂ and are cleared with rate constant 𝑐. Here the parameter 𝜏_𝑖 accounts for the time between the target cells of class 𝑖 are contacted by the virus particle and the contacting virus enters the cells. The recruitment of virus-producing cells at time 𝑡 is given by the number of cells that were newly infected cells at times 𝑡−𝜏_𝑖 and are still alive at time 𝑡. Also, 𝑚_𝑖 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 𝑡−𝜏_𝑖 to 𝑡 is 𝑒^{−𝑚𝑖𝜏𝑖}, 𝑖=1,2. The time between the virus has penetrated into a target cell of class 𝑖 and the emission of infectious (matures) virus particles is represented by 𝜔_𝑖. The probability of survival of an immature virus is given by 𝑒^{−𝑛𝑖𝜔𝑖}, where 𝑛_𝑖 is constant.

2.1. Initial Conditions

The initial conditions for system (2.1)–(2.5) take the form𝑥₁(𝜃)=𝜑₁(𝜃),𝜑𝑦₁(𝜃)=𝜑₂(𝜃),𝑥₂(𝜃)=𝜑₃(𝜃),𝑦₂(𝜃)=𝜑₄(𝜃),𝑣(𝜃)=𝜑₅(𝜃),_𝑖(𝜃)≥0,𝜃∈−max𝜏₁,𝜏₂,𝜔₁,𝜔₂,0,𝑖=1,…,5,(2.6) where (𝜑₁(𝜃),…,𝜑₅(𝜃))∈𝐶([−max{𝜏₁,𝜏₂,𝜔₁,𝜔₂},0],ℝ⁵⁺), the Banach space of continuous functions mapping the interval [−max{𝜏₁,𝜏₂,𝜔₁,𝜔₂},0]into ℝ⁵⁺, where ℝ⁵⁺={(𝑧₁,𝑧₂,…,𝑧₅)∶𝑧_𝑖≥0}.

By the fundamental theory of functional differential equations [32], system (2.1)–(2.5) has a unique solution(𝑥₁(𝑡),𝑦₁(𝑡),𝑥₂(𝑡),𝑦₂(𝑡),𝑣(𝑡)) 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 (𝑥₁(𝑡),𝑦₁(𝑡),𝑥₂(𝑡),𝑦₂(𝑡),𝑣(𝑡)) be any solution of (2.1)–(2.5) satisfying the initial conditions (2.6), then 𝑥₁(𝑡), 𝑦₁(𝑡), 𝑥₂(𝑡), 𝑦₂(𝑡),and 𝑣(𝑡) are all nonnegative for 𝑡≥0 and ultimately bounded.

Proof. From (2.1) and (2.3), we have 𝑥_𝑖(𝑡)=𝑥_𝑖(0)𝑒^{−∫𝑡0(𝑑𝑖+𝛽𝑖𝑣(𝜉))𝑑𝜉}+𝜆_𝑖^𝑡0𝑒^{−∫𝑡𝜂(𝑑𝑖+𝛽𝑖𝑣(𝜉))𝑑𝜉}𝑑𝜂,𝑖=1,2,(2.7) which indicates that 𝑥₁(𝑡)≥0, 𝑥₂(𝑡)≥0 for all 𝑡≥0. Now from (2.2), (2.4), and (2.5), we have 𝑦_𝑖(𝑡)=𝑦_𝑖(0)𝑒^{−𝑎𝑖𝑡}+𝑒^{−𝑚𝑖𝜏𝑖}^𝑡0𝑣(𝑡)=𝑣(0)𝑒_−𝑐𝑡𝛽_𝑖𝑥_𝑖𝜂−𝜏_𝑖𝑣𝜂−𝜏_𝑖𝑒^{−𝑎𝑖(𝑡−𝜂)}𝑑𝜂,𝑖=1,2,+^𝑡0𝑒^{−𝑛1𝜔1}𝑝₁𝑦₁𝜂−𝜔₁+𝑒^{−𝑛2𝜔2}𝑝₂𝑦₂𝜂−𝜔₂𝑒^{−𝑐(𝑡−𝜂)}𝑑𝜂,(2.8) confiming that 𝑦₁(𝑡)≥0, 𝑦₂(𝑡)≥0, 𝑣(𝑡)≥0 for all 𝑡∈[0,max{𝜏₁,𝜏₂,𝜔₁,𝜔₂}]. By a recursive argument, we obtain 𝑦₁(𝑡)≥0, 𝑦₂(𝑡)≥0, 𝑣(𝑡)≥0 for all 𝑡≥0.
To show the boundedness of the solutions, we let 𝑋₁(𝑡)=𝑒^{−𝑚1𝜏1}𝑥₁(𝑡−𝜏₁−𝜔₁)+𝑦₁(𝑡−𝜔₁) and 𝑋₂(𝑡)=𝑒^{−𝑚2𝜏2}𝑥₂(𝑡−𝜏₂−𝜔₂)+𝑦₂(𝑡−𝜔₂), then 𝑋𝑋̇̇₁(𝑡)≤𝜆₁𝑒^{−𝑚1𝜏1}−𝜎₁𝑋₁(𝑡),₂(𝑡)≤𝜆₂𝑒^{−𝑚2𝜏2}−𝜎₂𝑋₂(𝑡),(2.9) where 𝜎₁=min{𝑑₁,𝑎₁} and 𝜎₂=min{𝑑₂,𝑎₂}. Hence, limsup_𝑡→∞𝑋₁(𝑡)≤𝐿₁, and limsup_𝑡→∞𝑋₂(𝑡)≤𝐿₂, where 𝐿₁=𝜆₁𝑒^{−𝑚1𝜏1}/𝜎₁ and 𝐿₂=𝜆₂𝑒^{−𝑚2𝜏2}/𝜎₂. On the other hand, 𝑣̇(𝑡)≤𝑒^{−𝑛1𝜔1}𝑝₁𝐿₁+𝑒^{−𝑛2𝜔2}𝑝₂𝐿₂−𝑐𝑣,(2.10) then limsup_𝑡→∞𝑣(𝑡)≤𝐿₃, where 𝐿₃=(𝑒^{−𝑛1𝜔1}𝑝₁𝐿₁+𝑒^{−𝑛2𝜔2}𝑝₂𝐿₂)/𝑐. It follows that the solution (𝑥₁(𝑡),𝑦₁(𝑡),𝑥₂(𝑡),𝑦₂(𝑡),𝑣(𝑡)) 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𝑅₀=𝑒^{−(𝑚1𝜏1+𝑛1𝜔1)}𝑝₁𝛽₁𝑎₂𝑥⁰¹+𝑒^{−(𝑚2𝜏2+𝑛2𝜔2)}𝑝₂𝛽₂𝑎₁𝑥⁰²𝑎₁𝑎₂𝑐,(2.11) where 𝑥⁰¹=𝜆₁/𝑑₁ and 𝑥⁰²=𝜆₂/𝑑₂. We observe that 𝑅₀ can be written as𝑅₀=𝑅₁+𝑅₂,(2.12) where𝑅₁=𝑒^{−(𝑚1𝜏1+𝑛1𝜔1)}𝑝₁𝛽₁𝜆₁𝑎₁𝑑₁𝑐,𝑅²=𝑒^{−(𝑚2𝜏2+𝑛2𝜔2)}𝑝₂𝛽₂𝜆₂𝑎₂𝑑₂𝑐(2.13) 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 𝑅₀≤1, then system (2.1)–(2.5) has only one steady state 𝐸₀=(𝑥⁰¹,0,𝑥⁰²,0,0) which is called uninfected steady state, and if 𝑅₀>1, then system (2.1)–(2.5) has two steady states 𝐸₀ and infected steady state 𝐸₁=(𝑥^∗1,𝑦^∗1,𝑥^∗2,𝑦^∗2,𝑣^∗). The coordinates of the infected steady state are given by 𝑥^∗1=⎧⎪⎪⎨⎪⎪⎩𝛼⁵𝑐𝛼₁𝛼₅+𝛼₂𝛼₃−𝛼,if𝛼₄=0,₁𝛼₅+𝛼₂𝛼₃−𝛼₄𝑐+𝛼₁𝛼₅+𝛼₂𝛼₃−𝛼₄𝑐²+4𝛼₁𝛼₄𝛼₅𝑐2𝛼₁𝛼₄𝑥,if𝛼₄≠0,_∗2=⎧⎪⎪⎨⎪⎪⎩𝛼³𝑐𝛼₁𝛼₅+𝛼₂𝛼₃𝑐,if𝛼₄=0,𝛼₂+𝛼¹𝛼₅+𝛼₂𝛼₃−𝛼₄𝑐−𝛼₁𝛼₅+𝛼₂𝛼₃−𝛼₄𝑐²+4𝛼₁𝛼₄𝛼₅𝑐2𝛼₂𝛼₄𝑦,if𝛼₄≠0,_∗1=𝑑¹𝑎₁𝑒^𝑚1𝜏1𝑥𝑥^01∗1−1𝑥^∗1,𝑦^∗2=𝑑²𝑎₂𝑒^𝑚2𝜏2𝑥𝑥^02∗2−1𝑥^∗2,𝑣^∗=𝑑¹𝛽₁𝑥𝑥^01∗1−1,(2.14) where 𝛼₁=𝑒^{−(𝑛1𝜔1+𝑚1𝜏1)}𝑝₁𝛽₁𝑎₁,𝛼₂=𝑒^{−(𝑛2𝜔2+𝑚2𝜏2)}𝑝₂𝛽₂𝑎₂𝛼₄,𝛼₃=𝜆₂𝛽₁,=𝛽₁𝑑₂−𝛽₂𝑑₁,𝛼₅=𝜆₁𝛽₂.(2.15)

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: 𝑧=𝑧(𝑡), for any 𝑧∈{𝑥₁,𝑦₁,𝑥₂,𝑦₂,𝑣}. We also define a function 𝐻∶ℝ_>0→ℝ_≥0 as 𝐻(𝑧)=𝑧−1−ln𝑧.(2.16) It is clear that 𝐻(𝑧)≥0 for any 𝑧>0 and 𝐻 has the global minimum 𝐻(1)=0.

Theorem 2.2. (i) If 𝑅₀≤1, then 𝐸₀ is GAS for any 𝜏₁,𝜏₂,𝜔₁,𝜔₂≥0.
(ii) If 𝑅₀>1, then 𝐸₁ is GAS for any 𝜏₁,𝜏₂,𝜔₁,𝜔₂≥0.

Proof. (i) We consider a Lyapunov functional 𝑊₁=𝑊₁₁+𝑊₁₂+𝑊₁₃,(2.17) where 𝑊₁₁=𝑒^{−𝑚1𝜏1}𝑥⁰¹𝐻𝑥¹𝑥⁰¹+𝑦¹+𝛾𝑒^{−𝑚2𝜏2}𝑥⁰²𝐻𝑥²𝑥⁰²+𝑦²+𝑎𝑝¹₁𝑊𝑒^𝑛1𝜔1𝑣,₁₂=𝑒^{−𝑚1𝜏1}^𝜏1₀𝛽₁𝑥₁(𝑡−𝜃)𝑣(𝑡−𝜃)𝑑𝜃+𝛾𝑒^{−𝑚2𝜏2}^𝜏2₀𝑊𝛽₂𝑥₂(𝑡−𝜃)𝑣(𝑡−𝜃)𝑑𝜃,₁₃=𝑎₁^𝜔1₀𝑦₁(𝑡−𝜃)𝑑𝜃+𝛾𝑎₂^𝜔2₀𝑦₂(𝑡−𝜃)𝑑𝜃,(2.18) where 𝛾=(𝑝₂𝑎₁/𝑝₁𝑎₂)𝑒^{𝑛1𝜔1−𝑛2𝜔2}. We note that 𝑊₁ is defined, continuous, and positive definite for all (𝑥₁,𝑦₁,𝑥₂,𝑦₂,𝑣)>0. Also, the global minimum 𝑊₁=0 occurs at the uninfected steady state 𝐸₀.
The time derivatives of 𝑊₁₁ and 𝑊₁₂ are given by 𝑑𝑊₁₁𝑑𝑡=𝑒^{−𝑚1𝜏1}1−𝑥𝑥⁰¹₁𝑥̇¹+𝑦̇₁+𝛾𝑒^{−𝑚2𝜏2}1−𝑥𝑥⁰²₂𝑥̇²+𝑦̇₂+𝑎𝑝¹₁𝑑𝑊₁₂𝑒^𝑛1𝜔1𝑣,̇𝑑𝑡=𝑒^{−𝑚1𝜏1}^𝜏1₀𝑑𝑑𝑡𝛽¹𝑥₁(𝑡−𝜃)𝑣(𝑡−𝜃)𝑑𝜃+𝛾𝑒^{−𝑚2𝜏2}^𝜏2₀𝑑𝑑𝑡𝛽²=−𝑒_{−𝑚1𝜏1}𝑥₂(𝑡−𝜃)𝑣(𝑡−𝜃)𝑑𝜃^𝜏1₀𝑑𝑑𝜃𝛽¹𝑥₁(𝑡−𝜃)𝑣(𝑡−𝜃)𝑑𝜃−𝛾𝑒^{−𝑚2𝜏2}^𝜏2₀𝑑𝑑𝜃𝛽²=𝑒_{−𝑚1𝜏1}𝑥₂(𝑡−𝜃)𝑣(𝑡−𝜃)𝑑𝜃𝛽₁𝑥₁𝑣−𝛽₁𝑥₁𝑡−𝜏₁𝑣𝑡−𝜏₁+𝛾𝑒^{−𝑚2𝜏2}𝛽₂𝑥₂𝑣−𝛽₂𝑥₂𝑡−𝜏₂𝑣𝑡−𝜏₂.(2.19) Similarly, 𝑑𝑊₁₃/𝑑𝑡 is given by 𝑑𝑊₁₃𝑑𝑡=𝑎¹𝑦₁−𝑦₁𝑡−𝜔₁+𝛾𝑎₂𝑦₂−𝑦₂𝑡−𝜔₂.(2.20) It follows that 𝑑𝑊₁𝑑𝑡=𝑒^{−𝑚1𝜏1}1−𝑥𝑥⁰¹₁𝜆¹−𝑑₁𝑥₁−𝛽₁𝑥₁𝑣+𝑒^{−𝑚1𝜏1}𝛽₁𝑥₁𝑡−𝜏₁𝑣𝑡−𝜏₁−𝑎₁𝑦₁+𝛾𝑒^{−𝑚2𝜏2}1−𝑥𝑥⁰²₂𝜆²−𝑑₂𝑥₂−𝛽₂𝑥₂𝑣+𝑒^{−𝑚2𝜏2}𝛽₂𝑥₂𝑡−𝜏₂𝑣𝑡−𝜏₂−𝑎₂𝑦₂+𝑎₁𝑝₁+𝑒𝑒^𝑛1𝜔1_{−𝑚1𝜏1}𝑒^{−𝑛1𝜔1}𝑝₁𝑦₁𝑡−𝜔₁+𝑒^{−𝑛2𝜔2}𝑝₂𝑦₂𝑡−𝜔₂−𝑐𝑣+𝑎₁𝛽₁𝑥₁𝑣−𝛽₁𝑥₁𝑡−𝜏₁𝑣𝑡−𝜏₁+𝛾𝑒^{−𝑚2𝜏2}𝛽₂𝑥₂𝑣−𝛽₂𝑥₂𝑡−𝜏₂𝑣𝑡−𝜏₂=𝑒_{−𝑚1𝜏1}𝑦₁−𝑦₁𝑡−𝜔₁+𝛾𝑎₂𝑦₂−𝑦₂𝑡−𝜔₂𝜆₁2−𝑥𝑥⁰¹₁−𝑥¹𝑥⁰¹+𝑒^{−𝑚2𝜏2}𝛾𝜆₂2−𝑥𝑥⁰²₂−𝑥²𝑥⁰²+𝑒^{−𝑚1𝜏1}+𝑒_{−𝑚2𝜏2}𝛽₁𝑥⁰¹𝑣𝛾𝛽₂𝑥⁰²𝑣−𝑎¹𝑐𝑝₁=𝑒_{−𝑚1𝜏1}𝑒^𝑛1𝜔1𝑣𝜆₁2−𝑥𝑥⁰¹₁−𝑥¹𝑥⁰¹+𝑒^{−𝑚2𝜏2}𝛾𝜆₂2−𝑥𝑥⁰²₂−𝑥²𝑥⁰²+𝑎¹𝑐𝑒^𝑛1𝜔1𝑝₁𝑅₀−1𝑣.(2.21) 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 𝑅₀≤1, then 𝑑𝑊₁/𝑑𝑡≤0 for all 𝑥₁,𝑥₂,𝑣>0. By Theorem  5.3.1 in [32], the solutions of system (2.1)–(2.5) limit to 𝑀, the largest invariant subset of {𝑑𝑊₁/𝑑𝑡=0}. Clearly, it follows from (2.21) that 𝑑𝑊₁/𝑑𝑡=0 if and only if 𝑥₁=𝑥⁰¹, 𝑥₂=𝑥⁰², 𝑣=0. Noting that 𝑀 is invariant, for each element of 𝑀 we have 𝑣=0, 𝑣̇=0. From (2.5) we drive that 0=𝑣̇=𝑒^{−𝑛1𝜔1}𝑝₁𝑦₁𝑡−𝜔₁+𝑒^{−𝑛2𝜔2}𝑝₂𝑦₂𝑡−𝜔₂.(2.22) Since 𝑦₁(𝑡−𝜃)≥0 and 𝑦₂(𝑡−𝜃)≥0 for all 𝜃∈[0,max{𝜏₁,𝜏₂,𝜔₁,𝜔₂}], then 𝑒^{−𝑛1𝜔1}𝑝₁𝑦₁(𝑡−𝜔₁)+𝑒^{−𝑛2𝜔2}𝑝₂𝑦₂(𝑡−𝜔₂)=0 if and only if 𝑦₁(𝑡−𝜔₁)=𝑦₂(𝑡−𝜔₂)=0. Hence, 𝑑𝑊₁/𝑑𝑡=0 if and only if 𝑥₁=𝑥⁰¹, 𝑥₂=𝑥⁰², 𝑦₁=𝑦₂=𝑣=0. From LaSalle’s invariance principle, 𝐸₀ is GAS for any 𝜏₁,𝜏₂,𝜔₁,𝜔₂≥0.
(ii) Define a Lyapunov functional as 𝑊₂=𝑒^{−𝑚1𝜏1}𝑥^∗1𝐻𝑥¹𝑥^∗1+𝑦^∗1𝐻𝑦¹𝑦^∗1+𝛾𝑒^{−𝑚2𝜏2}𝑥^∗2𝐻𝑥²𝑥^∗2+𝑦^∗2𝐻𝑦²+𝑎𝑦^∗2₁𝑝₁𝑒^𝑛1𝜔1𝑣^∗𝐻𝑣𝑣^∗+𝑒^{−𝑚1𝜏1}𝛽₁𝑥^∗1𝑣^∗^𝜏1₀𝐻𝑥¹(𝑡−𝜃)𝑣(𝑡−𝜃)+𝑒_{−𝑚2𝜏2}𝑥^∗1𝑣^∗𝑑𝜃𝛾𝛽₂𝑥^∗2𝑣^∗^𝜏2₀𝐻𝑥²(𝑡−𝜃)𝑣(𝑡−𝜃)𝑥^∗2𝑣^∗𝑑𝜃+𝑎¹𝑦^∗1^𝜔1₀𝐻𝑦¹(𝑡−𝜃)+𝛾𝑎₂𝑦^∗1𝑑𝜃𝑦^∗2^𝜔2₀𝐻𝑦²(𝑡−𝜃)𝑦^∗2𝑑𝜃.(2.23) Differentiating with respect to time yields 𝑑𝑊₂𝑑𝑡=𝑒^{−𝑚1𝜏1}1−𝑥𝑥^∗1₁𝜆¹−𝑑₁𝑥₁−𝛽₁𝑥₁𝑣+1−𝑦^∗1𝑦₁𝑒^{−𝑚1𝜏1}+𝛾𝑒_{−𝑚2𝜏2}𝛽₁𝑥₁𝑡−𝜏₁𝑣𝑡−𝜏₁−𝑎₁𝑦₁1−𝑥𝑥^∗2₂𝜆²+1−𝑦_∗2−𝑑₂𝑥₂−𝛽₂𝑥₂𝑣𝑦₂𝑒^{−𝑚2𝜏2}+𝑎₁𝛽₂𝑥₂𝑡−𝜏₂𝑣𝑡−𝜏₂−𝑎₂𝑦₂𝑝₁𝑒^𝑛1𝜔11−𝑣𝑣^∗𝑒^{−𝑛1𝜔1}+𝑒_{−𝑚1𝜏1}𝑝₁𝑦₁𝑡−𝜔₁+𝑒^{−𝑛2𝜔2}𝑝₂𝑦₂𝑡−𝜔₂−𝑐𝑣𝛽¹𝑥₁𝑣−𝛽₁𝑥₁𝑡−𝜏₁𝑣𝑡−𝜏₁+𝛽₁𝑥^∗1𝑣^∗ln𝑥¹𝑡−𝜏₁𝑣𝑡−𝜏₁+𝛾𝑒_{−𝑚2𝜏2}𝑥₁𝑣𝛽²𝑥₂𝑣−𝛽₂𝑥₂𝑡−𝜏₂𝑣𝑡−𝜏₂+𝛽₂𝑥^∗2𝑣^∗ln𝑥²𝑡−𝜏₂𝑣𝑡−𝜏₂+𝑎₁𝑥₂𝑣𝑦¹−𝑦₁𝑡−𝜔₁+𝑦^∗1ln𝑦¹𝑡−𝜔₁+𝑎₂𝑦₁𝛾𝑦₂−𝑦₂𝑡−𝜔₂+𝑦^∗2ln𝑦²𝑡−𝜔₂=𝑒_{−𝑚1𝜏1}𝑦₂1−𝑥𝑥^∗1₁𝜆¹−𝑑₁𝑥₁+𝑒^{−𝑚1𝜏1}𝛽₁𝑥^∗1𝑣−𝑒^{−𝑚1𝜏1}𝑦^∗1𝛽₁𝑥₁𝑡−𝜏₁𝑣𝑡−𝜏₁𝑦₁+𝑎₁𝑦^∗1+𝛾𝑒^{−𝑚2𝜏2}1−𝑥𝑥^∗2₂𝜆²−𝑑₂𝑥₂+𝑒^{−𝑚2𝜏2}𝛽₂𝑥^∗2𝑣−𝑒^{−𝑚2𝜏2}𝑦^∗2𝛽₂𝑥₂𝑡−𝜏₂𝑣𝑡−𝜏₂𝑦₂+𝑎₂𝑦^∗2−𝑎¹𝑣^∗𝑦₁𝑡−𝜔₁𝑣−𝛾𝑎²𝑣^∗𝑦₂𝑡−𝜔₂𝑣−𝑎¹𝑐𝑝₁𝑒^𝑛1𝜔1𝑣+𝑎¹𝑐𝑝₁𝑒^𝑛1𝜔1𝑣^∗+𝑒^{−𝑚1𝜏1}𝛽₁𝑥^∗1𝑣^∗ln𝑥¹𝑡−𝜏₁𝑣𝑡−𝜏₁𝑥₁𝑣+𝛾𝑒^{−𝑚2𝜏2}𝛽₂𝑥^∗2𝑣^∗ln𝑥²𝑡−𝜏₂𝑣𝑡−𝜏₂+𝑎₁𝑥₂𝑣𝑦^∗1ln𝑦¹𝑡−𝜔₁𝑦₁+𝛾𝑎²𝑦^∗2ln𝑦²𝑡−𝜔₂𝑦₂.(2.24) Using the infected steady state 𝐸₁ conditions 𝜆₁=𝑑₁𝑎𝑥^∗1₂+𝛽₁𝑥^∗1𝑣^∗,𝜆₂=𝑑₂𝑥^∗2+𝛽₂𝑥^∗2𝑣^∗,𝑎₁𝑦^∗1𝑒^𝑚1𝜏1=𝛽₁𝑥^∗1𝑣^∗,𝑦^∗2𝑒^𝑚2𝜏2=𝛽₂𝑥^∗2𝑣^∗,𝑐𝑣^∗=𝑝₁𝑒^{−𝑛1𝜔1}𝑦^∗1+𝑝₂𝑒^{−𝑛2𝜔2}𝑦^∗2,(2.25) we obtain 𝑑𝑊₂𝑑𝑡=𝑒^{−𝑚1𝜏1}𝑑¹𝑥^∗1+𝑎₁𝑦^∗1𝑒^𝑚1𝜏1−𝑑₁𝑥₁−𝑥^∗1𝑥₁𝑑₁+𝑒_{−𝑚1𝜏1}𝑥^∗1+𝑎₁𝑦^∗1𝑒^𝑚1𝜏1−𝑑₁𝑥₁𝛽₁𝑥^∗1𝑣^∗𝑣𝑣^∗−𝑎¹𝑦^∗1𝑦^∗1𝑥₁𝑡−𝜏₁𝑣𝑡−𝜏₁𝑦₁𝑥^∗1𝑣^∗+2𝑎₁𝑦^∗1+𝛾𝑒^{−𝑚2𝜏2}𝑑₂𝑥^∗2+𝑎₂𝑦^∗2𝑒^𝑚2𝜏2−𝑑₂𝑥₂−𝑒^{−𝑚2𝜏2}𝑥^∗2𝑥₂𝑑₂+𝑒_{−𝑚2𝜏2}𝑥^∗2+𝑎₂𝑦^∗2𝑒^𝑚2𝜏2−𝑑₂𝑥₂𝛽₂𝑥^∗2𝑣^∗𝑣𝑣^∗−𝑎²𝑦^∗2𝑦^∗2𝑥₂𝑡−𝜏₂𝑣𝑡−𝜏₂𝑦₂𝑥^∗2𝑣^∗+2𝑎₂𝑦^∗2−𝑎₁𝑦^∗1𝑣^∗𝑦₁𝑡−𝜔₁𝑣𝑦^∗1−𝛾𝑎₂𝑦^∗2𝑣^∗𝑦₂𝑡−𝜔₂𝑣𝑦^∗2−𝑎¹𝑐𝑝₁𝑒^𝑛1𝜔1𝑣^∗𝑣+𝑎₁𝑣^∗𝑦^∗1ln𝑥¹𝑡−𝜏₁𝑣𝑡−𝜏₁𝑥₁𝑣+𝛾𝑎²𝑦^∗2ln𝑥²𝑡−𝜏₂𝑣𝑡−𝜏₂+𝑎₁𝑥₂𝑣𝑦^∗1ln𝑦¹𝑡−𝜔₁𝑦₁+𝛾𝑎²𝑦^∗2ln𝑦²𝑡−𝜔₂=𝑒_{−𝑚1𝜏1}𝑦₂𝑑₁𝑥^∗12−𝑥𝑥^1∗1−𝑥^∗1𝑥₁−𝑎¹𝑦^∗1𝑥^∗1𝑥₁−𝑎₁𝑦^∗1𝑦^∗1𝑥₁𝑡−𝜏₁𝑣𝑡−𝜏₁𝑦₁𝑥^∗1𝑣^∗−𝑎₁𝑦^∗1𝑣^∗𝑦₁𝑡−𝜔₁𝑣𝑦^∗1+3𝑎₁𝑦^∗1+𝑎₁𝑦^∗1ln𝑥¹𝑡−𝜏₁𝑣𝑡−𝜏₁𝑦₁𝑡−𝜔₁+𝛾𝑒_{−𝑚2𝜏2}𝑥₁𝑣𝑦₁𝑑₂𝑥^∗22−𝑥𝑥^2∗2−𝑥^∗2𝑥₂−𝑎²𝑦^∗2𝑥^∗2𝑥₂−𝑎₂𝑦^∗2𝑦^∗2𝑥₂𝑡−𝜏₂𝑣𝑡−𝜏₂𝑦₂𝑥^∗2𝑣^∗−𝑎₂𝑦^∗2𝑣^∗𝑦₂𝑡−𝜔₂𝑣𝑦^∗2+3𝑎₂𝑦^∗2+𝑎₂𝑦^∗2ln𝑥²𝑡−𝜏₂𝑣𝑡−𝜏₂𝑦₂𝑡−𝜔₂+𝑎𝑥₂𝑣𝑦₂₁𝑝₁𝑒^𝑛1𝜔1𝑒^{−𝑛1𝜔1}𝑝₁𝑦^∗1+𝑒^{−𝑛2𝜔2}𝑝₂𝑦^∗2−𝑐𝑣^∗𝑣𝑣^∗.(2.26) From (2.25), we see that the last term in (2.26) vanishes. Then, using the following equalities: ln𝑥¹𝑡−𝜏₁𝑣𝑡−𝜏₁𝑦₁𝑡−𝜔₁𝑥₁𝑣𝑦₁=ln𝑥𝑥^∗1₁+ln𝑦¹𝑡−𝜔₁𝑣^∗+ln𝑦_∗1𝑦^∗1𝑣𝑥₁𝑡−𝜏₁𝑣𝑡−𝜏₁ln𝑥₂𝑦₁𝑥^∗1𝑣^∗,𝑡−𝜏₂𝑣𝑡−𝜏₂𝑦₂𝑡−𝜔₂𝑥₂𝑣𝑦₂=ln𝑥𝑥^∗2₂+ln𝑦²𝑡−𝜔₂𝑣^∗+ln𝑦_∗2𝑦^∗2𝑣𝑥₂𝑡−𝜏₂𝑣𝑡−𝜏₂𝑦₂𝑥^∗2𝑣^∗,(2.27) we can rewrite (2.26) as 𝑑𝑊₂𝑑𝑡=𝑒^{−𝑚1𝜏1}𝑑₁𝑥^∗12−𝑥𝑥^1∗1−𝑥^∗1𝑥₁+𝑒^{−𝑚2𝜏2}𝛾𝑑₂𝑥^∗22−𝑥𝑥^2∗2−𝑥^∗2−𝑎₁𝑥₂𝑦^∗1𝐻𝑥𝑥^∗1₁+𝐻𝑦¹𝑡−𝜔₁𝑣^∗𝑦^∗1𝑣+𝐻𝑦^∗1𝑥₁𝑡−𝜏₁𝑣𝑡−𝜏₁−𝛾𝑎₂𝑦₁𝑥^∗1𝑣^∗𝑦^∗2𝐻𝑥𝑥^∗2₂+𝐻𝑦²𝑡−𝜔₂𝑣^∗𝑦^∗2𝑣+𝐻𝑦^∗2𝑥₂𝑡−𝜏₂𝑣𝑡−𝜏₂𝑦₂𝑥^∗2𝑣^∗.(2.28) 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 𝑥^∗1,𝑦^∗1,𝑥^∗2,𝑦^∗2,𝑣^∗>0, then 𝑑𝑊₂/𝑑𝑡≤0. By [32,Theorem  5.3.1], the solutions of system (2.1)–(2.5) limit to 𝑀, the largest invariant subset of {𝑑𝑊₂/𝑑𝑡=0}. It can be seen that 𝑑𝑊₂/𝑑𝑡=0 if and only if 𝑥₁=𝑥^∗1,𝑥₂=𝑥^∗2,𝑣=𝑣^∗, and 𝐻=0, that is, 𝑦₁𝑡−𝜔₁𝑣^∗𝑦^∗1𝑣=𝑦²𝑡−𝜔₂𝑣^∗𝑦^∗2𝑣=𝑦^∗1𝑥₁𝑡−𝜏₁𝑣𝑡−𝜏₁𝑦₁𝑥^∗1𝑣^∗=𝑦^∗2𝑥₂𝑡−𝜏₂𝑣𝑡−𝜏₂𝑦₂𝑥^∗2𝑣^∗=1.(2.29) If 𝑣=𝑣^∗, then from (2.29), we have 𝑦₁=𝑦^∗1 and 𝑦₂=𝑦^∗2, and hence 𝑑𝑊₂/𝑑𝑡 equal to zero at 𝐸₁. LaSalle’s invariance principle implies global stability of 𝐸₁.

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 with 𝑛 classes of target cells. Two types of discrete-time delays (𝜏_𝑖,𝜔_𝑖,𝑖=1,…,𝑛) 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 𝑛=1 and 𝜔₁=0, then the following model leads to the model presented in [11]. 𝑦̇_𝑖𝑥̇_𝑖=𝜆_𝑖−𝑑_𝑖𝑥_𝑖−𝛽_𝑖𝑥_𝑖𝑣,𝑖=1,…,𝑛,=𝑒^{−𝑚𝑖𝜏𝑖}𝛽𝑣_𝑖𝑥̇=_𝑖𝑡_𝑛−𝜏_𝑖𝑣𝑡−𝜏_𝑖−𝑎_𝑖𝑦_𝑖,𝑖=1,…,𝑛,_𝑖=1𝑒^{−𝑛𝑖𝜔𝑖}𝑝_𝑖𝑦_𝑖𝑡−𝜔_𝑖−𝑐𝑣,(3.1) where 𝑥_𝑖 and 𝑦_𝑖 represent the populations of the uninfected target cells and infected cells of class 𝑖, respectively, 𝑣 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𝑦𝑥_𝑗(𝜃)=𝜑_𝑗(𝜃),𝑗=1,…,𝑛,_𝑗(𝜃)=𝑣𝜙(𝜃)_𝑗+𝑛=(𝜃𝜑),_2𝑛+1𝑗=1,…,𝑛,𝜑(𝜃),_𝑗(𝜃)≥0,𝑗=1,…,2𝑛+1,𝜃∈−max𝜏₁,…,𝜏_𝑛,𝜔₁,…,𝜔_𝑛,0,(3.2) where (𝜑₁(𝜃),𝜑₂(𝜃),…,𝜑_2𝑛+1(𝜃))∈𝐶 and 𝐶=𝐶([−max{𝜏₁,…,𝜏_𝑛,𝜔₁,…,𝜔_𝑛},0],ℝ^+2𝑛+1) is the Banach space of continuous functions mapping the interval [−max{𝜏₁,…,𝜏_𝑛,𝜔₁,…,𝜔_𝑛},0] into ℝ^+2𝑛+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 𝐸₀=(𝑥⁰¹,…,𝑥^0𝑛,𝑦⁰¹,…,𝑦^0𝑛,𝑣⁰), where 𝑥^0𝑖=𝜆_𝑖/𝑑_𝑖, 𝑦^0𝑖=0, 𝑖=1,…,𝑛, and 𝑣⁰=0. The system can also have a positive infected steady state 𝐸₁(𝑥^∗1,…,𝑥^∗𝑛,𝑦^∗1,…,𝑦^∗𝑛,𝑣^∗). The coordinates of the infected steady state, if they exist, satisfy the equalities:𝑎𝜆_𝑖=𝑑_𝑖𝑥^∗𝑖+𝛽_𝑖𝑥^∗𝑖𝑣^∗,𝑖=1,…,𝑛,(3.3)_𝑖𝑦^∗𝑖=𝑒𝑐𝑣^{−𝑚𝑖𝜏𝑖}_∗𝛽_𝑖𝑥^∗𝑖𝑣^∗,𝑖=1,…,𝑛,(3.4)=^𝑛_𝑖=1𝑒^{−𝑛𝑖𝜔𝑖}𝑝_𝑖𝑦^∗𝑖.(3.5) The basic reproduction number of system (3.1) is given by𝑅₀=^𝑛_𝑖=1𝑅_𝑖=^𝑛_𝑖=1𝑒^{−(𝑚𝑖𝜏𝑖+𝑛𝑖𝜔𝑖)}𝛽_𝑖𝑝_𝑖𝜆_𝑖𝑎_𝑖𝑑_𝑖𝑐,(3.6) where 𝑅_𝑖 is the basic reproduction number for the dynamics of the interaction of the virus only with the target cells of class 𝑖.

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) If 𝑅₀≤1, then 𝐸₀ is GAS for any 𝜏_𝑖,𝜔_𝑖≥0, 𝑖=1,…,𝑛.
(ii) If 𝐸₁ exists, then it is GAS for any 𝜏_𝑖,𝜔_𝑖≥0, 𝑖=1,…,𝑛.

Proof. (i) Define a Lyapunov functional 𝑊₁ as follows: 𝑊₁=^𝑛_𝑖=1𝛾_𝑖𝑒^{−𝑚𝑖𝜏𝑖}𝑥^0𝑖𝐻𝑥^𝑖𝑥^0𝑖+𝑦^𝑖+𝑒^{−𝑚𝑖𝜏𝑖}𝛽_𝑖^𝜏𝑖₀𝑥_𝑖(𝑡−𝜃)𝑣(𝑡−𝜃)𝑑𝜃+𝑎_𝑖^𝜔𝑖₀+𝑎𝑦_𝑖(𝑡−𝜃)𝑑𝜃₁𝑝₁𝑒^𝑛1𝜔1𝑣,(3.7) where 𝛾_𝑖=(𝑎₁𝑝_𝑖/𝑎_𝑖𝑝₁)𝑒^{𝑛1𝜔1−𝑛𝑖𝜔𝑖}. The time derivative of 𝑊₁ along the solution of system (3.1) satisfies 𝑑𝑊₁𝑑𝑡=^𝑛_𝑖=1𝛾_𝑖𝑒^{−𝑚𝑖𝜏𝑖}1−𝑥𝑥^0𝑖_𝑖𝜆^𝑖−𝑑_𝑖𝑥_𝑖−𝛽_𝑖𝑥_𝑖𝑣+𝑒^{−𝑚𝑖𝜏𝑖}𝛽_𝑖𝑥_𝑖𝑡−𝜏_𝑖𝑣𝑡−𝜏_𝑖−𝑎_𝑖𝑦_𝑖+𝑎+𝑒^{−𝑚𝑖𝜏𝑖}𝛽_𝑖𝑥_𝑖𝑣−𝛽_𝑖𝑥_𝑖𝑡−𝜏_𝑖𝑣𝑡−𝜏_𝑖+𝑎_𝑖𝑦_𝑖−𝑦_𝑖𝑡−𝜔_𝑖₁𝑝₁𝑒^𝑛1𝜔1^𝑛_𝑖=1=𝑒^{−𝑛𝑖𝜔𝑖}𝑝_𝑖𝑦_𝑖𝑡−𝜔_𝑖−𝑐𝑣_𝑛_𝑖=1𝑒^{−𝑚𝑖𝜏𝑖}𝛾_𝑖𝜆_𝑖2−𝑥𝑥^𝑖0𝑖−𝑥^0𝑖𝑥_𝑖+𝑎¹𝑐𝑝₁𝑒^𝑛1𝜔1^𝑛_𝑖=1𝑒^{−(𝑚𝑖𝜏𝑖+𝑛𝑖𝜔𝑖)}𝑝_𝑖𝛽_𝑖𝑥^0𝑖=𝑎_𝑖𝑐−1_𝑛𝑣_𝑖=1𝑒^{−𝑚𝑖𝜏𝑖}𝛾_𝑖𝜆_𝑖2−𝑥𝑥^𝑖0𝑖−𝑥^0𝑖𝑥_𝑖+𝑎¹𝑐𝑝₁𝑒^𝑛1𝜔1𝑅₀−1𝑣.(3.8) 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 𝑅₀≤1, then 𝑑𝑊₁/𝑑𝑡≤0 for all 𝑥_𝑖,𝑦_𝑖,𝑣>0. Similar to the previous section, one can show that the maximal compact invariant set in {𝑑𝑊₁/𝑑𝑡=0} is the singleton {𝐸₀} when 𝑅₀≤1. The global stability of 𝐸₀ follows from LaSalle’s invariance principle.
To prove (ii), we consider the Lyapunov functional 𝑊₂=^𝑛_𝑖=1𝛾_𝑖𝑒^{−𝑚𝑖𝜏𝑖}𝑥^∗𝑖𝐻𝑥^𝑖𝑥^∗𝑖+𝑦^∗𝑖𝐻𝑦^𝑖𝑦^∗𝑖+𝑒^{−𝑚𝑖𝜏𝑖}𝛽_𝑖𝑥^∗𝑖𝑣^∗^𝜏𝑖₀𝐻𝑥^𝑖(𝑡−𝜃)𝑣(𝑡−𝜃)+𝑎_𝑖𝑥^∗𝑖𝑣^∗𝑑𝜃𝑦^∗𝑖^𝜔𝑖₀𝐻𝑦^𝑖(𝑡−𝜃)𝑦^∗𝑖𝑑𝜃+𝑎𝑝¹₁𝑒^𝑛1𝜔1𝑣^∗𝐻𝑣𝑣^∗.(3.9) Differentiating with respect to time yields 𝑑𝑊₂𝑑𝑡=^𝑛_𝑖=1𝛾_𝑖𝑒^{−𝑚𝑖𝜏𝑖}1−𝑥𝑥^∗𝑖_𝑖𝜆^𝑖−𝑑_𝑖𝑥_𝑖−𝛽_𝑖𝑥_𝑖𝑣+1−𝑦^∗𝑖𝑦_𝑖𝑒^{−𝑚𝑖𝜏𝑖}+𝑒_{−𝑚𝑖𝜏𝑖}𝛽_𝑖𝑥_𝑖𝑡−𝜏_𝑖𝑣𝑡−𝜏_𝑖−𝑎_𝑖𝑦_𝑖𝛽_𝑖𝑥_𝑖𝑣−𝛽_𝑖𝑥_𝑖𝑡−𝜏_𝑖𝑣𝑡−𝜏_𝑖+𝑒^{−𝑚𝑖𝜏𝑖}𝛽_𝑖𝑥^∗𝑖𝑣^∗ln𝑥^𝑖𝑡−𝜏_𝑖𝑣𝑡−𝜏_𝑖+𝑎_𝑖𝑥_𝑖𝑣𝑦_𝑖−𝑦_𝑖𝑡−𝜔_𝑖+𝑎_𝑖𝑦^∗𝑖ln𝑦^𝑖𝑡−𝜔_𝑖+𝑎𝑦_𝑖₁𝑒^𝑛1𝜔1𝑝₁1−𝑣𝑣^∗^𝑛_𝑖=1=𝑒^{−𝑛𝑖𝜔𝑖}𝑝_𝑖𝑦_𝑖𝑡−𝜔_𝑖−𝑐𝑣_𝑛_𝑖=1𝛾_𝑖𝑒^{−𝑚𝑖𝜏𝑖}𝜆^𝑖−𝑑_𝑖𝑥_𝑖−𝜆^𝑖𝑥^∗𝑖𝑥_𝑖+𝑑_𝑖𝑥^∗𝑖+𝛽_𝑖𝑥^∗𝑖𝑣−𝑒^{−𝑚𝑖𝜏𝑖}𝛽_𝑖𝑥_𝑖𝑡−𝜏_𝑖𝑣𝑡−𝜏_𝑖𝑦^∗𝑖𝑦_𝑖+𝑎_𝑖𝑦^∗𝑖+𝑒^{−𝑚𝑖𝜏𝑖}𝛽_𝑖𝑥^∗𝑖𝑣^∗ln𝑥^𝑖𝑡−𝜏_𝑖𝑣𝑡−𝜏_𝑖𝑥_𝑖𝑣+𝑎^𝑖𝑦^∗𝑖ln𝑦^𝑖𝑡−𝜔_𝑖−𝑎𝑦_𝑖₁𝑐𝑒^𝑛1𝜔1𝑝₁𝑣−𝑣^∗𝑣^𝑛_𝑖=1𝛾_𝑖𝑎_𝑖𝑦_𝑖𝑡−𝜔_𝑖+𝑎¹𝑐𝑒^𝑛1𝜔1𝑝₁𝑣^∗.(3.10) Using the infected steady state conditions (3.3)–(3.5), and the following equality: 𝑎₁𝑐𝑒^𝑛1𝜔1𝑝₁𝑣^∗=𝑎¹𝑒^𝑛1𝜔1𝑝_1𝑛𝑖=1𝑒^{−𝑛𝑖𝜔𝑖}𝑝_𝑖𝑦^∗𝑖=^𝑛_𝑖=1𝛾_𝑖𝑎_𝑖𝑦^∗𝑖,(3.11) we obtain 𝑑𝑊₂𝑑𝑡=^𝑛_𝑖=1𝛾_𝑖𝑒^{−𝑚𝑖𝜏𝑖}𝑑^𝑖𝑥^∗𝑖+𝑒^{𝑚𝑖𝜏𝑖}𝑎_𝑖𝑦^∗𝑖−𝑑_𝑖𝑥_𝑖−𝑥^∗𝑖𝑥_𝑖𝑑_𝑖𝑥^∗𝑖+𝑒^{𝑚𝑖𝜏𝑖}𝑎_𝑖𝑦^∗𝑖+𝑑_𝑖𝑥^∗𝑖−𝑎_𝑖𝑦^∗𝑖𝑦^∗𝑖𝑥_𝑖𝑡−𝜏_𝑖𝑣𝑡−𝜏_𝑖𝑦_𝑖𝑥^∗𝑖𝑣^∗+2𝑎_𝑖𝑦^∗𝑖−𝑎_𝑖𝑦^∗𝑖𝑣^∗𝑦_𝑖𝑡−𝜔_𝑖𝑣𝑦^∗𝑖+𝑎_𝑖𝑦^∗𝑖ln𝑥^𝑖𝑡−𝜏_𝑖𝑣𝑡−𝜏_𝑖𝑥_𝑖𝑣+𝑎^𝑖𝑦^∗𝑖ln𝑦^𝑖𝑡−𝜔_𝑖+𝑦_𝑖_𝑛_𝑖=1𝑒^{−𝑚𝑖𝜏𝑖}𝛾_𝑖𝛽_𝑖𝑥^∗𝑖−𝑎¹𝑐𝑒^𝑛1𝜔1=_𝑛𝑝₁𝑣_𝑖=1𝛾_𝑖𝑒^{−𝑚𝑖𝜏𝑖}𝑑_𝑖𝑥^∗𝑖2−𝑥𝑥^∗𝑖_𝑖−𝑥^𝑖𝑥^∗𝑖−𝑎^𝑖𝑦^∗𝑖𝑥^∗𝑖𝑥_𝑖−𝑎_𝑖𝑦^∗𝑖𝑣^∗𝑦_𝑖𝑡−𝜔_𝑖𝑣𝑦^∗𝑖−𝑎_𝑖𝑦^∗𝑖𝑦^∗𝑖𝑥_𝑖𝑡−𝜏_𝑖𝑣𝑡−𝜏_𝑖𝑦_𝑖𝑥^∗𝑖𝑣^∗+3𝑎_𝑖𝑦^∗𝑖+𝑎_𝑖𝑦^∗𝑖ln𝑥^𝑖𝑡−𝜏_𝑖𝑣𝑡−𝜏_𝑖𝑦_𝑖𝑡−𝜔_𝑖+𝑎𝑥_𝑖𝑣𝑦_𝑖₁𝑒^𝑛1𝜔1𝑝₁^𝑛_𝑖=1𝑒^{−𝑛𝑖𝜔𝑖}𝑝_𝑖𝑦^∗𝑖−𝑐𝑣^∗𝑣𝑣^∗.(3.12) From (3.5), we can see that the last term in (3.12) vanishes. Then, by using the following equality: ln𝑥^𝑖𝑡−𝜏_𝑖𝑣𝑡−𝜏_𝑖𝑦_𝑖𝑡−𝜔_𝑖𝑥_𝑖𝑣𝑦_𝑖=ln𝑥𝑥^∗𝑖_𝑖+ln𝑣^∗𝑦_𝑖𝑡−𝜔_𝑖𝑣𝑦^∗𝑖+ln𝑦^∗𝑖𝑥_𝑖𝑡−𝜏_𝑖𝑣𝑡−𝜏_𝑖𝑦_𝑖𝑥^∗𝑖𝑣^∗(3.1,3) we can rewrite (3.12) as 𝑑𝑊₂𝑑𝑡=^𝑛_𝑖=1𝛾_𝑖𝑒^{−𝑚𝑖𝜏𝑖}𝑑_𝑖𝑥^∗𝑖2−𝑥𝑥^∗𝑖_𝑖−𝑥^𝑖𝑥^∗𝑖−𝑎^𝑖𝑦^∗𝑖𝐻𝑥𝑥^∗𝑖_𝑖+𝐻𝑣^∗𝑦_𝑖𝑡−𝜔_𝑖+𝐻𝑦_∗𝑖𝑣𝑦^∗𝑖𝑥_𝑖𝑡−𝜏_𝑖𝑣𝑡−𝜏_𝑖𝑦_𝑖𝑥^∗𝑖𝑣^∗.(3.14) It is easy to see that if 𝑥^∗𝑖,𝑦^∗𝑖,𝑣^∗>0, 𝑖=1,…,𝑛, then 𝑑𝑊₂/𝑑𝑡≤0 for all (𝑥_𝑖,𝑦_𝑖,𝑣)>0 (the arithmetical mean is greater than or equal to the geometrical mean and 𝐻≥0). Clearly, the singleton {𝐸₁} is the only invariant set in {𝑑𝑊₂/𝑑𝑡=0}. LaSalle's invariance principle implies global stability of 𝐸₁.

4. Virus Dynamics Model with Saturation Infection Rate

In this section, we proposed a virus dynamics model which describes the interaction of the virus with 𝑛 classes of target cells taking into account the saturation infection rate and multiple intracellular delays: 𝑥̇_𝑖(𝑡)=𝜆_𝑖−𝑑_𝑖𝑥_𝑖(𝑡)−𝛽^𝑖𝑥_𝑖(𝑡)𝑣(𝑡)𝑦̇_𝑖1+𝛼_𝑖𝑣(𝑡),𝑖=1,…,𝑛,(𝑡)=𝑒^{−𝑚𝑡𝜏𝑖}𝛽_𝑖𝑥_𝑖𝑡−𝜏_𝑖𝑣𝑡−𝜏_𝑖1+𝛼_𝑖𝑣(𝑡)−𝑎^𝑖𝑣̇(𝑡)=_𝑛𝑦_𝑖(𝑡),𝑖=1,…,𝑛,_𝑖=1𝑒^{−𝑛𝑖𝜔𝑖}𝑝_𝑖𝑦_𝑖𝑡−𝜔_𝑖−𝑐𝑣(𝑡),(4.1) where 𝛼_𝑖,𝑖=1,…,𝑛 are positive constants. The variables and parameters of the model have the same definitions as given in Section 2. We mention that if 𝑛=1 and 𝜔₁=0, then model (4.1) leads to the model presented in [9], and if 𝑛=2 and 𝜔₁=𝜔₂=0, 𝛼₁=𝛼₂=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 𝐸₀=(𝑥⁰¹,…,𝑥^0𝑛,𝑦⁰¹,…,𝑦^0𝑛,𝑣⁰), where 𝑥^0𝑖=𝜆_𝑖/𝑑_𝑖, 𝑦^0𝑖=0, and 𝑣⁰=0. The system can also have a positive infected steady state 𝐸₁(𝑥^∗1,…,𝑥^∗𝑛,𝑦^∗1,…,𝑦^∗𝑛,𝑣^∗). The coordinates of the infected steady state, if they exist, satisfy the equalities:𝜆_𝑖=𝑑_𝑖𝑥^∗𝑖+𝛽^𝑖𝑥^∗𝑖𝑣^∗1+𝛼_𝑖𝑣^∗𝑎,𝑖=1,...,𝑛,_𝑖𝑦^∗𝑖=𝑒^{−𝑚𝑖𝜏𝑖}𝛽_𝑖𝑥^∗𝑖𝑣^∗1+𝛼_𝑖𝑣^∗𝑐𝑣_∗,𝑖=1,...,𝑛,=^𝑛_𝑖=1𝑒^{−𝑛𝑖𝜔𝑖}𝑝_𝑖𝑦^∗𝑖.(4.2) The basic reproduction number 𝑅₀ 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) If 𝑅₀≤1, then 𝐸₀ is GAS for any 𝜏_𝑖,𝜔_𝑖≥0, 𝑖=1,…,𝑛.
(ii) If 𝐸₁ exists, then it is GAS for any 𝜏_𝑖,𝜔_𝑖≥0, 𝑖=1,…,𝑛.

Proof. (i) Define a Lyapunov functional 𝑊₁ as follows: 𝑊₁=^𝑛_𝑖=1𝛾_𝑖𝑒^{−𝑚𝑖𝜏𝑖}𝑥^0𝑖𝐻𝑥^𝑖𝑥^0𝑖+𝑦^𝑖+𝑒^{−𝑚𝑖𝜏𝑖}𝛽_𝑖^𝜏𝑖₀𝑥_𝑖(𝑡−𝜃)𝑣(𝑡−𝜃)1+𝛼_𝑖𝑣(𝑡−𝜃)𝑑𝜃+𝑎^𝑖^𝜔𝑖₀+𝑎𝑦_𝑖(𝑡−𝜃)𝑑𝜃₁𝑒^𝑛1𝜔1𝑝₁𝑣,(4.3) The time derivative of 𝑊₁ along the trajectories of (4.1) satisfies 𝑑𝑊₁𝑑𝑡=^𝑛_𝑖=1𝛾_𝑖𝑒^{−𝑚𝑖𝜏𝑖}1−𝑥𝑥^0𝑖_𝑖𝜆^𝑖−𝑑_𝑖𝑥_𝑖−𝛽^𝑖𝑥_𝑖𝑣1+𝛼_𝑖𝑣+𝑒^{−𝑚𝑖𝜏𝑖}𝛽_𝑖𝑥_𝑖𝑡−𝜏_𝑖𝑣𝑡−𝜏_𝑖1+𝛼_𝑖𝑣𝑡−𝜏_𝑖−𝑎^𝑖𝑦_𝑖+𝑒^{−𝑚𝑖𝜏𝑖}𝛽^𝑖𝑥_𝑖𝑣1+𝛼_𝑖𝑣−𝛽^𝑖𝑥_𝑖𝑡−𝜏_𝑖𝑣𝑡−𝜏_𝑖1+𝛼_𝑖𝑣𝑡−𝜏_𝑖+𝑎^𝑖+𝑎𝑦_𝑖−𝑦_𝑖𝑡−𝜔_𝑖₁𝑒^𝑛1𝜔1𝑝₁^𝑛_𝑖=1=𝑒^{−𝑛𝑖𝜔𝑖}𝑝_𝑖𝑦_𝑖𝑡−𝜔_𝑖−𝑐𝑣_𝑛_𝑖=1𝛾_𝑖𝑒^{−𝑚𝑖𝜏𝑖}𝜆^𝑖−𝑑_𝑖𝑥_𝑖−𝜆_𝑖𝑥^0𝑖𝑥_𝑖+𝑑_𝑖𝑥^0𝑖+𝛽^𝑖𝑥^0𝑖𝑣1+𝛼_𝑖𝑣−𝑎¹𝑐𝑒^𝑛1𝜔1𝑝₁=_𝑛𝑣_𝑖=1𝛾_𝑖𝑒^{−𝑚𝑖𝜏𝑖}𝜆_𝑖2−𝑥𝑥^𝑖0𝑖−𝑥^0𝑖𝑥_𝑖−𝑎¹𝑐𝑒^𝑛1𝜔1𝑝₁𝑣+𝑎¹𝑐𝑒^𝑛1𝜔1𝑝_1𝑛𝑖=1𝑒^{−(𝑚𝑖𝜏𝑖+𝑛𝑖𝜔𝑖)}𝑝_𝑖𝛽_𝑖𝑥^0𝑖𝑣=_𝑛𝑎_𝑖𝑐1+𝛼_𝑖𝑣_𝑖=1𝛾_𝑖𝑒^{−𝑚𝑖𝜏𝑖}𝜆_𝑖2−𝑥𝑥^𝑖0𝑖−𝑥^0𝑖𝑥_𝑖−𝑎¹𝑐𝑒^𝑛1𝜔1𝑝₁𝑣+𝑎¹𝑐𝑒^𝑛1𝜔1𝑝_1𝑛𝑖=1𝑅_𝑖𝑣=_𝑛1+𝛼_𝑖𝑣_𝑖=1𝛾_𝑖𝑒^{−𝑚𝑖𝜏𝑖}𝜆_𝑖2−𝑥𝑥^𝑖0𝑖−𝑥^0𝑖𝑥_𝑖+𝑎¹𝑐𝑒^𝑛1𝜔1𝑝₁𝑅₀−1𝑣−𝑎¹𝑐𝑒^𝑛1𝜔1𝑝_1𝑛𝑖=1𝑅_𝑖𝛼_𝑖𝑣²1+𝛼_𝑖𝑣.(4.4) It is clear that if 𝑅₀≤1, then 𝑑𝑊₁/𝑑𝑡≤0 for all 𝑥_𝑖,𝑦_𝑖,𝑣>0, where equality occurs at 𝐸₀. The global stability of 𝐸₀ follows from LaSalle’s invariance principle.
To prove (ii), we consider the Lyapunov functional: 𝑊₂=^𝑛_𝑖=1𝛾_𝑖𝑒^{−𝑚𝑖𝜏𝑖}𝑥^∗𝑖𝐻𝑥^𝑖𝑥^∗𝑖+𝑦^∗𝑖𝐻𝑦^𝑖+𝑒_{−𝑚𝑖𝜏𝑖}𝑦^∗𝑖𝛽_𝑖𝑥^∗𝑖𝑣^∗1+𝛼_𝑖𝑣^∗^𝜏𝑖₀𝐻𝑥^𝑖(𝑡−𝜃)𝑣(𝑡−𝜃)1+𝛼_𝑖𝑣^∗+𝑎_𝑖𝑥^∗𝑖𝑣^∗1+𝛼_𝑖𝑣(𝑡−𝜃)𝑑𝜃𝑦^∗𝑖^𝜔𝑖₀𝐻𝑦^𝑖(𝑡−𝜃)𝑦^∗𝑖𝑑𝜃+𝑎¹𝑒^𝑛1𝜔1𝑝₁𝑣^∗𝐻𝑣𝑣^∗.(4.5) Differentiating with respect to time yields 𝑑𝑊₂𝑑𝑡=^𝑛_𝑖=1𝛾_𝑖𝑒^{−𝑚𝑖𝜏𝑖}1−𝑥𝑥^∗𝑖_𝑖𝜆^𝑖−𝑑_𝑖𝑥_𝑖−𝛽^𝑖𝑥_𝑖𝑣+1−𝑦_∗𝑖1+𝛼_𝑖𝑣𝑦_𝑖𝑒^{−𝑚𝑖𝜏𝑖}𝛽_𝑖𝑥_𝑖𝑡−𝜏_𝑖𝑣𝑡−𝜏_𝑖1+𝛼_𝑖𝑣𝑡−𝜏_𝑖−𝑎^𝑖𝑦_𝑖+𝑒_{−𝑚𝑖𝜏𝑖}𝛽^𝑖𝑥_𝑖𝑣1+𝛼_𝑖𝑣−𝛽^𝑖𝑥_𝑖𝑡−𝜏_𝑖𝑣𝑡−𝜏_𝑖+𝛽_𝑖1+𝛼_𝑖𝑣𝑡−𝜏_𝑖𝑥^∗𝑖𝑣^∗1+𝛼_𝑖𝑣^∗ln𝑥^𝑖𝑡−𝜏_𝑖𝑣𝑡−𝜏_𝑖1+𝛼_𝑖𝑣+𝑎_𝑖𝑥_𝑖𝑣1+𝛼_𝑖𝑣𝑡−𝜏_𝑖𝑦_𝑖−𝑦_𝑖𝑡−𝜔_𝑖+𝑎_𝑖𝑦^∗𝑖ln𝑦^𝑖𝑡−𝜔_𝑖+𝑎𝑦_𝑖₁𝑒^𝑛1𝜔1𝑝₁1−𝑣𝑣^∗^𝑛_𝑖=1=𝑒^{−𝑛𝑖𝜔𝑖}𝑝_𝑖𝑦_𝑖𝑡−𝜔_𝑖−𝑐𝑣_𝑛_𝑖=1𝛾_𝑖𝑒^{−𝑚𝑖𝜏𝑖}𝜆^𝑖−𝑑_𝑖𝑥_𝑖−𝜆^𝑖𝑥^∗𝑖𝑥_𝑖+𝑑_𝑖𝑥^∗𝑖+𝛽^𝑖𝑥^∗𝑖𝑣1+𝛼_𝑖𝑣−𝑒^{−𝑚𝑖𝜏𝑖}𝛽_𝑖𝑥_𝑖𝑡−𝜏_𝑖𝑣𝑡−𝜏_𝑖1+𝛼_𝑖𝑣𝑡−𝜏_𝑖𝑦𝑦^∗𝑖_𝑖+𝑎_𝑖𝑦^∗𝑖+𝑒^{−𝑚𝑖𝜏𝑖}𝛽_𝑖𝑥^∗𝑖𝑣^∗1+𝛼_𝑖𝑣^∗ln𝑥^𝑖𝑡−𝜏_𝑖𝑣𝑡−𝜏_𝑖1+𝛼_𝑖𝑣𝑥_𝑖𝑣1+𝛼_𝑖𝑣𝑡−𝜏_𝑖+𝑎^𝑖𝑦^∗𝑖ln𝑦^𝑖𝑡−𝜔_𝑖−𝑎𝑦_𝑖₁𝑐𝑒^𝑛1𝜔1𝑝₁𝑣−𝑎¹𝑒^𝑛1𝜔1𝑝₁𝑣^∗𝑣^𝑛_𝑖=1𝑒^{−𝑛𝑖𝜔𝑖}𝑝_𝑖𝑦_𝑖𝑡−𝜔_𝑖+𝑎¹𝑐𝑒^𝑛1𝜔1𝑝₁𝑣^∗.(4.6) Using the infected steady state conditions (4.2), and the following equality: 𝑎₁𝑐𝑒^𝑛1𝜔1𝑝₁𝑣=𝑎¹𝑐𝑒^𝑛1𝜔1𝑝₁𝑣^∗𝑣𝑣^∗=𝑎¹𝑒^𝑛1𝜔1𝑝₁𝑣𝑣^∗𝑛_𝑖=1𝑒^{−𝑛𝑖𝜔𝑖}𝑝_𝑖𝑦^∗𝑖=𝑣𝑣^∗𝑛_𝑖=1𝛾_𝑖𝑎_𝑖𝑦^∗𝑖,(4.7) we obtain 𝑑𝑊₂𝑑𝑡=^𝑛_𝑖=1𝛾_𝑖𝑒^{−𝑚𝑖𝜏𝑖}𝑑^𝑖𝑥^∗𝑖+𝑒^{𝑚𝑖𝜏𝑖}𝑎_𝑖𝑦^∗𝑖−𝑑_𝑖𝑥_𝑖−𝑥^∗𝑖𝑥_𝑖𝑑_𝑖𝑥^∗𝑖+𝑒^{𝑚𝑖𝜏𝑖}𝑎_𝑖𝑦^∗𝑖+𝑑_𝑖𝑥^∗𝑖+𝑎^𝑖𝑦^∗𝑖𝑣1+𝛼_𝑖𝑣^∗−𝑎_𝑖𝑣^∗1+𝛼_𝑖𝑣𝑦^∗𝑖𝑦^∗𝑖𝑥_𝑖𝑡−𝜏_𝑖𝑣𝑡−𝜏_𝑖1+𝛼_𝑖𝑣^∗𝑦_𝑖𝑥^∗𝑖𝑣^∗1+𝛼_𝑖𝑣𝑡−𝜏_𝑖+2𝑎^𝑖𝑦^∗𝑖−𝑎_𝑖𝑦^∗𝑖𝑣^∗𝑦_𝑖𝑡−𝜔_𝑖𝑣𝑦^∗𝑖+𝑎_𝑖𝑦^∗𝑖ln𝑥^𝑖𝑡−𝜏_𝑖𝑣𝑡−𝜏_𝑖𝑦_𝑖𝑡−𝜔_𝑖1+𝛼_𝑖𝑣𝑥_𝑖𝑣𝑦_𝑖1+𝛼_𝑖𝑣𝑡−𝜏_𝑖−𝑎^𝑖𝑦^∗𝑖𝑣=_𝑛𝑣^∗_𝑖=1𝛾_𝑖𝑒^{−𝑚𝑖𝜏𝑖}𝑑_𝑖𝑥^∗𝑖2−𝑥𝑥^∗𝑖_𝑖−𝑥^𝑖𝑥^∗𝑖−𝑎^𝑖𝑦^∗𝑖𝑥^∗𝑖𝑥_𝑖+3𝑎_𝑖𝑦^∗𝑖+𝑎_𝑖𝑦^∗𝑖𝑣1+𝛼^𝑖𝑣^∗−𝑎_𝑖𝑣^∗1+𝛼_𝑖𝑣−𝑣𝑣^∗𝑦^∗𝑖𝑦^∗𝑖𝑥_𝑖𝑡−𝜏_𝑖𝑣𝑡−𝜏_𝑖1+𝛼_𝑖𝑣^∗𝑦_𝑖𝑥^∗𝑖𝑣^∗1+𝛼_𝑖𝑣𝑡−𝜏_𝑖−𝑎^𝑖𝑦^∗𝑖𝑣^∗𝑦_𝑖𝑡−𝜔_𝑖𝑣𝑦^∗𝑖+𝑎_𝑖𝑦^∗𝑖ln𝑥^𝑖𝑡−𝜏_𝑖𝑣𝑡−𝜏_𝑖𝑦_𝑖𝑡−𝜔_𝑖1+𝛼_𝑖𝑣=𝑥_𝑖𝑣𝑦_𝑖1+𝛼_𝑖𝑣𝑡−𝜏_𝑖_𝑛_𝑖=1𝛾_𝑖𝑒^{−𝑚𝑖𝜏𝑖}𝑑_𝑖𝑥^∗𝑖2−𝑥𝑥^∗𝑖_𝑖−𝑥^𝑖𝑥^∗𝑖−𝑎^𝑖𝑦^∗𝑖𝑥^∗𝑖𝑥_𝑖+4𝑎_𝑖𝑦^∗𝑖−𝑎_𝑖𝑦^∗𝑖𝑦^∗𝑖𝑥_𝑖𝑡−𝜏_𝑖𝑣𝑡−𝜏_𝑖1+𝛼_𝑖𝑣^∗+𝑎_𝑖𝑦_𝑖𝑥^∗𝑖𝑣^∗1+𝛼_𝑖𝑣𝑡−𝜏_𝑖𝑦^∗𝑖−1+𝑣1+𝛼^𝑖𝑣^∗𝑣^∗1+𝛼_𝑖𝑣−𝑣𝑣^∗+1+𝛼^𝑖𝑣1+𝛼_𝑖𝑣^∗−1+𝛼^𝑖𝑣1+𝛼_𝑖𝑣^∗−𝑎^𝑖𝑦^∗𝑖𝑣^∗𝑦_𝑖𝑡−𝜔_𝑖𝑣𝑦^∗𝑖+𝑎_𝑖𝑦^∗𝑖ln𝑥^𝑖𝑡−𝜏_𝑖𝑣𝑡−𝜏_𝑖𝑦_𝑖𝑡−𝜔_𝑖1+𝛼_𝑖𝑣𝑥_𝑖𝑣𝑦_𝑖1+𝛼_𝑖𝑣𝑡−𝜏_𝑖.(4.8) Then using the following equalities: ln𝑥^𝑖𝑡−𝜏_𝑖𝑣𝑡−𝜏_𝑖𝑦_𝑖𝑡−𝜔_𝑖1+𝛼_𝑖𝑣𝑥_𝑖𝑣𝑦_𝑖1+𝛼_𝑖𝑣𝑡−𝜏_𝑖=ln𝑥𝑥^∗𝑖_𝑖+ln𝑣^∗𝑦_𝑖𝑡−𝜔_𝑖𝑣𝑦^∗𝑖+ln1+𝛼^𝑖𝑣+ln𝑦_∗𝑖1+𝛼_𝑖𝑣^∗𝑥_𝑖𝑡−𝜏_𝑖𝑣𝑡−𝜏_𝑖1+𝛼_𝑖𝑣^∗−1+𝑣1+𝛼_𝑖𝑦_𝑖𝑥^∗𝑖𝑣^∗1+𝛼_𝑖𝑣𝑡−𝜏_𝑖,𝑣^∗𝑣^∗1+𝛼_𝑖𝑣−𝑣𝑣^∗+1+𝛼^𝑖𝑣1+𝛼_𝑖𝑣^∗=−𝛼^𝑖𝑣−𝑣^∗²𝑣^∗1+𝛼_𝑖𝑣^∗1+𝛼_𝑖𝑣,(4.9) we obtain 𝑑𝑊₂𝑑𝑡=^𝑛_𝑖=1𝛾_𝑖𝑒^{−𝑚𝑖𝜏𝑖}𝑑_𝑖𝑥^∗𝑖2−𝑥𝑥^∗𝑖_𝑖−𝑥^𝑖𝑥^∗𝑖−𝑎^𝑖𝑦^∗𝑖𝛼_𝑖𝑣−𝑣^∗²−𝑎_𝑖𝑣^∗1+𝛼_𝑖𝑣^∗1+𝛼_𝑖𝑣𝑦^∗𝑖𝐻𝑥𝑥^∗𝑖_𝑖+𝐻𝑣^∗𝑦_𝑖𝑡−𝜔_𝑖𝑣𝑦^∗𝑖+𝐻1+𝛼^𝑖𝑣+𝐻𝑦_∗𝑖1+𝛼_𝑖𝑣^∗𝑥_𝑖𝑡−𝜏_𝑖𝑣𝑡−𝜏_𝑖1+𝛼_𝑖𝑣^∗𝑦_𝑖𝑥^∗𝑖𝑣^∗1+𝛼_𝑖𝑣𝑡−𝜏_𝑖.(4.10) It is easy to see that if 𝑥^∗𝑖,𝑦^∗𝑖,𝑣^∗>0,𝑖=1,…,𝑛, then 𝑑𝑊₂/𝑑𝑡≤0 for all (𝑥_𝑖,𝑦_𝑖,𝑣)>0. Clearly, the singleton {𝐸₁} is the only invariant set in {𝑑𝑊₂/𝑑𝑡=0}. LaSalle’s invariance principle implies global stability of 𝐸₁.

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 𝑛 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 𝑅₀ is less than unity, then the uninfected steady state is GAS and if 𝑅₀>1 (or the infected steady state exists) then the infected steady state is GAS for all time delays.