Hopf bifurcation in a CTL-inclusive HIV-1 infection model with two time delays

1.   Introduction

AIDS (acquired immunodeficiency syndrome) is a syndrome caused by HIV (human immunodeficiency virus). HIV is a lentivirus (a subetaoup of retrovirus). It infects vital cells in the human immune system, such as helper T cells (specifically CD4$^+$ T cells), macrophages, and dendritic cells ^[1]. When the number of CD4$^+$ T cells declines below a critical level, cell-mediated immunity is lost, and the body becomes progressively more susceptible to opportunistic infections, leading to the development of AIDS. Mathematical modeling has contributed a lot to the understanding of HIV infection (see, for example, the review by Perelson and Ribeiro ^[2] for within-host models).

In the simplest and earliest models of HIV infection, only the key players were taken into account. These models include uninfected target cells ($T$), productively infected cells ($T^*$), and free viruses ($V$). One such model is described by the following system of ordinary differential equations,

For more detail, we refer the readers to Ribeiro and Perelson ^[3]. Inspired by this model, researchers have proposed many other HIV models by considering, for example, different uninfected target cell growth and incidence, latently infected CD4$^+$ T cells, treatment, drug resistance, and immune response (to name a few, see ^{[4,5,6,7,8,9,10,11,12,13,14,15,16]}).

Time delay is commonly observed in many biological processes. For HIV infection, on the one hand, it roughly takes about 1 day for a newly infected cell to become productive and then to be able to produce new virus particles ^[17]. On the other hand, during CTL response, effector CTLs need time to recognize infected cells and destroy them. Herz et al. ^[18] were the first to introduce an intracellular delay to describe the time between the initial viral entry into a target cell and subsequent viral production. They obtained the effect of the delay on viral load change. Since then, delayed HIV models have attracted the attention of many researchers. See, for example, ^{[19,20,21,22,23,24,25,26]} and the references therein.

In this paper, motivated by the studies in ^[7,10,27], we propose and study the following delayed HIV model,

Here $T(t)$, $T^*(t)$, $V(t)$, and $E(t)$ represent the densities of uninfected CD4$^+$ T-cells, productively infected CD4$^+$ T-cells, free viruses, and immune effectors at time $t$, respectively. As in ^[7,28], $k_1 = ke^{-\alpha \tau_1}$, where $\alpha\in [d, \delta]$ is the death rate of infected cells before becoming productive. $\tau_1$ denotes the time delay between viral entry and viral production while $\tau_2$ stands for the time needed for the CTLs immune response to emerge to control viral replication. The interpretations of the parameters are summarized in Table 1, where their units and ranges will be given in Section 4. The logistic growth in target cells and natural growth of immune effectors combined is a new feature of Model (1.1). Our main focus is on the effects of delays, especially $\tau_2$, on the dynamics of (1.1).

The rest of the paper is organized as follows. In Section 2, we present some preliminary results of (1.1a), which include the positivity and boundedness of solution, the existence of steady states. Then we analyze the stability of steady states and possible Hopf bifurcation in Section 3. We conclude the paper with some numerical simulations to illustrate the main theoretical results.

2.   Preliminaries

The suitable phase space for (1.1) is $\mathcal{C} = \mathcal{C}_1\times \mathcal{C}_2\times \mathcal{C}_1\times R$, where $\mathcal{C}_i = C([-\tau_i, 0], R)$ is the Banach space of all continuous functions from $[-\tau_i, 0]$ to $R$ equipped with the supremum norm, $i = 1$, $2$. The norm on $\mathcal{C}$ is the usual product norm. The nonnegative cone of $\mathcal{C}_i$ is $\mathcal{C}_i^+ = C([-\tau_i, 0], R_+)$. To be biologically meaningful, in the sequel, the initial conditions of (1.1) will be always from $\mathcal{C}_+ = \mathcal{C}_1^+\times \mathcal{C}_2^+\times \mathcal{C}_1^+\times R_+$.

For each $\Phi = (\phi_1, \phi_2, \phi_3, \phi_4)\in \mathcal{C}_+$, by the standard theory of functional differential equations ^[29], Model (1.1) has a unique and global solution through it. For such a solution, we first claim that $T(t) > 0$ for $t > 0$. In fact, it is clear that there exists $t_0 > 0$ such that $T(t) > 0$ for $t\in (0, t_0)$. Suppose to the contrary that there exists $t_1 > t_0$ such that $T(t) > 0$ for $t\in (0, t_1)$ and $T(t_1) = 0$. Then by (1.1), $\frac{dT(t_1)}{dt} = s > 0$ and hence there exists $\varepsilon\in (0, t_1)$ such that $T(t) < 0$ for $t\in (t_1-\varepsilon, t_1)$, a contradiction. This proves the claim. Next, with step-by-step method we show that $T^*(t)\ge 0$ and $V(t)\ge 0$ for $t\ge 0$. Note that, for $t\ge 0$, by (1.1b) and (1.1c), we have

