Bifurcation analysis of HIV-1 infection model with cell-to-cell transmission and immune response delay.

A within-host viral infection model with both virus-to-cell and cell-to-cell transmissions and time delay in immune response is investigated. Mathematical analysis shows that delay may destabilize the infected steady state and lead to Hopf bifurcation. Moreover, the direction of the Hopf bifurcation and the stability of the periodic solutions are investigated by normal form and center manifold theory. Numerical simulations are done to explore the rich dynamics, including stability switches, Hopf bifurcations, and chaotic oscillations.


1.
Introduction. Human Immunodeficiency Virus (HIV) and Acquired Immune Deficiency Syndrome (AIDS) have spread in successive waves in various regions and kept being a serious threat to public health. HIV targets cells with CD4 receptors, including the CD4 + T-cells, and damages the body's immune system, leading to humoral and cellular immune function loss (the marker of the onset of AIDS), making the body susceptible to opportunistic infections. The earlier models of virus infection describe the interaction between virus and target cells by assuming that the infected cells produce virions instantaneously [1,2].
The early models of virus infections, given by ordinary differential equations (ODEs), ignore the time delays of the viral infection, production of subsequent virus particles, and activation of immune response. Ciupe et al. [3] have shown that allowing for time delays in the model better predicts viral load data when compared to models without delays. The introduction of delays make the models more realistic. A discrete delay was first introduced into HIV infection model by Herz et al. [4]. Various models of viral dynamics with discrete or distributed delays have generally been studied [5][6][7][8][9][10][11][12][13][14][15][16].
We noticed that most within-host virus models concentrate on the virus-to-cell transmission. In fact, the infection via cell-to-cell contact is found to be much more rapid and efficient than virus-to-cell transmission because it avoids several 344 JINHU XU AND YICANG ZHOU biophysical and kinetic barriers [17]. It has been reported that cell-to-cell spread of virus is favored over infections with cell-free virus inocula [18,19]. The data of Gummuluru et al. [20] support the hypothesis that cell-to-cell spread of HIV-1 is the predominant route of viral spread since viral replication in a system with rapid cell turnover kinetics depends on cell-to-cell transfer of virus. Cell-to-cell transmission has also been reported for many other infections, such as HCV [21][22][23], Epstein Barr Virus (EBV) [24], Herpes simplex virus type-1 (HSV-1) [25], and HTLV-1 [26]. The mechanisms cell-to-cell transmission mode were, however, not well understood until the recent description of the "virological synapses" (VSs) [27]. Cell-to-cell spread greatly influences pathogenesis, not only facilitates rapid viral dissemination but may also promote immune invasion and, thereby, influence the disease [28][29][30]. As far as cell-to-cell infection is concerned, much less has been done in mathematical modeling. Culshaw et al. [11] studied a delayed two-dimensional model of cell-tocell spread of HIV-1 in tissue cultures with logistic growth term for target cells, assuming that infection is spread directly from infected cells to healthy cells and neglecting the effects of free virus. Thereafter, Lai and Zou [31] studied a virus model with both virus-to-cell infection and cell-to-cell infection. These authors also considered a model which including both cell-to-cell infection and full logistic growth term for target cells [32], Here T (t), I(t) and V (t) represent the concentrations of susceptible CD4 + T cells (target cells), productively infected T cells and free virus particles at time t, respectively. Target cells are infected by free viral particles and infectious cells (productively infected cells) at rates β 1 T (t)V (t) and β 2 T (t)I(t), respectively. r, T max , γ, d 1 and d 2 represent the growth rate of a target cell, carrying capacity of target cells, the rate of free viral particles released by infected cells, the losing rate of productively infected cells and free viruses, respectively. α (α ≥ 1) is the limitation coefficient of infected cells imposed on the growth of target cells. The stability, persistence as well as Hopf bifurcation of model (1) have been investigated.
The immune response has not been considered in model (1) though antibodies, cytokines, natural killer cells, and T cells are essential components of a normal immune response to a virus. In most virus infections, cytotoxic T lymphocytes (CTLs) play a critical role in antiviral defense by attacking virus-infected cells. Indeed, in HIV infection, CTLs are the main host factors which determine viral load. The dynamics of HIV infection with CTL response has received much attention in the past decades [3,5,12,13,16,33,34]. For example, Ciupe et al. [3] considered the following delayed HIV model Here E(t) is the concentration of effector cells. The constant r is the growth rate of target cells and the growth is limited by a carrying capacity T max . Target cells are infected by free viral particles at rates kT (t)V (t). d 1 , d 3 , N , d 2 , p and d 4 represent the death rate of productive infected cells, the killing rate of infected cells by effector cells, the number of virions produced by an infected cell during its life span (burst size), the viral clearance rate and productive rate of the effector cells and the death rate of effector cells, respectively. The term I(t − τ ) accounts for the time needed to activate the CD8 + T cell response, where τ is a constant. The authors mainly focused on estimating the kinetic parameters of model (2) while the dynamics behavior of model (2) has not been studied. The cell-to-cell transmission has not been taken into consideration in model (2).
Motivated by [3,32], we consider the following model with initial conditions The constant s is the source of CD4 + T-cells from precursors, d is the natural death rate (d < r in general). The other parameters in model (3) have the same meaning with model (1) and (2).
The paper is organized as follows. In Section 2, we present some preliminaries. In Section 3, the dynamics behavior of infection-free steady state of model (3) is studied. Both the local stability of the infection steady state for model (3) and the conditions for the existence of Hopf bifurcation are presented. Furthermore, the properties of the Hopf bifurcation solutions have been investigated by applying normal form and center manifold theory. In Section 4, numerical simulations are carried out to show the rich and complex dynamics of model (3), such as Hopf bifurcation, stability switches phenomena and chaotic oscillations. Finally, a brief summary and discussions complete the paper.
2. Preliminaries. We denote by X = C([−τ, 0], R 4 + ) the Banach space of continuous functions mapping the interval [−τ, 0] into R 4 + equipped with the sup-norm. By the standard theory of functional differential equations [35] we know that for of model (3) with initial condition (4).
with initial condition (4). Then T (t), I(t), V (t), E(t) are positive for all t ≥ 0, and they are ultimately bounded. Moreover, there exists an η 0 > 0 such that lim inf Proof. At first, we prove that T (t) is positive for t ≥ 0. Otherwise, there exists a positive t 0 , such that T (t) > 0 for t ∈ [0, t 0 ) and T (t 0 ) = 0. By the first equation of model (3), we have T (t 0 ) = s > 0. T (t 0 ) = s > 0 implies that T (t) < 0 for t ∈ (t 0 − , t 0 ) and sufficiently small > 0. This contradicts T (t) > 0 for t ∈ [0, t 0 ). It follows that T (t) > 0 for t > 0. From the equation of (3) we have From those expressions and (4) we know that the solution of model (3) is positive for all t ≥ 0. Next, we show that the solution of model (3) is ultimately bounded. From the first equation of (3), we obtain From this inequality and the comparison principle we know that lim sup where K = s + rT max 4 and δ = min{d, d 1 } . Thus, we have lim sup t→∞ G ≤ K δ and I(t) is ultimately bounded. It follows from the third and fourth equations of (3), Therefore, we have lim sup t→∞ V ≤ N d1K d2δ and lim sup t→∞ E ≤ pK d4δ . That is V (t) and E(t) are ultimately bounded. Furthermore, from the first equation of model (3) we have, whereĨ andṼ are the upper bounds of I(t) and V (t) respectively. This shows that T (t) is uniformly bounded away from zero.
Model (3) has two steady states: the infection-free steady state P 0 = (T 0 , 0, 0, 0), and the infected steady state P * = (T * , I * , V * , E * ), where it is easy to validate that R 0 is the basic reproduction number of system (3). Biologically, R 0 represents the average number of secondary infections. In fact, the basic reproduction number R 0 includes two parts, we can rewrite R 0 as The first term is the average number of secondary infection caused by a virus, corresponding to virusto-cell infection mode; the second term is the average number of secondary infection caused by an infected cell, corresponding to cell-to-cell infection. We can see that the basic reproduction number R 0 which we have defined is larger than that given in existing models with only one infection mode. The basic reproduction number of the model neglecting either the virus-to-cell infection or cell-to-cell infection may underevaluate the spread risk.
If R 0 < 1, then there is only the infection-free steady state. From the expression I * , we know that the infected steady state exists if and only if β 1 d2 + β 2 T * < d 1 , thus there exists no infection steady state, i.e., only the infection-free steady state exists.
3. Dynamics analysis of model.

