Threshold dynamics of an HIV-1 model with both viral and cellular infections, cell-mediated and humoral immune responses

: Human speciﬁc immunity consists of two branches: humoral immunity and cellular immunity. To protect us from pathogens, cell-mediated and humoral immune responses work together to provide the strongest degree of e ﬃ cacy. In this paper, we propose an HIV-1 model with cell-mediated and humoral immune responses, in which both virus-to-cell infection and cell-to-cell transmission are considered. Five reproduction ratios, namely, immunity-inactivated reproduction ratio, cell-mediated immunity-activated reproduction ratio, humoral immunity-activated reproduction ratio, cell-mediated immunity-competed reproduction ratio and humoral immunity-competed reproduction ratio, are calculated and veriﬁed to be sharp thresholds determining the local and global properties of the virus model. Numerical simulations are carried out to illustrate the corresponding theoretical results and reveal the e ﬀ ects of some key parameters on viral dynamics.


Introduction
The human immunodeficiency virus (HIV) is a lentivirus that causes HIV infection and over time acquired immunodeficiency syndrome (AIDS). AIDS leads to progressive failure of the immune system, which allows life-threatening opportunistic infections and cancers to thrive. In the past decades, within host virus models have been investigated in some literatures, which helps us understand the biological interactions between viruses and host cells. Nowak et al. [20] designed a mathematical model including uninfected cells x(t), infected cells y(t) and free virus v(t) to describe the viral dynamics in HIV-1 infection: v(t) = ky(t) − uv(t), (1.1) where uninfected cells x(t) are produced at rate s and die at rate d; β 1 is the infection rate of virus-to-cell infection; a is the death rate of infected cells; k denotes the number of free virus particles produced by per infected cell; u is the remove rate of virus. System (1.1) has been further investigated by Perelson and Nelson [21] and Cangelosi et al. [1].
Faced with different virus infections, immunity system protects us against pathogens. Human specific immunity can be classified into cell-mediated immunity, for which the protective function is associated with cells and humoral immunity, where the protective function exists in the humor [2]. As for cell-mediated immunity, activated effector T cells can detect peptide antigens originating from various types of pathogens and remove virus-infected cells. Some HIV-1 infection models have been proposed to describe the virus dynamics with cell-mediated immune response (see, for example, [15,19,24,26,34]). While, in humoral immunity, matured B cells migrate from bone marrow to lymph nodes or other lymphatic organs, where they begin to eliminate pathogens [23]. There have been several works on virus models with humoral immune response (see, for example, [4,14,[28][29][30]). In [6], Fouts et al. pointed out that a guiding principle for HIV vaccine design has been that cellular and humoral immunities work together to provide the strongest degree of efficacy. In [33], Yan

and Wang considered both cell-mediated and humoral immune responses and put forward an HIV-1 infection model including both T cells and B cells, which only involves virus-to-cell infection mechanism.
It is mentioned in [17] that cell-to-cell transmission is a more potent and efficient means of virus propagation than the virus-to-cell infection mechanism. Cell-to-cell spread not only facilitates rapid viral dissemination but also reduce the effectiveness of neutralizing antibodies and viral inhibitors by immune evasion. In [25], Sigal et al. proved that cell-to-cell spread of HIV-1 does reduce the efficacy of antiretroviral therapy, since cell-to-cell transmission can cause multiple infections of target cells, which can in turn reduce the sensitivity to the antiretroviral drugs. In view of this, some mathematical analysis of virus models with cell-to-cell transmission has been performed. For instance, Li and Wang [13] dealt with the global dynamics of an HIV infection model which incorporated direct cell-to-cell transmission. Meanwhile, Lai and Zou [11,12] studied the effect of cell-to-cell transfer of HIV-1 on the virus dynamics.
It was assumed in system (1.1) that the infection process is governed by the mass-action principle, namely, the infection rate per host or per virus is a constant. In [22], Regoes et al. illustrated that the infection rate is often found to be a sigmoidal rather than a linear function of the parasite dose to which it is exposed, and presented a dose-dependent infection rate (v/ID 50 ) κ /[1 + (v/ID 50 ) κ ], where ID 50 denotes the infectious dose at which 50% of the hosts are infected and κ measures the slope of the sigmoidal curve at ID 50 . In [10], Huang et al. indicated that the bilinear incidence rate is insufficient to describe the infection process in detail and proposed a class of nonlinear incidence. Besides, to place the model on more sound biological grounds, Xu [31] and Elaiw et al. [5] incorporated a saturation incidence β 1 v(t)/(1 + αv(t)) to replace the mass-action infection rate.
Motivated by the works of Fouts et al. [6], Yan and Wang [33], Sigal et al. [25] and Regoes et al. [22], in the present paper, we are concerned with the effects of cell-to-cell transmission, saturation incidence, both cell-mediated and humoral immune responses on the global dynamics of HIV-1 infection model. To this end, we consider the following delay differential equations: where x(t), y(t), v(t), z(t), w(t) denote the concentration of uninfected cells, infected cells, virus, T cells and B cells at time t, respectively, and other parameters are described in Table 1. A simple schematic diagram for the virus infection corresponding to system (1.2) is shown in Figure 1. The initial condition for system (1.2) takes the form It can be proved by the fundamental theory of functional differential equations [7] that system (1.2) has a unique solution (x(t), y(t), v(t), z(t), w(t)) satisfying the initial condition (1.3). It is easy to show that all solutions of system (1.2) with initial condition (1.3) are defined on [0, +∞) and remain positive for all t ≥ 0.  This paper is organized as follows. In Section 2, we calculate the reproduction ratios to system (1.2) and establish the existence of feasible equilibria. In Section 3, the local asymptotic stability of each of feasible equilibria is studied. In Section 4, we investigate the global asymptotic stability of each of feasible equilibria. In Section 5, we present numerical simulations to illustrate the theoretical results and study the effects of cell-to-cell transmission, viral production rate, death rate of infected cells and viral removal rate on viral dynamics, respectively. Besides, we perform a sensitivity analysis of reproduction ratios. The paper ends with a conclusion in Section 6.

Reproduction ratios and feasible equilibria
Clearly, system (1.2) always has an infection-free equilibrium E 0 (s/d, 0, 0, 0, 0). Denote where R 0 is called immunity-inactivated reproduction ratio of system (1.2), which represents the number of newly infected cells produced by one infected cell during its lifespan [3]. It is easy to show that if R 0 > 1, system (1.2) has an immunity-inactivated equilibrium E 1 (x 1 , y 1 , v 1 , 0, 0), where and in which, R 1 is called cell-mediated immunity-activated reproduction ratio, which denotes the average number of T cells activated by infectious cells when virus infection is successful and humoral immune response has not been established. If R 1 > 1, in addition to E 0 and E 1 , system (1.2) has a cell-mediated immunity-activated equilibrium E 2 (x 2 , y 2 , v 2 , z 2 , 0), where We further denote in which R 2 is called humoral immunity-activated reproduction ratio, which denotes the average number of B cells activated by viruses when virus infection is successful and cell-mediated response has not been established. When R 2 > 1, system (1.2) has a humoral immunity-activated equilibrium where w 3 is the positive real root of the following quadratic equation: Denote where R 3 is called humoral immunity-competed reproduction ratio and represents the average number of B cells activated by viruses under the condition that cell-mediated immune response has been established, while, R 4 is called cell-mediated immunity-competed reproduction ratio and represents the average number of T cells activated by infectious cells under the condition that humoral immune response has been established. If R 3 > 1 and R 4 > 1, system (1.2) has an immunity-activated equilibrium E * (x * , y * , v * , z * , w * ), where in which cell-mediated and humoral immune responses take effect simultaneously.

Local asymptotic stability
In this section, we are concerned with the local asymptotic stability of each of feasible equilibria to system (1.2) by analyzing the distribution of roots of corresponding characteristic equations .
1) It is clear that (3.1) has negative real roots λ = −b 1 , λ = −b 2 , λ = −d and other roots are determined by the following equation: Substituting R 0 and R 02 into (3.2) yields Now, we claim that all roots of (3.3) have negative real parts. Otherwise, there exists a root λ 1 = Reλ 1 + iImλ 1 with Reλ 1 ≥ 0. In this case, if R 0 < 1, it is easy to see that It follows that which contradicts to (3.3). Therefore, if R 0 < 1, all roots of (3.1) have negative real parts and E 0 is locally asymptotically stable. If R 0 > 1, we denote the left side of (3.2) by G(λ): where G(0) = au(1 − R 0 ) < 0 and G(λ) → ∞ as λ → ∞. Noting that G(λ) is a continuous function in respect to λ, if R 0 > 1, Eq. (3.1) has a positive real root, then E 0 is unstable. (3.5) Note that in which where It is clear that (3.5) has negative real roots λ = c 1 y 1 − b 1 and λ = c 2 v 1 − b 2 , and other roots are determined by the following equation: For the sake of contradiction, let λ 2 = Reλ 2 + iImλ 2 with Reλ 2 ≥ 0. In this case, it is easy to see that Direct calculation shows that which contradicts to (3.8). Thus, if R 0 > 1, R 1 < 1 and R 2 < 1, all roots of Eq. (3.5) have negative real parts, and E 1 is locally asymptotically stable.
and other roots are determined by the following equation: Similarly, we claim that all roots of (3.10) have negative real parts. Otherwise, there exists a root λ 3 = Reλ 3 + iImλ 3 with Reλ 3 ≥ 0. In this case, it is obvious that It follows that which contradicts to (3.10). Hence, if R 1 > 1 and R 3 < 1, all roots of Eq. (3.9) have negative real parts, and E 2 is locally asymptotically stable. (3.11) Note that It is obvious that (3.11) has negative real root λ = c 1 y 3 − b 1 , and other roots are determined by the following equation: (3.13) Similarly, we claim that all roots of (3.13) have negative real parts. If not, there exists a root λ 4 = Reλ 4 + iImλ 4 with Reλ 4 ≥ 0. In this case, it is easy to see that Direct calculation yields which contradicts to (3.13). Therefore, if R 2 > 1 and R 4 < 1, all roots of Eq. (3.11) have negative real parts, and E 3 is locally asymptotically stable. (3.14) Similarly, we claim that all roots of (3.14) have negative real parts. Otherwise, there exists a root λ 5 = Reλ 5 + iImλ 5 with Reλ 5 ≥ 0. In this case, it is clear that Direct calculation shows that which contradicts to (3.14). Therefore, if R 3 > 1 and R 4 > 1, all roots of Eq. (3.14) have negative real parts, and E * is locally asymptotically stable.