and

respectively. It follows from (2.1) that $T^*(t)\ge 0$ for $t\in [0, \tau_1]$. This, combined with (2.2), gives $V(t)\ge 0$ for $t\in [0, \tau_1]$, which together with (2.1) yields $T^*(t)\ge 0$ for $t\in [0, 2\tau_1]$. In turn from (2.2) we have $V(t)\ge 0$ for $t\in [0, 2\tau_1]$. Continuing this way gives the desired result. Finally, from (1.1d), we get

for $t\ge 0$ and hence $E(t)\ge 0$ for $t\ge 0$. Therefore, the solution of (1.1a) with initial condition in $\mathcal{C}_+$ is nonnegative.

Next, we consider the boundedness of solutions. Firstly, we obtain from (1.1a) that

for $t\ge 0$. It follows that

where

is the unique positive zero of $s-dT+rT(1-\frac{T}{T_{\max}})$. Moreover, if $T(0)\le T_0$ then $T(t)\le T_0$ for $t\ge 0$. Secondly, consider the Lyapunov functional

The derivative of $L_1$ along solutions of (1.1) is

where $d_1 = \min\{\delta, d\}$ and $M_0 = \frac{rT_{\max}+4s}{4}$ ($> 0$). Then $\limsup\limits_{t\to\infty}L_1(t)\le \frac{M_0}{d_1}$. In particular, $\limsup\limits_{t\to\infty}T^*(t)\le \frac{k_1M_0}{kd_1}$. Finally, this combined with (1.1c) and (1.1d) immediately gives

Lastly, we study the lower boundedness of $T$. For any $\varepsilon > 0$, there exists $t_2 > 0$ such that $V(t)\le \frac{N\delta k_1M_0}{ckd_1}+\varepsilon$ for $t\ge t_2$. This, together with (1.1a), gives us

Hence, as $\varepsilon$ is arbitrary, we get

To summarize, we have shown the following result.

Proposition 1. The solutions of (1.1) with initial conditions in $\mathcal{C}_+$ are nonnegative and bounded. Moreover, the region

is a positively invariant and attracting subset of (1.1) in $\mathcal{C}_+$.

In the remaining of this section, we consider the steady states of (1.1). Note that a steady state is a solution of the following system of algebraic equations,

It follows from (2.3c) that $V = \frac{N\delta T^*}{c}$. Sustituting it into (2.3b) gives

Then $T^* = 0$ or $T = \frac {c(\delta+d_xE)}{N\delta k_1}$. When $T^* = 0$, we get the infection-free steady state $P_0 = (T_0, 0, 0, E_0)$, where $E_0 = \frac{\lambda_E}{d_E}$. Now assume $T = \frac {c(\delta+d_xE)}{N\delta k_1}$. Combining it with $E = \frac{\lambda_E+pT^*}{d_E}$ obtained from (2.3d), we can get after a little computation that

Then

Substituting it into (2.3a) yields

where

Note that $G$ always has a positive zero and it only has one positive zero. However, for infected steady states, we have $T > \frac {c(\delta+\frac{d_x\lambda_E}{d_E})}{N\delta k_1}$ from (2.4), or equivalently, $G\big(\frac {c(\delta+\frac{d_x\lambda_E}{d_E})}{N\delta k_1}\big) > 0$ or $\frac {c(\delta+\frac{d_x\lambda_E}{d_E})}{N\delta k_1} < T_0$. Thus there is an infected steady state if and only if $R_0 > 1$, where

In summary we have obtained the following result.

Theorem 2.1. (i) If $R_0\le 1$ then (1.1) only has the infection-free steady state $P_0$.

(ii) If $R_0 > 1$ then, besides $P_0$, (1.1) also has a unique infected steady state $P_1 = (T_1, T^{\ast}_1, V_1, E_1)$, where

Note that, in epidemiology, $R_0$ is called the basic reproduction number, whose expression can also be derived by the procedure in ^[30].

3.   Stability and bifurcation

3.1.   Global stability of $P_0$

We start with the local stability of the infection-free steady state $P_0$.

Theorem 3.1. (i) If $R_0 < 1$, then the infection-free steady state $P_0$ of (1.1) is locally asymptotically stable.

(ii) If $R_0 > 1$, then $P_0$ is unstable.

Proof. The characteristic equation at $P_0$ is

Clearly, $-d_E$ and $-(d-r+\frac{2rT_0}{T_{\max}}) = -\frac{s}{T_0}-\frac{rT_0}{T_{\max}}$ are eigenvalues and both are negative. The other eigenvalues are roots of the following transcendental equation,

Noting

we can rewrite (3.1) as

or equivalently

(ⅰ) Assume $R_0 < 1$. We claim that all roots of (3.3) have negative real parts. Otherwise, (3.3) has a root $\lambda = \sigma+\omega i$ with $\sigma\ge 0$ and $\sigma^2+\omega^2 > 0$ since $0$ is not a root by $R_0 > 1$. Taking moduli of both sides of (3.3) gives

