Global Dynamics of Virus Infection Model with Antibody Immune Response and Distributed Delays

Abstract

We present qualitative behavior of virus infection model with antibody immune response. The incidence rate of infection is given by saturation functional response. Two types of distributed delays are incorporated into the model to account for the time delay between the time when uninfected cells are contacted by the virus particle and the time when emission of infectious (matures) virus particles. Using the method of Lyapunov functional, we have established that the global stability of the steady states of the model is determined by two threshold numbers, the basic reproduction number 𝑅₀ and antibody immune response reproduction number 𝑅₁. We have proven that if 𝑅₀≤1, then the uninfected steady state is globally asymptotically stable (GAS), if 𝑅₁≤1<𝑅₀, then the infected steady state without antibody immune response is GAS, and if 𝑅₁>1, then the infected steady state with antibody immune response is GAS.

1. Introduction

In the past ten years there has been a growing interest in modeling viral infections for the study and characterization of host infection dynamics. The mathematical models, based on biological interactions, present a framework which can be used to obtain new insights into the viral dynamics and to interpret experimental data. Many authors have formulated mathematical models to describe the population dynamics of several viruses such as, human immunodeficiency virus (HIV) [1–14], hepatitis B virus (HBV) [15, 16], and hepatitis C virus (HCV) [17, 18]. During viral infections, the host immune system reacts with antigen-specific immune response. In particular, B cells play a critical role in antiviral defense by attacking free virus particles by making antibodies to clear antigens circulating in blood and lymph. The antibody immune response is more effective than the cell-mediated immune in some diseases like in malaria infection [19].

Mathematical models for virus dynamics with the antibody immune response have been developed in [20–26]. The basic virus dynamics model with antibody immune response was introduced by Murase et al. [21] as 𝑥̇(𝑡)(𝑡)(𝑡)𝑣(𝑡),𝑦(̇𝑡)(𝑡)𝑣(𝑡)(𝑡),𝑣(̇𝑡)(𝑡)(𝑡)(𝑡)𝑧(𝑡),𝑧(̇𝑡)(𝑡)𝑧(𝑡)(𝑡),=𝜆−𝑑𝑥−𝛽𝑥=𝛽𝑥−𝛿𝑦=𝑁𝛿𝑦−𝑐𝑣−𝑞𝑣=𝑟𝑣−𝜇𝑧 (1) where 𝑥(𝑡), 𝑦(𝑡), 𝑣(𝑡), and 𝑧(𝑡) represent the populations of uninfected cells, infected cells, virus, and B cells at time 𝑡, respectively; 𝜆 and 𝑑 are the recruited rate and death rate constants of uninfected cells, respectively; 𝛽 is the infection rate constant; 𝑁 is the number of free virus produced during the average infected cell life span; 𝛿 is the death rate constant of infected cells; 𝑐 is the clearance rate constant of the virus; 𝑟 and 𝜇 are the recruited rate and death rate constants of B cells; and 𝑞 is the B cells neutralizion rate. Model (1) is based on the assumption that the infection could occur and that the viruses are produced from infected cells instantaneously, once the uninfected cells are contacted by the virus particles. Other accurate models incorporate the delay between the time viral entry into the uninfected cell and the time the production of new virus particles, modeled with discrete time delay or distributed time delay using functional differential equations (see e.g., [7–14, 27]). In [7–14, 27], the viral infection models are presented without taking into consideration the antibody immune response. In [25, 28], global stability of viral infection models with antibody immune response and discrete delays has been studied.

In model (1) the infection rate is assumed to be bilinear in 𝑥 and 𝑣 however, this bilinear incidence rate associated with the mass action principle is insufficient to describe the infection process in detail [29, 30]. For example, a less than linear response in 𝑣 could occur due to saturation at high virus concentration, where the infectious fraction is high so that exposure is very likely. Thus, it is reasonable to assume that the infection rate is given by saturation functional response [31]. In [26], a virus infection model with antibody immune response and with saturation incidence rate has been considered. However, the time delay was not considered.

In this paper, we assume that the infection rate is given by saturation functional response. We incorporate two types of distributed delays into the model to account for the time delay between the time when uninfected cells are contacted by the virus particle and the time of emission of infectious (matures) virus particles. The global stability of this model is established using Lyapunov functionals, which are similar in nature to those used in [32]. We prove that the global dynamics of this model is determined by the basic reproduction number 𝑅₀ and antibody immune response reproduction number 𝑅₁. If 𝑅₀≤1, then the uninfected steady state is globally asymptotically stable (GAS), if 𝑅₁≤1<𝑅₀, then the infected steady state without antibody immune response is GAS, and if 𝑅₁>1, then the infected steady state with antibody immune response is GAS.

2. The Model

