Stability and Hopf bifurcation of an HIV infection model with two time delays

1.   Introduction

It is all known that acquired immunodeficiency syndrome (AIDS), an incurable chronic infectious disease, is caused by the human immunodeficiency virus (HIV). The dynamics of HIV in a host can be broadly divided into three stages ^[1]. First of all, after the human body infection with HIV, CD4$^{+}$T cells in the body decrease dramatically and the viral load reaches a sharp peak, which is called the acute stage. Secondly, the body is subjected to a prolonged asymptomatic phase in which the viral load slowly increases as the body's immune system works normally. Finally, the body's immune system is disrupted by HIV, resulting in death from various illnesses. In order to gain a clearer understanding of the disease and develop various drug treatments against it, several HIV models have been proposed ^[1,2]. The classical model of viral infection includes the interaction between four variables, namely $x$, $y$, $v$ and $z$, which stand for the density of healthy target cells, infected cells, free virus particles, and cytotoxic T lymphocyte (CTL), respectively ^[3,4,5]. The specific ordinary differential equations form is as follows

where new target cells are produced at a constant rate $\Lambda$ and die at a rate $d$. And $\alpha$ indicates the transmission rate of viral infection. The mortality of infected cells is $a$. CTLs fight against infection at a rate $p$. $k$ denotes the rate at which infected cell produces virus particles. $u$ is the rate at which virus particles are removed. $f(x, y, z)$ describes the rate of immune response activated, and $b$ is the decay rate of the specific CTL.

Based on the biological phenomenon that the average lifetime of virus particles is usually significantly shorter than that of infected cells, Arnaout et al. ^[6] gave a quasi-steady-state hypothesis. The assumption shows that the viral load reaches a quasi-equilibrium level relatively quickly compared to the slow change in the level of infected cells. According to the third equation $\dot{v}(t) = 0$ in the model (1.1), $v(t) = k y(t)/u$ holds, that is, the concentration of free virus is simply proportional to the concentration of infected cells. Thus, model (1.1) can be modified to

where the rates of cell infection and viral multiplication are both denoted by $\beta = \alpha k/u$. Model (1.2) is also considered an HIV infection model involving cell-to-cell transmission ^[7,8]. Due to the complexity of the human immune system ^{[3,9,10,11,12,13]}, several types of proliferation of the immune response function $f(x, y, z)$ have been proposed: (i) $f(x, y, z) = cy(t)$, CTLs are generated only by the stimulation of the levels of infected cells; (ii) $f(x, y, z) = cy(t)z(t)$, the proliferation of CTLs is caused by the interaction between infected cells and their own cells; (iii) $f(x, y, z) = cx(t)y(t)z(t)$, the activation and proliferation of CTLs are dependent on the levels of the three cells mentioned above, i.e., CTLs, infected cells and healthy CD4$^{+}$T target cells. Here when HIV invades the body, it targets the CD4$^{+}$T cells, often referred to as "helper" T cells. These cells can be considered "messengers", or the command centers of the immune system. They signal other immune cells that an invader is to be fought ^{[14,15,16,17]}. Taking form (iii) and assuming that HIV evolves toward higher replication rates, Huang et al. ^[18] proposed and discussed a possible mechanism that enables HIV to escape immune control. Therefore, an HIV infection model with form (iii) may lead to more meaningful and realistic results.

Indeed, when HIV invades the body, it takes time for both the infection and the immune response to occur. Time delays are often used to explain various biological transitions. The intracellular delay describes the latent period between the time when target cells are infected. Since the complexity and uncertainty of the principle by which CTLs are activated by infected cells, the immune response delay is introduced. Mathematical models of HIV infection accounting for time delays have been extensively explored ^{[18,19,20,21,22,23,24,25,26]}. Combining data from clinical experiments, Herz et al. ^[9] first used intracellular delay to describe the time between the initial viral entry into a target cell and subsequent viral production. Wang et al. ^[27] analyzed the dynamic properties of a three-dimensional delay differential equation by choosing $f = cy(t-\tau)$. Huang et al. ^[20] discussed CTLs immune response delay in two different situations $f = cy(t-\tau)z(t)$ and $f = cy(t-\tau)z(t-\tau)$. Extending to a four-dimensional system, Zhu et al. ^[17] took $f = cx(t-\tau)y(t-\tau)z(t-\tau)$ as an example to study the influence of time delay on the system.

