Viral infection dynamics with mitosis, intracellular delays and immune response

1.   Introduction

Epidemic diseases induced by viral infection have been widely studied in recent decades. Great efforts have been made in mathematical modeling and analysis of viral dynamics ^[1,2,3,4]. These models are described by differential equations involving three compartments: uninfected target cells, infected cells, and free viruses. Some works introduce the models incorporating one more compartment related to human immunity ^[5,6,7]. Immune response in viral dissemination is indispensable for controlling or even eradicating infectious diseases. The immune cells like Cytotoxic T Lymphocyte cells (CTLs) attack and eliminate the infected cells in antiviral defense. Mathematical modeling with immune response can provide us with a more comprehensive understanding of viral kinetics and more effective treatment strategies to contain viral infection.

Virus-to-cell transmission and cell-to-cell transmission are the two main infection modes for within-host viral infection dynamics ^[8,9]. In virus-to-cell infection, free virions infect uninfected target cells. In cell-to-cell infection, viral particles are transferred directly from an infected source cell to a susceptible target cell through the formation of virological synapses ^[10,11]. Prior works on dynamical behavior of human immunodeficiency virus (HIV) and hepatitis B virus (HBV) only considered the virus-to-cell infection routine ^[12,13,14]. However, the cell-to-cell transmission also plays a critical role in viral infection which should not be neglected. Great efforts for examining the mechanism of cell-to-cell infection in cell cultures or infected individuals have also been done in many works ^[8,9,15,16]. Sigal et al. ^[8] showed that the cell-to-cell spread of HIV can still occur even with the presence of antiretroviral therapy. Through experimental-mathematical investigation, Iwami et al. ^[9] demonstrated that cell-to-cell infection mode may account for more than a half of the viral infections. Lai and Zou ^[15] revealed that both two infection modes contribute to the value of the basic reproduction number. Global threshold dynamical analyses were established in ^[16] for a within-host viral infection model incorporating both virus-to-cell and cell-to-cell transmissions.

The biological perspective shows that viral infection or immune response is not instantaneous in vivo, and the time delay is indispensable to account for a series of processes. The infected cells need some time to become active and generate viral particles after initial infection; the newly produced virions may go through a maturation time to acquire infectivity; and antigenic stimulation generating immune cells also requires a time interval ^[17,18,19]. Therefore, it is necessary and important to consider time delays when modeling viral spread and immune response. Chen et al. ^[20] developed an HIV infection model with cellular delay and immune delay, and assumed that the immune cells are produced at a linear rate. They observed that immune delay and cellular delay in viral infection lead to stability switches at the infected equilibrium under certain conditions. A general viral infection model in ^[21] incorporating two infection modes was proposed to examine the global stability of viral dynamics, and provide some general global stability results applicable to immune system dynamics. It also has been shown that the intracellular delays during the processes of viral infection and viral production will lead to periodic oscillations in within-host models only with the right kind of target-cell dynamics, and the time delay during immune response can induce sustained oscillations ^[5,22,23].

We denote $T(t), I(t), V(t)$ and $Z(t)$ as the concentrations of uninfected target cells, actively infected target cells, mature viruses, and virus-specific CTLs at time $t$, respectively. The viral infection model with mitosis of uninfected target cells, two infection modes, and CTL immune response can be described by the following system of delay differential equations

Here, uninfected target cells are assumed to be produced at a constant rate $\lambda$ and die at a per capita rate $d$. The mitosis of uninfected target cells is described by the logistic term $r T(t) (1-T(t)/T_m)$, where $r$ is the intrinsic mitosis rate and $T_m$ is the carrying capacity. Infection of target cells by virus-to-cell transmission and cell-to-cell transmission are assumed to occur at the rates $\beta_1 T(t) V(t)$ and $\beta_2 T(t) I(t)$, respectively. The infected cells produce virions at a rate $kI$ and are cleared by CTLs at the rate of $pIZ$. The CTLs are recruited at a rate $qIZ$. Parameters $\mu_1$, $\mu_2$, and $\mu_3$ are the per capita death rates of infected cells, virions, and CTLs, respectively. Parameters $s_1, s_2$ and $s_3$ denote the death rates of inactively infected cells, immature viruses, and inactively virus-specific CTLs, respectively. The delays during the processes of viral infection, viral production, and CTLs recruitment are $\tau_1$, $\tau_2$ and $\tau_3$, respectively. The term $e^{-s_1 \tau_1}$ accounts for the survival probability of infected cells that are infected at time $t$ and become active at $\tau_1$ time units past the infection. The term $e^{-s_2 \tau_2}$ describes the survival probability that start budding from activated infected cells at time $t$ and become free mature viruses at $\tau_2$ time later. The term $e^{-s_3 \tau_3}$ represents the survival rate of virus-specific CTLs during the delay between cell encounters and subsequent recruitment. All parameters are assumed to be positive.

In our study, we investigate the impact of the intrinsic mitosis rate, and the intracellular delays during the processes of viral infection, viral production, and CTLs recruitment on viral dynamics in vivo incorporating both virus-to-cell and cell-to-cell transmissions. The rest of this paper is organized as follows. In Section 2, we present some preliminary results concerning the well-posedness of model (1.1), the existence of equilibria, and the basic reproduction numbers. Local and global stability analysis of equilibria is established in Section 3. In Section 4, we investigate the stability switches and local and global Hopf bifurcations at the positive equilibrium. In Section 5, we carry out some numerical simulations to examine the applicability of our theoretical results and show rich viral dynamics. A simulation of two-parameter bifurcation analysis is further given to explore the joint impacts on viral dynamics for the interplay between nonlinear target-cell dynamics and the CTLs recruitment delay. Finally, we give a summary and discussion in Section 6.

2.   Preliminary results and the existence of equilibria

To establish the well-posedness of system (1.1), we choose the phase space $\mathcal{C}\times \mathcal{C}\times \mathcal{C}\times \mathcal{C}$, where the Banach space $\mathcal{C} = C\left([-\tau, 0], \mathbb{R}\right)$ is equipped with the supremum norm and $\tau = \max\{\tau_1, \tau_2, \tau_3\}$. As usual, $\phi_t\in \mathcal{C}$ is defined by $\phi_t(\theta) = \phi(t+\theta)$ for $\theta\in [-\tau, 0]$. For biological applications, the initial condition of system (1.1) is given as

where $\mathcal{C}^{+} = \mathcal{C}\left([-\tau, 0], \mathbb{R}_{+}\right)$ is the nonnegative cone of $\mathcal{C}$. The existence and uniqueness of the solution of model (1.1) follow from the standard theory of functional differential equations^[24]. For simplicity, we denote

and

Using a similar argument as that in the proof of [5,Lemma 2.1] or [25,Proposition 2.1], we can obtain the well-posedness of solutions of system (1.1).

Lemma 2.1. The solutions $(T(t), I(t), V(t), Z(t))$ of system (1.1) with initial conditions in $X$ are nonnegative, and the region

is positively invariant and absorbing in $X$, that is, all solutions ultimately enter $\Gamma$.

System (1.1) always has an infection-free equilibrium (IFE) $E_0 = (\overline{T}, 0, 0, 0)$. Based on the linearized model at the IFE and the method of calculating the basic reproduction number in ^[26], we define the basic reproduction number for viral infection of (1.1) as

where $R_0^1$ and $R_0^2$ represent the average numbers of infected cells generated by virus-to-cell infection and cell-to-cell infection, respectively. Besides $E_0$, model (1.1) may admit an immune-inactivated equilibrium (IIE) $E_1 = (T_1, I_1, V_1, 0)$ or an immune-activated equilibrium (IAE) $E_2 = (T_2, I_2, V_2, Z_2)$, where all variables are positive reading

and

To explore the existence of equilibria, we define

Direct calculation yields that $T_2 < \overline{T}$, thus $R_1 < R_0$. Considering the linearized equation of (1.1) at $E_1$, we then obtain the basic reproduction number for immune response as

Note that $R_0$ can be rewritten as $R_0 = \overline{T}/T_1$. Thus we have $I_1 > 0$ if and only if $R_0 > 1$. It follows from the expression of $R_1$ that $Z_2 = \mu_1 (R_1-1)/p$, which indicates that $Z_2 > 0$ if and only if $R_1 > 1$. By the equilibria equations of $E_1$ and $E_2$, we have $\text{Sign}(T_2-T_1) = \text{Sign}(I_1-I_2) = \text{Sign}(V_1-V_2)$. Since $R_1$ can be rewritten as $R_1 = T_2/T_1$, it then implies $\text{Sign}(T_2-T_1) = \text{Sign}(R_1-1)$. Therefore, together with (2.5), we can derive the following result.