Global asymptotic stability
In this section, we study the global stability of each of feasible equilibria to system (1.2) by suitable Lyapunov functionals and LaSalle's invariance principle. First, we discuss the boundedness of solutions.
Proof. Let (x(t), y(t), v(t), z(t), w(t)) be any solution of system (1.2) with initial condition (1.3). Denote Calculating the derivatives of B 1 (t) and B 2 (t) in respect to t yieldṡ Therefore, the following set is positively invariant set for system (1.2): It is easy to see that x(t), y(t), v(t), z(t) and w(t) are bounded in the invariant set Ω.
Next, define a function g(x) = x − 1 − lnx, which will be used in Lyapunov functionals of this section.

Direct calculation yieldṡ
It follows thatV 1 (t) ≤ 0 with equality holding if and only if Noting that if R 0 < 1, E 0 is locally asymptotically stable, thus we obtain the global asymptotic stability of E 0 from LaSalle's invariance principle.

(t) ≤ 0 with equality holding if and only if
Noting that if R 1 > 1 and R 3 < 1, E 2 is locally asymptotically stable, thus we obtain the global asymptotic stability of E 2 from LaSalle's invariance principle.

Numerical simulations
In this section, we want to illustrate the theoretical results for system (1.2) by numerical simulations. Besides, we investigate the effects of cell-to-cell transmission, viral production rate, death rate of infected cells and viral remove rate on viral dynamics. Furthermore, sensitivity analysis is used to quantify the range of variables in reproduction ratios and identify the key factors giving rise to reproduction ratios, which can be helpful to design treatment strategies and provide insights on evaluating effective antiviral drug therapies. Following [18,26,27,32], we choose appropriate parameters and simulate each of feasible equilibria, respectively.
Case 1: Corresponding parameters are listed in Case 1 of Table 2. The immunity-inactivated reproduction ratio is calculated as R 0 = 0.5640 < 1. From Theorem 3.1, we derive that infection-free equilibrium E 0 is locally asymptotically stable, which is illustrated in Figure 2.
Case 2: Corresponding parameters are listed in Case 2 of Table 2. By simple computing, we obtain that R 0 = 1.0217 > 1, R 1 = 0.9635 < 1 and R 2 = 0.9625 < 1. From Theorem 3.2, we derive that immunity-inactivated equilibrium E 1 is locally asymptotically stable, which is in accord with Figure 3.
Case 3: Corresponding parameters are listed in Case 3 of Table 2. Similarly, we obtain that R 1 = 5.1140 > 1 and R 3 = 0.2459 < 1. From Theorem 3.3, we derive that cell-mediated immunity-activated equilibrium E 2 is locally asymptotically stable, which is illustrated in Figure 4.
Case 4: Corresponding parameters are listed in Case 4 of Table 2. Likewise, we obtain that R 2 = 14.1830 > 1 and R 4 = 0.2108 < 1. From Theorem 3.4, we derive that humoral immunity-activated equilibrium E 3 is locally asymptotically stable, which is in keeping with in Figure 5.
Case 5: Corresponding parameters are listed in Case 5 of Table 2. By calculation, we obtain that R 3 = 24.5895 > 1 and R 4 = 3.7395 > 1. From Theorem 3.5, we derive that immunity-activated equilibrium E * is locally asymptotically stable, which is consistent with observation in Figure 6.  Figure 2. The temporal solutions of x(t), y(t), v(t), z(t) and w(t) versus t of system (1.2) where R 0 = 0.5640 < 1.