In this section we propose a delay mathematical model of viral infection with saturation functional response which describes the interaction of the virus with uninfected and infected cells, taking into account the effect of antibody immune response. Consider 𝑥̇(𝑡)(𝑡)−(𝑡)𝑣(𝑡)=𝜆−𝑑𝑥𝛽𝑥(𝑡),𝑦(̇𝑡)1+𝛼𝑣=𝛽^ℎ₀𝑓(𝜏)𝑒^−𝑚𝜏𝑥(𝑡−𝜏)𝑣(𝑡−𝜏)()(𝑡),𝑣(̇𝑡)1+𝛼𝑣𝑡−𝜏𝑑𝜏−𝛿𝑦=𝑁𝛿^𝜔₀𝑔(𝜏)𝑒^−𝑛𝜏𝑦()(𝑡)(𝑡)𝑧(𝑡),𝑧(̇𝑡)(𝑡)𝑧(𝑡)(𝑡),𝑡−𝜏𝑑𝜏−𝑐𝑣−𝑞𝑣=𝑟𝑣−𝜇𝑧 (2)(3)(4)(5) where 𝛼 is a positive constant, and all the variables and parameters of the model have the same meanings as given in model (1). To account for the time lag between viral contacting a target cell and the production of new virus particles, two distributed intracellular delays are introduced. It is assumed that the uninfected cells that are contacted by the virus particles at time 𝑡−𝜏 become infected cells at time 𝑡, where 𝜏 is a random variable with a probability distribution 𝑓(𝜏) over the interval [0,ℎ] and ℎ is limit superior of this delay. The factor 𝑒^−𝑚𝜏 accounts for the probability of surviving the time period of delay, where 𝑚 is the death rate of infected cells but not yet virus producer cells. On the other hand, it is assumed that a cell infected at time 𝑡−𝜏 starts to yield new infectious virus at time 𝑡 where 𝜏 is distributed according to a probability distribution 𝑔(𝜏) over the interval [0,𝜔] and 𝜔 is limit superior of this delay. The factor 𝑒^−𝑛𝜏 accounts for the probability of surviving during the time period of delay, where 𝑛 is constant. All the parameters are supposed to be positive.

The probability distribution functions 𝑓(𝜏) and 𝑔(𝜏) are assumed to satisfy 𝑓(𝜏)>0 and 𝑔(𝜏)>0, and ^ℎ₀𝑓(𝜏)𝑑𝜏=^𝜔₀𝑔(𝜏)𝑑𝜏=1,_ℎ0𝑓(𝑢)𝑒^𝑠𝑢𝑑𝑢<∞,^𝜔₀𝑔(𝑢)𝑒^𝑠𝑢𝑑𝑢<∞, (6) where 𝑠 is a positive number. Then 0<^ℎ₀𝑓(𝜏)𝑒^−𝑚𝜏𝑑𝜏≤1,for𝑚≥0,0<^𝜔₀𝑔(𝜏)𝑒^−𝑛𝜏𝑑𝜏≤1,for𝑛≥0. (7) Let 𝐹=∫^ℎ₀𝑓(𝜏)𝑒^−𝑚𝜏𝑑𝜏 and 𝐺=∫^𝜔₀𝑔(𝜏)𝑒^−𝑛𝜏𝑑𝜏. The initial conditions for system (2)–(5) take the form 𝑥(𝜃)=𝜑₁(𝜃),𝑦(𝜃)=𝜑₂(𝜃),𝑣(𝜃)=𝜑₃(𝜃),𝑧(𝜃)=𝜑₄(𝜃),𝜑_𝑗(𝜃)𝜑≥0,𝜃∈−𝜌,0,𝑗=1,…,4,_𝑗(0)>0,𝑗=1,…,4, (8) where 𝜌=max{ℎ,𝜔},(𝜑₁(𝜃),𝜑₂(𝜃),…,𝜑₄(𝜃))∈𝐶([−𝜌,0],ℝ⁴₊), where 𝐶([−𝜌,0],ℝ⁴₊) is the Banach space of continuous functions mapping the interval [−𝜌,0]into ℝ⁴₊. By the fundamental theory of functional differential equations [33], system (2)–(5) has a unique solution satisfying the initial conditions (8).

2.1. Nonnegativity and Boundedness of Solutions

In the following, we establish the nonnegativity and boundedness of solutions of (2)–(5) with initial conditions (8).

Proposition 1. Let (𝑥(𝑡),𝑦(𝑡),𝑣(𝑡),𝑧(𝑡)) be any solution of (2)–(5) satisfying the initial conditions (8); then 𝑥(𝑡),𝑦(𝑡),𝑣(𝑡), and 𝑧(𝑡) are all nonnegative for 𝑡≥0 and ultimately bounded.

