Bifurcation Analysis of a Delayed Infection Model with General Incidence Function

In this paper, an infection model with delay and general incidence function is formulated and analyzed. Theoretical results reveal that positive equilibrium may lose its stability, and Hopf bifurcation occurs when choosing delay as the bifurcation parameter. The direction of Hopf bifurcation and the stability of the periodic solutions are also discussed. Furthermore, to illustrate the numerous changes in the local stability and instability of the positive equilibrium, we conduct numerical simulations by using four different types of functional incidence, i.e., bilinear incidence, saturation incidence, Beddington–DeAngelis response, and Hattaf–Yousfi response. Rich dynamics of the model, such as Hopf bifurcations and chaotic solutions, are presented numerically.


Introduction
Mathematical modeling has been proven to be valuable in exploring mechanisms and dynamical behaviors of the viral infection process. e analysis of these models can provide insights into developing effective control strategies for infections and evaluating antiviral therapies. Based on the basic virus infection model introduced by Nowak et al. [1], various models have been formulated by many authors to describe the dynamics of virus population qualitatively and quantitatively, and lots of interesting phenomena have been discussed, such as the global stability of the equilibria, bifurcations, periodic oscillations, limit cycles, and the effects of time delays.
For models of virus infections, it is observed that the time delays introduced into them can predict the viral infection process better when compared to models without delays, such as the time between viral entry into a target cell and the production of subsequent viral particles, the time necessary for the newly produced virus to become mature and then infectious particles, and that needed to activate the immune response, all of which cannot be ignored when describing the interaction of them [2,3]. In the field of virus dynamics, many models with discrete or distributed delays have been studied, showing that rich dynamics including stability switches, Hopf bifurcation, and chaotic oscillations can be generated by the mechanisms of time delays. Specifically, stability switches can occur due to a delayed immune response and exponentially decayed delay-dependent parameters [4][5][6][7][8][9], which leads to the switches in the stability of equilibria because of the change in the value of parameters related with these delays.
It is known that the function forms of the incidence rate of the infection have a crucial role in the modeling of the virus dynamics, which are important in determining qualitative behavior of the proposed models and in giving a reasonable description of the dynamics. Huang et al. [10] and Zhang and Xu [11] formulated a virus dynamical model with Beddington-DeAngelis (BD) infection rate βxυ/ (1 + a 1 x + a 2 υ), where x represents the concentration of uninfected cells and υ represents that of virus, with a 1 , a 2 ≥ 0 being constants. is function is similar to the well-known Holling type II functional response, and the term a 2 υ in the denominator reflects the mutual interference between viruses. In [12], Zhou and Cui used a Crowley-Martin (CM) function response of the form βxυ/(1 + a 1 x+ a 2 υ + a 1 a 2 xυ), with a 1 , a 2 , a 3 ≥ 0, which are constants. It can be seen that in BD and CM functional responses, when a 1 > 0, a 2 � 0, these two functions are simplified to the Holling type II functional response. And when a 1 � 0, a 2 > 0, they express a saturation response. Moreover, when a 1 � 0, a 2 � 0, they are the mass action process (or Holling type I functional response). Other more general incidence rates have also been proposed. McCluskey and Yang used function f(x, υ) as incidence rate under some biologically motivated assumptions to study the global stability of a virus model [13]. Hattaf et al. in [14,15] used an incidence rate f(x, y, υ) that covers a variety of incidence functions such as BD response when a 3 � 0 and CM response when a 3 � a 1 a 2 if f(x, y, υ) � βx/(1+ a 1 x + a 2 υ + a 3 xυ), which will be called Hattaf-Yousfi response in the following, with y being infected cells.
In this paper, we will study the influence of general incidence function and time delay on the dynamical behaviors of an infection model. Motivated by the works of [7,13], we will focus on the dynamics of a delayed model that incorporates the immune response in antiviral defense, time delay in activating the immune response for the virus, and two general incidence functions for the transmission of virus-to-cell and cell-to-cell, respectively. e main aim of this paper is to study stability switches of the positive equilibrium, which are generated by the mechanism of time delay, implying the existence of changes in local stability and the instability of this equilibrium. By using the characteristic equation with delay-dependent parameters, normal theory, and center manifold theorem, we can indicate the existence of pure imaginary roots and then Hopf bifurcation. Furthermore, to verify the theoretical results, numerical simulations are performed for four different forms of incidence functions, including bilinear incidence, saturation incidence, Beddington-DeAngelis response, and Hattaf-Yousfi functional response, respectively. e organization of the paper is as follows. In the next section, we introduce the virus dynamic model. e boundedness and nonnegativity of solutions of the model as well as the existence and uniqueness of equilibria are discussed in Section 3. In Section 4, both the stabilities of two equilibria and the conditions for the existence of Hopf bifurcation are studied. Furthermore, numerical simulations are presented in Section 5. A brief summary is given in Section 6, and the properties of the Hopf bifurcation solutions have been provided in Appendix.