This is impossible as the right side of the above is $> 1$ and $R_0 < 1$. This proves the claim and hence $P_0$ is locally asymptotically stable if $R_0 < 1$.

(ⅱ) Assume $R_0 > 1$. In this case, (3.2) has a positive root. In fact, this follows from the Intermediate Value Theorem and

Therefore, $P_0$ is unstable if $R_0 > 1$. This completes the proof.

In fact, the local stability of $P_0$ implies its global stability.

Theorem 3.2. If $R_0 < 1$, then the infection-free steady state $P_0$ of (1.1) is globally asymptotically stable.

Proof. Define the Lyapunov functional

Then the time derivative of $W_0$ along solutions of (1.1) is

Moreover, $\frac{dW_0}{dt} = 0$ if and only if $T^*(t) = 0$ and $V(t)(T(t)-T_0) = 0$. Then one can see that the largest invariant subset of $\{\frac{dW_0}{dt} = 0\}$ is $\{P_0\}$. By the Lyapunov-LaSalle invariance principle (see [29,Theorem 5.3.1] or [31,Theorem 3.4.7]) and Theorem 3.1, we conclude that if $R_0 < 1$ then $P_0$ is globally asymptotically stable.

3.2.   Stability of $P_1$ and bifurcation analysis

Recall that $P_1$ exists only when $R_0 > 1$, which implies that necessarily $\frac{N\delta k_1T_0}{c(\delta+d_x\frac {\lambda_E}{d_E})} > 1$ as $k_1 = ke^{-\alpha \tau_1}$. The purpose of this paper is to consider the effects of delays on the dynamics. As a result, in the sequel of this section, we always assume that $\frac{N\delta k_1T_0}{c(\delta+d_x\frac {\lambda_E}{d_E})} > 1$ and denote

Then $R_0 > 1$ is equivalent to $\tau_1 < \hat{\tau}_1$.

The characteristic equation at $P_1$ is

In the following, we follow the arguments in ^[23] to first show that $P_1$ is locally stable for $\tau_1\in [0, \hat{\tau}_1)$ and $\tau_2 = 0$. Then for given $\tau_1\in [0, \hat{\tau}_1)$, we discuss the possible bifurcations.

Theorem 3.3. Suppose $\tau_1\in [0, \hat{\tau}_1)$, $\tau_2 = 0$, and $0\le r < \frac {d}{1-\frac {T_1}{T_{\max}}}$. Then the infected steady state $P_1$ is locally asymptotically stable.

Proof. When $\tau_2 = 0$, the characteristic equation (3.4) reduces to

We will prove that all roots of (3.5) have negative real parts in three steps.

Firstly, we show that (3.5) has no roots on the imaginary axis with contradictory arguments. Let $\lambda = i\omega_0$ with $\omega_0\ge 0$ be a root of (3.5). Then

Note that the modulus of the right hand side of (3.6) is

However, since

and

it follows that the modulus of the left hand side of (3.6) is strictly larger than that of its right hand side, a contradiction. Thus we have verified that (3.5) has no roots on the imaginary axis.

Secondly, we show that (3.5) has no nonnegative real roots. Again, by contradiction, assume that (3.5) has a root $\lambda_0\geq0$ and we know that $e^{-\lambda_0\tau_1}\in(e^{-\lambda_0\hat{\tau}_1}, 1]$. Noting

we get from (3.5) that

But

as $N\delta k_1T_1 = c(\delta+d_xE_1)$, which contradicts with (3.7). This proved that (3.5) has no nonnegative real root.

Finally, we claim that there exists $\eta_0 > 0$ such that all roots of (3.5) have negative real parts when $1 < R_0 < 1+\eta_0$. If this is not true, then there exists a sequence of values for the parameters where $R_0 \ (>1) \ \to 1$ such that for each set of values for the parameters there exists a pair of conjugate roots for (3.5) with positive real parts (which follows from the results just proved above). Note that roots of (3.5) having nonnegative real parts are uniformly bounded. Without loss of generality, we suppose that the sequence of the conjugate roots converges to $\alpha_0\pm i \beta_0$, otherwise just consider a subsequence. Then $\alpha_0\ge 0$ and $\alpha_0\pm i\beta_0$ are roots of the characteristic equation of the infection-free steady state $P_0$ when $R_0 = 1$. However, when $R_0 = 1$, this characteristic equation has no roots with nonnegative real parts except the simple root $0$. Then $\alpha_0 = \beta_0 = 0$, which implies that $0$ is a root of at least multiplicity $2$ of this characteristic equation, a contradiction. This proves the claim.

Now the proof is done by noting the fact that the roots of (3.5) depend continuously on the parameters.

Theorem 3.4. If $\tau_1 = \tau_2 = 0$ and $0\leq r < \frac{d}{1-\frac{T_1}{T_{\max}}}$ holds, then the infected steady state $P_1$ is global asymptotically stable.