Lemma 2.2. If the IIE $E_1$ and IAE $E_2$ of system (1.1) exist, then we have $\mathit{\text{Sign}}(T_2-T_1) = \mathit{\text{Sign}}(I_1-I_2) = \mathit{\text{Sign}}(V_1-V_2) = \mathit{\text{Sign}}(R_1-1) = \mathit{\text{Sign}}(R_{IM}-1)$.

In view of Lemma 2.2, $R_1$ and $R_{IM}$ are equivalent while comparing them to $1$. Since $R_{IM}$ is more biologically relevant, we will employ $R_{IM}$ as the threshold parameter in the following discussion. To summarize, we can state the existence and uniqueness of the equilibria for (1.1).

Theorem 2.1.

$({\rm{i}})$ If $R_0\le 1$, then $E_0 = (\overline{T}, 0, 0, 0)$ is the unique equilibrium for system (1.1).

$({\rm{ii}})$ If $R_{IM}\leq 1 < R_0$, then, besides $E_0$, there is a unique IIE $E_1 = (T_1, I_1, V_1, 0)$, and no IAE.

$({\rm{iii}})$ If $R_{IM} > 1$, then, besides $E_0$ and $E_1$, there is a unique IAE $E_2 = (T_2, I_2, V_2, Z_2)$.

3.   Threshold dynamics

3.1.   Global stability of the IFE $E_0$

We first investigate the global stability of the IFE $E_0$ and obtain the result as follows.

Theorem 3.1. If $R_0\leq 1$, then the IFE $E_0$ of system (1.1) is globally asymptotically stable in $X$, while $E_0$ is unstable if $R_0 > 1$.

Proof. The characteristic equation associated with the linearization of (1.1) at $E_0$ is

where $n'(\overline{T}) = r-d-2 r \overline{T}/T_m$ and

Two eigenvalues are $-\mu_3$ and $n'(\overline{T}) < 0$, and all other eigenvalues are determined by $F_0(\xi) = 0$, which is equivalent to

If $R_0 < 1$, then $F_0(0) = \mu_1 \mu_2(1-R_0) > 0$. Thus $0$ is not the eigenvalue. Suppose $a+bi$ to be a root of $F_0(\xi) = 0$ with $a\geq 0$ and $a^2+b^2 > 0$. Then it follows from (3.2) that

This is a contradiction and hence the IFE $E_0$ is locally asymptotically stable if $R_0 < 1$.

If $R_0 > 1$, then $F_0(0) < 0$ and $\lim\limits_{\xi\to \infty} F_0(\xi) = \infty$. Hence, there exists at least one positive eigenvalue and then $E_0$ is unstable if $R_0 > 1$.

If $R_0 = 1$, then $F_0(0) = 0$ and $0$ is a simple eigenvalue. Using a similar argument as above we can show that all other eigenvalues have negative real parts. Now, we examine the local stability of $E_0$ by using the center manifold theory and the normal forms. Let $U = \{\xi\in\mathbb{C}, \xi \ \text{is an eigenvalue of} \ (3.1) \ \text{with} \ \text{Re}\xi = 0\}$. Then $U = \{0\}$ if $R_0 = 1$, and (1.1) satisfies the nonresonance condition relative to $U$. Let $u(t) = (u_1(t), u_2(t), u_3(t), u_4(t))^T = (\overline{T}-T(t), I(t), V(t), Z(t))^T$. By applying the standard notation in delay differential equations $u_t(\theta) = u(t+\theta)$, we rewrite system (1.1) as an abstract ODE

where $A$ is a linear operator defined as $(A \phi)(\theta) = \phi'(\theta)$ for $\theta\in [-\tau, 0)$ with

and $R$ is a nonlinear operator defined as $(R(\phi))(\theta) = 0$ for $\theta\in [-\tau, 0)$ and

for $\phi = (\phi_1, \phi_2, \phi_3, \phi_4)^T\in \mathcal{C}^4$. For $\psi = (\psi_1, \psi_2, \psi_3, \psi_4)\in \left(C([0, \tau], \mathbb{R})\right)^4$, we define a bilinear inner product

We set $\varphi = (1, \varphi_2, \varphi_3, 0)^T$ and $\psi = (0, 1, \psi_3, 0)$ to be the right and left eigenvectors of the linear operator $A$ with respect to the eigenvalue $0$, respectively, where

We have the decomposition $u_t = z \varphi+y$ such that $\langle \psi, y \rangle = 0$. It implies that $\dot{u}_t = \dot{z} \varphi+\dot{y}$ and $\langle \psi, \dot{y} \rangle = 0$. Together with the facts that $A\varphi = 0$ and $\langle \psi, Ay \rangle = 0$, we have

If the initial value is a small perturbation of $E_0$, then $z$ is also small with a positive initial value $z(0)$ and $y = O(z^2)$. Through Taylor expansion, the normal form of (3.3) at the origin follows

where $\hat{a} = (\varphi_2+\varphi_3 \psi_3+\tau_1 \overline{T}(\beta_1 \varphi_3+\beta_2 \varphi_2) e^{-s_1 \tau_1}+k \tau_2 \varphi_2 \psi_3 e^{-s_2 \tau_2})^{-1} < 0$. By (3.4), the zero solution of (3.5) is locally asymptotically stable, which proves the local asymptotic stability of $E_0$ for (1.1) if $R_0 = 1$.

To prove the global stability of $E_0$ when $R_0\le1$, we construct the following Lyapunov functional $L_0: \Gamma\to \mathbb{R}$ to show the global attractivity of $E_0$.

Calculating the time derivative of $L_0$ along a solution of (1.1) yields

Moreover, the maximal compact invariant set in $\{L_0' = 0\}$ is the singleton $\{E_0\}$. By the LaSalle invariance principle [24,Theorem 5.3.1], $E_0$ is globally attractive in $\Gamma$. Since $\Gamma$ is absorbing in $X$, we conclude that $E_0$ is globally attractive in $X$. Furthermore, the above result combined with the local stability implies that the IFE $E_0$ is globally asymptotically stable in $X$ if $R_0\le1$.

3.2.   Global stability of the IIE $E_1$

To investigate the global stability of $E_1$, we first claim the uniform persistence result of the model (1.1) by using a similar argument as that in the proof of [21,Lemma 4.2] or [27,Theorem 4.1].

Lemma 3.1. Assume that $R_0 > 1$, then there exists an $\eta > 0$ such that $\liminf\limits_{t\to \infty} T(t)\geq \eta$, $\liminf\limits_{t\to \infty} I(t)\geq \eta$ and $\liminf\limits_{t\to \infty} V(t)\geq \eta$ for any solution of (1.1) with initial condition in $X_1$, where

It has been shown in ^[5,25] that the intrinsic mitosis rate $r$ may induce sustained oscillations through Hopf bifurcation. To acquire the global convergence of the IIE $E_1$, we assume that $T_m(r-d) < r T_1$. We now establish global stability of the IIE $E_1$ if $R_{IM}\leq 1 < R_0$.

Theorem 3.2. If $R_{IM}\leq 1 < R_0$ and $T_m(r-d) < r T_1$, then the IIE $E_1$ of system (1.1) is globally asymptotically stable in $X_1$, while $E_1$ is unstable if $R_{IM} > 1$.

Proof. The characteristic equation of the linearized system of (1.1) at $E_1 = (T_1, I_1, V_1, 0)$ is

where $F_1(\xi) = \xi+\mu_3-q I_1 e^{-s_3 \tau_3} e^{-\xi \tau_3}$ and

It follows from Lemma 6 in ^[28] that all roots of $F_1(\xi) = 0$ have negative real parts if and only if $R_{IM} < 1$; there exists at least one positive root if $R_{IM} > 1$; and $0$ is a simple root and all other roots have negative real parts if $R_{IM} = 1$.

Since $F_2(0) = \mu_1 \mu_2(\beta_1 V_1+\beta_2 I_1) > 0$, then $0$ is not the root of $F_2(\xi) = 0$. We now claim that all roots of $F_2(\xi) = 0$ have negative real parts. Otherwise, suppose $a+bi$ to be a root of $F_2(\xi) = 0$ with $a\geq 0$ and $a^2+b^2 > 0$. We rewrite $F_2(\xi) = 0$ as $F_{2L}(\xi) = F_{2R}(\xi)$ with

Note that $T_m(r-d) < r T_1$ indicates $n'(T_1) < 0$. Then we have $|F_{2L}(\xi)| > 1$ and $|F_{2R}(\xi)|\leq 1$. This is a contradiction. Hence, we obtain that all eigenvalues of (3.6) have negative real parts if and only if $R_{IM} < 1$; there exists at least one positive root if $R_{IM} > 1$; and $0$ is a simple eigenvalue and all other eigenvalues have negative real parts if $R_{IM} = 1$. Using a similar manner as that in the proof of Theorem 3.1, we derive the normal form of (1.1) at $E_1$ when $R_{IM} = 1$ as what follows:

Thus, $E_1$ is locally asymptotically stable if $R_{IM}\leq 1 < R_0$ and $T_m(r-d) < r T_1$ hold, while $E_1$ is unstable if $R_{IM} > 1$.

To establish the global stability of $E_1$ if $R_{IM}\leq 1 < R_0$, we introduce a Lyapunov functional $L_1: X_1\to \mathbb{R}$ as

where $\phi = (\phi_1, \phi_2, \phi_3, \phi_4)\in X_1$ and $u(\theta) = \theta-1-\ln\theta$. It follows from Lemma 3.1 that $L_1$ is well-defined in $X_1$. Direct calculation yields

The condition $T_m(r-d) < r T_1$ implies that $(n(T)-n(T_1))(1-T_1/T)\leq 0$ for any $T > 0$. Then $L_1'\leq 0$ if $R_{IM}\leq1$, and the largest compact invariant set in $\{L_1' = 0\}$ is the singleton $\{E_1\}$. By the LaSalle invariance principle ^[24] and the local asymptotical stability, we obtain that $E_1$ is globally asymptotically stable in $X_1$ if $R_{IM}\leq 1 < R_0$ and that $T_m(r-d) < r T_1$.