Model Formulation
To account for the possible effect of the latent period in viral infection, we introduce a constant time delay τ to represent the time needed to activate the immune response. In this model, two transmissions of virus-to-cell and cell-to-cell have been both taken into consideration, and we hope to employ the function of incidence rate given by more general form. Consequently, the recruitment of infected cells is given by two functions, i.e., f 1 (u, v)v and f 2 (u, w)w, which were used in [16][17][18], here u, w, and v denote the concentrations of uninfected cells, infected cells, and free virus particles. z is used to denote the concentration of CTL cells. Assume that the uninfected cells are produced at a constant rate s. d 1 , d 2 , d 3 , and d 4 represent the death rates of uninfected cells, infected cells, free virus, and CTL cells. p denotes the killing rate of infected cells by CTL cells. k and c are the production rates of free viral particles and the effector cells, respectively. en, we construct the following model: 3,4), and all constant parameters are positive.
of model (1) with initial condition (2) follow from the standard theory of functional differential equations [19]. Furthermore, X is positively invariant for (1). In the following, we discuss the boundedness [20] and nonnegativity of solutions of model (1) as well as the existence of equilibria.
e solutions Y(t, φ) of model (1) with initial condition (2) are nonnegative and ultimately uniformly bounded for all t ≥ 0.
us, it is valid that u(t) > 0 for all t ≥ 0. In order to show w(t) > 0, v(t) > 0, and z(t) > 0 for all t ≥ 0, assume that there exists t 1 > 0 such that min w(t 1 ), v(t 1 ), z(t 1 ) � 0 for the first time. Firstly, if w(t 1 ) � 0, then from the second equation of (1); we have and then there exists ε 1 small enough such that it is obvious that v(t 1 ) � 0 or z(t 1 ) � 0 contradicts with the above expression about v(t) and z(t) at t � t 1 . erefore, the solution Y(t, φ) with initial condition (2) is nonnegative for all t > 0.
For convenience, we denote Define the basic reproductive number as From biological viewpoint, R 0 can be divided into two parts, with R 1 0 measuring the average number of secondary infected generation caused by an existing free virus and R 2 0 being that caused by an infected cell, which give the basic reproduction number corresponding to virus-to-cell and cell-to-cell infections, respectively.
It is easy to obtain v * � (k/d 3 )w * and z * � (c/d 4 )w * . By the first and second equations in (1), we have which gives In order to get the infected steady state, it is necessary that From the equation of w, it follows that Denote the left side of (11) as G(w * ), and then the positive equilibrium of model (1) are given by G(w * ) � 0 for w * ∈ (0, w]. Noting that G(0) � 0 and for u � 0 when w � w, and it is sufficient to show that G ′ (0) > 0 if there is a positive equilibrium. In fact, when and then G(w * ) is positive for sufficiently small w * ; therefore, there exists an infected steady state E * . Now we examine the derivative of and thus, we have

Computational and Mathematical Methods in Medicine
denoting the sum of negative terms, and thus function G is strictly decreasing at each of its zeros. Suppose there exists more than one infected steady state, then there must exist an equilibrium which contradicts with the above fact. erefore, there is only one positive equilibrium E * when the basic reproductive number R 0 > 1, and it is easy to verify that there exists no infection equilibrium when R 0 < 1. en the following result is obtained.
When R 0 > 1, by using the persistence theory in [21], it can be shown that system (1) is uniformly persistent.

Stability and Hopf Bifurcation
In this section, we will study the stability of the two equilibria and the conditions for the existence of Hopf bifurcation.

4.1.
Stability of E 0 . Firstly, we linearize model (1) at the steady states to discuss its local stability. e Jacobian matrix leads to the following characteristic equation: and thus, the linearization of system (1) at E 0 can be expressed by It is clear that (21) has two negative real roots λ 1 � −d 1 and λ 2 � −d 4 . e stability of E 0 is determined completely by using the following equation: For the case R 0 < 1, the roots of (22) have only negative real parts. en the infection-free equilibrium E 0 is locally asymptotically stable. Otherwise, when R 0 > 1, equation (22) has at least one root with positive real part, which implies that E 0 is unstable. (1) is globally asymptotically stable, and it is unstable

Proof. For a continuous and bounded function
en for any solution of (1), we have By using the fluctuation lemma in [22], we know that there is a sequence t n with t n ⟶ ∞ such that w t n ⟶ w ∞ , and w ′ t n ⟶ 0, n ⟶ ∞.