Proof. We define a Lyapunov functional

Then the time derivative of $W_1(t)$ along solutions of system (1.1) is

It is clear that $\frac{dW_1}{dt} = 0$ if and only if $(T, T^*, V, E) = P_1$. By the Laypunov asymptotic stability theorem, we conclude that $P_1$ of system (1.1) is globally asymptotically stable.

In the following, we study the effects of delays by fixing $\tau_1 = \tau_1^0\in [0, \hat{\tau}_1)$ and changing $\tau_2$. Then the characteristic equation (3.4) can be rewritten as

where

Since

we know that $\lambda = 0$ is not a root of (3.8). Therefore, for stability changes of $P_1$ to occur, we first look for $\tau_2$ where (3.8) has a pair of conjugate roots $\lambda = \pm i\omega(\tau_1^0, \tau_2)$ with $\omega(\tau_1^0, \tau_2) > 0$. Substitute $\lambda = i\omega(\tau_1^0, \tau_2)$ into (3.8) and then separate the real and imaginary parts to obtain

where

It follows that

Using $\sin^2\omega\tau_2+\cos^2\omega\tau_2 = 1$, we see that $\omega(\tau_1^0, \tau_2)$ satisfies $F(\omega, \tau_1^0) = 0$, where

with

Therefore, $\omega(\tau_1^0, \tau_2)$ is independent of $\tau_2$. Denote

which is a finite set. If $I_{\tau_1^0} = \emptyset$ then $P_1$ is stable for $\tau_1 = \tau_1^0$ and $\tau_2\ge 0$. Now, suppose $I_{\tau_1^0}\neq \emptyset$. For example, this is true if $A_0+B_0 < C_0$ since

and $A_0+B_0+C_0 > 0$.

Assume $I_{\tau_1^0} = \{\omega_1(\tau_1^0), \ldots, \omega_{j(\tau_1^0)}(\tau_1^0)\}$. For $j\in \mathbb{N}(\tau_1^0) = \{1, \ldots, j(\tau_1^0)\}$, choose the unique angle $\theta_j(\tau_1^0) \in [0, 2\pi)$ such that

Now, define

Then the characteristic equation (3.8) at $\tau_2 = \tau_{2j}^n$ has a pair of conjugate eigenvalues $\lambda = \pm i\omega_j(\tau_1^0)$ for $j\in \mathbb{N}(\tau_1^0)$ and $n\in \mathbb{N}$. The following result comes from [32,Theorem 2.2].

Theorem 3.5. Let $\tau_1^0\in [0, \hat{\tau}_1)$. Then the following two statements are true.

${\rm (i)}$ If $I_{\tau_1^0} = \emptyset$, then $P_1$ is locallay asymptotically stable for $\tau_1 = \tau_1^0$ and $\tau_2\ge 0$.

${\rm (ii)}$ If $I_{\tau_1^0}\neq \emptyset$, then a pair of simple conjugate pure imaginary roots $\lambda(\tau_{2j}^n(\tau_1^0)) = \pm i\omega_j(\tau_1^0)$ of (3.8) exists at $\tau_2 = \tau_{2j}^n(\tau_1^0)$ for $j\in \mathbb{N}(\tau_1^0)$ and $n\in \mathbb{N}$, which crosses the imaginary axis from left to right if $\delta(\tau_{2j}^n(\tau_1^0)) > 0$ and crosses the imaginary axis from right to left if $\delta(\tau_{2j}^n(\tau_1^0)) < 0$, where

By Theorem 3.5, for any $\tau_1\in [0, \hat{\tau_1})$, there exists a $\tau_2^*(\tau_1)\in (0, \infty]$ such that $P_1$ is locally asymptotically stable for $\tau < \tau_2^*(\tau_1)$.

In general, it is hard to determine whether $I_{\tau_1}$ is empty or not. Moreover, if $I_{\tau_1}\neq \emptyset$ and has more than one element, then Theorem 3.5 indicates that there may be stability switches for $P_1$. To get a clear picture of it, we consider the case where $\tau_1^0 = 0$. Moreover, $F(\omega, 0)$ in (3.10) reduces to $h(\omega^2)$, where

with