Motivated by the works in ^[17,20,21], we suppose that the specific CTL proliferates at a rate $cxyz$. Then two-time delays are considered in the model (1.2). To be specific, the first delay $\tau_1$ characterizes the intracellular latency for cell-to-cell infection and the second delay $\tau_2$ describes the time lag in the activation of CTLs induced by infected cells. The following delay differential equations are investigated.

The dynamics of the delay-induced system convey more intricate scenarios than the dynamics portrayed in the ODE system without delay. This paper aims to discuss the effects of time delays on the local dynamics of model (1.3), such as the effects of time delay on the steady-state and the periodic solution.

The article is arranged as follows: In Section 2, the positivity of the solution and equilibria of the model (1.3) are performed. In Section 3, we give the local stability analysis and Hopf bifurcation results. In Section 4, the normal form near the Hopf bifurcation is derived by a series of calculations. We exemplify the obtained analytical results by numerical simulations in Section 5. In Section 6, we conclude with the obtained analysis results.

2.   Positivity of solution and equilibrium

From a biological model with practical implications, the non-negative initial conditions of model (1.3) need satisfy

where $\tau = \max\{\tau_1, \tau_2\}$, $\phi = (\phi_1, \phi_2, \phi_3)^{\top}\in \mathbb{C}$ and $\mathbb{C}$ is the Banach space $\mathbb{C}([-\tau, 0], \mathbb{R}_{+}^{3})$ of all continuous functions that map from $[-\tau, 0]$ into $\mathbb{R}_{+}^{3}$. By the fundamental mathematical theory ^[28], the uniqueness and existence of the solution on $[0, +\infty)$ of model (1.3) with conditions (2.1) hold. In addition, we reach the following conclusion.

Theorem 2.1. Let $S(t, \phi) = (x(t, \phi), y(t, \phi), z(t, \phi))^{\top}$ be a solution of model (1.3) satisfying the initial conditions (2.1). Then $S(t, \phi)$ is non-negative.

Proof. Using the method in ^[29,30], we can rewrite model (1.3) as

where $S_t(\theta) = S(t + \theta)$ for $\theta\in [-\tau, 0]$ and

for $\phi = (\phi_1, \phi_2, \phi_3)^{\top}\in \mathbb{C}$. Clearly for any $\phi\in \mathbb{C}$, $\phi_i(0) = 0$, $i = 1, 2, 3$, we can prove $\mathcal{F}_i(\phi) \geq 0$. Then, on the basis of Theorem 2.1 in ^[29], it can be deduced that $S(t, \phi) \geq 0$ for all $t\geq 0$.

For the convenience of discussion, we apply the following notation in ^[18] to model (1.3)

Here, $R_0$ and $Q_0$ are denoted as the basic reproduction number of infected cells and immune response. The dynamic properties of model (1.3) without time delays are mainly determined by these two threshold parameters, such as the existence and stability of equilibria. Below we give the results already discussed in ^[18].

Lemma 2.1. For model (1.3) with $\tau_1 = 0$ and $\tau_2 = 0$.

$\textrm{(i)}$ There is always an infection-free equilibrium $E_0 = (x_0, y_0, z_0)$, where

When $R_0 < 1$, $E_0$ is the unique equilibrium and stable.

$\textrm{(ii)}$ When $R_0 > 1$, the system also has an immunity-absent equilibrium $E_1 = (x_1, y_1, z_1)$, which is stable if $Q_0 < 1$, where

$\textrm{(iii)}$ When $Q_0 > 1$, an immunity-present equilibrium $E^{*} = (x^{*}, y^{*}, z^{*})$ exists and is stable, where

3.   Local stability and Hopf bifurcation

This section mainly focuses on the effects of time delays on the local stability of the model (1.3) at the equilibrium. By linearizing model (1.3), we obtain the characteristic equation at the equilibrium $E_*(x_*, y_*, z_*)$. The corresponding characteristic equation is as follows

From the above expression (3.1), the stability results for equilibria $E_0$ and $E_1$ can be given as follows.

Theorem 3.1. Consider model (1.3) for any $\tau_1 \geq 0$ and $\tau_2 \geq 0$.