(25)
Substituting the sequence t n into the first equation of (1) and taking the limit gives and then we have u ∞ ≤ (λ/d 1 ) � u 0 . A similar argument to the second and third equations of model (1) yields and By the assumption of (4), combining equality (28) into  (27) gives which leads to d 2 w ∞ ≤ d 2 R 0 w ∞ . Noticing that w ∞ is nonnegative since it is the supremum of the function w(t), then there are possible cases of w ∞ > 0 or w ∞ � 0. If w ∞ > 0, then we have R 0 ≥ 1, which gives a contradiction. erefore, w ∞ � 0 is valid. From (28), we have v ∞ � 0. When considering the fourth equation about z(t) in model (1), we obtain z ∞ � 0. As for u(t), we can find that lim t⟶∞ u(t) � s/d 1 by applying the limit theory to the first equation in system (1), and thus, the proof is completed.
In the following calculation, we will use the following equilibrium equations: en, we have Linearized model (1) at the infection equilibrium E * and the characteristic equation can be expressed as and B 2 � pw * c, It can be seen that A i (i � 0, 1, 2, 3) and B j (j � 0, 1, 2) are all positive. When the delay τ � 0, A i (i � 0, 1, 2, 3) and B j (j � 0, 1, 2) are independent with τ and the characteristic equation (33) can be reduced to Computational and Mathematical Methods in Medicine 5 By using the Routh-Hurwitz criterion, we know that all solutions of (35) have negative real parts if and only if the following conditions are satisfied: Let Theorem 4. If (37) holds and τ � 0, then the infection equilibrium E * is globally asymptotically stable when R 0 > 1.
However, if the delay τ ≠ 0, it may destabilize the infected steady state and lead to Hopf bifurcation, which will be discussed in the following subsection.

Hopf Bifurcation.
It is shown that all roots of (33) locate in the left side of the imaginary axis if τ � 0 and (36) is satisfied, which leads to local stability of the infected equilibrium E * . When τ increases from 0 to variant positive delays, it is possible that the roots of (35) pass through the imaginary axis and enter the right side in the complex plane. erefore, the stability switch may occur, and it is necessary to study the transcendental equation (33) when λ � iω, which is a critical value under small perturbation. In the following analysis, we will study the occurrence of any possible stability switching in this case.
Substituting λ � iω into (33) and separating the real and imaginary parts give