3.1.
Stability of infection-free steady states P 0 . We linearize the model at steady states of model (3) to study the local stability. The characteristic equation is We have the following result for the infection-free steady state.
Theorem 3.1. The infection-free steady state P 0 of model (3) is locally asymptotically stable when R 0 < 1 and unstable when R 0 > 1.
Proof. At the infection-free steady state P 0 , the characteristic equation becomes There are two negative real roots: The inequality R 0 < 1 implies that d 2 + d 1 − β 2 T 0 > 0, and all the roots of (6) have negative real part. Then the infection-free steady state P 0 is locally asymptotically stable. If R 0 > 1, then (6) has at least one root with positive real part. Thus, the infection-free steady state P 0 is unstable.
Proof. For a continuous and bounded function f (t), we define We claim that From the fluctuation lemma [36], the second and third equations of model (3), we know that there is a sequence t n with t n → ∞ such that Those two inequalities lead to I ∞ is nonnegative since it is the supremum of the function I(t). If I ∞ > 0, then the inequality in (10) yields which is contradiction with R 0 < 1. The possible case is I ∞ = 0, which implies lim t→∞ I(t) = 0. From the inequality (9) and I ∞ = 0, we have V ∞ = 0, which implies 3.2. Stability of infected steady state P * and Hopf bifurcation. In this section, we investigate the stability of the infected steady state and the existence of Hopf bifurcations. The infected steady state P * (T * , I * , V * , E * ) satisfies The characteristic equation at the infected steady state P * is where a i > 0 (i = 0, 1, 2, 3), b i > 0 (i = 0, 1, 2), and When τ = 0, the corresponding characteristic equation becomes By Routh-Hurwitz criterion we know that all solutions of (12) have negative real parts if and only if The stability is given in the following theorem.
The root of (11) depends on τ continuously [38]. All roots of (11) locate in the left side of the imaginary axis if τ = 0 since the endemic equilibrium P * is stable. A root of (11) may pass through the imaginary axis and enter the right side when τ increases. λ = iω is the critical case since a root may enter the right side or the left side under small perturbation when it locates on the imaginary axis. After substituting λ = iω into (11) and separating the real and imaginary parts, we have The equations of (14) lead to where F (λ, τ ) = 0 has a purely imaginary root iω is equivalent to that G(z) = 0 has a positive real root z.
(i) All roots of (11) have negative real parts for τ ∈ [0, τ 0 ) if any one of the following conditions holds: (ii) All roots of (11) have negative real parts for τ ≥ 0 if the conditions in (i) are not satisfied.
According to the Hopf bifurcation theorem for functional differential equations [40, Theorem 1.1 in Chapter 11] and together with Lemmas 3.4, 3.5 and 3.6 we have following result. Theorem 3.7. Let τ 0 , z 0 be defined by (16). Suppose that (13) holds.
(i) If the conditions of (a)-(c) of Lemma 3.5 are not satisfied, then infected steady state P * is asymptotically stable for all τ ≥ 0. In this subsection, we study the direction and stability of the Hopf bifurcation by using the normal theory and the center manifold theorem [41]. We always assume that model (3) undergoes Hopf bifurcation at the steady state P * = (T * , I * , V * , E * ) for τ =τ = τ The detailed calculation of µ 2 , β 2 and T 2 is given in Appendix A. For the parameter values Data 1 given in Table.1. It is easy to see that R 0 = 0.7514 < 1, from Theorem 3.2 shows that the infection-free steady state P 0 is globally asymptotically stable for any τ ≥ 0 (see Fig. 1) Under the condition (13) the infected steady state P * is locally asymptotically stable independent of the size of the delay, though the time delay does cause transient oscillations in all components. For the parameter values given in the last column of Table 1, we can compute that R 0 = 78.05 > 1 and the infected steady state is P * = (37.71854, 16.3522, 1362.6799, 0.5053). The positive real roots of (15) are z 1 = 0.3563, z 2 = 0.1039, and the pure imaginary roots of (11) are λ 1 = iω 1 and λ 2 = iω 2 with ω 1 = 0.5969 > ω 2 = 0.3223. Furthermore, we have G (z 1 ) > 0 and G (z 2 ) < 0. From the transversal condition (17) and [15], we can have following results. (a) At τ (1) j , j = 0, 1, 2, · · · , a pair of characteristic roots of (11) crosses the imaginary axis from left to the right. (b) At τ (2) j , j = 0, 1, 2, · · · , a pair of characteristic roots of (11) crosses the imaginary axis from right to the left.
j−1 . From (c) we know that there exists an integer k such that τ The results of (a)-(c) imply the stability switch as τ increases: a pair of characteristic roots will cross the imaginary axis to the right at τ Furthermore, there exists a k = 1 such that τ 1 . The infected steady state P * is stable for τ < 2.5901, unstable for τ ∈ (2.5901, 10.2828), stable for τ ∈ (10.2828, 13.1167), and unstable for τ > 13.1167, which is presented in Fig. 2 and the corresponding stability and bifurcation is shown in Fig. 3 (left). The horizontal axis is the delay τ , and the vertical axis is the virus V . For τ < 2.5901 and τ ∈ (10.2828, 13.1167), there is a line in Fig. 3 (left), which is given V = V * = 1362.6799, the infected steady state P * is locally asymptotically stable. For τ ∈ (2.5901, 10.2828), when τ cross τ   Fig. 4 shows that when τ = 31, chaotic motions occurs. Furthermore, from Fig. 5, when τ becomes large, say τ = 49, the infected steady state P * is unstable, and the system trajectory exhibits a transient seemingly chaotic solution for a longer time (see the small figures in Fig. 5) then involves into a final nonchaotic state such as a quasi-periodic solution. Thus, immune response delay has an effect on the control of the disease.
From Fig. 6, which shows that though the components I and V of the infection steady state are only slightly changed, the time needed for the system converge to the steady state is much shorter as β 2 increases, i.e. the system will spend a shorter time to reach the steady state for high value of β 2 . Moreover, the figures in Fig.  7 together with the Fig. 3 (left one) shows that the stable intervals is enlarged as β 2 increases, though the amplitude of the periodic solutions is smaller. And when β 2 = 0.01, the Fig. 7 (right) shows that multiple stability switches can occurs. Hence, neglecting cell-to-cell transmission (β 2 ) may lose some dynamics behavior. As for parameter s, which is the source of new health target cells from precursors. Fig. 8, shows that the components I and V of the infection steady state increases as s increases. This is because high value of s increases the pool of susceptible target cells. Moreover, comparing with Fig. 9 and Fig. 3 (left) we can see that the amplitude of the periodic solutions increases as s increases, but the stable interval decreases. Furthermore, numerical simulation shows that there is a period-doubling solution (see Fig. 10). From Fig. 11, we can see that the component I and V of the infection steady state decreases as d increases. Comparing with Fig. 3 (left) and Fig. 12, we can see that the amplitude of the periodic solutions decreases as d increases, whereas the length of the stable interval increases. Hence, both the recruitment rate s and the death rate d of the target cells do have some impact on the dynamics of the model.
In general, the existence of logistic term may lead to rich dynamics for a model, especially for a model without delay, logistic term may cause Hopf bifurcation. In order to show the impact of the delay on the dynamics of the model without the effect of logistic term. Then, we give the bifurcation diagram when the system is in absence of logistic term, i.e., r = 0 in Fig. 3 (right), which implies that stability switch and Hopf bifurcation still exist when there is no logistic growth term for the system. Thus, we can see that stability switch, Hopf bifurcation and chaotic oscillation exist in both cases. We can claim that when take immune responses into consideration, time delay may be the main factor for periodic oscillations. Furthermore, the two figures in Fig. 3 show that the stable intervals for the system in absence of logistic growth term is much larger than the system with logistic growth term and the oscillation interval will be enlarged as r increases, though the existence of logistic growth term may not change the main dynamics behavior of the system.
Simulations are also done to show the impact of logistic growth term r on the dynamics of the model (see Fig. 13 and Fig. 14). The simulation shows that the infected steady state P * may be stable or unstable, and the model may has periodic solutions, or chaotic motions, depending on r. From Fig. 13, we see that only Hopf bifurcation occurs and no chaotic motions when τ , say τ = 2, stay in a stable interval; while both hopf bifurcation and chaotic motions exist when τ , say τ = 5, stay in an unstable interval, and the corresponding bifurcation diagrams are given in Fig. 14, respectively. Moreover, Fig. 14 (left) shows that at the left end of the r range, though there exists an interval for which the infected steady state P * is asymptotically stable, the viral load increases, then Hopf bifurcation occurs and the amplitude of the bifurcating periodic solutions increase and then decrease as r increase. Thus, the simulation shows that the logistic growth term also plays an important role on the dynamics of the model. Our results suggest that both immune delay and logistic growth term are responsible for rich dynamics of the model.   5. Summary and discussion. In this paper, we extend the previous work to a more realistic delayed model including cell-to-cell transmission. The basic reproduction number R 0 include two parts: cell-to-cell infection and virus-to-cell infection. It is easy to see that the basic reproduction number will be underestimated for models neglecting the cell-to-cell infection or virus-to-cell infection. Mathematical analysis gives the conditions for the existence of the equilibria and shows the influence of the time delay on the stability of equilibrium states. It is proved that the local stability of the uninfected steady state is independent of the size of the delay. Furthermore, the global stability of the uninfected steady state P 0 is obtained if R 0 is less than one by applying the fluctuation lemma. Our results show that      increasing the delay can destabilize the infected steady state by leading to a Hopf bifurcation and periodic solutions. The bifurcation direction and the stability of the periodic solutions are investigated by using normal form and center manifold. The theoretical analysis shows the importance of time delay on HIV dynamics. Numerical simulations show that both cell-to-cell transmission and time delay τ have an impact on the dynamics of the model, and rich dynamics can occur for large τ in the realistic parameter space. Compared to the earlier studies [3,32], our analysis shows that the introduction of the immune delay not only destabilize the stability of the infected steady state, leading to a Hopf bifurcation and periodic solutions, but also stability switch occurs as time delay τ increases, which has not been observed in [3,16,32]. Chaotic oscillations were observed for large τ . We also  find that the viral load may be destabilized into oscillations with the increase of the logistic growth rate r for T-cells. Moreover, numerical simulations also show that the oscillation interval will be enlarged as r increases. The system can occurs multiple stability switches for high value of cell-to-cell transmission β 2 , and the stable intervals is larger though the amplitude is smaller. Furthermore, numerical simulations show that larger β 2 makes the system converge to the steady state more easily. Hence, our results suggest that the logistic growth rate for T-cells r and the immune response delay τ and the cell-to-cell transmission are responsible for the complex dynamics.
Although the delay of immune response and cell-to-cell transmission considered in this paper is a good way to improve the viral dynamic model, the intracellular delays should also be taken into consideration for more realistic models. The dynamical analysis of the epidemic models with multiple delays will be more complex and bigger challenge in the future.