$\textrm{(i)}$ The infection-free equilibrium $E_0$ is locally asymptotically stable if $R_0 < 1$ and unstable if $R_0 > 1$. When $R_0 = 1$, model (1.3) undergoes a fold bifurcation at $E_0$.

$\textrm{(ii)}$ When the immunity-absent equilibrium $E_1$ exists, it is locally asymptotically stable if $Q_0 < 1$ and unstable if $Q_0 > 1$.

Proof. (i) At $E_0 = (x_0, y_0, z_0)$, the characteristic equation (3.1) turns into

and it has roots $\lambda = -b < 0$, $\lambda = -d < 0$ and remaining root will be executed from

We assume that $\lambda$ has non-negative real parts, then take the modulus of both sides of Eq (3.3) so that the left-hand side becomes

but the right-hand side of Eq (3.3) when $R_0 < 1$ becomes

There is a contradiction. Hence $\lambda$ has no non-negative real parts. Since $\lambda = 0$ is not the root of Eq (3.3) when $R_0 < 1$, then all the roots of the characteristic equation of $E_0$ have negative real parts. Therefore, $E_0$ is locally asymptotically stable if $R_0 < 1$. When $R_0 = 1$, the characteristic equation (3.2) has a zero root $\lambda = 0$, and

Therefore, model (1.3) undergoes a fold bifurcation at $E_0$. Besides, Eq (3.3) is equivalent to

Thus, $h(0) = a-\beta x_0 < 0$ when $R_0 > 1$, and $\lim\limits_{\lambda \to +\infty} h(\lambda) = +\infty$. This yields that Eq (3.3) has at least one positive root, which implies that $E_0$ is unstable if $R_0 > 1$.

(ii) At $E_1 = (x_1, y_1, z_1)$, the characteristic equation (3.1) turns into

Obviously, it has an eigenvalue $\lambda = cx_1y_1-b = b(Q_0-1)$. If $Q_0 > 1$, $\lambda > 0$ holds which indicates that $E_1$ is unstable. If $Q_0 < 1$, we have $\lambda < 0$. Next, we analyze the transcendental equation

Clearly, $\lambda = 0$ is not the root of Eq (3.4). Using the same method as above, suppose the real part of $\lambda$ is positive. Then

however,

The occurrence of the contradiction illustrates that the real parts of $\lambda$ is negative. Thus, $E_1$ is locally asymptotically stable when $R_0 > 1$ and $Q_0 < 1$.

In the following contents, the effect of $\tau_1$ and $\tau_2$ on the existence of local Hopf bifurcation of model (1.3) in the immunity-present equilibrium $E^{*}$ will be described in detail. To simplify the analysis, let $X = x(t)- x^{*}$, $Y = y(t)-y^{*}$, $Z = z(t)-z^{*}$, model (1.3) turns into

where

The linearized part is given separately,

Then, related characteristic equation of model (3.5) can be expressed as

where $A_1 = -a_{11}-a_{22}$, $A_2 = a_{11}a_{22}$, $A_3 = -a_{12}a_{23}a_{31}$, $B_1 = -b_{22}$, $B_2 = a_{11}b_{22}-a_{12}b_{21}$, $C_1 = -a_{23}c_{32}$, $C_2 = a_{11}a_{23}c_{32}$. For $\tau_1, \tau_2 \geq 0$, we will classify four cases to discuss the stability of $E^{*}$ respectively.

${\bf{Case}}$ Ⅰ: $\tau_1 = \tau_2 = 0$.

Equation (3.6) becomes

where $A_{11} = d+\beta y^{*} > 0$, $A_{12} = \beta^{2}x^{*} y^{*}+bpz^{*} > 0$, $A_{13} = bdpz^{*} > 0$. Since $A_{11}A_{12}-A_{13} = d\beta^{2}x^{*} y^{*}+\beta y^{*}(\beta^{2}x^{*} y^{*}+bpz^{*}) > 0$ holds, it is proved that the real parts of all roots of Eq (3.7) are positive by the Routh-Hurwitz criterion. Then the immunity-present equilibrium $E^{*}$ is locally asymptotically stable.

${\bf{Case}}$ Ⅱ: $\tau_1 > 0$ and $\tau_2 = 0$.