⎧ ⎨ ⎩ (40)
Squaring and adding the two equalities lead to where Note that equation F(ω) � 0 has at least one positive root when M 0 < 0, i.e., A 0 < B 0 . And let z � ω 2 , then the characteristic equation (33) has a purely imaginary root iω which is equivalent to F(z) � 0 has a positive real root z. Take the transformation y � z + (M 3 /4); then for F ′ (z) � 0, it is equivalent to y 3 + m 1 y + m 0 � 0 with m 1 � (M 2 /2) − (3M Then for j � 0, 1, 2, 3, . . ., we can obtain the following expression of delay τ: From this lemma, we can see that the infection equilibrium is asymptotically stable for all τ ≥ 0 if the conditions (c 1 ) ∼ (c 3 ) are not satisfied. Otherwise, if one of these three conditions is satisfied, then the infection equilibrium is asymptotically stable for τ ∈ [0, min i,j τ i j ), and a Hopf bifurcation can occur at this equilibrium when τ � τ * with τ * being the critical value of τ i j .
Theorem 5. Suppose (36) holds, when τ * � τ i j (correspondingly ω * � ω i for some i � 1, 2, 3, 4), then the characteristic equation (35) admits a pair of simple conjugate pure imaginary roots λ � iω * and λ � −iω * , which crosses the imaginary axis from left to right (from right to left) if Proof. Differentiating the characteristic equation (33) with respect to delay τ and arranging it can give the expression of (dλ/dτ) −1 as Meanwhile, noticing the fact and then some computation leads to which gives the desired result. When model (1) undergoes stability switch at E * for τ � τ * , the conditions and theorem for the direction and stability of Hopf bifurcation can be discussed and proven by the normal theory and center manifold theorem in [24], which are presented in Appendix.

Numerical Simulation
In model (1), the incidence rates for transmission of virusto-cell and cell-to-cell are taken as the general form. By choosing four specific types of functions, i.e., bilinear incidence, saturation incidence, Beddington-DeAngelis response, and Hattaf-Yousfi response, we conduct numerical simulations to investigate the complex dynamics that model (1) can have, when taking time delay τ as the bifurcation parameter. e values of parameters, d 3 � 2.4, d 4 � 1.618, p � 0.812, and k � 200, remain the same in the following different cases.

Case 1.
Consider the functions of f 1 (u, v) � β 1 u and f 2 (u, w) � β 2 u, then the incidence rates have the bilinear forms of β 1 uv and β 2 uw. We take the other parameters as c � 0.05, d 1 � 0.01, d 2 � 0.4, β 1 � 0.00025, and β 2 � 0.00065. Figure 1 presents the different results for varying producing rate of uninfected cells, plotted with s � 2 (Figure 1(a)) and s � 10 (Figure 1(b)), which shows the effect of parameter s on dynamical behavior of the model in this case. In this figure, we can see that there exists periodic solutions bifurcated from the infection equilibrium. e maximal and minimal values of v(t) are denoted by the two curves, and the line implies the local stability of E * . It is clearly shown that the value of v(t) increases with that of s. Figure 2 shows that the infection equilibrium E * is asymptotically stable when τ � 15 and a periodic    solution exists when τ � 20, which corresponds to the bifurcation diagram in Figure 1(b). Case 2. e type of saturation incidence is used in this case, i.e., f 1 (u, v) � β 1 u/(1 + av) and f 2 (u, w) � β 2 u/(1 + bw); then the dynamical behaviors of the model are simulated in Figures 3 and 4, with the parameters taken as c � 1, d 1 � 0.01, and d 2 � 0.03. In Figures 3(a), we set s � 10, β 1 � 0.002, β 2 � 0.003, a � 1, and b � 1, and in Figure 3   Computational and Mathematical Methods in Medicine v * � 198.0481, and the infected steady state E * is locally asymptotically stable when τ < 4.7403. More complicated dynamics for model (1) are caused by the delay perturbation in Figure 3(b). A period-three solution occurs for τ � 35 and τ � 45, which are displayed clearly by time series of v(t) for the time interval [5000, 5500] and their corresponding phase portraits of v(t) and u(t) plotted in Figure 4, respectively. Case 3. Beddington-DeAngelis response is used in this case; here, f 1 (u, v) � β 1 u/(1 + a 1 u + a 2 v) and f 2 (u, w) � β 2 u/ (1 + a 1 u + a 2 w). By this type of incidence function, chaotic motions also occur when increasing τ from 0 to 50 with s � 20, c � 1, d 1 � 0.1, d 2 � 0.3, β 1 � 0.02, and β 2 � 0.03. When altering parameter values of a 1 and a 2 , it can lead to the different bifurcation diagrams, which can be seen in Figures 5  and 6. We take a 1 � 0.08, a 2 � 0.5 and a 1 � 0.01, a 2 � 0.2, respectively, and then the simulation results in Figure 5 show the effect of these two parameters on the system. From Figure 5, we can see that model (1) Figure 7(b) shows that system (1) can have solution of period three. Figure 8 gives the time series of v(t) and phase portrait of v(t) and u(t) for τ � 49 to illustrate the solution of system (1), which corresponds to Figure 7(b).

Discussion
To model the mechanisms of infectious diseases mathematically and explore dynamical behaviors of infection processes, some types of functional incidence have been used in [10,[12][13][14][15], which play an important role in determining qualitative behaviors of the proposed models and in giving reasonable descriptions of the dynamics. General incidence rate has been introduced into many epidemic models with an aim to include different situations as much as possible. In our model, we incorporate the effect of time delay needed to activate the immune response for the virus and two general incidence functions for the transmission of virus-to-cell and cell-to-cell. e basic reproductive number, which works as a crucial threshold and determines the dynamics of the model, has been defined as the sum of two parts related with infections of virus-to-cell and cell-to-cell, respectively.

Computational and Mathematical Methods in Medicine
In this paper, we are concerned with the stabilities of two equilibria and the existence of Hopf bifurcation through which the positive equilibrium loses its stability and periodic solutions occur. In the analysis, the time delay is chosen as the bifurcation parameter, which can destabilize the positive equilibrium when it increases. By using the characteristic equation with delay-dependent parameters, normal theory, and center manifold theorem, the existence of pure imaginary roots is verified and then the system experiences Hopf bifurcation. To illustrate the dynamical behavior of stability switches, simulations are conducted to show the process numerically when taking general incidence function as four specific types of bilinear incidence, saturation incidence, Beddington-DeAngelis response, and Hattaf-Yousfi response. In each case, the dynamical behavior is simulated when increasing τ from 0 to 50, and two sets of parameters are taken to compare the bifurcation results. Furthermore, for some fixed values of τ, we also present the plot of time series, phase portrait, or solution orbit to show the complicated dynamics. e effect of delay perturbation and different forms of incidence can be seen in the figures.
It should be noted that when other factors are taken into consideration for more realistic mechanism, such as the delays describing the intracellular latency for virus-to-cell infection and cell-to-cell infection, or varying producing rate of the uninfected cells instead of constant, the model can be extended reasonably and the dynamical analysis will become more challenging than the present one, which we may discuss in future work.