Effect of cell-to-cell transmission
In order to investigate the effect of cell-to-cell transmission, we carry out some numerical simulations to show the contribution of cell-to-cell transmission during the whole infection. First, we let β 2 as zero to compare the virus infection without cell-to-cell transmission with the infection which has both transmissions. Figure 7 (β 2 = 0, β 2 = 4.7 × 10 −7 ) shows that cell-to-cell transmission is of benefit to HIV-1 transmission and the time to reach the peak level of virus is shorter. Then, we increase β 2 to study the change of the peak level of infected cells and virus, and the time to reach the peak level. Figure 7 (β 2 = 4.7 × 10 −7 , β 2 = 4.7 × 10 −6 , β 2 = 4.7 × 10 −5 ) shows that infected cells and virus reach the peak level more quickly as β 2 increases, meanwhile, the peak level become larger as β 2 increases, too. Therefore, cell-to-cell transmission plays an important role in the whole virus infection.

Effect of viral production rate
Viral production rate also has great influence on the dynamical behavior of the model. We set the viral production rate k as 11.349, 34.047, 68.094. In Figure 8, we observe that the time to reach the peak level of infected cells and virus becomes shorter as k increases, which means that larger viral production rate contributes to the viral infection. Meanwhile, T cells and B cells increase more quickly as k increases, especially, larger viral production rate can stimulate more B cells. Hence, the peak level of infected cells and virus decreases as k increases. In terms of the prevention and treatment of HIV, it implies that antiretroviral therapies, such as, reverse transcriptase inhibitors and protease inhibitors are effective methods to decrease k, namely, to inhibit virus reproduction.