Equation (3.6) becomes

where $A_{21} = A_1$, $A_{22} = A_2+C_1$, $A_{23} = A_3+C_2$. The stability of equilibrium may be broken by time delay ^[31]. We study how the stability of $E^{*}$ varies from bifurcation parameter $\tau_1$ and results in oscillation. Suppose Eq (3.8) has a pair of pure imaginary roots $\lambda = \pm i\omega_1$ ($\omega_1 > 0$). Solving Eq (3.8) with $\lambda = i\omega_1$ and separating the real and imaginary parts, we get

Squaring the two equations of (3.9) and then adding them together gives

where $D_{21} = A_{21}^{2}-B_1^{2}-2A_{22}$, $D_{22} = A_{22}^{2}-2A_{21}A_{23}-B_2^2$, $D_{23} = A_{23}^2$. Let $h_1 = \omega_1^{2}$, then

and

Then, the roots of $f_1^{\prime}(h_1) = 0$ is

Because of $D_{23} = A_{23}^{2} > 0$, we understand that $f_1^{\prime}(h_1) = 0$ has no positive roots when $\Delta = D_{21}^{2} -3D_{22} < 0$. If hypothesis ${\bf{(H_1)}}$

holds, then $f_1(\omega_{1}^{2}) = 0$ has two positive roots $\omega_{11}^{2}$ and $\omega_{12}^{2}$. Let $\omega_{11}^{2} < \omega_{12}^{2}$ and then $f_1^{\prime}(\omega_{11}^{2}) < 0$ and $f_1^{\prime}(\omega_{12}^{2}) > 0$ (see ^[32]). From Eqs (3.9), we can deduce

Let $\tau_{10} = \min_{k = 1, 2}\{\tau_{1k}^{(0)}\}$. Differentiating both sides of characteristic equation (3.8) with respect to $\tau_1$, we obtain

So we know

Therefore, we assume ${\bf{(H_2)}}$

Theorem 3.2. For model (1.3), when $\tau_1 > 0$ and $\tau_2 = 0$, the following conclusions hold:

$\textrm{(i)}$ if $\Delta < 0$, then $E^{*}$ is locally asymptotically stable for all $\tau_1$,

$\textrm{(ii)}$ if hypothesis ${\bf{(H_1, H_2)}}$ are true, $E^{*}$ is locally asymptotically stable for $\tau_1 \in [0, \tau_{10})$, Hopf bifurcation occurs when $\tau_1 = \tau_{10}$.

${\bf{Case}}$ Ⅲ: $\tau_1 = 0$ and $\tau_2 > 0$.

Equation (3.6) becomes

where $A_{31} = A_1+B_1$, $A_{32} = A_2+B_2$, $A_{33} = A_3$. We do a similar analysis with Case Ⅱ. Defining a purely imaginary root of Eq (3.10) as $\lambda = i\omega_2 (\omega_2 > 0)$ and substitute it into Eq (3.10), by conventional calculation we get

where $D_{31} = A_{31}^{2}-2A_{32}$, $D_{32} = A_{32}^{2}-2A_{31}A_{33}-C_1^2$, $D_{33} = A_{33}^2 -C_2^{2}$. Let $h_2 = \omega_2^{2}$, then

Notice that $D_{33} = A_{3}^2 -C_2^{2} < 0$, which shows that Eq (3.11) has at least one positive root.

In general, suppose that $h_{21}$, $h_{22}$ and $h_{23}$ are the positive roots of (3.11), then $\omega_{2k} = \sqrt{h_{2k}}$, $k = 1, 2, 3$. In the same way, we obtain

Let $\tau_{20} = \min_{k = 1, 2, 3}\{\tau_{2k}^{(0)}\}$, we have

So

Thus, making the hypothesis ${\bf{(H_3)}}$

Theorem 3.3. For model (1.3), when $\tau_1 = 0$ and $\tau_2 > 0$, if hypothesis ${\bf{(H_3)}}$ is true, $E^{*}$ is locally asymptotically stable for $\tau_2 \in [0, \tau_{20})$, Hopf bifurcation occurs when $\tau_2 = \tau_{20}$.

${\bf{Case}}$ Ⅳ: $\tau_1 > 0$ and $\tau_2 > 0$.

