Global dynamics of a delay virus model with recruitment and saturation effects of immune responses.

In this paper, we formulate a virus dynamics model with the recruitment of immune responses, saturation effects and an intracellular time delay. With the help of uniform persistence theory and Lyapunov method, we show that the global stability of the model is totally determined by the basic reproductive number R0. Furthermore, we analyze the effects of the recruitment of immune responses on virus infection by numerical simulation. The results show ignoring the recruitment of immune responses will result in overestimation of the basic reproductive number and the severity of viral infection.


1.
Introduction. In recent years, virus dynamics attracts more and more attentions of researchers and plays a crucial role in many diseases research, including AIDS, hepatitis and influenza. Many mathematical models have provided insights into virus infection and dynamics, as well as on how an infection can be managed, reduced or even eradicated ( [3], [4], [7], [15], [17], [27], [38], [43], [44]). Since the basic three-dimensional viral infection model was proposed by Nowak et al. [21], Perelson et al. [26], Perelson and Nelson [25], Nowak and May [20], many people have established different within-host infection model, which help us to better understand virus infection and various drug therapy strategies by mathematical analysis, numerical simulations and clinical data ( [13], [19], [22], [28], [29]). Note that immune responses play a critical part in the process of viral infections. Concretely, cytotoxic T lymphocyte (CTL) cells can attack infected cells, and antibody cells can neutralize viruses. To better understand the role of the immune function during virus infection, Wodarz proposed the following model with both CTL and antibody immune responses [41], where a dot denotes the differentiation with respect to time t, T (t), I(t), V (t), C(t) and A(t) are the concentrations of healthy cells, infected cells, free virus, CTL cells and antibody cells at time t, respectively. d 1 , d 2 , d 3 , d 4 , d 5 are the death rates of healthy cell, infected cell, free virus, CTL cells and antibody cells, respectively. λ represents a constant production of the healthy cells. The term βT V represents the rate for the heathy T cells to be infected by virus. Furthermore, infected cells are killed by CTL cells at a rate pIC. Free virus are produced by infected cells at a rate rd 2 I, and are killed by antibody responses at a rate qAV . CTL immune responses are activated at a rate proportional to the abundance of CTL cells and infected cells, k 1 IC. Antibody responses are produced at a rate proportional to the abundance of antibodies and free virus, k 2 AV . Biologically, all parameters are positive.
After that, some researchers have taken into account the effect of immune responses including CTL responses or antibody responses ( [24], [35], [36], [37], [39]). Some other researchers have incorporated the effect of CTL responses and intracellular delays ( [11], [16], [18], [32], [45]). Concretely, the global dynamics of (1) with and without intracellular time delay is given in [24] and [42], respectively. Note that model (1) assumes that CTL and antibody responses are produced at bilinear rates. However, De Boer [5] pointed out that the bilinear rates cannot model several immune responses that are together controlling a chronic infection. In [5], De Boer has proposed an immune response function with the saturation. Incorporating the saturation effects of immune responses and the delay, [12] also obtained the global stability of the model, which is totally determined by the corresponding reproductive numbers. These results preclude the complicated behaviors such as the backward bifurcations and Hopf bifurcations which may be induced by saturation factors and time delay.
Note also that most of models assume CTL responses are activated by infected cells/antigenic stimulation, and antibody responses are activated by virus in these studies. However, as pointed out by Nowak and May [20], CTL responses have another function of self-regulating, i.e., the CTL responses are triggered by encountering foreign antigen and then adopts a constant level which is independent of the concentration of virions or infected cell. Bocharvor et al. have provided evidence the export of precursor CTL cells from the thymus [2]. Pang and Cui et al. have studied the export of specific precursor CTL cells from the thymus in [23], but they didn't considered intracellular time delay and antibody responses. Similarly, Wang and Wang have considered that neutralizing antibodies are produced at a constant rate after the injection [37], but they didn't take into account the effect of CTL responses and intracellular time delay.
Motivated by the above studies, we will formulate and analyze a virus dynamics model with the recruitment of immune responses, saturation effects of immune responses and an intracellular time delay, which can be described by the following functional differential equations: Here, we use λ 1 to describe the export of specific precursor CTL cells from the thymus, λ 2 to describe the recruitment rate of antibody responses. The production of infected cells is delayed in such a way that the production of new virus at time t depends on the population of virus and infected cells at a previous time t − τ , and only a fraction of e −sτ can survive after the interval τ , where 1/s is the average lifetime of infected cells without reproduction. h 1 , h 2 are saturation constants. Other parameters are same as that in model (1).
The main aim of the present paper is to explore the effects of the recruitment of immune responses on virus infection. The organization of this paper is as follows. In the next section, some preliminary analyzes of the model (2) will be given. Stability of all equilibria are given in Section 3. In Section 4, some numerical simulations are given to explain the effects of λ 1 , λ 2 , i.e., the term of recruitment of immune responses. Lastly, some brief conclusions are given in Section 5.