Effect of death rate of infected cells and viral remove rate
Usually, the death rate of infected cells is larger than the death rate of uninfected cells due to the fact that HIV infection can kill more host cells. We present some numerical simulations to study the effect of death rate of infected cells on the dynamical behavior of the model. We can observe from Figure 9 that, infected cells and virus increase more slowly as a increases, which indicates that increasing the death rate of infected cells can slow down the virus infection. Humoral immunity is mainly used to clear virus in our humor, so the viral remove rate has an effect on viral infection as well. Figure 10 implies that as the viral remove rate increases, infected cells and virus increase more slowly, which has the similar results to a. In the clinic treatment of HIV, promoting body's immune capacity contributes to increasing the death rate of infected cells and viral remove rate.

Sensitivity analysis
Sensitivity analysis is used to quantify the range of variables in reproduction ratios and to identify the key factors giving rise to reproduction ratios. In [9,16], Latin hypercube sampling (LHS) is found to be a more efficient statistical sampling technique which has been introduced to the field of disease modelling.
LHS allows an un-biased estimate of the reproduction ratios, with the advantage that it requires fewer samples than simple random sampling to achieve the same accuracy. For each parameter of     reproduction ratios, a probability density function is defined based on experimental data and stratified into N equiprobable serial intervals. A single value is then selected randomly from every interval and this is done for every parameter. In this way, an parameter value from each sampling interval is used only once in the analysis but the entire parameter space is equitably sampled in an efficient manner. Distributions of the reproduction ratios can then be derived directly by running the model N times with each of the sampled parameter sets.
In terms of the prevention and treatment of HIV, we pay more attention to antiretroviral therapies, which is directly related to viral production rate and viral remove rate. Figure 11 shows the scatter plots of R 0 , R 1 and R 2 in respect to k and u, which implies that k is a positively correlative variable with R 0 and R 2 , while u is a negatively correlative variable. As for R 1 , we find that the correlation between k and R 1 or u and R 1 is not clear.
In [16], Marino et al. mentioned that Partial Rank Correlation Coefficients (PRCCs) provide a measure of the strength of a linear association between the parameters and the reproduction ratios. Furthermore, PRCCs are useful for identifying the most important parameters. The positive or negative of PRCCs respectively denote the positive or negative correlation with the reproduction ratios, and the sizes of PRCCs measure the strength of the correlation. First, we investigate the immunity-inactivated reproduction ratio R 0 , as we can see in Figure 12, β 1 and k are positively correlative variables with R 0 while others are negatively correlative variables. In order of correlative strength, it goes: β 1 , d, a, k, u and β 2 . Similarly, we obtain the PRCCs of R 1 and R 2 (see Figure 13). Specially, we observe that k and u is weakly correlative in respect to R 1 , which accords with the scatter plots of R 1 .