In this case, $\delta(\tau_{2j}^n(0))$ in Theorem 3.5 is $\mathrm{sign}(h'(\omega^2_j(0)))$. As a result, for Hopf bifurcation to occur, we only focus on the situations where $h(z)$ defined by (3.12) has simple positive real zeros.

Though a simple calculation gives $b_3 = d_E^2+(c+\delta+d_xE_1)^2+(\frac{s}{T_1}+\frac{rT_1}{T_{\max}})^2 > 0$, we cannot easily get the signs of the other coefficients. By Descartes' rule of sign, the polynomial $h(z)$ has at most three positive real zeros. In fact, Yan and Li ^[33] have obtained the conditions on the existence of at least one positive real zero for $h(z)$. To cite the result, let

where $\rho$ is one of the cubic roots of the complex number $-\frac{q}{2}+\sqrt{\Delta}$ and $\epsilon = \frac{-1+\sqrt{3}i}{2}$. Note that when $\Delta > 0$, $z_1^*$ is the only real zero of $h'(z)$; when $\Delta = 0$, $z_1 = -\frac{b_3}{4}-2\sqrt[3]{\frac{q}{2}}$, $z_2 = z_3 = -\frac{b_3}{4}+\sqrt[3]{\frac{q}{2}}$ are the three real zeros of $h'(z)$; when $\Delta < 0$, $-\frac{b_3}{4}+2\mathrm {Re}\{\rho\}$, $-\frac{b_3}{4}+2\mathrm{Re}\{\rho\epsilon\}$, and $-\frac{b_3}{4}+2\mathrm{Re}\{\rho\bar{\epsilon}\}$ are the three real zeros of $h'(z)$ and we arrange them as $\hat{z}_1 < \hat{z}_2 < z_3^*$.

Proposition 3.1 ([33], Lemma 2.1]).

${\rm (i)}$ If $b_0 < 0$, then $h(z)$ has at least one positive real zero. ${\rm (ii)}$ If $b_0\ge 0$, then $h(z)$ has no positive real zero if one of the following conditions holds.

${\rm (ii-1)}$ $\Delta > 0$ and $z_1^* < 0$;

${\rm (ii-2)}$ $\Delta = 0$ and $z_2^* < 0$;

${\rm (ii-3)}$ $\Delta < 0$ and $z_3^* < 0$.

$\rm (iii)$ If $b_0\ge 0$, then $h(z)$ has at least one positive real zero if one of the following conditions holds.

${\rm (iii-1)}$ $\Delta>0$, $z_1^*>0$ and $h(z_1^*)<0$;

${\rm (iii-2)}$ $\Delta = 0$, $z_2^*>0$ and $h(z_2^*)<0$;

${\rm (iii-3)}$ $\Delta<0$, $z_3^*>0$ and $h(z^*_3)<0$.

When $\tau_1 = 0$, by Proposition 3.1 we can have the following result.

Theorem 3.6. Assume $\tau_1 = 0$ and one of the conditions in statement (ii) of Proposition 3.1 holds. Then the infected steady state $P_1$ is locally asymptotically stable for all $\tau_2\ge 0$.

In the following result, we characterize the situations where $h(z)$ has simple positive real zeros, which is not difficult to see by considering the possible graphs of $h(z)$ and $h'(z)$. Recall that $h(z)$ can have at most three positive real zeros.

Proposition 3.2. For the polynomial $h(z)$ defined by (3.12), the following results hold.

${\rm (i)}$ $h(z)$ has one simple positive zero and no other positive zeros if and only if (H₁): one of the following conditions hold.

${\rm (i-1)}$ $\Delta\ge 0$ and $b_0<0$.

${\rm (i-2)}$ $\Delta>0$, $b_0 = 0$, and $z_1^*>0$.

${\rm (i-3)}$ $\Delta = 0$, $b_0 = 0$ and ($z_2^* = z_1>0$ or $z_2^* = z_2>z_1>0$).

${\rm (i-4)}$ $\Delta<0$, $b_0<0$ and $\hat{z}_2\le 0$.

${\rm (i-5)}$ $\Delta<0$, $b_0<0$, $\hat{z}_2>0$ and $h(\hat{z}_2)<0$.

${\rm (i-6)}$ $\Delta<0$, $b_0<0$, $\hat{z}_2>0$, $h(\hat{z}_2)>0$ and $h(z_3^*)>0$.

${\rm (i-7)}$ $\Delta<0$, $b_0 = 0$ and $\hat{z}_2<0<z_3^*$.

${\rm (i-8)}$ $\Delta<0$, $b_0 = 0$, $\hat{z}_1>0$, $h(\hat{z}_2)>0$ and $h(z_3^*)>0$.

${\rm (i-9)}$ $\Delta<0$, $b_0 = 0$, $\hat{z}_1>0$ and $h(\hat{z}_2)<0$.

$\rm (ii)$ $h(z)$ has two simple positive zeros and no other positive zeros if and only if (H₂) : one of the following conditions holds.

${\rm (ii-1)}$ $\Delta>0$, $b_0>0$, $z_1^*>0$ and $h(z_1^*)<0$.

${\rm (ii-2)}$ $\Delta = 0$, $b_0>0$, $z_2^* = z_1>0$ and $h(z_2^*)<0$.

${\rm (ii-3)}$ $\Delta = 0$, $b_0>0$, $z_2^* = z_2>z_1>0$ and $h(z_1)<0$.

${\rm (ii-4)}$ $\Delta<0$, $b_0 = 0$, $\hat{z}_1\le 0<\hat{z}_2$ and $h(z_3^*)<0$.

${\rm (ii-5)}$ $\Delta<0$, $b_0>0$, $\hat{z}_2\le0<z_3^*$ and $h(z_3^*)<0$.