3.3.   Global stability of the IAE $E_2$

Similarly to Lemma 3.1, we have the following persistence result of model (1.1).

Lemma 3.2. If $R_{IM} > 1$, then there exists an $\epsilon > 0$ such that $\liminf\limits_{t\to \infty} T(t)\geq \epsilon$, $\liminf\limits_{t\to \infty} I(t)\geq \epsilon$, $\liminf\limits_{t\to \infty} V(t)\geq \epsilon$ and $\liminf\limits_{t\to \infty} Z(t)\geq \epsilon$ for any solution of (1.1) with initial condition in $X_2$, where

The existing works ^[20,23] have revealed that the CTLs recruitment delay may generate complicated dynamical behavior such as stability switches and sustained oscillations in viral infection models. To explore the global convergence of the IAE $E_2$, we assume that $\tau_3 = 0$ and $T_m(r-d) < r T_2$.

Theorem 3.3. Consider model (1.1) with $\tau_3 = 0$. If $R_{IM} > 1$ and that $T_m(r-d) < r T_2$ holds, then the unique IAE $E_2$ is globally asymptotically stable in $X_2$.

Proof. The characteristic equation of the linearized system of (1.1) with $\tau_3 = 0$ at $E_2$ reads

which is equivalent to $F_{3L}(\xi) = F_{3R}(\xi)$, where

Assumption $T_m(r-d) < r T_2$ indicates that $n'(T_2) = -d+r-2 r T_2/T_m < 0$. Thus we have $F_3(0) = p \mu_2 \mu_3 Z_2(\beta_1 V_2+\beta_2 I_2-n'(T_2)) > 0$, which implies that $0$ is not an eigenvalue. Suppose that $a+bi$ with $a\geq0$ and $a^2+b^2 > 0$ is an eigenvalue of (3.7). Note that $p Z_2 = \mu_1 (R_1-1)$, then $|F_{3L}(\xi)| > R_1$ and $|F_{3R}(\xi)|\leq R_1$. This contradicts $|F_{3L}(\xi)| = |F_{3R}(\xi)|$ and thus all eigenvalues of (3.7) have negative real parts, which implies that $E_2$ is locally asymptotically stable if $R_{IM} > 1$ and $\tau_3 = 0$ provided that $T_m(r-d) < r T_2$.

We now define the Lyapunov functional $L_2: X_2 \to \mathbb{R}$ as

where $\phi = (\phi_1, \phi_2, \phi_3, \phi_4)\in X_2$. Then the time derivative of $L_2$ along solutions of (1.1) with $\tau_3 = 0$ is

Assumption $T_m(r-d) < r T_2$ implies that $(n(T)-n(T_2))(1-T_2/T)\leq 0$ for any $T > 0$. Hence, $L_2'\leq 0$ if $R_{IM} > 1$, and the largest compact invariant set in $\{L_2' = 0\}$ is the singleton $\{E_2\}$. The LaSalle invariance principle and a similar argument as that in the proof of Theorem 3.2 show that the unique IAE $E_2$ of system (1.1) with $\tau_3 = 0$ is globally asymptotically stable in $X_2$ if $R_{IM} > 1$ and that $T_m(r-d) < r T_2$ holds.

4.   Stability of the IAE $E_2$ and bifurcation analysis

Theorem 3.3 shows the global asymptotic stability of the unique IAE $E_2$ of model (1.1) without CTLs recruitment delay when $R_{IM} > 1$. In this section, we will investigate the stability of $E_2$ and identify parameter regions in which the time delay can destabilize $E_2$ and lead to Hopf bifurcation. Throughout this section, we assume that $R_{IM} > 1$, which guarantees the existence of the IAE $E_2$.

4.1.   Local Hopf bifurcation analyses

To analyze the effect of the CTLs recruitment delay on the dynamics of (1.1), we choose $\tau_3$ as the bifurcation parameter and fix $\tau_1$ and $\tau_2$ as constants, to explore the stability of $E_2$ and identify parameter regions in which the CTLs recruitment delay may induce sustained oscillations through Hopf bifurcation. The method to deal with the case with positive constants $\tau_1, \tau_2$ and the case $\tau_1 = \tau_2 = 0$ are similar. For simplicity, we set $\tau_1 = \tau_2 = 0$. Linearizing system (1.1) at the IAE $E_2$ yields the following characteristic equation

where

Here $\Lambda = \lambda/T_2+r T_2/T_m$. In the absence of delays, that is, $\tau_1 = \tau_2 = \tau_3 = 0$, the characteristic equation (4.1) reduces to

where

Since $R_{IM} > 1$, which indicates $R_1 > 1$ by Lemma 2.2, then $p_i > 0$ for $i = 0, 1, 2$. By the Routh-Hurwitz criterion, all eigenvalues of (4.2) have negative real parts if and only if $q_1: = p_2 p_3-p_1 > 0$ and $q_2: = p_1 q_1-p_3^2 p_0 > 0$. Denote $k_0 = k\beta_1 T_2/\mu_2$, then direct calculation yields

Clearly, $q_1, q_2 > 0$ if $\mu_2 < \beta_2 T_2$. To summarize, we obtain the stability of $E_2$ for system (1.1) without delay.

Lemma 4.1. Consider model (1.1) with $\tau_1 = \tau_2 = \tau_3 = 0$. Assume that $R_{IM} > 1$, then the unique IAE $E_2$ is locally asymptotically stable if and only if $q_1 > 0$ and $q_2 > 0$. In particular, $E_2$ is locally asymptotically stable if $\mu_2 < \beta_2 T_2$.

In view of $a_0+b_0 = \Lambda \mu_1\mu_2\mu_3(R_1-1) > 0$ if $R_{IM} > 1$, $0$ is not an eigenvalue of (4.1) for any $\tau_3\ge0$. Then we investigate the existence of purely imaginary eigenvalues of (4.1) for $\tau_3 > 0$. Substituting $\xi = i \omega$ ($\omega > 0$) into (4.1) and separating the real and imaginary parts, we obtain

By squaring and adding both the above equations, $\pm i \omega$ are a pair of purely imaginary eigenvalues of (4.1) only if $\omega$ is a positive root of the polynomial $F$ defined by

where

Denote $\tau_{3, max}$ as the largest value of $\tau_3$ such that $E_2$ exists, that is $0\le\tau_3 < \tau_{3, max}$. This is true if and only if $R_{IM} > 1$, which is

Define

If $I = \emptyset$, then $E_2$ is locally asymptotically stable for any $\tau_3 \in [0, \tau_{3, max})$. If $I\neq\emptyset$, then let $\omega = \omega(\tau_3) > 0$ be the positive root of $F(\omega(\tau_3), \tau_3) = 0$ and define $\theta(\tau_3)\in [0, 2 \pi)$ such that $\sin \theta(\tau_3) = g_1(\tau_3)$ and $\cos \theta(\tau_3) = g_2(\tau_3)$. Thus it follows

Following the method in ^[29], we define the functions

Hence, $\pm i \omega(\tau_3^*)$ are purely imaginary eigenvalues of (4.1) if and only if $S_n(\tau_3^*) = 0$ for some $n\in \mathbb{N}$. According to [29,Theorem 2.2], we have the following lemma concerning the transversality condition.