Choosing $\tau_2$ as the bifurcation parameter, we take the root of Eq (3.6) to be $\lambda = i\omega_2^{*}(\omega_2^{*} > 0)$. According to the same computing process, it is easy to get

Based on the sum of squares of above equations, we have

where

${\bf{(H_4)}}$: there are finite positive roots $\omega_{2k}^{*}, k = 1, 2, ..., l_1$ for Eq (3.12). Then the critical value is shown as

where

Let $\tau_{20}^{*} = \min\{{\tau_{2k}^{*}}^{(0)}\}$, differentiate equation (3.6) concerning $\tau_2$,

So

where

Thus, there exists a hypothesis ${\bf{(H_5)}}$

Theorem 3.4. For model (1.3), when $\tau_1 > 0$, $\tau_2 > 0$, if the hypothesis ${\bf{(H_4, H_5)}}$ are true, then $E^{*}$ is locally asymptotically stable for $\tau_2 \in [0, \tau_{20}^{*})$, Hopf bifurcation occurs when $\tau_2 = \tau_{20}^{*}$.

4.   Direction and stability of Hopf bifurcation

From the analysis in Section 3, sufficient conditions for model (1.3) to undergo Hopf bifurcation at $E^{*}$ have been obtained. Then we will study the bifurcation properties when $\tau_1 > 0$ and $\tau_2 = \tau_{20}^{*}$ by using the normal form method and center manifold theorem ^[33].

For the sake of discussion, suppose $\tau_1 < \tau_{20}^{*}$, where $\tau_2 \in [0, \tau_{20}^{*})$. Rescaling the time by $t = s\tau_2$, let $\tau_2 = \tau_{20}^{*}+\mu$, $X(s\tau_2) = \bar{X}(s)$, $Y(s\tau_2) = \bar{Y}(s)$, $Z(s\tau_2) = \bar{Z}(s)$. In general, redefining $\bar{X}(s), \bar{Y}(s), \bar{Z}(s)$ as $X(t), Y(t), Z(t)$ and rewriting model (1.3), we get a FDE in $\mathbb{C}([-1, 0], \mathbb{R}^{3})$

where $u(t) = (X(t), Y(t), Z(t))^{\top} \in \mathbb{R}^{3}$. $L_{\mu}(\phi): \mathbb{C} \rightarrow \mathbb{R}^{3}$ and $f\left(\mu, u_{t}\right)$ are described respectively as

where

and

where

By the Riesz representation theorem, there exists a function $\eta(\theta, \mu)$ of bounded variation for $\theta \in[-1, 0]$ such that

for $\phi \in \mathbb{C}\left([-1, 0], \mathbb{R}^{3}\right)$. In fact, we can choose

For $\phi \in \mathbb{C}\left([-1, 0], \mathbb{R}^{3}\right)$, define

and

Then Eq (4.1) turns into

For $\theta \in[-1, 0]$, $\psi \in \mathbb{C}^{1}\left([-1, 0], (\mathbb{R}^{3})^{*}\right)$, define a operator

and a bilinear inner product

where $\eta(\theta) = \eta(\theta, 0)$, $\mathcal{A}(0)$ and $\mathcal{A}^{*}$ are adjoint operators. The eigenvalues of $\mathcal{A}(0)$ are known to be $\pm i\omega_{20}^{*}\tau_{20}^{*}$ by previous discussion which are also eigenvalues of $\mathcal{A}^{*}$. Assume that the eigenvectors of $\mathcal{A}(0)$ and $\mathcal{A}^{*}$ corresponding to the eigenvalues $i\omega_{20}^{*}\tau_{20}^{*}$ and $-i\omega_{20}^{*}\tau_{20}^{*}$ are $q(\theta) = (1, q_2, q_3)^{\top}e^{i \omega_{20}^{*} \tau_{20}^{*}\theta}, \theta \in [-1, 0]$ and $q^{*}(s) = D(1, q_2^{*}, q_3^{*})^{\top}e^{i \omega_{20}^{*}\tau_{20}^{*} s}, s \in [0, 1]$, respectively, such that

Thus, we can figure out

From $\langle q^{*}(s), q(\theta)\rangle = 1$, $\langle q^{*}(s), \bar{q}(\theta)\rangle = 0$, we have