Proof. First, we prove that 𝑥(𝑡)>0 for all 𝑡≥0. Assume that 𝑥(𝑡) loses its non-negativity on some local existence interval [0,ℓ] for some constant ℓ, and let 𝑡₁∈[0,ℓ] be such that 𝑥(𝑡₁)=0. From (2) we have 𝑥̇(𝑡₁)=𝜆>0. Hence, 𝑥(𝑡)>0 for some 𝑡∈(𝑡₁,𝑡₁+𝜀), where 𝜀>0 is sufficiently small. This leads to contradiction, and hence 𝑥(𝑡)>0 for all 𝑡≥0. Now from (3), (4), and (5) we have 𝑦(𝑡)=𝑦(0)𝑒^−𝛿𝑡+𝛽^𝑡₀𝑒^{−𝛿(𝑡−𝜂)}×^ℎ₀𝑓(𝜏)𝑒^−𝑚𝜏𝑥𝜂−𝜏𝑣𝜂−𝜏𝑣(𝑡)(0)𝑒1+𝛼𝑣𝜂−𝜏𝑑𝜏𝑑𝜂,=𝑣^{−∫𝑡0(𝑐+𝑞𝑧(𝜉))𝑑𝜉}+𝑁𝛿^𝑡₀𝑒^{−∫𝑡𝜂(𝑐+𝑞𝑧(𝜉))𝑑𝜉}×^𝜔₀𝑔(𝜏)𝑒^−𝑛𝜏𝑦𝑧(𝑡)(0)𝑒𝜂−𝜏𝑑𝜏𝑑𝜂,=𝑧^{−∫𝑡0(𝜇−𝑟𝑣(𝜉))𝑑𝜉}, (9) confirming that 𝑦(𝑡)≥0,𝑣(𝑡)≥0, and 𝑧(𝑡)≥0 for all 𝑡∈[0,𝜌]. By a recursive argument, we obtain 𝑦(𝑡)≥0,𝑣(𝑡)≥0, and 𝑧(𝑡)≥0 for all 𝑡≥0.
Next we show the boundedness of the solutions. From (2) we have 𝑥̇(𝑡)≤𝜆−𝑑𝑥(𝑡). This implies limsup_𝑡→∞𝑥(𝑡)≤𝜆/𝑑. Let 𝑋(𝑡)=∫^ℎ₀𝑓(𝜏)𝑒^−𝑚𝜏𝑥(𝑡−𝜏)𝑑𝜏+𝑦(𝑡); then 𝑋̇(𝑡)=^ℎ₀𝑓(𝜏)𝑒^−𝑚𝜏()−()𝑣()𝜆−𝑑𝑥𝑡−𝜏𝛽𝑥𝑡−𝜏𝑡−𝜏()+1+𝛼𝑣𝑡−𝜏𝑑𝜏_ℎ0𝑓(𝜏)𝑒^−𝑚𝜏𝛽𝑥(𝑡−𝜏)𝑣(𝑡−𝜏)()(𝑡)1+𝛼𝑣𝑡−𝜏𝑑𝜏−𝛿𝑦=𝜆^ℎ₀𝑓(𝜏)𝑒^−𝑚𝜏𝑑𝜏−𝑑×^ℎ₀𝑓(𝜏)𝑒^−𝑚𝜏𝑥(𝑡−𝜏)𝑑𝜏−𝛿𝑦(𝑡)≤𝜆^ℎ₀𝑓(𝜏)𝑒^−𝑚𝜏𝑑𝜏−𝜎₁×^ℎ₀𝑓(𝜏)𝑒^−𝑚𝜏𝑥(𝑡−𝜏)𝑑𝜏+𝑦(𝑡)≤𝜆−𝜎₁𝑋(𝑡), (10) where 𝜎₁=min{𝑑,𝛿}. Hence, limsup_𝑡→∞𝑋(𝑡)≤𝐿₁, where 𝐿₁=𝜆/𝜎₁. Since ∫^ℎ₀𝑓(𝜏)𝑒^−𝑚𝜏𝑥(𝑡−𝜏)𝑑𝜏>0, then limsup_𝑡→∞𝑦(𝑡)≤𝐿₁. On the other hand, let 𝑍(𝑡)=𝑣(𝑡)+(𝑞/𝑟)𝑧(𝑡); then 𝑍̇(𝑡)=𝑁𝛿^𝜔₀𝑔(𝜏)𝑒^−𝑛𝜏𝑦(𝑡−𝜏)𝑑𝜏−𝑐𝑣(𝑡)−𝑞𝜇𝑟𝑧(𝑡)≤𝑁𝛿𝐿₁^𝜔₀𝑔(𝜏)𝑒^−𝑛𝜏𝑑𝜏−𝜎₂𝑣(𝑡)+𝑞𝑟𝑧(𝑡)=𝑁𝛿𝐿₁^𝜔₀𝑔(𝜏)𝑒^−𝑛𝜏𝑑𝜏−𝜎₂𝑍(𝑡)≤𝑁𝛿𝐿₁−𝜎₂𝑍(𝑡), (11) where 𝜎₂=min{𝑐,𝜇}. Hence, limsup_𝑡→∞𝑍(𝑡)≤𝐿₂, where 𝐿₂=𝑁𝛿𝐿₁/𝜎₂. Since 𝑣(𝑡)≥0, and 𝑣(𝑡)≥0, then limsup_𝑡→∞𝑣(𝑡)≤𝐿₂ and limsup_𝑡→∞𝑧(𝑡)≤𝐿₂.
Therefore, 𝑥(𝑡),𝑦(𝑡),𝑣(𝑡), and 𝑧(𝑡) are ultimately bounded.

2.2. Steady States

We define the basic reproduction number for system (2)–(5) as 𝑅₀=𝑁𝐹𝐺𝛽𝜆,𝑐𝑑 (12) and the antibody immune response reproduction number 𝑅₁=𝑅₀.1+𝜇𝛼𝑑+𝜇𝛽/𝑟𝑑 (13) Clearly 𝑅₁<𝑅₀. It can be seen that system (2)–(5) has an uninfected steady state 𝐸₀=(𝑥₀,0,0,0), where 𝑥₀=𝜆/𝑑. In addition to 𝐸₀, the system has an infected steady state without immune response 𝐸₁(𝑥₁,𝑦₁,𝑣₁,0) and infected steady state with immune response 𝐸₂(𝑥₂,𝑦₂,𝑣₂,𝑧₂), where 𝑥₁=𝑥₀1+𝛼𝑣₁𝑅₀,𝑦₁=𝑐𝑣𝑁𝛿𝐺¹,𝑣₁=𝑑𝑅𝛼𝑑+𝛽⁰,𝑥−1₂=𝜆𝑑+𝛽𝑣₂/1+𝛼𝑣₂,𝑦₂=𝐹𝛽𝜆𝑣₂𝛿𝑑1+𝛼𝑣₂+𝛽𝑣₂,𝑣₂=𝜇𝑟,𝑧₂=𝑐𝑞𝑅₁−1. (14) From the above we have the following:(i)if 𝑅₀>1, then there exists a positive steady state 𝐸₁(𝑥₁,𝑦₁,𝑣₁,0);(ii)if 𝑅₁>1, then there exists a positive steady state 𝐸₂(𝑥₂,𝑦₂,𝑣₂,𝑧₂).