2.
Preliminary analyses of the model. In this section, we will first prove the positivity and boundedness of solutions, and then derive the expression of the basic reproduction number for model (2).

Positivity and boundedness of solutions. Let
The initial functions for model (2) are provided with φ ∈ X + = C([−τ, 0], R 5 + ). Proposition 1. Under the above initial conditions, all solutions of model (2) are nonnegative. In particular, the solution (T (t), Proof. We first verify that T (t) is positive in the existence interval of the solution. Suppose not. Then there is t 1 > 0 such that T (t 1 ) = 0 and T (t) > 0, t ∈ [0, t 1 ). Indeed, if T (0) = 0, we haveṪ (0) = λ > 0. Thus, T (t) > 0 for small positive t. Evidently, this remains valid if T (0) > 0. As a result, the existence of t 1 follows if the claim is not true. Note thatṪ (t 1 ) = λ > 0. Thus, there is a sufficiently small ε > 0, such that T (t) < 0 for all t ∈ (t 1 − ε, t 1 ). We get a contradiction because T (t) > 0, t ∈ [0, t 1 ). Hence, T (t) is positive in the existence interval of the solution. In the same way, we obtain C(t) and A(t) are positive in the existence interval of the solution. Then, we verify that I(t) and V (t) are positive. In the same way, we assume that there is a first time t 2 such that V (t 2 ) = 0, from the third equation of model (2) we haveV (t 2 ) = rd 2 I(t 2 ). By solving the second equation of model (2), we obtain It follows thatV (t 2 ) > 0, hence V (t) > 0. Furthermore, From the above expression of I(t) solution, we get I(t) > 0.
It follows easily that I(t) ≥ 0, V (t) ≥ 0 in the existence interval of the solution if the initial functions are in X + , and Proposition 2. All solutions of model (2) in X + are ultimately bounded.
Proof. Set Calculating the derivative of L along the solution of (2), we geṫ Since It follows that the nonnegative solutions of (2) exist on [0, ∞) and are ultimately bounded. Moreover, From the first equation of model (2), we geṫ It follows that Then, From the third equation of model (2) and (4), we have Further, let The dynamics of model (2) can be analyzed in the following bounded feasible region 2.2. The basic reproductive number. Based on the concept of the basic reproductive number for an epidemic disease presented in [6,35], we know the basic reproductive number R 0 of virus is the expected number of viruses that one virion gives rise to in an totally uninfected cell population during its lifetime.
From model (2), it is clear that healthy cells, CTL cells and antibody cells will stabilize to λ/d 1 , λ 1 /d 4 and λ 2 /d 5 if there is not infection, i.e., I(t) = V (t) = 0. In this case, suppose that one virion is introduced, it can produce a maximum amount of P 1 (ϕ 1 ) = β λ d1 e −sτ ϕ 1 infections during its mean lifetime of ϕ 1 = 1/(d 3 + q λ2 d5 ). In addition, an infected cell has an average lifetime of ϕ 2 = 1/(d 2 + p λ1 d4 ), and hence, it can averagely generate P 2 (ϕ 2 ) = rd 2 ϕ 2 virus. Therefore, the basic reproduction number of virus for model (2) can be defined as Based on the above expression, we know that there are inverse proportional relationship between the basic reproduction number of virus (R 0 ) and the recruitment rate of immune responses (λ 1 and λ 2 ). Thus, R 0 will decreases along with λ 1 , λ 2 increasing, which means that ignoring the effects of recruitment rate of immune responses will overestimate the basic reproduction number of virus.
3. Stability of the equilibria. In this section, we first discuss the existence of infection-free equilibrium, and then analyze its stability. Besides, using the uniform persistence theory, we obtain the existence of an endemic equilibrium. After that, the stability of an endemic equilibrium was proved by constructing Lyapunov functional.
3.1. Infection-free equilibrium. Apparently, there is always an infection-free equilibrium in system (2): Next, we discuss the stability of the infection-free equilibrium E 0 .
Theorem 3.1. When R 0 < 1, the infection-free equilibrium E 0 is globally asymptotically stable in region Γ.
Proof. First we define a Lyapunov functional L 0 by Calculating the time derivative of L 0 along the solution of system (2), we obtaiṅ Since the geometric mean is less than or equal to the arithmetical mean, it follows from R 0 < 1 thatL 0 ≤ 0. Let It is easy to show that E 0 = (T 0 , 0, 0, C 0 , A 0 ) is the largest invariant set in D 0 . By the Lyapunov-LaSalle invariance principle [8], E 0 is globally asymptotically stable.
3.2. Uniform persistence. In order to obtain the the existence of an endemic equilibrium, in this subsection, we investigate the uniform persistence of (2). We first introduce a preliminary theory. Let X be a metric space and Φ be a semiflow on X. Suppose that X 0 is an open set in X, X 0 ⊂ X, X 0 ∩ X 0 = ∅, and X 0 ∪ X 0 = X. Define M ∂ = {x ∈ X 0 : Φ t (x) ∈ X 0 , t ≥ 0}, which may be empty. A continuous p : X → [0, ∞) satisfying condition: p(Φ t (x)) > 0 for t > 0 if either p(x) = 0 and x ∈ X 0 or if p(x) > 0, will be called a generalized distance function for Φ.
Then there exists δ > 0 such that for any compact chain transitive set L with L ⊂ M i for i = 1, ..k, there holds min x∈L p(x) > δ.
Let ω(ψ) be the omega limit set of the orbit Φ(t) through ψ ∈ X + . Note that system (2) has an unique boundary equilibrium It then follows from the result in [14] that Then {E 0 } is isolated and acyclic covering and the conditions (i), (ii) and (iii) of Lemma 3.2 are satisfied. Since Thus, there is sufficiently small σ such that Since the non-diagonal elements of (8) are positive, and from (6), we obtain |A σ | < 0. By Perron-Frobrniuss Theorem, we can obtain the maximum eigenvalue α > 0 of A σ , and it has an eigenvector u = (u 1 , u 2 ), u 1 > 0, u 2 > 0 , then choose sufficiently small l such that I * (t 0 ) > lu 1 , V * (t 0 ) > lu 2 . Now consider the following auxiliary system Note I(t), V (t) are solutions of (9) with initial condition I(t 0 ) = lu 1 , V (t 0 ) = lu 2 .
From the Theorem 3.1, we are easy to get that E 0 is unstable if R 0 > 1, and by Proposition 2 and the uniformly persistent of system (2), we can obtain that if R 0 > 1, system (2) exists at least one endemic equilibrium E 1 = (T 1 , I 1 , V 1 , C 1 , A 1 ).
3.3. The endemic equilibrium. Now, we discuss the stability of the endemic equilibrium E 1 .
Theorem 3.4. When R 0 > 1, the endemic equilibrium E 1 is globally asymptotically stable in region Γ.
Proof. Set Define a Lyapunov functional L 1 by Calculating the time derivative of L 1 along the solution of system (2), we obtaiṅ Since we havė Since the geometric mean is less than or equal to the arithmetical mean and 1 − x + ln x ≤ 0 for any x > 0, it follows thatL 1 ≤ 0 . Let It is easy to verify thatL 1 (t) = 0 if and only if Thus, T (t) = T 1 andṪ (t) = λ − d 1 T 1 − βT 1 V (t) = 0. As a result, we have V (t) = V 1 , and then I(t) = I 1 . From the second equation and the third equation of model (2), we have which implies C(t) = C 1 , A(t) = A 1 . Therefore, the largest invariant set in D 1 is E 1 . Thus, when R 0 > 1, all positive solutions converge to E 1 by the LaSalle invariance principle [8].
4. Numerical simulations. In this section, we implement numerical simulations to explore the effects of the recruitment of immune responses (λ 1 and λ 2 ) on the infected cells (I 1 ) and virus load (V 1 ) at the endemic equilibrium E 1 . The parameter values are chosen from literatures ( [2], [26], [33], [32], [41], [38], [40], [46]). Especially, according to [9], we choose the range of λ 1 , λ 2 from 0 to 1. The all parameter values are shown in Table 1. Figure 1 illustrates that I 1 , V 1 decrease along with λ 1 , λ 2 increasing in which implies that ignoring the recruitment of immune responses will overestimate the severity of the infection.  Varied Recruitment rate of antibody [9] 5. Conclusions. In this paper, the global dynamics of a within-host model with immune responses and intracellular time delay has been studied. By the method of Lyapunov functional and persistence theory, we obtain the global stability of the model (2)  Considering the basic reproductive number of virus R 0 = R(τ ) = λβrd 2 e −sτ d 1 (d 2 + p λ1 d4 )(d 3 + q λ2 d5 ) as a function of τ , we can find that it is decreasing in τ and it tends to 0 if the time delay tends to ∞. Furthermore, comparing with the previous studies ( [11], [12], [16], [23], [24], [39], [42], [45]), we find that the expression of R 0 for model (2) is different, i.e., it includes the parameters λ 1 , λ 2 which reflect the recruitment of immune responses. This implies that ignoring the recruitment of immune responses will result in overestimation of the reproductive number. Numerical simulations also show that the part of V 1 , I 1 at steady state will decrease along with λ 1 , λ 2 increasing (see Fig.1), which mean that the recruitment of immune responses play a significant role in eradication of diseases. To sum up, we can conclude that ignoring the recruitment of immune responses will overestimate the infection degree and the severity of disease. These may provide a new insight for developing antiviral drug therapy strategies, which is to increase the recruitment rate of immune responses.