${\rm (ii-6)}$ $\Delta<0$, $b_0>0$, $\hat{z}_1\le 0<\hat{z}_2$ and $h(z_3^*)<0$.

${\rm (ii-7)}$ $\Delta<0$, $b_0>0$, $\hat{z}_1>0$, $h(\hat{z}_1)>0$ and $h(z_3^*)<0$.

${\rm (ii-8)}$ $\Delta<0$, $b_0>0$, $\hat{z}_1>0$, $h(\hat{z}_1)<0$ and $h(\hat{z}_2)<0$.

${\rm (ii-9)}$ $\Delta<0$, $b_0>0$, $\hat{z}_1>0$, $h(\hat{z}_1)<0$, $h(\hat{z}_2)>0$ and $h(z_3^*)>0$.

$\rm (iii)$ $h(z)$ has three simple positive zeros and no other positive zeros if and only if (H₃): one of the following conditions holds.

${\rm (iii-1)}$ $\Delta<0$, $b_0<0$, $\hat{z}_2>0$, $h(\hat{z}_2)>0$, and $h(z_3^*)<0$.

${\rm (iii-2)}$ $\Delta<0$, $b_0 = 0$, $\hat{z}_1>0$, $h(\hat{z}_2)>0$ and $h(z_3^*)<0$.

If (H₁) holds, let $\tilde{z}>0$ be the unique simple positive zero of $h(z)$ and denote $\tilde{\omega} = \sqrt{\tilde{z}}$. Solving (3.11) to obtain the unique $\tilde{\theta}\in [0, 2\pi)$. Define

As $h'(\tilde{z}) > 0$, we have $\delta (\tau_2^n) = 1$ and hence the following result holds.

Theorem 3.7. Assume that $\tau_1 = 0$ and assumption (H₁) holds. Then there exists a sequence $0 < \tau_2^0 < \tau_2^1 < \tau_2^2 < \cdots$ such that $P_1$ is locally asymptotically stable for $\tau_2\in [0, \tau_2^0)$ and unstable for $\tau > \tau_2^0$, and system (1.1) undergoes a Hopf bifurcation at $P_1$ when $\tau_2 = \tau_2^n$ for $n\in \mathbb{N}$.

Now, assume (H₂) holds. Let $\tilde{z}_2 < \tilde{z}_1$ be the only positive real zeros of $h(z)$, which are also simple. Similarly as for the case of (H₁), we can get two increasing positive sequences $\{\tau_{21}^n\}$ and $\{\tau_{22}^n\}$, associated with $\tilde{z}_1$ and $\tilde{z}_2$, respectively. Since $h'(\tilde{z}_1) > 0 > h'(\tilde{z}_2)$, we easily see that $\tau_{21}^0 < \tau_{22}^0$. Since $\tilde{\omega}_1 = \sqrt{\tilde{z}_1} > \sqrt{\tilde{z}_2} = \tilde{\omega}_2$, we have $\frac{2\pi}{\tilde{\omega}_1} < \frac{2\pi}{\tilde{\omega}_2}$. Thus we define

Such $k$ exists due to $\tau_{21}^{l+1}-\tau_{22}^l\to-\infty$ as $l\to\infty$. Then the first few Hopf bifurcation values can be ordered as

Theorem 3.8. Assume $\tau_1 = 0$ and (H₂) holds. Given $n\in \mathbb{N}$ and $j\in \{1, 2\}$, system (1.1) undergoes Hopf bifurcation at $\tau_2 = \tau_{2j}^n$ if $\tau_{2j}^n\neq \tau_{2(3-j)}^l$ for all $l\in \mathbb{N}$. Furthermore, the stability of $P_1$ switches off (namely, it becomes unstable) when $\tau_2$ crosses $\tau_{21}^0$, $\ldots$, $\tau_{21}^k$ and switches on (namely, it becomes stable) when $\tau_2$ crosses $\tau_{22}^0$, $\ldots$, $\tau_{22}^{k-1}$. In other words, $P_1$ is stable when $\tau_2\in [0, \tau_{21}^0)\cup (\tau_{22}^0, \tau_{21}^1)\cup\cdots \cup (\tau_{22}^{k-1}, \tau_{21}^k)$ and unstable when $\tau_2\in (\tau_{21}^0, \tau_{22}^0)\cup \cdots\cup (\tau_{21}^{k-1}, \tau_{22}^{k-1})\cup (\tau_{21}^k, \infty)$.

We mention that we can study the global continuation of Hopf bifurcation in Theorem 3.8 as in Li and Shu ^[34]. We believe that the Hopf bifurcation branches are bounded and each joins a pair of $\tau_{21}^n$ and $\tau_{22}^n$ for $n\in\mathbb{N}$. As a result, for $\tau_2\in (\tau_{21}^{k+1}, \tau_{22}^k)$, there will be two stable periodic orbits.