Lemma 4.2. Assume that $R_{IM} > 1$ and $I\neq\emptyset$. If $S_n(\tau_3) = 0$ has a positive root $\tau_3^*\in I$ for some $n\in \mathbb{N}$, then (4.1) has a pair of simple purely imaginary roots $\pm i \omega(\tau_3^*)$ with $\omega(\tau_3^*) > 0$ at $\tau_3^*$ and

Furthermore, this pair of simple purely imaginary eigenvalues $\pm i \omega(\tau_3^*)$ cross the imaginary axis from left to right (resp. from right to left) at $\tau_3^*$ provided that $\mathit{\text{Sign}} (\mathit{\text{Re}} \xi'(\tau_3^*)) > 0$ (resp. $\mathit{\text{Sign}} (\mathit{\text{Re}} \xi'(\tau_3^*)) < 0$).

If $\sup_{\tau_3\in I} S_0(\tau_3) < 0$, $S_n(\tau_3)$ has no zeros in $[0, \tau_{3, max})$ for all nonnegative integer $n$. This excludes the existence of purely imaginary eigenvalues and thus implies that the stability of $E_2$ does not change for $\tau_3\in[0, \tau_{3, max})$. If $\sup_{\tau_3\in I} S_0(\tau_3) = 0$, $S_0(\tau_3)$ has a unique zero of even multiplicity in $[0, \tau_{3, max})$, denoted by $\tau_3^*$, and $S_0'(\tau_3^*) = 0$. Lemma 4.2 implies that the transversality condition at $\tau_3^*$ is not satisfied and all eigenvalues can not cross the imaginary axis. Thus, the stability of $E_2$ should be the same for $\tau_3\in[0, \tau_{3, max})$.

If $\sup_{\tau_3\in I} S_0(\tau_3) > 0$ and $q_1, q_2 > 0$, it then follows from Lemma 4.1 that $S_0(0) < 0$. Then $S_0(\tau_3)$ has at least one zero, which satisfies the transversality condition. Denote

Then Hopf bifurcation occurs at $\tau_3 = \tau_3^H$, and the stability of $E_2$ changes when $\tau_3$ crosses $\tau_3^H$.

Theorem 4.1. Consider model (1.1) with $\tau_1 = \tau_2 = 0$. Assume that $R_{IM} > 1$, and $q_1 > 0, q_2 > 0$ hold. Let $S_n(\tau_3)$ be defined in (4.8).

$({\rm{i}})$ If $\sup_{\tau_3\in I} S_0(\tau_3)\leq 0$, then $E_2$ is locally asymptotically stable for $\tau_3\in [0, \tau_{3, max})$.

$({\rm{ii}})$ If $\sup_{\tau_3\in I} S_0(\tau_3) > 0$, then the model undergoes a Hopf bifurcation at $E_2$ when $\tau_3 = \tau_3^H$, $E_2$ is locally asymptotically stable for $\tau_3\in[0, \tau_3^H)$, and becomes unstable for $\tau_3\in(\tau_3^H, \tau_3^H+\zeta)$ for some $\zeta > 0$.

In general, the equation $F(\omega, \tau_3) = 0$ could have multiple positive roots. It then follows from Lemma 4.2 that stability switches of $E_2$ may occur. To figure out the existence of simple positive zeros of $F(\omega, \tau_3) = 0$, we denote $z = \omega^2$ and rewrite (4.4) as

Clearly, $F'(z) = 4 z^3+3 c_3 z^2+2 c_2 z+c_1$. Denote $n_1 = c_2/2-3 c_3^2/16$, $n_2 = c_3^3/32-c_2 c_3/8+c_1/4$ and

It is known from Cardano's formulae for the cubic algebra equation that the existence of the real (imaginary) root of $F'(z) = 0$ is determined by the sign of $\Delta$. If $\Delta > 0$, then $F'(z) = 0$ has a unique real root

If $\Delta = 0$, then $F'(z) = 0$ has three real roots

If $\Delta < 0$, then $F'(z) = 0$ has three unequal real roots, denoted by $z_1 < z_2 < z_3$, which are

where $\alpha = (-n_2/2+\sqrt{\Delta})^{1/3}$. $F'(z) = 0$ has a unique simple real root when $\Delta\geq 0$, which implies that $F(z)$ achieves its local minimum at $z_0$. Here, $z_0 = z_1^0$ if $\Delta > 0$ and $z_0 = z_1^*$ if $\Delta = 0$. We now count the number of simple positive zeros for $F(z)$ in the following lemma.

Lemma 4.3. Let $F(z)$ be given as in (4.9) with general real coefficients.

$({\rm{i}})$ $F(z)$ does not have any positive zero if and only if $({\bf H_0})$: one of the following conditions holds.

$({\rm{i}}.1)$ $\Delta < 0$, $c_0\geq0$ and $z_3\leq0$;

$({\rm{i}}.2)$ $\Delta < 0$, $c_0\geq0$, $z_1\leq0 < z_3$ and $F(z_3) > 0$;

$({\rm{i}}.3)$ $\Delta < 0$, $c_0 > 0$, $z_1 > 0$ and $\min\{F(z_1), F(z_3)\} > 0$;

$({\rm{i}}.4)$ $\Delta\geq 0$, $c_0\geq0$ and $z_0\leq0$;

$({\rm{i}}.5)$ $\Delta\geq 0$, $c_0 > 0$, $z_0 > 0$ and $F(z_0) > 0$.

$({\rm{ii}})$ $F(z)$ has a unique simple positive zero and no other positive zeros if and only if $({\bf H_1})$: one of the following conditions holds.

$({\rm{ii}}.1)$ $\Delta < 0$, $c_0 < 0$ and $z_2\le0$;

$({\rm{ii}}.2)$ $\Delta < 0$, $c_0 = 0$ and $z_2 < 0 < z_3$;

$({\rm{ii}}.3)$ $\Delta < 0$, $c_0 < 0$, $z_2 > 0$ and $F(z_2) F(z_3) > 0$;

$({\rm{ii}}.4)$ $\Delta < 0$, $c_0 = 0$, $z_1 > 0$ and $F(z_2) F(z_3) > 0$;

$({\rm{ii}}.5)$ $\Delta\geq 0$ and $c_0 < 0$;

$({\rm{ii}}.6)$ $\Delta\geq 0$, $c_0 = 0$ and $z_0 > 0$.

$({\rm{iii}})$ $F(z)$ has two simple positive zeros and no other positive zeros if and only if $({\bf H_2})$: one of the following conditions holds.

$({\rm{iii}}.1)$ $\Delta < 0$, $c_0 > 0$, $z_1\leq0 < z_3$ and $F(z_3) < 0$;

$({\rm{iii}}.2)$ $\Delta < 0$, $c_0 = 0$, $z_1\leq0 < z_2$ and $F(z_3) < 0$;

$({\rm{iii}}.3)$ $\Delta < 0$, $c_0 > 0$, $z_1 > 0$ and $F(z_1) F(z_2) F(z_3) < 0$;

$({\rm{iii}}.4)$ $\Delta\geq 0$, $c_0 > 0$, $z_0 > 0$ and $F(z_0) < 0$.

$({\rm{iv}})$ $F(z)$ has three simple positive zeros and no other positive zeros if and only if $({\bf H_3})$: one of the following conditions holds.

$({\rm{iv}}.1)$ $\Delta < 0$, $c_0 < 0$, $z_2 > 0$ and $F(z_2) F(z_3) < 0$;

$({\rm{iv}}.2)$ $\Delta < 0$, $c_0 = 0$, $z_1 > 0$ and $F(z_2) F(z_3) < 0$.

$({\rm{v}})$ $F(z)$ has exactly four simple positive zeros if and only if $({\bf H_4})$: $\Delta < 0$, $c_0 > 0$, $z_1 > 0$, $F(z_1) < 0$ and $F(z_2) F(z_3) < 0$ holds.

We now explore the dynamics of model (1.1) under the conditions (ⅰ)–(ⅲ) in Lemma 4.3. The last two cases $({\rm{iv}})$ and $({\rm{v}})$ in Lemma 4.3 are extremely complicated to analyze and we omit here. For case $({\rm{i}})$ in Lemma 4.3, we set

Theorem 4.2. Consider model (1.1) with $\tau_1 = \tau_2 = 0$. Assume that $R_{IM} > 1$, $q_1 > 0$, $q_2 > 0$, and $0\in I_0$, then the IAE $E_2$ is locally asymptotically stable for all $\tau_3\in [0, \sup\mathcal{I}_0)$, where $\mathcal{I}_0$ is the maximal connected subinterval of $I_0$ and $0\in\mathcal{I}_0$.

Next we focus on the case $({\rm{ii}})$ in Lemma 4.3. If $({\bf H_1})$ is satisfied, then $F(\omega, \tau_3)$ in (4.4) admits exactly one simple positive zero, denoted by $\overline{\omega}$. Let