3. Global Stability

In this section, we prove the global stability of the steady states of system (2)–(5) employing the method of Lyapunov functional which is used in [32] for SIR epidemic model with distributed delay. Next we will use the following notation: 𝑏=𝑏(𝑡) for any 𝑏∈{𝑥,𝑦,𝑣,𝑧}. We also define a function 𝐻∶(0,∞)→[0,∞) as 𝐻(𝑏)=𝑏−1−ln(𝑏). (15) It is clear that 𝐻(𝑏)≥0 for any 𝑏>0 and 𝐻 has the global minimum 𝐻(1)=0.

Theorem 2. If 𝑅₀≤1, then 𝐸₀ is GAS.

Proof. Define a Lyapunov functional 𝑊₀ as follows: 𝑊₀=𝑁𝐹𝐺𝑥₀𝐻𝑥𝑥₀+1𝐹𝛽𝑦+𝐹×^ℎ₀𝑓(𝜏)𝑒^−𝑚𝜏×^𝜏₀𝑥(𝑡−𝜃)𝑣(𝑡−𝜃)()+𝛿1+𝛼𝑣𝑡−𝜃𝑑𝜃𝑑𝜏𝐹𝐺^𝜔₀𝑔(𝜏)𝑒^{−𝑛𝜏𝜏}₀𝑦()𝑞𝑡−𝜃𝑑𝜃𝑑𝜏+𝑣+𝑟𝑧. (16) The time derivative of 𝑊₀ along the trajectories of (2)–(5) satisfies 𝑑𝑊₀𝑥𝑑𝑡=𝑁𝐹𝐺1−⁰𝑥𝜆−𝑑𝑥−𝛽𝑥𝑣+𝛽1+𝛼𝑣𝐹^ℎ₀𝑓(𝜏)𝑒^−𝑚𝜏𝑥(𝑡−𝜏)𝑣(𝑡−𝜏)()𝛿1+𝛼𝑣𝑡−𝜏𝑑𝜏−𝐹𝑦+𝛽𝐹^ℎ₀𝑓(𝜏)𝑒^−𝑚𝜏×𝑥𝑣−𝑥()𝑣()1+𝛼𝑣𝑡−𝜏𝑡−𝜏()+𝛿1+𝛼𝑣𝑡−𝜏𝑑𝜏𝐹𝐺^𝜔₀𝑔(𝜏)𝑒^−𝑛𝜏𝑦−𝑦(𝑡−𝜏)𝑑𝜏+𝑁𝛿^𝜔₀𝑔(𝜏)𝑒^−𝑛𝜏𝑦(𝑡−𝜏)𝑑𝜏−𝑐𝑣−𝑞𝑣𝑧+𝑞𝑣𝑧−𝑞𝜇𝑟𝑧=𝑁𝐹𝐺−𝑑𝑥−𝑥₀²𝑥+𝛽𝑥₀𝑣1+𝛼𝑣−𝑐𝑣−𝑞𝜇𝑟𝑧=−𝑁𝐹𝐺𝑑𝑥−𝑥₀²𝑥𝑅−𝑐𝑣+𝑐𝑣⁰(1+𝛼𝑣−𝛼𝑣)−1+𝛼𝑣𝑞𝜇𝑟𝑧=−𝑁𝐹𝐺𝑑𝑥𝑥−𝑥₀²+𝑐𝑣𝑅₀−𝑅−1₀𝑐𝛼𝑣²−1+𝛼𝑣𝑞𝜇𝑟𝑧. (17) If 𝑅₀≤1, then 𝑑𝑊₀/𝑑𝑡≤0 for all 𝑥,𝑣,𝑧>0. By Theorem 5.3.1 in [33], the solutions of system (2)–(5) limit to 𝑀, the largest invariant subset of {𝑑𝑊₀/𝑑𝑡=0}. Clearly, it follows from (17) that 𝑑𝑊₀/𝑑𝑡=0 if and only if 𝑥=𝑥₀, 𝑣=0, and 𝑧=0. Noting that 𝑀 is invariant for each element of 𝑀, we have 𝑣=0, and 𝑧=0, and then 𝑣̇=0. From (4) we derive that 0=𝑣̇=𝑁𝛿^𝜔₀𝑔(𝜏)𝑒^−𝑛𝜏𝑦(𝑡−𝜏)𝑑𝜏. (18) This yields 𝑦=0. Hence, 𝑑𝑊₀/𝑑𝑡=0 if and only if 𝑥=𝑥₀, 𝑦=0, 𝑣=0, and 𝑧=0. From LaSalle’s Invariance Principle, 𝐸₀ is GAS.

Theorem 3. If 𝑅₁≤1<𝑅₀, then 𝐸₁ is GAS.