Hence, $\bar{D}^{-1} = 1+q_2\bar{q}_2^{*}+q_3 \bar{q}_3^{*}+\bar{q}_2^{*}(b_{21}+b_{22}q_2)\tau_{1}e^{-i \omega_{20}^{*} \tau_{1}}+c_{32}q_2\bar{q}_3^{*}\tau_{20}^{*}e^{-i \omega_{20}^{*} \tau_{20}^{*}}$.

Next, we compute the center manifold $C_{0}$ at $\mu = 0$. Let $u_{t}$ be the solution of Eq (4.1), define

On the center manifold $C_{0}$, we have

where $m, \bar{m}$ are local coordinates for center manifold $C_{0}$ in the direction of $q^{*}$ and $\bar{q}^{*}$. For the solution $u_{t} \in C_{0}$ of (4.1), there exists $\langle\psi, \mathcal{A} \phi \rangle = \langle\mathcal{A}^{*}\psi, \phi\rangle$ when $\mu = 0$, we get

The above equation is equivalent to

where

then

Defining

we have

Combining (4.3) with (4.4), we obtain

where

Comparing the coefficient with Eq (4.3), we have

with

where $G_{1}$ and $G_{2}$ are governed by

and

After the above analysis, the explicit expression of Eq (4.5) can be evaluated. Then, it is easy to acquire the following critical values:

By the classical bifurcation theorem ^[33], we can state the following theorem:

Theorem 4.1. $\textrm{(i)}$ $\mu_{2}$ determines the direction of Hopf bifurcation, if $\mu_{2} > 0 (<0)$, then Hopf bifurcation is supercritical (subcritical), and the bifurcating periodic orbits of model (1.3) at $E^{*}$ exist for $\tau_2 > \tau_{20}^{*}$,

$\textrm{(ii)}$ $\beta_{2}$ determines the stability of the bifurcating periodic orbits, if $\beta_{2} < 0 (>0)$, then the bifurcating periodic orbits are stable $($unstable$)$,

$\textrm{(iii)}$ $T_{2}$ determines the period of bifurcating periodic solution, if $T_{2} > 0 (<0)$, the period of bifurcating periodic solution increases $($decreases$)$.

5.   Numerical simulations

In this section, we perform some numerical simulations of model (1.3) to validate our analytical results with some fixed parameters. For the set of parameter values: $\Lambda = 1.35, \beta = 0.25, d = 0.25, a = 0.45, p = 0.02, c = 0.15, b = 0.45$, we can calculate $R_0 = 3 > 1$ and $Q_0 = 1.2 > 1$. Then model (1.3) has a unique equilibrium $E^{*} = (2.4, 1.25, 7.5)$.

Following the calculation process shown in ^[34,35], we can denote $\Omega$ as the set with complex conjugate roots. Then we establish

where $P_0(\lambda) = \lambda^{3}+A_1\lambda^{2}+A_2\lambda+A_3$, $P_1(\lambda) = B_1\lambda^{2}+B_2\lambda$, $P_2(\lambda) = C_1\lambda+C_2$ and $\bar{P_1}$ is the conjugate of $P_1$. Through simulations, there exist two positive roots for $F (\omega) = 0$, namely $\omega_{-} \approx 0.0527238$ and $\omega_{+} \approx 0.1435193$ (see Figure 1(a)). Then, we obtain $\Omega = [0.0527238, 0.1435193]$. Figure 1(b)) illustrates the stability switching curves $\mathcal{T}$ in the crossing set $\Omega$.

In Case Ⅰ, $E^{*}$ is locally asymptotically stable for $\tau_1 = \tau_2 = 0$ (see Figure 2).

In Case Ⅱ, fixed $\tau_2 = 0$, then $\Delta = -0.2739636 < 0$, the condition that Hopf bifurcation appears is not satisfied. As shown in Figure 1(b)), no curve intersects the horizontal axis, that is, the length of the $\tau_1$ does not change the dynamical properties of the model (1.3) and $E^{*}$ is stable for all $\tau_1\geq 0$. Below, we take $\tau_1 = 10$ and $\tau_1 = 30$ as examples respectively to give the solution curves of the model (1.3) (see Figure 3).