For $\tau_3\in I_1$, we have $F'_\omega(\omega(\tau_3)) > 0$. It then follows from Lemma 4.2 that $\pm i \overline{\omega}(\tau_3)$ are a pair of simple purely imaginary eigenvalues of (4.1) if and only if $S_n(\tau_3) = 0$ has a positive root $\tau_3\in I_1$ for some $n\in \mathbb{N}$, and

If $I_1\neq\emptyset$, we denote $\bar\tau_{3}: = \sup I_1$. Then $F(\omega, \tau_3)$ has either two simple positive zeros or no positive zero when $\tau_3 = \bar{\tau}_3+\varepsilon$ with sufficiently small $\varepsilon > 0$. The former case will be studied later. In the latter case, we have $\lim\limits_{\tau_3\to \bar\tau_{3}} \overline{\omega}(\tau_3) = 0$. This, together with (4.7), implies that $\lim\limits_{\tau_3\to \bar\tau_{3}^{-}} \theta(\tau_3) = \pi$ and $\lim\limits_{\tau_3\to \bar\tau_{3}^{-}} S_n(\tau_3) = -(2n+1)\pi < 0$ for all $n\in \mathbb{N}$. One can easily observe that $S_{n+1}(\tau_3) < S_n(\tau_3)$ for all $n\in \mathbb{N}$. In view of Lemma 4.1, $S_n(0) < 0$ for all $\tau_3\in I_1$ and $n\in \mathbb{N}$ provided that $0\in I_1$ and $q_1, q_2 > 0$. Thus, the function $S_n(\tau_3)$ has at least two zeros in $I_1$ provided that $\sup\limits_{\tau_3\in I_1} S_n(\tau_3) > 0$ for $n\in \mathbb{N}$. To investigate the stability of $E_2$ and the existence of Hopf bifurcations, we assume that

{$({\bf A_1})$} $I_1 = [0, \bar\tau_3)$, $\lim\limits_{\tau_3\to \bar\tau_{3}} \overline{\omega}(\tau_3) = 0$, $\sup\limits_{\tau_3\in I_1} S_0(\tau_3) > 0$, and $S_n(\tau_3)$ has at most two zeros (counting multiplicity) for any $n\in\mathbb{N}$.

Assumption $({\bf A_1})$ implies that there exists a positive integer $K$ such that, for all integer $n\in[0, K-1]$, $S_n(\tau_3)$ has two simple zeros, denoted by $\tau_3^n < \tau_3^{2K-1-n}$. Since $S_n(\tau_3)$ is strictly decreasing in $n$, we have $0 < \tau_3^0 < \tau_3^1 < \cdots < \tau_3^{2K-1} < \bar{\tau}_3$. Obviously, we have $S'_n(\tau_3^n) > 0$ and $S'_n(\tau_3^{2K-1-n}) < 0$ for each integer $0\leq n\leq K-1$. It then follows from (4.12) that a pair of purely imaginary eigenvalues $\pm i \overline{\omega}(\tau_3)$ of (4.1) cross the imaginary axis from left to right (resp. right to left) at $\tau_3^n$ (resp. $\tau_3^{2K-1-n}$). Hence, system (1.1) with $\tau_1 = \tau_2 = 0$ undergoes a Hopf bifurcation at $E_2$ when $\tau_3 = \tau_3^i$ with integer $0\leq i\leq 2K-1$. Let $P_i$ be the period of periodic solution bifurcated at $\tau_3^i$. Applying the Hopf bifurcation theorem for delay differential equations ^[24], we have

for integer $0\le i\le K-1$. Thus, $P_0 > \tau_3^0$, $P_{2K-1} > \tau_3^{2K-1}$, ${\tau_3^n/(n+1)} < P_n\leq{\tau_3^n/n}$ and $\tau_3^{2K-n-1}/(n+1) < P_{2K-n-1}\le \tau_3^{2K-n-1}/n$ for integer $1\le n\le K-1$. To summarize, we have the following results.

Theorem 4.3. Consider model (1.1) with $\tau_1 = \tau_2 = 0$, denote $\bar\tau_{3}: = \sup I_1$. Assume that $R_{IM} > 1$, $q_1 > 0, q_2 > 0$ and $({\bf A_1})$ hold, then there exist exactly $2K$ local Hopf bifurcation values, namely, $\tau_3^0 < \tau_3^1 < \cdots < \tau_3^{2K-1}$ such that the model undergoes a Hopf bifurcation at the IAE $E_2$ when $\tau_3 = \tau_3^i$ for integer $0\leq i\leq 2K-1$. $E_2$ is locally asymptotically stable for $\tau_3\in [0, \tau_3^0)\cup (\tau_3^{2K-1}, \bar{\tau}_3)$ and unstable for $\tau_3\in (\tau_3^0, \tau_3^{2K-1})$. Moreover, for all integers $0\leq i\leq 2K-1$, when $\tau_3$ is sufficiently close to $\tau_3^i$, the period $P_i$ of periodic solution bifurcated at $\tau_3^i$ satisfies $P_0 > \tau_3^0$, $P_{2K-1} > \tau_3^{2K-1}$, ${\tau_3^n/(n+1)} < P_n\leq{\tau_3^n/n}$ and $\tau_3^{2K-n-1}/(n+1) < P_{2K-n-1}\le \tau_3^{2K-n-1}/n$ for integer $1\le n\le K-1$.

We now study the case $({\rm{iii}})$ in Lemma 4.3. If $({\bf H_2})$ holds, then $F(\omega, \tau_3)$ in (4.4) has exactly two simple positive zeros, denoted by $\tilde{\omega}_- < \tilde{\omega}_+$. Set

and $\tilde{\theta}_\pm$ is defined in (4.7) when $\omega = \tilde{\omega}_\pm$. From (4.8), we define

The relation $\tilde{\omega}_- < \tilde{\omega}_+$ implies that $F'_\omega(\tilde\omega_-(\tau_3)) < 0$ and $F'_\omega(\tilde\omega_+(\tau_3)) > 0$ for $\tau_3\in I_2$. If $S^+_n(\tau_3)$ (resp. $S^-_n(\tau_3)$) has a positive zero $\tau_{3+}$ (resp. $\tau_{3-}$) for some $n\in \mathbb{N}$, then (4.1) admits a pair of simple purely imaginary eigenvalues $\pm i \tilde{\omega}_+(\tau_{3+})$ (resp. $\pm i \tilde{\omega}_-(\tau_{3-})$), and

If $I_2\neq\emptyset$, we denote $\tilde\tau_{3}: = \sup I_2$. Then there exist two cases for the existence of positive zeros of $F(\omega, \tau_3)$ when $\tau_3 = \tilde{\tau}_3+\varepsilon$ with sufficiently small $\varepsilon > 0$. One reads that one more positive real root occurs besides $\tilde{\omega}_\pm$, which can be converted to case $({\rm{iv}})$ in Lemma 4.3. Another case follows that $\tilde{\omega}_\pm$ collide together at one zero with multiplicity two, which implies that $\lim\limits_{\tau_3\to \tilde\tau_{3}}\tilde{\omega}_+(\tau_{3}) = \lim\limits_{\tau_3\to \tilde\tau_{3}}\tilde{\omega}_-(\tau_{3})$, then $\lim\limits_{\tau_3\to \tilde\tau_{3}}S_n^+(\tau_{3}) = \lim\limits_{\tau_3\to \tilde\tau_{3}}S_n^-(\tau_{3})$. For simplicity, we assume that

$({\bf A_2})$ $I_2 = [0, \tilde\tau_{3})$; $\lim\limits_{\tau_3\to \tilde\tau_{3}}\tilde{\omega}_+(\tau_{3}) = \lim\limits_{\tau_3\to \tilde\tau_{3}}\tilde{\omega}_-(\tau_{3})$; $S_0^-(\tau_3) < S_0^+(\tau_3)$; $\sup\limits_{\tau_3\in I_2} S_0^+(\tau_3) > 0$, and $S_n^{\pm}(\tau_3)$ has at most two zeros (counting multiplicity) for any $n\in\mathbb{N}$.

Lemma 4.4. Assume that $R_{IM} > 1$, $q_1 > 0, q_2 > 0$ and $({\bf A_2})$ hold. Then for any $n\in \mathbb{N}$, $S_n^{\pm}(0) < 0$; $S_n^-(\tau_3) < S_n^+(\tau_3)$; $S_n^{\pm}(\tau_3)$ is strictly decreasing in $n$; $\lim\limits_{\tau_3\to \tilde\tau_{3}}S_n^+(\tau_{3}) = \lim\limits_{\tau_3\to \tilde\tau_{3}}S_n^-(\tau_{3})$, where $\tilde\tau_{3}: = \sup I_2$.