Proof. We construct the following Lyapunov functional: 𝑊₁=𝑁𝐹𝐺𝑥₁𝐻𝑥𝑥₁+1𝐹𝑦₁𝐻𝑦𝑦₁+1𝐹𝛽𝑥₁𝑣₁1+𝛼𝑣₁×^ℎ₀𝑓(𝜏)𝑒^−𝑚𝜏×^𝜏₀𝐻𝑥(𝑡−𝜃)𝑣(𝑡−𝜃)1+𝛼𝑣₁𝑥₁𝑣₁((+1+𝛼𝑣𝑡−𝜃))𝑑𝜃𝑑𝜏𝛿𝑦₁𝐹𝐺^𝜔₀𝑔(𝜏)𝑒^{−𝑛𝜏𝜏}₀𝐻𝑦(𝑡−𝜃)𝑦₁𝑑𝜃𝑑𝜏+𝑣₁𝐻𝑣𝑣₁+𝑞𝑟𝑧. (19) The time derivative of 𝑊₁ along the trajectories of (2)–(5) is given by 𝑑𝑊₁𝑥𝑑𝑡=𝑁𝐹𝐺1−¹𝑥𝜆−𝑑𝑥−𝛽𝑥𝑣+11+𝛼𝑣𝐹𝑦1−¹𝑦×𝛽^ℎ₀𝑓(𝜏)𝑒^−𝑚𝜏𝑥(𝑡−𝜏)𝑣(𝑡−𝜏)()+𝛽1+𝛼𝑣𝑡−𝜏𝑑𝜏−𝛿𝑦𝐹^ℎ₀𝑓(𝜏)𝑒^−𝑚𝜏×𝑥𝑣−𝑥()𝑣()1+𝛼𝑣𝑡−𝜏𝑡−𝜏()+𝑥1+𝛼𝑣𝑡−𝜏¹𝑣₁1+𝛼𝑣₁𝑥()𝑣()×ln𝑡−𝜏𝑡−𝜏)(1+𝛼𝑣((+𝛿𝑥𝑣1+𝛼𝑣𝑡−𝜏))𝑑𝜏𝐹𝐺^𝜔₀𝑔(𝜏)𝑒^−𝑛𝜏×𝑦−𝑦(𝑡−𝜏)+𝑦₁𝑦()ln𝑡−𝜏𝑦+𝑣𝑑𝜏1−¹𝑣𝑁𝛿^𝜔₀𝑔(𝜏)𝑒^−𝑛𝜏𝑦(𝑡−𝜏)𝑑𝜏−𝑐𝑣−𝑞𝑧𝑣+𝑞𝑣𝑧−𝑞𝜇𝑟𝑧. (20) Using the steady state conditions for 𝐸₁: 𝜆=𝑑𝑥₁+𝛽𝑥₁𝑣₁1+𝛼𝑣₁,𝛿𝑦₁=𝐹𝛽𝑥₁𝑣₁1+𝛼𝑣₁,𝑐𝑣₁=𝑁𝛿𝐺𝑦₁, (21) we have 𝑑𝑊₁𝑑𝑡=𝑁𝐹𝐺−𝑑𝑥−𝑥₁²𝑥+𝛽𝑥₁𝑣₁1+𝛼𝑣₁−𝛽𝑥₁𝑣₁1+𝛼𝑣₁𝑥₁𝑥+𝛽𝑥₁𝑣−11+𝛼𝑣𝐹𝛽𝑥₁𝑣₁1+𝛼𝑣₁×^ℎ₀𝑓(𝜏)𝑒^−𝑚𝜏×𝑦₁𝑥(𝑡−𝜏)𝑣(𝑡−𝜏)1+𝛼𝑣₁𝑦𝑥₁𝑣₁((+𝛿1+𝛼𝑣𝑡−𝜏))𝑑𝜏𝐹𝑦₁+1𝐹𝛽𝑥₁𝑣₁1+𝛼𝑣₁×^ℎ₀𝑓(𝜏)𝑒^−𝑚𝜏𝑥()𝑣()×ln𝑡−𝜏𝑡−𝜏)(1+𝛼𝑣((+𝑥𝑣1+𝛼𝑣𝑡−𝜏))𝑑𝜏𝛿𝑦₁𝐹𝐺^𝜔₀𝑔(𝜏)𝑒^−𝑛𝜏𝑦()ln𝑡−𝜏𝑦𝑑𝜏+𝑁𝛿𝑦₁^𝜔₀𝑔(𝜏)𝑒^−𝑛𝜏𝑣₁𝑦(𝑡−𝜏)𝑣𝑦₁𝑑𝜏−𝑐𝑣+𝑐𝑣₁+𝑞𝑣₁𝑧−𝑞𝜇𝑟𝑧. (22) Using the following equalities: 𝑥()𝑣()ln𝑡−𝜏𝑡−𝜏)(1+𝛼𝑣((𝑦𝑥𝑣1+𝛼𝑣𝑡−𝜏))=ln¹𝑥(𝑡−𝜏)𝑣(𝑡−𝜏)1+𝛼𝑣₁𝑦𝑥₁𝑣₁((𝑥1+𝛼𝑣𝑡−𝜏))+ln¹𝑥𝑣+ln¹𝑦𝑣𝑦₁+ln1+𝛼𝑣1+𝛼𝑣₁,𝑦()ln𝑡−𝜏𝑦=ln𝑣𝑦₁𝑣₁𝑦𝑣+ln¹𝑦(𝑡−𝜏)𝑣𝑦₁, (23) then collecting terms of (22), we obtain 𝑑𝑊₁𝑑𝑡=𝑁𝐹𝐺−𝑑𝑥−𝑥₁²𝑥−𝛿𝑦₁𝐹𝑥₁𝑥𝑥−1−ln¹𝑥+𝛿𝑦₁𝐹𝑣1+𝛼𝑣₁𝑣₁()−1+𝛼𝑣𝛿𝑦₁𝐹²×^ℎ₀𝑓(𝜏)𝑒^−𝑚𝜏×𝑦₁𝑥(𝑡−𝜏)𝑣(𝑡−𝜏)1+𝛼𝑣₁𝑦𝑥₁𝑣₁((𝑦1+𝛼𝑣𝑡−𝜏))−1−ln¹𝑥(𝑡−𝜏)𝑣(𝑡−𝜏)1+𝛼𝑣₁𝑦𝑥₁𝑣₁((−1+𝛼𝑣𝑡−𝜏))𝑑𝜏𝛿𝑦₁𝐹𝐺^𝑤₀𝑔(𝜏)𝑒^−𝑛𝜏×𝑣₁𝑦(𝑡−𝜏)𝑣𝑦₁𝑣−1−ln¹𝑦(𝑡−𝜏)𝑣𝑦₁+𝑑𝜏𝛿𝑦₁𝐹ln1+𝛼𝑣1+𝛼𝑣₁−𝑣𝑣₁+𝑞𝑧𝑣₁−𝜇𝑟=𝑁𝐹𝐺−𝑑𝑥−𝑥₁²𝑥+𝛿𝑦₁𝐹×𝑣−1+1+𝛼𝑣₁𝑣₁()−𝑣1+𝛼𝑣𝑣₁+1+𝛼𝑣1+𝛼𝑣₁−𝛿𝑦₁𝐹𝐻𝑥₁𝑥−𝛿𝑦₁𝐹²×^ℎ₀𝑓(𝜏)𝑒^−𝑚𝜏𝑦×𝐻¹𝑥(𝑡−𝜏)𝑣(𝑡−𝜏)1+𝛼𝑣₁𝑦𝑥₁𝑣₁((−𝛿1+𝛼𝑣𝑡−𝜏))𝑑𝜏𝑦𝐹𝐺^1𝜔₀𝑔(𝜏)𝑒^−𝑛𝜏𝐻𝑣₁𝑦(𝑡−𝜏)𝑣𝑦₁−𝑑𝜏𝛿𝑦₁𝐹𝐻1+𝛼𝑣1+𝛼𝑣₁+𝑞𝑧𝑣₁−𝜇𝑟=𝑁𝐹𝐺−𝑑𝑥−𝑥₁²𝑥−𝛿𝑦₁𝐹𝛼𝑣−𝑣₁²𝑣₁(1+𝛼𝑣)1+𝛼𝑣₁−𝛿𝑦₁𝐹𝐻𝑥₁𝑥−𝛿𝑦₁𝐹²×^ℎ₀𝑓(𝜏)𝑒^−𝑚𝜏𝑦×𝐻¹𝑥(𝑡−𝜏)𝑣(𝑡−𝜏)1+𝛼𝑣₁𝑦𝑥₁𝑣₁((−𝛿1+𝛼𝑣𝑡−𝜏))𝑑𝜏𝑦𝐹𝐺^1𝜔₀𝑔(𝜏)𝑒^−𝑛𝜏𝐻𝑣₁𝑦(𝑡−𝜏)𝑣𝑦₁−𝑑𝜏𝛿𝑦₁𝐹𝐻1+𝛼𝑣1+𝛼𝑣₁+𝑞𝑑𝛼𝑑+𝛽1+𝜇𝛼𝑑+𝜇𝛽𝑅𝑟𝑑¹−1𝑧. (24) Thus, if 𝑅₀>1, then 𝑥₁,𝑦₁, and 𝑣₁>0, and hence, if 𝑅₁≤1, then 𝑑𝑊₁/𝑑𝑡≤0 for all 𝑥,𝑦,𝑣>0. By Theorem 5.3.1 in [33], the solutions of system (2)–(5) limit to 𝑀, the largest invariant subset of {𝑑𝑊₁/𝑑𝑡=0}. It can be seen that 𝑑𝑊₁/𝑑𝑡=0 if and only if 𝑥=𝑥₁,𝑣=𝑣₁, and 𝐻=0; that is, 𝑦₁𝑥(𝑡−𝜏)𝑣(𝑡−𝜏)1+𝛼𝑣₁𝑦𝑥₁𝑣₁((=𝑣1+𝛼𝑣𝑡−𝜏))¹𝑦(𝑡−𝜏)𝑣𝑦₁.=1foralmostall𝜏∈0,𝜌 (25) From (25), if 𝑣=𝑣₁, then 𝑦=𝑦₁, and hence 𝑑𝑊₁/𝑑𝑡 equals zero at 𝐸₁. LaSalle’s Invariance Principle implies global stability of 𝐸₁.

Theorem 4. If 𝑅₁>1, then 𝐸₂ is GAS.

Proof. We construct the following Lyapunov functional: 𝑊₂=𝑁𝐹𝐺𝑥₂𝐻𝑥𝑥₂+1𝐹𝑦₂𝐻𝑦𝑦₂+1𝐹𝛽𝑥₂𝑣₂1+𝛼𝑣₂×^ℎ₀𝑓(𝜏)𝑒^−𝑚𝜏×^𝜏₀𝐻𝑥(𝑡−𝜃)𝑣(𝑡−𝜃)1+𝛼𝑣₂𝑥₂𝑣₂((+1+𝛼𝑣𝑡−𝜃))𝑑𝜃𝑑𝜏𝛿𝑦₂𝐹𝐺^𝜔₀𝑔(𝜏)𝑒^{−𝑛𝜏𝜏}₀𝐻𝑦(𝑡−𝜃)𝑦₂𝑑𝜃𝑑𝜏+𝑣₂𝐻𝑣𝑣₂+𝑞𝑟𝑧₂𝐻𝑧𝑧₂. (26) The time derivative of 𝑊₂ along the trajectories of (2)–(5) is given by 𝑑𝑊₂𝑥𝑑𝑡=𝑁𝐹𝐺1−²𝑥𝜆−𝑑𝑥−𝛽𝑥𝑣+11+𝛼𝑣𝐹𝑦1−²𝑦×𝛽^ℎ₀𝑓(𝜏)𝑒^−𝑚𝜏𝑥(𝑡−𝜏)𝑣(𝑡−𝜏)()+𝛽1+𝛼𝑣𝑡−𝜏𝑑𝜏−𝛿𝑦𝐹^ℎ₀𝑓(𝜏)𝑒^−𝑚𝜏×𝑥𝑣−𝑥()𝑣()1+𝛼𝑣𝑡−𝜏𝑡−𝜏()+𝑥1+𝛼𝑣𝑡−𝜏²𝑣₂1+𝛼𝑣₂𝑥()𝑣()×ln𝑡−𝜏𝑡−𝜏)(1+𝛼𝑣((+𝛿𝑥𝑣1+𝛼𝑣𝑡−𝜏))𝑑𝜏𝐹𝐺^𝜔₀𝑔(𝜏)𝑒^−𝑛𝜏×𝑦−𝑦(𝑡−𝜏)+𝑦₂𝑦()ln𝑡−𝜏𝑦+𝑣𝑑𝜏1−²𝑣×𝑁𝛿^𝜔₀𝑔(𝜏)𝑒^−𝑛𝜏𝑦()+𝑧𝑡−𝜏𝑑𝜏−𝑐𝑣−𝑞𝑧𝑣1−²𝑧𝑞𝑣𝑧−𝑞𝜇𝑟𝑧. (27) Using the steady state conditions for 𝐸₂: 𝜆=𝑑𝑥₂𝑥+𝛽²𝑣₂1+𝛼𝑣₂,𝑐𝑣₂𝛿=𝑁𝐹𝐺𝐹𝑦₂−𝑞𝑣₂𝑧₂,𝜇=𝑟𝑣₂, (28) we obtain 𝑑𝑊₂𝑑𝑡=𝑁𝐹𝐺−𝑑𝑥−𝑥₂²𝑥+𝛽𝑥₂𝑣₂1+𝛼𝑣₂−𝛽𝑥₂𝑣₂1+𝛼𝑣₂𝑥₂𝑥+𝛽𝑥₂𝑣−11+𝛼𝑣𝐹𝛽𝑥₂𝑣₂1+𝛼𝑣₂×^ℎ₀𝑓(𝜏)𝑒^−𝑚𝜏×𝑦₂𝑥(𝑡−𝜏)𝑣(𝑡−𝜏)1+𝛼𝑣₂𝑦𝑥₂𝑣₂((+𝛿1+𝛼𝑣𝑡−𝜏))𝑑𝜏𝐹𝑦₂+1𝐹𝛽𝑥₂𝑣₂1+𝛼𝑣₂×^ℎ₀𝑓(𝜏)𝑒^−𝑚𝜏𝑥()𝑣()×ln𝑡−𝜏𝑡−𝜏)(1+𝛼𝑣((+𝑥𝑣1+𝛼𝑣𝑡−𝜏))𝑑𝜏𝛿𝑦₂𝐹𝐺^𝜔₀𝑔(𝜏)𝑒^−𝑛𝜏𝑦()ln𝑡−𝜏𝑦𝑑𝜏−𝑁𝛿𝑦₂^𝜔₀𝑔(𝜏)𝑒^−𝑛𝜏𝑣₂𝑦(𝑡−𝜏)𝑣𝑦₂𝑑𝜏−𝑐𝑣+𝑐𝑣₂+𝑞𝑣₂𝑧−𝑞𝑣𝑧₂−𝑞𝜇𝑟𝑧+𝑞𝜇𝑟𝑧₂. (29) Using the following equalities: 𝑐𝑣=𝑐𝑣₂𝑣𝑣₂𝛿=𝑁𝐹𝐺𝐹𝑦₂𝑣𝑣₂−𝑞𝑣𝑧₂,𝑥()𝑣()ln𝑡−𝜏𝑡−𝜏)(1+𝛼𝑣((𝑦𝑥𝑣1+𝛼𝑣𝑡−𝜏))=ln²𝑥(𝑡−𝜏)𝑣(𝑡−𝜏)1+𝛼𝑣₂𝑦𝑥₂𝑣₂((𝑥1+𝛼𝑣𝑡−𝜏))+ln²𝑥𝑣+ln²𝑦𝑣𝑦₂+ln1+𝛼𝑣1+𝛼𝑣₂,𝑦()ln𝑡−𝜏𝑦=ln𝑣𝑦₂𝑣₂𝑦𝑣+ln²𝑦(𝑡−𝜏)𝑣𝑦₂, (30) and collecting terms of (29), we obtain 𝑑𝑊₂𝑑𝑡=𝑁𝐹𝐺−𝑑𝑥−𝑥₂²𝑥−𝛿𝑦₁𝐹𝑥₂𝑥𝑥−1−ln²𝑥+𝛿𝑦₂𝐹𝑣1+𝛼𝑣₂𝑣₂()−1+𝛼𝑣𝛿𝑦₂𝐹²×^ℎ₀𝑓(𝜏)𝑒^−𝑚𝜏×𝑦₂𝑥(𝑡−𝜏)𝑣(𝑡−𝜏)1+𝛼𝑣₂𝑦𝑥₂𝑣₂((𝑦1+𝛼𝑣𝑡−𝜏))−1−ln²𝑥(𝑡−𝜏)𝑣(𝑡−𝜏)1+𝛼𝑣₂𝑦𝑥₂𝑣₂((−1+𝛼𝑣𝑡−𝜏))𝑑𝜏𝛿𝑦₂𝐹𝐺^𝜔₀𝑔(𝜏)𝑒^−𝑛𝜏𝑣₂𝑦(𝑡−𝜏)𝑣𝑦₂𝑣−1−ln²𝑦(𝑡−𝜏)𝑣𝑦₂+𝑑𝜏𝛿𝑦₂𝐹ln1+𝛼𝑣1+𝛼𝑣₂−𝑣𝑣₂=𝑁𝐹𝐺−𝑑𝑥−𝑥₂²𝑥−𝛿𝑦₂𝐹𝛼𝑣−𝑣₂²𝑣₂(1+𝛼𝑣)1+𝛼𝑣₂−𝛿𝑦₂𝐹𝐻𝑥₂𝑥−𝛿𝑦₂𝐹²×^ℎ₀𝑓(𝜏)𝑒^−𝑚𝜏𝑦×𝐻²𝑥(𝑡−𝜏)𝑣(𝑡−𝜏)1+𝛼𝑣₂𝑦𝑥₂𝑣₂((−𝛿1+𝛼𝑣𝑡−𝜏))𝑑𝜏𝑦𝐹𝐺^2𝜔₀𝑔(𝜏)𝑒^−𝑛𝜏𝐻𝑣₂𝑦(𝑡−𝜏)𝑣𝑦₂−𝑑𝜏𝛿𝑦₂𝐹𝐻1+𝛼𝑣1+𝛼𝑣₂. (31) Thus, if 𝑅₁>1, then 𝑥₂,𝑦₂,𝑣₂, and 𝑧₂>0 and 𝑑𝑊₂/𝑑𝑡≤0. By Theorem 5.3.1 in [33], the solutions of system (2)–(5) limit to 𝑀, the largest invariant subset of {𝑑𝑊₂/𝑑𝑡=0}. It can be seen that 𝑑𝑊₂/𝑑𝑡=0 if and only if 𝑥=𝑥₂,𝑣=𝑣₂, and 𝐻=0; that is, 𝑦₂𝑥(𝑡−𝜏)𝑣(𝑡−𝜏)1+𝛼𝑣₂𝑦𝑥₂𝑣₂((=𝑣1+𝛼𝑣𝑡−𝜏))²𝑦(𝑡−𝜏)𝑣𝑦₂=1+𝛼𝑣1+𝛼𝑣₂.=1foralmostall𝜏∈0,𝜌 (32) Then from (32) 𝑦=𝑦₂, and hence 𝑑𝑊₂/𝑑𝑡=0 at the steady state 𝐸₂. LaSalle’s Invariance Principle implies global stability of 𝐸₂.

4. Conclusion

In this paper, we have proposed a virus infection model which describes the interaction of the virus with the uninfected and infected cells taking into account the antibody immune response. The infection rate is given by saturation functional response. Two types of distributed time delays have been incorporated into the model to describe the time needed for infection of uninfected cell and virus replication. Using the method of Lyapunov functional, we have established that the global dynamics of the model is determined by two threshold parameters 𝑅₀ and 𝑅₁. The basic reproduction number viral infection 𝑅₀ determines whether a chronic infection can be established, and the antibody immune response reproduction number 𝑅₁ determines whether a persistent antibody immune response can be established. We have proven that if 𝑅₀≤1, then the uninfected steady state is GAS, and the viruses are cleared. If 𝑅₁≤1<𝑅₀, then the infected steady state without antibody immune response is GAS, and the infection becomes chronic but without persistent antibody immune response. If 𝑅₁>1, then the infected steady state with antibody immune response is GAS, and the infection is chronic with persistent antibody immune response. We note that the effect of the time delay appears in the parameters 𝐹 and 𝐺. Since 0<𝐹≤1 and 0<𝐺≤1, then the intracellular delay can reduce the parameters 𝑅₀ and 𝑅₁. As a consequence, ignoring the delay will produce overestimation of 𝑅₀.

Conflict of Interests

The authors declare that there is no conflict of interests regarding the publication of this paper.

Acknowledgments

This article was funded by the Deanship of Scientific Research (DSR), King Abdulaziz University, Jeddah, Saudi Arabia. The authors, therefore, acknowledge with thanks DSR technical and financial support. The authors are also grateful to Professor Victor Kozyakin and to the anonymous reviewers for constructive suggestions and valuable comments, which improved the quality of the article.