Also in Theorem 3.8, we exclude the situation where $\tau_{21}^n = \tau_{22}^l$ for some $n$, $l\in \mathbb{N}$. If this happens, then $n > l$ since $\tau_{21}^0 < \tau_{22}^0$ and $\frac{2\pi}{\tilde{\omega}_1} < \frac{2\pi}{\tilde{\omega}_2}$. In this critical situation, as $\tau_2$ crosses this common critical value, two pairs of purely imaginary eigenvalues $\pm i\tilde{\omega}_1$ and $\pm i\tilde{\omega}_2$ appear and all other eigenvalues have nonzero real parts. Therefore, a double Hopf bifurcation occurs.

Theorem 3.9. Assume $\tau_1 = 0$ and (H₂) holds. If there exist integers $n > l\ge 0$ such that $\tau_{21}^n = \tau_{22}^l = \tau_{20}$, then (1.1) undergoes a double Hopf bifurcation at $P_1$ when $\tau_2 = \tau_{20}$.

When (H₃) holds, we can similarly get three sequences of critical values for $\tau_2$. Similar results as those in Theorem 3.8 and Theorem 3.9 can be obtained. Moreover, the global Hopf bifurcation associated with the third sequence is unbounded (one can refer to Li et al. ^[35] for a similar discussion).

4.   Numerical simulations

In this paper, we rigorously analyzed an HIV-infection model with CTL-immune response and two time delays. The model incorporates a logistic growth term for the target cell growth and a natural resource for the immune effectors. The basic reproduction number $R_0$ played an important in the infection dynamics. If $R_0 < 1$ then the infection-free steady state is globally asymptotically stable. Note that $R_0$ explicitly depends on $\tau_1$. It follows that if $\tau_1$ is large enough then the virus will be cleared. We emphasize that this has not been noted in most existing study. Of course, in the real situation, this delay between viral entry and subsequent viral production usually is not very big. This leads to the complicated dynamics when $R_0 > 1$. In this case, the unique infected steady state could be stable or unstable, depending on the parameter values. In particular, we focused on the effects of time delays. Theoretical results indicate that there can be Hopf bifurcation, double Hopf bifurcation and stability switches.

We conclude this paper with some numerical simulations to illustrate the above mentioned main results. The ranges of the parameters except the delays are summarized in Table 2.

For simplicity, we use the same initial condition $(T_0, T_0^*, V_0, E_0) = (100, 0, 10^{-2}, 0.6)$ in all simulations.

First, we take $s = 10$, $d = 0.01$, $r = 0.03$, $T_{\max} = 1500$, $\alpha = 0.02$, $\delta = 0.3$, $d_x = 0.01$, $N = 21$, $c = 3$, $\lambda_E = 1$, $p = 0.3$, $d_E = 0.1$, $\tau_1 = 1.2$ and $k = 2.4\times10^{-5}$. Then $R_0 = 0.1680 < 1$. By Theorem 3.2, the infection-free steady state $P_0 = (1366, 0, 0, 10)$ is globally asymptotically stable (see Figure 1 with $\tau_2 = 4$).

Next, we only change $k$ to $k = 2.4\times10^{-3}$ and keep the others as above. In this case, we have $\hat{\tau}_1 = 142.2801$. Take $\tau_1 = 1.2\in [0, \hat{\tau}_1)$. Then $R_0 = 16.8038 > 1$. Through numerical calculations, we get the first few critical values $\tau_{21}^0 = 2.9940$ and $\tau_{21}^1 = 32.9360$ associated with $\omega_1 = 0.2098$ and $\tau_{22}^0 = 28.8271$ associated with $\omega_2 = 0.1066$. By Theorem 3.5, the infected steady state $P_1 = (171.7615, 14.8383, 31.1604, 54.5149)$ is locally asymptotically stable for $\tau_2 < \tau_{21}^0$ (see Figure 2 with $\tau_2 = 2$).

In fact, numerical simulations indicates that there is Hopf bifurcation for $\tau_2\in (\tau_{21}^0, \tau_{22}^0)$ (see Figure 3 with $\tau_2 = 5$) and there is a stability switch at $\tau_2 = \tau_{22}^0$, that is, $P_1$ is locally asymptotically stable for $\tau_2\in (\tau_{22}^0, \tau_{21}^1)$ (see Figure 4 with $\tau_2 = 32$).

In the following we illustrate this more clearly with the special case where $\tau_1 = 0$. We distinguish three cases.

