Global Stability and Bifurcation Analysis of a Virus Infection Model with Nonlinear Incidence and Multiple Delays

: In order to investigate the impact of general nonlinear incidence, cellular infection, and multiple time delays on the dynamical behaviors of a virus infection model, a within-host model describing the virus infection is formulated and studied by taking these factors into account in a single model. Qualitative analysis of the global properties of the equilibria is carried out by utilizing the methods of Lyapunov functionals. The existence and properties of local and global Hopf bifurcations are discussed by regarding immune delay as the bifurcation parameter via the normal form, center manifold theory, and global Hopf bifurcation theorem. This work reveals that the immune delay is mainly responsible for the existence of the Hopf bifurcation and rich dynamics rather than the intracellular delays, and the general nonlinear incidence does not change the global stability of the equilibria. Moreover, ignoring the cell-to-cell infection may underevaluate the infection risk. Numerical simulations are carried out for three kinds of incidence function forms to show the rich dynamics of the model. The bifurcation diagrams and the identiﬁcation of the stability region show that increasing the immune delay can destabilize the immunity-activated equilibrium and induce a Hopf bifurcation, stability switches, and oscillation solutions. The obtained results are a generalization of some existing models.


Introduction
Human immunodeficiency virus (HIV) is regarded as one of the major threats to the health of human society and is an important research topic in the field of public health. In recent years, an increasing number of scholars have investigated the dynamic behaviors of HIV infection models by incorporating the different factors that impact the infection procedure. The analysis and understanding of the dynamical behaviors of HIV in the host by modeling HIV infection plays an important role in exploring the mechanism of virus infection. The classical virus infection model is composed of three components: uninfected cells, infected T-cells, and the virus [1]. As the research progresses, improved models have been proposed by many researchers (see [2][3][4][5][6][7][8][9][10][11][12][13][14][15][16][17] and the references therein).
During the process of HIV infection, it takes some time for the initial virus to enter the target cells and for the subsequent viral latency, as well as for infected CD4+ T cells to release infectious free virus particles. Moreover, because of the presence of latently infected cells, HIV cannot be completely eradicated and can be reactivated, continuing to replicate even after antiviral therapy. Thus, this may be an important reason for explaining the failure to eradicate HIV virus infection. In addition, time is needed during the activation procedure of latent cells; thus, there exists a time delay for latent cells to be activated and converted to infected cells [15,18]. Motivated by the above facts, dynamical models with intracellular delays and latency have been investigated (see [11,[14][15][16][19][20][21] and the references therein).
As we know, the main two primary immunity modes are humoral and cellular immunity, which are dominated by B cells and cytotoxic T lymphocytes (CTLs), respectively. In virus infections, cellular immunity is reduced by CTLs attacking the infected cells, whereas the B-cell immune response is prevalent during viral infections by attacking the viruses. Both modes are regarded as an important path of eliminating or controlling viral gain when HIV invades the human body. The investigation of the two modes of immune response using viral infection models has received attention from many researchers (see [2,4,5,10,19,[22][23][24] and the references therein). However, it is not clear which immune response mode is the most effective one. The existing literature on malaria infection shows that the humoral immune response is more effective compared to the cellular immune response [25]. Thus, the humoral immunity response is considered in the current study. Admittedly, an effective immune response requires the combination of humoral and cellular immunity due to the fact that the humoral immune response alone may not eradicate the infection [26].
In addition, viral stimulation of antigens also requires some time to generate an effective humoral immune response. Studies have been conducted that analyzed the effect of humoral immune delay on the equilibrium stability of viral infection models (see [6,16,27,28] and the references therein). The results have shown that the dynamics of the models incorporating immune delays become more complex, and the existence of time delays can destabilize the steady state and result in Hopf bifurcation, or even chaos solutions, in the corresponding models. However, what the dynamics will be when taking both the intracellular delay and immune response delay into account in a single model remains unknown. The answer will be addressed in this study.
Recently, the literature has implied that the HIV infection mode within the host has another important mechanism, i.e., cell-to-cell infection, in addition to virus-to-cell infection. It has been found that cell-to-cell infection is a more potent and efficient way of virus propagation than the virus-to-cell infection mode [34][35][36][37][38]. Motivated by this fact, models incorporating cell-to-cell infection have been proposed and studied by many researchers (see [5,20,[39][40][41][42][43][44][45][46][47] and the references therein). Note that only virus-to-cell infection was considered and the B-cell immune response in the host was ignored in Model (1). However, the B-cell immune response plays a key role in the immune response by detecting and eliminating HIV virus particles during infection. Thus, motivated by [6,14,16], we consider a delayed virus infection model incorporating general nonlinear incidence and humoral immunity in line with Model (1). In addition, two infection mechanisms are taken into consideration in the model. Therefore, the model takes the form: where B(t) denotes the concentrations of the B cells at time t. Infected cells in the host stimulate B-cell production at a rate of qVB, and free virus particles are cleared by antibodies at a rate of pVB. d 5 is the mortality rate of the B cells, β 1 is the rate of virus-to-cell infection, β 2 is the rate of cell-to-cell infection, τ 3 is the time delay for latent infected cells to become active infected cells, and τ 5 is the time delay of the activation of the B-cell immune system. The other parameters have the same meanings as in Model (1) (see [14]). Here, the incidences are assumed to be the nonlinear forms β 1 T f (V) and β 2 Tg(I), respectively, and f (V) and g(I) satisfy the following properties In line with (3), one can obtain Assume that Model (2) satisfies the following initial condition where (ψ 1 (θ), ψ 2 (θ), ψ 3 (θ), ψ 4 (θ), ψ 5 (θ)) ∈ C([−τ, 0], R 5 + ). The rest of the paper is organized as follows. In Section 2, we present some basic results, including the existence of equilibria and the positivity and boundedness of the solutions for Model (5). In Section 3, both the global stability of all equilibria of Model (5) and the existence of local and global Hopf bifurcations are investigated. Moreover, the properties of the Hopf bifurcation solutions are analyzed by applying the normal form and center manifold theory. Some numerical simulations are carried out in Section 4. The summary in Section 5 concludes the paper.

Preliminary Results
Before analyzing the dynamical behaviors of the model, we first present some preliminary results. According to the results presented in [48], the solution of Model (2) with initial conditions (5) is non-negative. In the following, we show the boundedness of the solution. From the first equation in Model (2), a simple calculation yields 2kq B(t + τ 4 + τ 5 ). The derivation of P(t) yields Thus, lim sup t→+∞ P(t) ≤ λ d , which implies that T(t), L(t), I(t), V(t), and B(t) are bounded. Based on the above analysis, we have the following result. Theorem 1. The solutions (T(t), L(t), I(t), V(t), B(t)) of Model (2) with initial conditions (5) are non-negative and ultimately bounded for all t ≥ 0.
Based on the above discussion, it is reasonable to suppose that there exists a positive constant M > 0 such that T ≤ T 0 , L, I, V, B ≤ M for large t. Thus, in the following, we analyze the dynamic behaviors of Model (2) in a bounded feasible region, which is given by In the following, we show the existence of equilibria for Model (2), including the infection-free equilibrium E 0 = (T 0 , 0, 0, 0, 0) (i.e., no infection exists), immunity-inactivated equilibrium E 1 = (T 1 , L 1 , I 1 , V 1 , 0) (i.e., the immune response has not been activated), and immunity-activated equilibrium E 2 = (T 2 , L 2 , I 2 , V 2 , B 2 ) (i.e., immune response has been activated and coexists with the virus). For convenience, we denote D 1 = αη This is the only biologically meaningful equilibrium if which is the basic reproduction number. It follows from the expression of the basic reproduction number R 0 that neglecting either the virus-to-cell infection or the cell-to-cell infection may underevaluate the infection risk. Any other equilibrium (2) is determined using the following equations.
When B 1 = 0 and V 1 = 0, it follows from (6) that In order to have T 1 > 0 and I 1 > 0 at equilibrium, we must have I 1 ∈ (0, λD 1 d 3 ). From the first equation in (7), we have T 1 = , and then substituting T 1 into the first equation in (6) yields =: H(I).
For all I > 0, it follows from (4) that Further, from (3), we have This implies that the immunity-inactivated equilibrium E 1 = (T 1 , L 1 , When B 2 = 0, a simple calculation from (6) implies that where ). In order to have T 2 > 0 and I 2 > 0, we must have . From the first equation in (8), we have then, substituting T 2 into the first equation in (6) leads to For all I > 0, from (4), we can obtain < 0. Then, there exists a unique I 2 ∈ (0, λD 1 d 3 ) such that F(I 2 ) = 0. This implies that the immunity-activated equilibrium E 2 = (T 2 , L 2 , I 2 , V 2 , B 2 ) exists only if R 1 > 1.
Before the proof of the global stability of the immunity-inactivated equilibrium E 1 , we present the following result.
kq > 0. and the inequalities R 1 = qV 1 d 5 = kqI 1 d 4 d 5 > 1 are equivalent to the inequality I 1 > I. From the monotonicity of function H(I), it follows that R 1 > 1 is equivalent to H( I) > 0 for H(I 1 ) = 0.
Proof. We employ a particular function, From the equilibrium conditions of E 1 , we have which leads tȯ It follows from (3) that Moreover, from Lemma 1, we haveV 2 ≤ 0, with equality holding if and only if T = T 1 , L = L 1 , I = I 1 , V = V 1 , B = B 1 , which implies that the largest compact invariant set in X ∈ Γ|V 2 = 0 is the singleton {E 1 }. Hence, the globally asymptotic stability of the unique immunity-inactivated equilibrium E 1 follows from LaSalle's invariance principle [49]. This completes the proof.
By calculating the derivative of V 3 along the solutions of Model (2) and using the following equilibrium condition at E 2 It follows from (3) that Therefore, the globally asymptotic stability of the unique immunity-activated equilibrium E 2 follows from LaSalle's invariance principle [49]. This completes the proof.

Bifurcation Analysis at E
Then, the characteristic equation in Model (11) at E 2 is given by Calculating the corresponding determinant gives where ).

Direction of Hopf Bifurcations
Based on the above discussion, a Hopf bifurcation occurs when τ 5 = τ 0 . In this section, we study the direction of the Hopf bifurcation using the normal theory and center manifold theorem [50]. We always assume that Model (2) undergoes a Hopf bifurcation at 5k , and iω * is the corresponding purely imaginary root of the characteristic equation at E 2 for τ 5 = τ * . The conditions for the direction and stability of the Hopf bifurcation are summarized in the following theorem. Theorem 6. (i) The direction of the Hopf bifurcation is determined by the sign of µ 2 , i.e., it is a supercritical bifurcation when µ 2 > 0 and a subcritical bifurcation when µ 2 < 0. (ii) The stability of the bifurcated periodic solution is determined by β 2 , i.e., the periodic solution is stable when β 2 < 0 and unstable when β 2 > 0. (iii) The period of bifurcated periodic solutions is determined by T 2 , i.e., the period increases when T 2 > 0 and decreases when T 2 < 0.
The detailed calculations of µ 2 ,β 2 , andT 2 are given in Appendix A.

Global Hopf Bifurcation
In the following, we explore the global Hopf bifurcation for R 1 > 1 using the global Hopf bifurcation theorem of functional equations [51]. Throughout this section, we assume that (16) has a unique positive root. The following lemma is obvious from the definition of the positive invariant set in Model (2).
It follows from Theorem 4 that E 2 is globally asymptotically stable when τ 5 = 0. As a special case of Model (2), Model (11) has no nonconstant periodic solution for τ 5 = 0. Thus, the following lemma holds.  The following conclusion shows that Model (11) cannot have a τ 5 -periodic solution. Proof. Let U(t) = (T(t), L(t), I(t), V(t), B(t)) be a nontrivial τ 5 -periodic solution of Model (11). Then, U(t + τ 5 ) = (T(t + τ 5 ), L(t + τ 5 ), I(t + τ 5 ), V(t + τ 5 ), B(t + τ 5 )) = U(t) is also a τ 5periodic solution of the following corresponding model It follows from Theorem 4 that Model (19) has no periodic solution when R 1 > 1. Thus, this contradiction completes the proof of the lemma. (11) is equivalent to the following functional differential equations where The mapping F : X × R + × R + → R 5 + is completely continuous. By restricting F to the subspace of the constant functions in X, we obtain the mapping F = F Then, it can be shown that the assumptions (A1)-(A3) in [51] hold. A point ( W, τ, T) is called a stationary solution of (20) T ) = 0 for some integer m. A center ( W, τ, T ) is said to be isolated if it is the only center in some neighborhood of ( W, τ, T ), and it has finite characteristic values of the form im 2π T . Let is the mth crossing number of (E 2 , τ (j) 5k , 2π ω k ) and m is an integer.
Proof. It follows from [51] that E 2 , τ (j) 5k , 2π ω k is an isolated center, and there is a smooth , which verifies the assumption (A4) in [51] for m = 1. Moreover, if we put = 0 that the crossing number and then, we have 5k . For a contradiction, we assume that the projection of E 2 , τ (j) 5k , 2π ω k onto the τ-space is bounded, which implies that there exists a τ * > 0 such that the projection of E 2 , τ (j) 5k , 2π ω k onto the τ-space is included in the interval (0, τ * ). Noting that 2π ω k < τ (j) 5k and from Lemma 3, we have T < τ * for (W, τ, T ) ∈ E 2 , τ ω k is bounded. This yields a contradiction and completes the proof.

Numerical Simulations
In this section, we demonstrate the above theoretical results through numerical simulations by choosing three different forms of the general incidence function in Model (2). Most of the parameter values come from [6,17,22,31].
As shown in Figure 6 (left), when τ i = 0 (i = 1, 2, 3, 4) and by choosing τ 5 as the bifurcation parameter, the bifurcation diagram with respect to τ 5 shows that the immunityactivated equilibrium is stable for a smaller immune delay (τ 5 < τ 0 ) and unstable for a larger immune delay, leading to irregular oscillations. However, the theoretical results obtained using Theorems 1-4 show that the intracellular delays τ 1 , τ 2 , τ 3 , and τ 4 cannot change the global stability of the equilibria. Nevertheless, how the intracellular delays impact the dynamical behaviors of Model (2) also needs to be considered by regarding τ 5 as the bifurcation parameter and varying τ i (i = 1, 2, 3, 4).By comparing Figure 6 (left) and Figure 7 (left column), it can be seen that the presence of τ 1 has little impact on the dynamics of Model (2) when τ 1 increases. However, in Model (2), the phenomenon of stability switches occurs when τ 2 appears, except for the Hopf bifurcation and irregular oscillations, as shown in Figure 7 (right column). This indicates that the stability and instability of immunity-activated equilibrium alternate a finite number of times and finally become unstable. A similar phenomenon occurs in the presence of τ 3 (left column) and τ 4 (right column), as shown in Figure 8. The numerical results show that the immune delay (τ 5 ) combined with the intracellular delays τ 2 and τ 4 can lead to more rich and complex dynamical behaviors of Model (2) compared to τ 1 and τ 3 . In order to more clearly compare the effect of intracellular delays τ i (i = 1, 2, 3, 4) combined with the immune delay τ 5 , we carried out a stability analysis by varying the values of the intracellular delays 0.25 ≤ τ i (i = 1, 2, 3, 4) ≤ 15 and the immune delay 0.25 ≤ τ 5 ≤ 50. The stability region was obtained numerically and the corresponding results are shown in Figure 9. We can see that the combination of τ 2 and τ 4 with τ 5 leads to much richer dynamics. The stability or stability switches are important characteristics from a biological perspective, as they determine whether the eradication of infection is possible in the case of stability regions with the implementation of effective controlling measures. In contrast, it will be difficult to eradicate the infection in unstable regions. Furthermore, in order to investigate the impact of the coexistence of all intracellular delays τ i > 0(i = 1, 2, 3, 4) and the immune delay τ 5 on the dynamics of Model (2), we fixed τ 1 = 2, τ 2 = 3, τ 3 = 6, and τ 4 = 9. The bifurcation diagram shows that stability switches still exist, except for the Hopf bifurcation and complex oscillations in Model (2), as shown in Figure 6 (right).    1+0.1V(t) and g(I(t)) = 1+0.4I(t) . Clearly, conditions (3) and (4) are true. When τ 5 > 0, τ i = 0 (i = 1, 2, 3, 4) using the same parameter values as those in Figure 3 and choosing τ 5 as the bifurcation parameter. The numerical results indicate that only one positive root of Q(z) = 0 exists, i.e., z = 0.1037, which implies that ω = 0.3220 and the critical value τ 0 = 2.9134. Figure 10 shows that the immunity-activated equilibrium E 2 is locally asymptotically stable when τ 5 = 1 < τ 0 = 2.9134, and a periodic solution bifurcates from an unstable E 2 when τ 5 = 3 > τ 0 = 2.9134, as shown in Figure 11. The corresponding bifurcation diagram with respect to τ 5 in the absence of intracellular delays τ i = 0, (i = 1, 2, 3, 4) shows that a Hopf bifurcation and periodic oscillations occur in Model (2), as shown in Figure 12 (left). Moreover, the bifurcation diagram with respect to τ 5 with intracellular delays τ i > 0 (i = 1, 2, 3, 4) shows that a Hopf bifurcation and periodic solutions also exist, and we can see that the stable interval is longer for τ 1 = 2, τ 2 = 1, τ 3 = 3, and τ 4 = 1, as shown in Figure 12 (right). This implies that different combinations of intracellular delays may be a type of infection control.  Figure 10. The immunity−activated equilibrium E 2 is asymptotically stable when R 1 = 4.5504 > 1,  Figure 11. When R 1 = 4.5504 > 1, τ 5 = 3 > τ 0 , τ i = 0 (i = 1, 2, 3, 4), a periodic solution is bifurcated from the unstable immunity−activated equilibrium E 2 .

Figure 12.
On the left is the bifurcation diagram of Model (2) with respect to τ 5 when τ i = 0 (i = 1, 2, 3, 4). On the right is the bifurcation diagram of Model (2) with respect to τ 5 when τ 1 = 2, τ 2 = 1, τ 3 = 3, τ 4 = 1 are fixed. We can see that the stable interval is longer on the right compared to the left, which implies that different combinations of intracellular delays may be a type of infection control.

Figure 15.
On the left is the bifurcation diagram of Model (2) with respect to τ 5 when τ i = 0 (i = 1, 2, 3, 4). On the right is the bifurcation diagram of Model (2) with respect to τ 5 when τ 1 = 2, τ 2 = 3, τ 3 = 6, τ 4 = 9 are fixed. We can see that different oscillations occur, and the stable interval on the right is smaller compared to that on the left.

Conclusions
This paper investigated the dynamical behaviors of a class of viral infection models with cell-to-cell infection, general nonlinear incidence, and multiple delays. The aim of this work was to study the effect of intracellular delays and the immune delay on the dynamics of the model. The threshold dynamics of the equilibria were obtained by constructing Lyapunov functionals. It was found that the infection-free equilibrium E 0 is globally asymptotically stable when R 0 < 1, which means that the infection dies out. In addition, the immunity-inactivated equilibrium E 1 is globally asymptotically stable when R 1 < 1 < R 0 , which implies that the virus load is insufficient to activate humoral immunity. Moreover, the immunity-activated equilibrium E 2 is globally asymptotically stable when R 1 > 1 in the absence of an immune delay, which means that the virus can coexist with immune cells and reach a balance within the host. By choosing the immune delay as the bifurcation parameter, we found that Model (2) generated a Hopf bifurcation, and periodic oscillations and stability switches occurred, which indicates that the immune delay can destabilize the immunity-activated equilibrium. The properties of the Hopf bifurcation were investigated by applying the normal theory and center manifold theorem. Moreover, the global existence of the Hopf bifurcation was studied. Finally, numerical simulations were carried out to show how the delays impact the dynamics of Model (2).
The obtained results imply that both the intracellular delays and immune delay are responsible for the rich dynamics of the model. However, it follows from the bifurcation diagrams in Figures 6-8, 12 and 15 that the immune delay is the main factor for the existence of the Hopf bifurcation and dominates over the intracellular delays in this viral infection model. This implies that the immune system itself has very complicated procedures during virus infection. Moreover, from a biological perspective, stability and stability switches are important, as they determine whether the eradication of infection may be possible in the case of stable regions, whereas it will be difficult in unstable regions. Thus, the stable regions can provide some insights into infection control. In addition, from a mathematical perspective, we can also provide some theoretical suggestions. For example, it was observed in the bifurcation diagrams that the immune delay may be the main reason for the bifurcation, and the model easily reached a stable state for a small immune delay. This indicates that developing effective drugs to decrease the time for activating the immune response system can contribute to infection control. Furthermore, the expression of the basic reproduction number shows that the infection risk will be underevaluated when ignoring either the cell-to-cell infection or virus-to-cell infection. Thus, developing effective drugs for preventing cell-to-cell infection should be considered, except for virus-to-cell infection.
The theoretical analysis of this virus infection model, which incorporates cellular infection, latency, immune response, general nonlinear incidence, and multiple delays, is rare. The investigation of the global stabilities of the corresponding equilibria by applying the Lyapunov method is a generalization of some existing models. Note that only the humoral immune response has been taken into account in this model. However, both the CTL immune response and humoral immune response are the two main immune mechanisms. Also, it takes some time to build up the CTL immune response, so time delays exist during the activation of the CTL immune response. Thus, exploring the dynamics of a model that incorporates both kinds of immune responses and their corresponding immune delays is an interesting project. We leave this for future work.

Data Availability Statement:
No data were used to support this study.

Conflicts of Interest:
The authors declare no conflict of interest.