In Case Ⅲ, fixed $\tau_1 = 0$, we get $\omega_{20} = 0.0185994$, $\tau_{20} = 6.3676306$ and $f_2^{\prime}(\omega_{20}^{2}) = 0.0531887 > 0$, then the transversality condition for Hopf bifurcation is satisfied. Figure 4 illustrates when $\tau_2 = 5 < \tau_{20}$, $E^{*}$ is locally asymptotically stable. Figure 5 illustrates when $\tau_2 = 7 > \tau_{20}$, Hopf bifurcation occurs and period orbit bifurcated from $E^{*}$.

In Case Ⅳ, fixed $\tau_1 = 6.2$, the crucial values for Hopf bifurcation are $\omega_{20}^{*} = 0.0946198$, $\tau_{20}^{*} = 6.3793047$ and $I_{41}+I_{42}+I_{43} = 0.0011718 > 0$, then the transversality condition for Hopf bifurcation is satisfied. So when $\tau_2 = 5 < \tau_{20}^{*}$, $E^{*}$ is locally asymptotically stable; when $\tau_2 = 7 > \tau_{20}^{*}$, Hopf bifurcation occurs and period orbit bifurcated from $E^{*}$ (see Figures 6 and 7). Furthermore, we compute $\mu_2 = -93.4286070 < 0$, $\beta_2 = 1.6297392 > 0$, $T_2 = -1.8245185 < 0$. So the Hopf bifurcation is subcritical, the bifurcating periodic solution is unstable and the period of the bifurcating periodic solution decreases. Besides, Figure 8 illustrates an irregular periodic oscillation of model (1.3) when $\tau_1 = 86.81$ and $\tau_2 = 6.9$.

6.   Conclusions

In this paper, we have studied an HIV infection model with intracellular delay and immune response delay. The intracellular delay $\tau_1$ describes the intracellular latency for cell-to-cell infection. Since the antigenic stimulation generating CTLs may require a time lag, we assume that CTLs produced at time $t$ depends on the number of infected cells at time $t-\tau_2$ and uninfected cells and CTLs at time $t$. For time delays $\tau_1\geq 0$ and $\tau_2\geq 0$, when $R_0 < 1$, the infection-free equilibrium $E_0$ is locally asymptotically stable; when $R_0 > 1$ and $Q_0 < 1$, the immunity-absent equilibrium $E_1$ is locally asymptotically stable. That is, time delays $\tau_1$ and $\tau_2$ have no effect on the stability of the infection-free equilibrium $E_0$ and immunity-absent equilibrium $E_1$. If $Q_0 > 1$, model (1.3) has an immunity-present equilibrium $E^{*}$. In Case Ⅰ, we know that without time delays, $E^{*}$ is stable. In Case Ⅱ, by numerical simulation, we find that $E^{*}$ is locally asymptotically stable for $\tau_1 > 0$. In Case Ⅲ, only one delay $\tau_2$ exists, when $\tau_2$ is sufficiently small, $E^{*}$ is stable, but the periodic solution is bifurcated from $E^{*}$ when the delay crosses the critical value. In Case Ⅳ, two delays coexist, restricting time delay $\tau_1$, we conclude that when $\tau_2$ is within a certain range, $E^{*}$ is stable, but the periodic solution is bifurcated from $E^{*}$ when the delay crosses the critical value. The direction and stability are discussed by the center manifold and normal form in Case Ⅳ. By comparing the four cases, it is concluded that only considering intracellular delay $\tau_1$ will not change the dynamic behavior of the model (1.3), but introducing immune response delay $\tau_2$ will break the stability of the positive equilibrium of the model (1.3) and cause population oscillations.

In fact, depending on the biological phenomena, it is also very interesting to introduce spatial heterogeneity in the model (1.3) or consider a nonautonomous model with drug therapy. Some recent research work suggests that more complex or meaningful dynamical properties may arise for partial differential equations ^[36,37,38]. Therefore, we will expand model (1.3) according to the actual situation and leave it for our further research work.

Acknowledgments

We thank the editor and the referees for their exceptionally helpful comments during the review process. This work is partially supported by the Fundamental Research Funds for the Central Universities, China University of Geosciences (Wuhan) (CUGSX01).

Conflict of interest

The authors declare there is no conflict of interest.