Case 1: (ⅱ) of Proposition 3.1 holds. We take $s = 5$, $d = 0.2$, $r = 0.03$, $T_{\max} = 1500$, $\alpha = 0.2$, $\delta = 0.3$, $d_x = 0.01$, $N = 2800$, $c = 15$, $\lambda_E = 1$, $p = 0.3$, $d_E = 0.1$, $\tau_1 = 0$ and $k = 2.4\times10^{-3}$. Then $R_0 = 9.8484 > 1$ and system (1.1) has the unique infected steady state $P_1 = (4.5277, 6.9510,389.2547, 30.8529)$. In this case, $\Delta = 4.1507\times10^{6} > 0$, $b_0 = 0.6078 > 0$, and $z_1^* = -1.9553\times10^{-2} < 0$. This means that (ⅱ-1) of Proposition 3.1 holds. It follows from Theorem 3.6 that the infected steady state $P_1$ is locally asymptotically stable for all $\tau_2\ge 0$ (see Figure 5 with $\tau_2 = 1$).

Case 2: Assumption (H₁) holds. We choose the parameter values $s = 10$, $d = 0.01$, $r = 0.25$, $T_{\max} = 1500$, $\alpha = 0.02$, $\delta = 0.3$, $d_x = 0.01$, $N = 21$, $c = 3$, $\lambda_E = 1$, $p = 0.3$, $d_E = 0.1$ and $k = 2.4\times10^{-4}$. Then $R_0 = 1.8655 > 1$ and system (1.1) has the unique infected steady state $P_1 = (1448.2000, 10.9961, 23.0918, 42.9882)$. Note that $\Delta = -0.2247 < 0$, $b_0 = -6.0220\times10^{-4} < 0$, and $\hat{z}_2 = -0.0442 < 0$, namely, assumption (H₁) (ⅰ-4) holds. We can get $\omega = 0.1519$ and $\tau_2^j = 4.3175+\frac{2j\pi}{\omega}$ for $j\in\mathbb{N}$. It follows from Theorem 3.7 that the infected steady state $P_1$ is locally asymptotically stable for $\tau_2\in [0, \tau_2^0)$ (see Figure 6 with $\tau_2 = 4 < \tau_2^0$) and unstable for $\tau > \tau_2^0$. Moreover, system (1.1) undergoes Hopf bifurcation at $\tau_2 = \tau_2^j$ for $j\in\mathbb{N}$. Figure 7 supports this with $\tau_2 = 5 > \tau_2^0$. Figure 8 provides the bifurcation diagram.

Case 3: Assumption (H₂) holds. This time we replace $k$ with $k = 2.4\times10^{-3}$, and $r$ with $r = 0.03$, and keep the other parameter values as in Case 2. It follows that $R_0 = 17.2119 > 1$ and system (1.1) has the unique infected steady state $P_1 = (168.9145, 15.0443, 31.5930, 55.1329)$. Moreover, $\Delta = -0.4441 < 0$, $b_0 = 3.0321\times10^{-4} > 0$, $\hat{z}_1 = -11.1936 < 0 < \hat{z}_2 = 0.0018 > 0$, and $h(z_3^*) = -9.3650\times10^{-4} < 0$. Therefore, assumption (H₂) (ⅱ-6) holds. In this case, we have $\tilde{\omega}_1 = 0.2845$, $\tilde{\omega}_2 = 0.1401$, $\tau_{21}^{j} = 2.0033+\frac{2j\pi}{\tilde{\omega}_1}$, and $\tau_{22}^j = 22.3237+\frac{2 j\pi}{\tilde{\omega}_2}$ for $j\in \mathbb{N}$. Then the first few Hopf bifurcation values are ordered as $\tau_{21}^0 < \tau_{22}^0 < \tau_{21}^1 < \tau_{21}^2 < \tau_{22}^1 < \cdots$. By Theorem 3.8, the infected steady state $P_1$ is locally asymptotically stable for $\tau_2\in[0, \tau_{21}^0)\cup(\tau_{22}^0, \tau_{21}^1)$ (see Figure 9 and Figure 10 with $\tau_2 = 1.5 < \tau_{22}^0$ and $\tau_2 = 24\in (\tau_{22}^0, \tau_{21}^1)$, respectively) and is unstable for $\tau_2\in(\tau_{21}^0, \tau_{22}^0)\cup(\tau_{21}^1, \infty)$. Thus there are stability switches for $P_1$. Moreover, there are supercritical Hopf bifurcation at $\tau_2 = \tau_{21}^j$ and subcritical Hopf bifurcation at $\tau_2 = \tau_{22}^j$ (see Figures 11 and 12). The numerical simulations also strongly indicate that the global Hopf bifurcation branches are bounded and each branch connects a pair $\tau_{21}^j$ and $\tau_{22}^j$, which we will confirm in a future work.

Acknowledgment

The authors are grateful to the anonymous reviewers for their constructive comments. JW is supported by the Nanhu Scholar Program for Young Scholars of Xinyang Normal University. CQ is supported by Scientific Research Foundation of Graduate School of Xinyang Normal University (No. 2018KYJJ12). XW is supported by the NSFC (No. 11771374), the Nanhu Scholar Program for Young Scholars of Xinyang Normal University, the Program for Science and Technology Innovation Talents in Universities of Henan Province (17HASTIT011). YC is supported by NSERC of Canada.

Conflict of interest

All authors declare no conflicts of interest in this paper.