Denote $\widetilde{S}_n: = \lim\limits_{\tau_3\to \tilde\tau_{3}}S_n^+(\tau_{3})$. In the following, we assume that $R_{IM} > 1$, $q_1 > 0, q_2 > 0$ and $({\bf A_2})$ hold. Then Lemma 4.4 indicates that if $\widetilde{S}_n > 0$ for some $n\in \mathbb{N}$, then $S_n^{\pm}(\tau_3)$ have exactly one simple zero $\tau_{3\pm}^n$, $S'_n(\tau_{3\pm}^n) > 0$, and $\tau_{3+}^n < \tau_{3-}^n$ due to Lemma 4.1 and (4.15); if $\widetilde{S}_n < 0$ for some $n\in \mathbb{N}$, then each $S_n^{\pm}(\tau_3)$ has either nonzero or two simple zeros, denoted by $\tau_{3\pm}^{n1} < \tau_{3\pm}^{n2}$, $S'_n(\tau_{3\pm}^{n1}) > 0$ and $S'_n(\tau_{3\pm}^{n2}) < 0$. Thus, the total number of all simple zeros of $S_n^{\pm}(\tau_3)$ for all $n\in\mathbb{N}$ is even provided that $\widetilde{S}_m\neq0$ for any $m\in \mathbb{N}$. Then all $S_n^+(\tau_3)$ ($S_n^-(\tau_3)$) have $K_1$ ($K_2$) simple zeros for all integer $n\in \mathbb{N}$, and list them in an increasing order as $0 < \tau_{3+}^0 < \tau_{3+}^1 < \cdots < \tau_{3+}^{K_1-1} < \tilde{\tau}_3$ ($0 < \tau_{3-}^0 < \tau_{3-}^1 < \cdots < \tau_{3-}^{K_2-1} < \tilde{\tau}_3$). Clearly, $\tau_{3+}^0 < \tau_{3-}^0$, $K_1\ge K_2$ and $K_1+K_2$ is an even number. Now, we consider the set of all zeros of $S_n^{\pm}(\tau_3)$. If a value appears more than once in the set, then there are at least two pairs of purely imaginary roots and thus the condition of Hopf bifurcation is violated. For this reason, we only keep the values that appear exactly once in the set and rearrange them in an increasing order, denoted by

Based on the above analysis, using a similar method as that in [32,Theorem 4.9], we have the following Hopf bifurcation and multiple stability switches theorem.

Theorem 4.4. Consider model (1.1) with $\tau_1 = \tau_2 = 0$. Assume that $R_{IM} > 1$, $q_1 > 0, q_2 > 0$ and $({\bf A_2})$ hold, $\widetilde{S}_n\neq0$ for any $n\in \mathbb{N}$, $\tau_3^j$ are defined in (4.16). Then there exist exactly $2K$ Hopf bifurcation values, namely, $0 < \tau_3^0 < \tau_3^1 < \cdots < \tau_3^{2K-1} < \tilde{\tau}_3$ such that the model undergoes a Hopf bifurcation at $E_2$ when $\tau_3 = \tau_3^j$ for $0\le j\le 2K-1$. $E_2$ is locally asymptotically stable for $\tau_3\in [0, \tau_{3}^0)\cup (\tau_3^{2K-1}, \tilde{\tau}_3)$, and there are three cases on the stability switches of $E_2$ for $\tau_3\in (\tau_{3}^0, \tau_3^{2K-1})$:

$({\rm{i}})$ $E_2$ is unstable for $\tau_3\in (\tau_{3}^0, \tau_3^{2K-1})$;

$({\rm{ii}})$ there exist $l+2$ stability switches for an even integer $2\le l\le K$: $\tau_3^0, \tau_3^1, \cdots, \tau_3^l$ and $\tau_3^{2K-1}$, namely, $E_2$ is locally asymptotically stable for $\tau_3\in (\tau_{3}^1, \tau_{3}^2)\cup\cdots\cup (\tau_{3}^{l-1}, \tau_{3}^l)$, and unstable for $\tau_3\in(\tau_{3}^0, \tau_{3}^1) \cup(\tau_{3}^2, \tau_{3}^3)\cup\cdots\cup(\tau_3^{l}, \tau_3^{2K-1})$;

$({\rm{iii}})$ all Hopf bifurcation values $\tau_3^j \ (0\le j\le 2K-1)$ are stability switches, namely, $E_2$ is locally asymptotically stable for $\tau_3\in (\tau_{3}^1, \tau_{3}^2)\cup\cdots\cup(\tau_3^{2K-3}, \tau_3^{2K-2})$, and unstable for $\tau_3\in(\tau_{3}^0, \tau_{3}^1) \cup(\tau_{3}^2, \tau_{3}^3)\cup\cdots\cup(\tau_3^{2K-2}, \tau_3^{2K-1})$.

4.2.   Global Hopf bifurcation analyses

Theorems 4.3 and 4.4 give the sufficient conditions on the existence of periodic solutions bifurcated at the IAE $E_2$ when the CTLs recruitment delay $\tau_3$ is near the local Hopf bifurcation values, denoted by $\tau_3^*$. In this subsection, we show the global existence of the bifurcating periodic solutions by using the global Hopf bifurcation theorem for delay differential equations ^[31]. Throughout this subsection, we assume that $R_{IM} > 1$, $q_1 > 0, q_2 > 0$ and $({\bf A_1})$, that is, the conditions in Theorem 4.1 hold to analyze the global Hopf branches.

Let $u(t) = (T(\tau_3 t), I(\tau_3 t), V(\tau_3 t), Z(\tau_3 t))^{T}$. Model (1.1) with $\tau_1 = \tau_2 = 0$ can be rewritten as the following functional differential equation

where $u_t(\theta) = u(t+\theta)$ for $\theta\in[-1, 0]$, and $u_t\in \mathcal{X}: = \mathbb{R}_+\times \mathcal{C}_0\times \mathbb{R}_+\times \mathcal{C}_0$ with $\mathcal{C}_0: = C([-1, 0], \mathbb{R})$. $I_1$ is defined in (4.11), and

for $\phi = (\phi_1, \phi_2, \phi_3, \phi_4)^{T}\in \mathcal{X}$. Restricting $F$ to the subspace of $\mathcal{X}$ consisting of all constant functions with $\mathbb{R}_+^4$, we obtain a map

for $u = (u_1, u_2, u_3, u_4)^{T}$. Clearly, $\widehat{F}$ is twice continuously differentiable, thus assumption (A1) in ^[31] is satisfied. Denote the set of stationary solutions of system (4.17) by

It follows from Theorem 2.1 that ${\bf N}(F) = \{(E_0, \tau_3, P), (E_1, \tau_3, P), (E_2, \tau_3, P)\}$ if $R_{IM} > 1$. For any stationary solution $(\hat{u}, \tau_3, P)$, the characteristic equation of (4.17) is

where $DF(\hat{u}, \tau_3, P)$ is the Jacobian matrix of $F$ at the stationary solutions. Clearly, $0$ is not the eigenvalue of any stationary solution of (4.17) when $R_{IM} > 1$, which implies that statement (A2) in ^[31] holds. By the expression of $F$, it can be checked easily that (A3) in ^[31] is satisfied.

A stationary solution $(\hat{u}, \tau_3, P)$ of (4.17) is called a center if $\det(\Delta_{\hat{u}}(i m \frac{2 \pi}{\widehat{P}})) = 0$ for some positive integer $m$. If there exist no other centers in some neighborhoods of $(\hat{u}, \tau_3, P)$ and only finite purely imaginary eigenvalues of the form $i m \frac{2 \pi}{\widehat{P}}$, then the center is isolated. By Theorem 4.3, the stationary solution $(E_2, \tau_3^j, 2 \pi/(\omega_j\tau_3^j))$ is an isolated center of (4.17) for each integer $0\leq j\leq 2K-1$, where $\omega_j = \overline{\omega}(\tau_3^j)$ is the unique positive zero of $F(\omega, \tau_3)$. In addition, only one purely imaginary characteristic value of the form $i m \frac{2 \pi}{\widehat{P}}$ exists with $m = 1$ and $P = 2 \pi/(\omega_j \tau_3^j)$. Thus the set of all positive integers $m$ is the singleton $\{1\}$. By the transversality condition (4.12), the crossing number $\gamma(E_2, \tau_3^j, 2 \pi/(\omega_j\tau_3^j))$ follows

which implies that condition (A4) in ^[31] is satisfied. We now define a closed subset $\Sigma(F)$ of $\mathcal{X}\times I_1\times \mathbb{R}_+$ by