Conclusion
In this paper, we have considered an HIV-1 infection model to describe cell-to-cell transmission, saturation incidence, both cell-mediated and humoral immune responses. By a complete mathematical analysis, the threshold dynamics of the model is established and it can be fully determined by reproduction ratios. If R 0 < 1, the infection-free equilibrium E 0 is locally and globally asymptotically stable; if R 0 > 1, R 1 < 1 and R 2 < 1, the immunity-inactivated equilibrium E 1 is locally and globally  Figure 11. Scatter plots of R 0 , R 1 and R 2 in respect to k and u.  Figure 13. Tornado plots of PRCCs in regard to R 1 and R 2 . asymptotically stable; if R 1 > 1 and R 3 < 1, the cell-mediated immunity-activated equilibrium E 2 is locally and globally asymptotically stable; if R 2 > 1 and R 4 < 1, the humoral immunity-activated equilibrium E 3 is locally and globally asymptotically stable; if R 3 > 1 and R 4 > 1, the immunity-activated equilibrium E * is locally and globally asymptotically stable.
Numerical simulations vividly illustrate our main results of stability analysis for system (1.2). Besides, we have investigated the effects of cell-to-cell transmission, viral production rate, death rate of infected cells and viral remove rate on viral dynamics. It is worth mentioning that as the infection rate of cell-to-cell transmission β 2 increases, virus load rises quickly and largely, which implies that cell-to-cell transmission facilitates virus spread. Furthermore, we perform a sensitivity analysis of reproduction ratios, which implies some useful consequences on the prevention and treatment of HIV-1.
It is easy to see that immunity-inactivated reproduction ratios R 0 is the sum of the reproduction ratio determined by virus-to-cell infection, R 01 , and that determined by cell-to-cell transmission, R 02 . In other words, immunity-inactivated reproduction ratio R 0 becomes larger when the model includes cell-to-cell transmission. Meanwhile, we find that our research includes some existing work. When β 2 = 0 and α = 0, our virus model is similar to the model in [33] and the immunity-inactivated reproduction ratio R 0 reduces to R 01 . Based on the model in [33], Wang et al. [27] consider nonlinear incidence and continuous intracellular delay, which is similar to our model with β 2 = 0 only. Besides, when we only consider one of the immune responses, our model reduces to the models in [14] and [26].