We denote by ${\bf C}(E_2, \tau_3^j, 2 \pi/(\omega_j\tau_3^j))$ the connected component of $(E_2, \tau_3^j, 2 \pi/(\omega_j\tau_3^j))$ in $\Sigma(F)$ for each integer $0\leq j\leq 2K-1$. It follows from Theorem 4.3 that ${\bf C}(E_2, \tau_3^j, 2 \pi/(\omega_j\tau_3^j))$ is a nonempty subset of $\Sigma(F)$. We now show the boundedness of the periodic solutions of system $(4.17)$.

Lemma 4.5. Assume that $R_{IM} > 1$, then all nonnegative periodic solutions of (4.17) are uniformly bounded, namely, $\epsilon\leq T(t), I(t), V(t), Z(t)\leq M$ for all $t\in \mathbb{R}_+$, where $M = \max\{\overline{T}, \overline{I}, \overline{V}, \overline{Z}\}$, and $\epsilon > 0$ is defined in Lemma 3.2.

Proof. We claim that $T(t)\leq M$ for all $t\in \mathbb{R}_+$. Otherwise, if there exists $t_1\in \mathbb{R}_+$ such that $T(t_1) > M$, then $\lim\limits_{n\rightarrow \infty}T(t_1+nP) = T(t_1) > M$, where $P$ is the period of the periodic solution $(T(t), I(t), V(t), Z(t))$. This contradicts the fact that $\liminf\limits_{t\to \infty} T(t)\leq \overline{T}\leq M$ in Lemma 2.1. Hence, $M$ is a uniform upper bound of $T(t)$. Similarly, we can prove that $I(t), V(t), Z(t)$ have a uniform upper bound $M$ from Lemma 2.1, and $T(t), I(t), V(t), Z(t)$ have a uniform lower bound $\epsilon$ from Lemma 3.2. This completes the proof.

Lemma 4.5 shows that the projection of ${\bf C}(E_2, \tau_3^j, 2 \pi/(\omega_j\tau_3^j))$ onto $\mathcal{X}$ is bounded for any integer $0\leq j\leq 2K-1$. Our next lemma excludes the existence of periodic solutions of (4.17) of period 1.

Lemma 4.6. Assume that $R_{IM} > 1$ and that $T_m(r-d) < r T_2$ holds, then system (4.17) has no nontrival periodic solution of period 1.

Proof. Assume to the contrary that $u(t)$ is a nontrivial periodic solution of (4.17) with period 1, that is $u(t-1) = u(t)$, which is equivalent to

Then $(T(t), I(t), V(t), Z(t))$ satisfies the following system

Theorem 3.3 indicates that the above system has no nontrivial periodic solutions. This is a contradiction and the proof is complete.

Note that system (4.17) has no periodic solutions of period 1 in Lemma 4.6, thus no periodic solutions of period $1/n$ or $1/(n+1)$ for any positive integer $n$. It follows from Theorem 4.3 and Lemma 4.6 that the period $P_j$ of a periodic solution on the connected component ${\bf C}(E_2, \tau_3^n, 2 \pi/(\omega_n \tau_3^n))$ satisfies $P_0 > 1$, $P_{2K-1} > 1$ and

where $K$ is defined in Theorem 4.3. As we shall see later in the numerical simulation, it seems hard to find an upper bound for the periods $P_0$ and $P_{2K-1}$. In what follows, we will restrict our investigation to the set

and assume that $J\neq\emptyset$. Especially, we will discuss the global continuation of periodic solutions bifurcated from the point $(E_2, \tau_3^j)$ for $1\le j\le 2K-2$ as the bifurcation parameter $\tau_3$ varies. Thus, the projection of ${\bf C}(E_2, \tau_3^j, 2 \pi/(\omega_j\tau_3^j))$ onto the $P$-space is bounded for any integer $1\le j\le 2K-2$. Clearly, $\tau_3\in I_1$ is a bounded interval, and Lemma 4.5 implies that the projection of ${\bf C}(E_2, \tau_3^j, 2 \pi/(\omega_j\tau_3^j))$ onto $\mathcal{X}$ is bounded for all integer $0\leq j\leq 2K-1$. Therefore, ${\bf C}(E_2, \tau_3^j, 2 \pi/(\omega_j\tau_3^j))$ is bounded in $\mathbb{R}_+^4\times I_1\times \mathbb{R}_+$.

The periodic solutions are all bounded away from zero by Lemma 4.5. Thus there is no need to consider the boundary equilibrium. We define ${\bf N}_1(F) = \{(E_2, \tau_3, P), \ (\tau_3, P)\in I_1\times \mathbb{R}_+\}$. By using the global Hopf bifurcation theorem [31,Theorem 3.3], we have $\mathcal{E}: = {\bf C}(E_2, \tau_3^j, 2 \pi/(\omega_j\tau_3^j))\bigcap {\bf N}_1(F)$ is finite and

Based on the above discussion, using a similar method as that in the proof [32,Theorem 5.4], we can now present the result describing the global continuation of Hopf bifurcation as follows.

Theorem 4.5. Assume that $R_{IM} > 1$, $T_m(r-d) < r T_2$, and $J: = \{\tau_3^j: 1\le j\le 2K-2\}\neq\emptyset$, where $\tau_3^j$ and $K$ are defined in Theorem 4.3. Then for (4.17), we have the following results.

$({\rm{i}})$ All global Hopf branches ${\bf C}(E_2, \tau_3^j, 2 \pi/(\omega_j\tau_3^j))$ are bounded for any $1\le j\le 2K-2$.

$({\rm{ii}})$ Two global Hopf branches ${\bf C}(E_2, \tau_3^n, 2 \pi/(\omega_n \tau_3^n))$ and ${\bf C}(E_2, \tau_3^{2K-n-1}, 2 \pi/(\omega_{2K-n-1} \tau_3^{2K-n-1}))$ coincide with each other and thus connect a pair of Hopf bifurcation values $\tau_3^n$ and $\tau_3^{2K-n-1}$ for any $1\le n\le K-1$.

$({\rm{iii}})$ For any $1\le n\le K-1$, there exists at least one periodic solution for each $\tau_3\in (\tau_3^n, \tau_3^{2K-n-1})$ with period in $(1/(n+1), 1/n)$.

$({\rm{iv}})$ For any $1\le i, j\le 2K-2$, ${\bf C}(E_2, \tau_3^j, 2 \pi/(\omega_j\tau_3^j))\cap {\bf C}(E_2, \tau_3^i, 2 \pi/(\omega_i\tau_3^i)) = \emptyset$ if $j\neq 2K-i-1$.

In Theorem 4.5, we have assumed that $T_m(r-d) < r T_2$. This condition is essential in finding a uniform upper bound for the periods of periodic solutions in a global Hopf branch. If $T_m(r-d)\ge r T_2$, we will not be able to find an upper bound for the periods of periodic solutions. However, as we see later in the numerical exploration, we still observe and thus conjecture that the global Hopf branches are bounded.

By using a similar method as that in the proof of Theorem 4.5, we can also analyze the global continuation of Hopf bifurcation when the conditions in Theorem 4.4 hold. We will numerically show the complex dynamics of system (1.1) for the case (ⅰ)–(ⅲ) of Theorem 4.4.

5.   Numerical exploration

In this section, we apply numerical explorations to illustrate the dynamical behavior of system (1.1). We choose $\tau_3$ as the bifurcation parameter and set the parameter values as follows.

It is readily seen that the conditions in Theorem 4.3 are satisfied. Direct calculation gives $\bar{\tau}_3\approx27.059 < \tau_{3, max} = 39.283$, and $R_{IM} > 1$ if and only if $0\le\tau_3 < \tau_{3, max}$. There exist exactly $8 \ (\text{with} \ K = 4)$ local Hopf bifurcation values

see also in Figure 1. Theorem 4.3 implies that $E_2$ is locally asymptotically stable for $\tau_3\in[0, \tau_3^0)\cup (\tau_3^7, \bar{\tau}_3)$, and unstable for $\tau_3\in (\tau_3^0, \tau_3^7)$.

In Figure 2(a), by using the Matlab package DDE-BIFTOOL, we plot the global Hopf branches ${\bf C}(E_2, \tau_3^j, 2 \pi/(\omega_j\tau_3^j))$ bifurcated near each bifurcation point $\tau_3^j$ for integer $0\le j\le 7$. It is observed that the two branches ${\bf C}(E_2, \tau_3^k, 2 \pi/(\omega_k\tau_3^k))$ and ${\bf C}(E_2, \tau_3^{7-k}, 2 \pi/(\omega_{7-k}\tau_3^{7-k}))$ are bounded and connected for any integer $k = 1, 2, 3$, which coincides with Theorem 4.5. A periodic solution is unstable if and only if its principal Floquet multiplier is larger than one ^[24]. We numerically calculate the associated principal Floquet multipliers for periodic solutions on each global Hopf branch, see Figure 2(b). To further investigate the dynamics when $\tau_3 > \bar{\tau}_3$, we depict the bifurcation diagram in Figure 2(c), (d). We observe that chaotic behavior occurs approximately for $\tau_3\in(14.4, 21.2)\cup(24.6, 25)$, and the stable equilibrium $E_2$ and a stable periodic solution coexist for $\tau_3\in(29.6, 31.8)$.

We now numerically explore the dynamical behavior of system (1.1) for Theorem 4.4, see Figure 3. To explain the dynamics specifically, we list the local Hopf bifurcation values following an increasing order and derive the corresponding stability results of $E_2$ as follows.

$({\rm{i}})$ In Figure 3(a), all Hopf bifurcation values are $\tau_{3+}^0 < \tau_{3+}^1 < \tau_{3-}^0 < \tau_{3+}^2 < \tau_{3-}^1 < \tau_{3+}^3 < \tau_{3+}^4 < \tau_{3+}^5$, $E_2$ is stable for $\tau_3\in [0, \tau_{3+}^0)\cup(\tau_{3+}^5, \tilde{\tau}_3)$ and unstable for $\tau_3\in (\tau_{3+}^0, \tau_{3+}^5)$. There are two stability switches $\tau_{3+}^0$ and $\tau_{3+}^5$.

$({\rm{ii}})$ In Figure 3(c), all Hopf bifurcation values are $\tau_{3+}^0 < \tau_{3-}^0 < \tau_{3+}^1 < \tau_{3-}^1 < \tau_{3+}^2 < \tau_{3-}^2 < \tau_{3+}^3 < \tau_{3-}^3$, $E_2$ is stable for $\tau_3\in [0, \tau_{3+}^0)\cup(\tau_{3-}^0, \tau_{3+}^1)\cup(\tau_{3-}^1, \tau_{3+}^2)\cup (\tau_{3-}^2, \tau_{3+}^3)\cup(\tau_{3-}^3, \tilde{\tau}_3)$ and unstable for $\tau_3\in (\tau_{3+}^0, \tau_{3-}^0)\cup(\tau_{3+}^1, \tau_{3-}^1)\cup(\tau_{3+}^2, \tau_{3-}^2) \cup(\tau_{3+}^3, \tau_{3-}^3)$. Here, all Hopf bifurcation values are stability switches.

$({\rm{iii}})$ In Figure 3(e), all Hopf bifurcation values are $\tau_{3+}^0 < \tau_{3-}^0 < \tau_{3+}^1 < \tau_{3+}^2 < \tau_{3-}^1 < \tau_{3+}^3 < \tau_{3+}^4 < \tau_{3-}^2$, $E_2$ is stable for $\tau_3\in [0, \tau_{3+}^0)\cup(\tau_{3-}^0, \tau_{3+}^1)\cup(\tau_{3-}^2, \tilde{\tau}_3)$ and unstable for $\tau_3\in (\tau_{3+}^0, \tau_{3-}^0)\cup(\tau_{3+}^1, \tau_{3-}^2)$. There are four stability switches $\tau_{3+}^0, \tau_{3-}^0, \tau_{3+}^1$ and $\tau_{3-}^2$.

The above results are following Theorem 4.4. An interesting behavior where two stable periodic solutions coexist for (1.1), depicted in Figure 4. To explore the effect of $\tau_3$ on the characteristics of periodic solutions, we choose the Hopf branches in Figure 3(d) to depict the phase orbits of $I(t), V(t)$ and $Z(t)$, see Figure 5.

To understand joint effects of the intrinsic mitosis rate of the uninfected target cells $r$ and the CTLs recruitment delay $\tau_3$ on viral dynamics in vivo, we numerically carry out two-parameter bifurcation analysis of model (1.1) with bifurcation parameters $r$ and $\tau_3$. According to the biological significance ^[30], we restrict $r\in[0, 3]$ and the other parameters are selected the same as those in Figure 3(e) except $s_3 = 0.6$. As shown in Figure 6, we observe that both $r$ and $\tau_3$ can destabilize the IAE $E_2$ and cause Hopf bifurcations. However, they do behave differently. $\tau_3$ can cause Hopf bifurcations only when $r$ is sufficiently large, and $r$ can cause Hopf bifurcations only when $\tau_3$ is neither too large nor too small.

It follows from Theorem 3.3 that the intracellular delay of viral infection $\tau_1$ and viral production delay $\tau_2$ will not lead to Hopf bifurcations or periodic oscillations when $r\in [0, T_m d/(T_m-T_2))$. To measure the oscillatory phenomenon of model (1.1), we choose the parameters in Figure 6 with $r = 2$ such that $r > T_m d/(T_m-T_2)$ to examine the impact of these two delays on in vivo viral infections. As shown in Figure 7, both $\tau_1$ and $\tau_2$ can destabilize the stability of $E_2$ and generate sustained oscillations as well as chaotic solutions, $E_2$ does not exist for all large $\tau_1$ regardless of the initial state. Then all solutions of (1.1) converge to the IFE $E_0$ when $\tau_3 > 17$ in Figure 7(a), (b), which implies that prolonging the intracellular delay in the process of viral infection is effective to suppress viral transmission.

6.   Summary and discussion

Because there exist time delays for viral infection, viral production, and CTLs recruitment, we proposed a delayed viral infection model with mitosis in the target-cell dynamics, two infection modes (virus-to-cell transmission and cell-to-cell transmission), and immune response. We have investigated whether these time delays and a mitotic term in the target-cell dynamics can independently lead to periodic oscillations, or more specifically, whether intracellular delays can lead to periodic oscillations without mitosis in the target-cell dynamics.

It is shown that this model admits three types of equilibria: the infection-free equilibrium (IFE), the immune-inactivated equilibrium (IIE), and the immune-activated equilibrium (IAE). The dynamics of our model are shown to be determined by two critical values: the basic reproduction number for infection $R_0$, and the basic reproduction number for immune response $R_{IM}$. More precisely, we have proved the following: (ⅰ) the IFE $E_0$ is globally asymptotically stable if $R_0\le1$, which means that the viral particles are cleared; (ⅱ) the IIE $E_1$ is globally asymptotically stable if $R_{IM}\le1 < R_0$ provided that $T_m(r-d) < r T_1$ holds, that is, the infection becomes chronic with no sustained immune responses; (ⅲ) the IAE $E_2$ is globally asymptotically stable for model (1.1) in the absent of the CTLs recruitment delay if $R_{IM} > 1$ provided that $T_m(r-d) < r T_2$ is satisfied, which indicates that both virus infection and CTL response will be permanent. We have also shown that both the virus-to-cell and cell-to-cell infection modes contribute to the basic reproduction number $R_0$, which implies that ignoring the cell-to-cell transmission will produce an underestimation of $R_0$.

In the case $R_{IM} > 1$, we conduct bifurcation analysis using the CTLs recruitment delay $\tau_3$ as the bifurcation parameter. We obtain the sufficient and necessary conditions for the local stability of the IAE $E_2$. If the condition is violated, we give the existence of local Hopf bifurcation and stability switch values for the CTLs recruitment delay. These results imply that the CTLs recruitment delay $\tau_3$ can stabilize or destabilize $E_2$ through Hopf bifurcation. We further apply the global Hopf bifurcation theorem to demonstrate the global continuity of each Hopf branch. It is proved that each global Hopf branch except the first and the last Hopf branches is bounded and it connects two Hopf bifurcation points. Moreover, multiple periodic solutions may coexist in the overlapped intervals of the Hopf branches. By calculating the principal Floquet multipliers, we obtain the stability of periodic orbits on the Hopf bifurcation branches and find some intervals on which multiple stable periodic solutions coexist.

From the numerical simulation, the CTLs recruitment delay $\tau_3$ can induce multiple stability switches, the coexistence of multiple stable limit cycles, and chaotic solutions, among other rich dynamical behaviors. By using DDE-BIFTOOL to plot the global Hopf bifurcation diagram, we also observe that the first and the last Hopf branches may not connect. A detailed theoretical study of this interesting phenomenon should enrich the global Hopf bifurcation theory, thus we leave it as a future work. We further numerically carry out a two-parameter bifurcation analysis. Our results show that while both $r$ and $\tau_3$ have a strong impact on the viral dynamics, the ways they exert their impact are quite different.

Acknowledgments

H. Shu is partially supported by the National Natural Science Foundation of China (No. 11971285), and the Fundamental Research Funds for the Central Universities (No. GK202201002). P. Jiang is partially supported by the National Natural Science Foundation of China (No. 72274119).

Conflict of interest

The authors declare there is no conflict of interest.