Effect of seasonality on the dynamics of an imitation-based vaccination model with public health intervention

1.   Introduction

According to clinical data, drug treatment against persistent human infections sometimes fails to consistently eradicate the infection from the host ^[1,2,3], such as human immunodeficiency virus (HIV), hepatitis B virus (HBV), and hepatitis C virus (HCV). For all these viruses, drug treatment is not effective, and for HIV, lifelong therapy is commonly required to control viral replication ^[4,5]. The main reason is that these viruses are able to suppress and impair immune responses, resulting in persistent infection ^[6,7]. An alternative strategy is to use drug therapies to boost virus-specific immune response and then induce sustained viral suppression in the absence of drugs. Antiviral therapy aimed at boosting specific immune responses has attracted more and more attentions recently. In 2003, Komarova et al. ^[6] used mathematical models to study immune response dynamics during therapy in the context of immunosuppressive infections. They showed that a single phase of antiviral drug therapy is also possible to establish sustained immunity. They assumed that the degree of immune expansion depends on virus load and the response inhibits virus growth. This assumption is natural for any branch of the adaptive immune system, such as CD4 T cells, CD8 T cells, or antibodies. Denote by $y(t)$ and $z(t)$, respectively, the virus population and immune cell population at time $t$. Komarova et al. ^[6] proposed the following system of ordinary differential equations:

The virus population grows at a rate described by the logistic function. The viral replication rate is a linear decreasing function of viral loads with a maximum value of $r$ and it vanishes at a viral load $K$. The treatment is modeled as a reduction in $r$. Immune cells kill virus at a rate $pyz$. $b$ and $\mu$ are the natural decay rates of viral particles and immune cells, respectively. Immune expansion is modeled by the growth function ${cyz\over 1+dy}$. When the virus load is low, the level of immune response is simply proportional to both the viral load and the immune cells. The immune response saturates when the virus load is sufficiently high. Immune cells are inhibited by the virus at a rate $qyz$. Komarova et al. presented a simple relationship between the timing of therapy and efficacy of the drugs. In the presence of strong viral suppression, they showed that therapy should be stopped relatively early because a longer duration of treatment may lead to a failure. On the other hand, in the presence of weaker viral suppression, stopping treatment too early is detrimental, and therapy has to be continued beyond a time threshold. Essentially, model (1.1) has two stable equilibria: a virus dominant equilibrium with no sustained immunity and an immune control equilibrium with sustained immunity. This bistability allows a solution from the basin of the attraction of the virus dominant equilibrium to be lifted to that of the immune control equilibrium via a single phase of therapy.

Model (1.1) assumes that the immune response is instantaneous. Nevertheless, it has been established that the immune response process involves a sequence of events such as antigenic activation, selection, and proliferation of the immune cells ^[8]. Moreover, oscillatory viral loads and immune cells were observed from clinical data ^[9], while the oscillatory behavior is not exhibited in model (1.1). A natural question is what is the cause of the oscillations. In 2014, Shu et al. ^[10] incorporated the time lag during the immune response process into model (1.1). By studying the dynamics between an immunosuppressive infection and antiviral immune response, they demonstrated that the time lag is the main factor causing oscillations. Here, $\tau$ is denoted as the time lag for the immune system to trigger a sequence of events. The activation rate of immune cells at time $t$ depends on the virus load and the number of immune cells at time $t-\tau$. The above assumptions lead to the following system:

It was shown that the delayed antiviral immune response may induce sustained periodic oscillations, transient oscillations and even sustained aperiodic oscillations (chaos), and the model can admit two stable equilibria and can also allow a stable equilibrium to coexist with a stable periodic solution.

In the aforementioned work, the mortality of the population during the time lag has been ignored. Consideration of the survival probability during the time lag requires an additional delay-dependent multiplier in the nonlinear delayed feedback term. The death of the free viral particles depends on the mature immune cells instead of all the newly generated immune cells. We denote $z(t)$ as the population of mature immune cells at time $t$, and consider the following delay differential equation with a delay-dependent coefficient:

Here, $f(y)$ describes the virus population growth function; $\delta > 0$ is the death rate of the immature immune cells during the time lag; $e^{-\delta\tau}$ represents the survival probability of the immature immune cells surviving $\tau$ time units before becoming mature; $h(y)z$ denotes the immune cell growth function; the other parameters have the same meaning as in model (1.2). In addition to more general virus and immune cell growth functions, model (1.3) differs from model (1.2) in the sense that it introduces an additional delay-dependent survival probability $e^{-\delta\tau}$. As we shall see later, this quantity will make a significant difference in the bifurcation analysis. For instance, model (1.3) does not possess any positive equilibrium for sufficiently large $\tau$ and its global Hopf branches always have bounded $\tau$ components. However, for model (1.2), the existence condition of positive equilibrium is independent of $\tau$ and the global Hopf branches always have unbounded $\tau$ components. Furthermore, there exists only one stability switch of positive equilibrium for model (1.2), but the stability of positive equilibrium for our model (1.3) may switch multiple times.

Actually, the model (1.3) can be derived from a stage structured population model for $u(t, a)$ as below

with the stage-specific decay rate of immune cells, $\mu(a)$, inhibition rate of immune cells, $q(a)$,

where $u(t, a)$ is the population of immune cells at age $a$ and time $t$. Note that the population of mature immune cells is $z(t) = \int_{\tau}^{\infty}u(t, a) da$. We integrate along characteristic lines to find

Note that $u(t, \infty) = 0$ and $u(t-\tau, 0) = h(y(t-\tau))z(t-\tau)$. The equation for $z(t)$ follows.

Throughout this paper, we assume that the virus population growth function $f(y)$ and the immune expansion rate $h(y)$ satisfy the following conditions:

(H1) $f(y)$ is continuously differentiable, $f(0) = 0$ and $f'(0) > b$; there exists $\bar{y} > 0$ such that $f(\bar{y}) = 0$, and $(y-\bar{y})f(y) < 0$ for positive $y\neq\bar{y}$; $f''(y) < 0$ for $y\in[0, \bar{y}]$.

(H2) $h(y)$ is continuously differentiable, $h(0) = 0$ and $h'(0) > q$; $h'(y) > 0$ and $h''(y) < 0$ for all $y\geq0$.

The assumptions (H1)-(H2) are biological conditions for infection dynamics. For instance, $f(0) = 0$ means no new infection in the absence of virus; $f'(0) > b$ indicates that the intrinsic growth rate is greater than the decay rate. Similarly, $h(0) = 0$ implies that immune cells cannot reproduce without the presence of virus; $h'(0) > q$ guarantees that the initial expansion rate is greater than the inhibition rate.

The rest of this paper is organized as follows. We first present some preliminary results concerning the well-posedness of model (1.3), and describe global dynamics of the ordinary differential equation system (2.4) in Section 2. Local and global stability analysis of equilibria as well as local and global Hopf bifurcation analysis are given in Section 3. In Section 4, numerical simulations based on the bifurcation analysis are reported and discussed. Finally, we give a brief summary in Section 5.

2.   Preliminaries

To establish the well-posedness of model (1.3), we choose the phase space $\mathcal{C}\times\mathcal{C}$, where $\mathcal{C}$ is the Banach space of continuous functions on $[-\tau, 0]$ defined by $\mathcal{C}: = \{\phi:[-\tau, 0]\rightarrow \mathbb{R} \ \text{is continuous} \}$, and the norm is defined by $\|\phi\| = \sup\limits_{-\tau\leq\theta\leq0}\phi(\theta)$. The nonnegative cone of $\mathcal{C}$ is denoted by $\mathcal{C^{+}}$. 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 (1.3) is given as

The existence and uniqueness of the solution of model (1.3) follows from the standard theory of functional differential equations ^[11]. Using the same manner as in proving [12,Proposition 2.1], we can easily show that the solution of (1.3) with initial conditions in $X$ is nonnegative, which implies that $X$ is positively invariant under the solution map of (1.3).

Next, we shall prove that all solutions of (1.3) are ultimately bounded. It follows from the first equation in (1.3) that $y'(t)\leq f(y)$, which yields $\limsup\limits_{t\rightarrow\infty}y(t)\leq\bar y$. By the nonnegativity of the solution of (1.3) and $(\textbf{H}_{1})$-$(\textbf{H}_{2})$, we obtain

where $M_f = \max\limits_{y\geq0}f(y)$ and $\gamma = \min\{b, \mu\}$. Thus,

which implies that $y(t)$ and $z(t)$ are ultimately bounded in $X$. To summarize, we obtain the following proposition.

Proposition 2.1. The solutions of model (1.3) with the initial conditions in $X$ are nonnegative, and the region

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

Clearly, model (1.3) always has a trivial equilibrium $E_0 = (0, 0)$. Since $f'(0) > b$, there also exists a boundary equilibrium $E_b = (\tilde{y}, 0)$, where $\tilde{y}$ is the unique positive root of $f(y) = by$. We call $E_b$ the virus dominant equilibrium (VDE). Assume that $E^* = (y^*, z^*)$ is a positive equilibrium with $y^* > 0$ and $z^* > 0$, then it must satisfy the following equilibrium equations

Obviously, $z^* > 0$ if and only if $y^* < \tilde{y}$. Hence, $E^*$ exists if and only if $y^*$ is a positive root of

with $y^* < \tilde{y}$. It is easily to calculate that $\text{sign}\left({h(y)\over y+\mu/q}\right)' = \text{sign} \ g(y)$, where $g(y) = h'(y)(y+{\mu\over q})-h(y)$. Note that $g(0) > 0, g(\infty) < 0$ and $g'(y) = h''(y)(y+\mu/q) < 0$. There exists a unique $y_c > 0$ such that $g(y_c) = 0$, and $(y-y_c)g(y) < 0$ for $y\neq y_c$. Therefore, ${h(y)\over y+\mu/q}$ is strictly increasing on $[0, y_c)$, and strictly decreasing on $(y_c, \infty)$. We denote

If $s < qe^{\delta\tau}$, (2.2) has no positive root. If $s\ge qe^{\delta\tau}$, (2.2) has exactly two positive roots (counting multiplicity) $y_1^*$ and $y_2^*$ such that $y_1^*\le y_c\le y_2^*$. These two roots coincide (i.e., $y_1^* = y_2^* = y_c$) if and only if $s = qe^{\delta\tau}$. Throughout this paper, we make the following assumption to ensure the existence of positive roots for (2.2).

(H₃) $s > q$ and $\tau\le\tau_{max}$, where $\tau_{max} = \frac{1}{\delta}\ln\frac s q$.

If ($\textbf{H}_{3}$) is violated, there does not exist any positive equilibrium for our model and the boundary equilibrium is locally asymptotically stable. Summarizing the above analysis, we obtain the following result.

Proposition 2.2. Assume that (H₁)-(H₃) are satisfied.

(i) If $y_1^*\geq \tilde{y}$ holds, then there are two equilibria: $E_{0} = (0, 0)$ and $E_{b} = (\tilde{y}, 0)$.

(ii) If $y_1^* < \tilde{y}\leq y_2^* \ \mathit{\text{or}} \ \ y_1^* = y_2^* < \tilde{y}$ holds, then there are three equilibria: $E_{0},$ $E_{b}$ and $E_1^* = (y_1^*, z_1^*)$.

(iii) If $y_1^* < y_2^* < \tilde{y}$ holds, then there are four equilibria: $E_{0},$ $E_b,$ $E_1^*$ and $E_2^* = (y_2^*, z_2^*).$

When $\tau = 0$, model (1.3) reduces to an ODE system which generalizes (1.1)

We next provide a complete description about its global dynamics.

Lemma 2.3. System (2.4) has no closed orbits.

Proof. Using an argument similar to Proposition 2.1, we can show that

is positively invariant and absorbing in $\mathbb{R}_+^2$. Thus, it suffices to show that System (2.4) has no closed orbits in $\Gamma_0$.

Let $(P(y, z), Q(y, z))$ denote the vector field of (2.4), then we have

Note that ($\textbf{H}_{1}$) implied that $(\frac{f(y)}{y})' < 0$ for all $y\leq \bar y$. Thus, $\frac{\partial(BP)}{\partial y}+\frac{\partial(BQ)}{\partial z} < 0$ in $\Gamma_0$. The nonexistence of closed orbits for (2.4) follows from the classical Dulac-Bendixson criterion. This ends the proof.

We linearize system (2.4) at $E_0 = (0, 0)$ and calculate that the two eigenvalues are $f'(0)-b > 0$ and $-\mu < 0$. Thus, $E_0$ is a saddle point. For the linearized system at $E_b = (\tilde{y}, 0)$, the two eigenvalues are

From the proof of Proposition 2.2, we know $h(\tilde{y})-q\tilde{y}-\mu < 0$ if $\tilde{y} < y_1^*$ or $\tilde{y} > y_2^*$, and $h(\tilde{y})-q\tilde{y}-\mu > 0$ provided that $y_1^* < \tilde{y} < y_2^*$. This shows that $E_b$ is stable if $y_1^* > \tilde{y}$ or $y_2^* < \tilde{y}$, and $E_b$ is a saddle point if $y_1^* < \tilde{y} < y_2^*$. For the linearized system at $E_1^* = (y_1^*, z_1^*)$, the associated characteristic equation is

It follows from (2.1) and ($\textbf{H}_{1}$) that $-f'(y_1^*)+pz_1^*+b = -f'(y_1^*)+f(y_1^*)/y_1^* > 0.$ This, together with $h'(y_1^*)-q > 0$, implies that all eigenvalues must have negative real parts, and hence $E_1^*$ is asymptotically stable. At $E_2^* = (y_2^*, z_2^*)$, the two eigenvalues $\lambda_1$ and $\lambda_2$ satisfy

This implies that one eigenvalue must be positive and the other negative. Thus, $E_2^*$ is a saddle point whenever it exists. Summarizing the above analysis, we obtain the following global stability results on model (2.4).

Theorem 2.4. Assume that ($\textbf{H}_{1}$)-($\textbf{H}_{2}$) are satisfied and $s > q$.

(i) If $y_1^* > \tilde{y}$ holds, then $E_0$ is a saddle point and $E_b$ attracts all solutions with initial conditions in $\{(y(t), z(t))\in \mathbb{R}_+^2: y(0) > 0, z(0)\geq0\}$.

(ii) If $y_1^* < \tilde{y} < y_2^*$ holds, then both $E_0$ and $E_b$ are saddle points and $E_1^*$ attracts all solutions with initial conditions in $\{(y(t), z(t))\in \mathbb{R}_+^2: y(0) > 0, z(0) > 0\}$.

(iii) If $y_1^* < y_2^* < \tilde{y}$ holds, then both $E_0$ and $E_2^*$ are saddle points, and for any initial condition $(y(0), z(0))$ with $y(0) > 0$, the solution of (2.4) approaches either $E_b$ or $E_1^*$, the stable manifold of $E_2^*$ separate the basins of attraction of $E_b$ and $E_1^*$.

3.   Dynamics of (1.3)

In this section, we study the dynamics of model (1.3) with $\tau > 0$.

3.1.   Global stability of the boundary equilibria

By linearizing (1.3) at $E_0 = (0, 0)$, we find two eigenvalues $f'(0)-b > 0$ and $-\mu < 0$. Hence, $E_0$ is a saddle point. The characteristic equation associated with the linearization of model (1.3) at $E_b$ is

One eigenvalue is $\lambda_{1} = f'(\tilde{y})-b < 0$. Hence $E_b$ is locally asymptotically stable if all zeros of

have negative real parts. By [13,Lemma 6], there exists at least one positive eigenvalue if $q\tilde{y}+\mu < e^{-\delta\tau}h(\tilde{y})$, or equivalently, $y_1^* < \tilde{y} < y_2^*$. On the other hand, if $q\tilde{y}+\mu > e^{-\delta\tau}h(\tilde{y})$, that is, either $\tilde{y} < y_1^*$ or $\tilde{y} > y_2^*$, then all eigenvalues have negative real parts. For the critical case $q\tilde{y}+\mu = e^{-\delta\tau}h(\tilde{y})$; namely, either $\tilde{y} = y_1^*$ or $\tilde{y} = y_2^*$, one eigenvalue is $0$, and all other eigenvalues have negative real parts. By using the method of normal forms, we obtain that $E_b$ is locally asymptotically stable. To summarize, we have the following result on the local stability of $E_0$ and $E_b$.

Theorem 3.1. Consider model (1.3) under the assumptions ($\textbf{H}_{1}$)-($\textbf{H}_{3}$). $E_0$ is a saddle point. If either $\tilde{y}\le y_1^*$ or $\tilde{y}\ge y_2^*$, then $E_b$ is locally asymptotically stable; while $E_b$ is unstable if $y_1^* < \tilde{y} < y_2^*$.

To establish the global stability of $E_b$, we first show that the infection is persistent.

Lemma 3.2. $\liminf\limits_{t\rightarrow\infty}y(t) > 0$ for any solution of (1.3) with the initial condition in $X_1 = \{(\phi, \psi)\in X: \phi(0) > 0, \psi(0)\geq0\}$.

Proof. We claim that $\limsup\limits_{t\rightarrow\infty}y(t) > 0$. If not, then $y(t)\to0$ as $t\to\infty$. It then follows from (1.3) that $z(t)\to0$ as $t\to\infty$. This contradicts to the fact that $E_0$ is unstable. Thus, $y(t)$ is weakly persistent. By using [14,Theorem 3.4.6] and Proposition 2.1, there exists a global attractor for the solution semiflow of (1.3). This, together with [15,Theorem 2.3], implies that $y(t)$ is actually strongly persistent; that is, $\liminf\limits_{t\rightarrow\infty}y(t) > 0$.

Our next result is concerned with global stability of $E_b$.

Theorem 3.3. Assume that ($\textbf{H}_{1}$)-($\textbf{H}_{3}$) hold. Then $E_b$ of model (1.3) is globally asymptotically stable in the region $X_1$ provided that $\tau\geq\tau_b$, where $\tau_b: = {1\over \delta}\ln{h'(0)\tilde{y}\over q\tilde{y}+\mu}$. Moreover, the condition $\tau\geq\tau_b$ implies that $\tilde{y}\leq y_1^*$.

Proof. By Proposition 2.1 and Theorem 3.1, it suffices to show that $E_b$ is globally attractive in $X_1\cap\Gamma$, which is a positively invariant set of (1.3). Let $(y(t), z(t))$ be a solution of (1.3) with the initial condition in $X_1$, it follows from Lemma 3.2 that $y(t) > 0$ for $t > 0$. Motivated by the earlier work in ^[16], we construct a Lyapunov functional $L: X_1\cap\Gamma\rightarrow\mathbb{R}$ by

where $c$ is a constant to be determined later. Calculating the time derivative of $L$ along the solution, and using $b = f(\tilde{y})/\tilde{y}$, $h(y)\leq h'(0)y$ for all $y\geq0$, we obtain

Assumption ($\textbf{H}_{1}$) implies that $(\frac{f(y)}{y}-\frac{f(\tilde{y})}{\tilde{y}})(y-\tilde{y})\leq0$ for all $y\leq\bar y$. Since $h'(0)e^{-\delta\tau}-q > 0$, the condition $\tau\geq\tau_b$ is the same as ${p\tilde y\over\mu}\leq {p\over e^{-\delta\tau}h'(0)-q}$. We may choose a constant $c\in[{p\tilde y\over\mu}, {p\over e^{-\delta\tau}h'(0)-q}]$ such that $dL/dt\leq0$. Moreover, $dL/dt = 0$ if and only if $y(t) = \tilde y$ and $z(t) = 0$. Thus the maximal compact invariant set in $\{dL/dt = 0\}$ is the singleton $\{E_b\}$. By the LaSalle invariance principle ^[11], $E_b$ is globally attractive in $X_1\cap\Gamma$. Since $X_1\cap\Gamma$ is absorbing in $X_1$, we conclude that $E_b$ is globally attractive in $X_1$. By Theorem 3.2, $E_b$ is globally asymptotically stable in $X_1$ if $\tau\geq\tau_b$.

If $\tau\geq\tau_b$, then $e^{-\delta\tau}h'(0)\tilde y-q\tilde y-\mu\leq0 = e^{-\delta\tau}h(y_1^*)-qy_1^*-\mu \le e^{-\delta\tau}h'(0)y_1^*-qy_1^*-\mu$. Since $e^{-\delta\tau}h'(0) > q$, we obtain $\tilde y\le y_1^*$.

Here, we remark that if $y_1^* = y_2^* < \tilde{y}$, then two positive equilibria $E_1^*$ and $E_2^*$ coincide, and $E_b$ is locally asymptotically stable.

3.2.   Stability of the positive equilibria and Hopf bifurcation

In this subsection, we investigate the stability of $E_1^*$ and identify parameter regions in which the time delay can destabilize $E_1^*$, lead to Hopf bifurcation and induce sustained oscillations.

The characteristic equation associated with the linearization of system (1.3) at equilibrium $E_i^* \ (i = 1, 2)$ is

where

We first investigate the stability of $E_1^*$, regarding $\tau$ as the bifurcation parameter. In the proof of Theorem 2.4, we have demonstrated that, when $\tau = 0$, all eigenvalues of (3.1) with $i = 1$ lie to the left of the imaginary axis and $E_1^*$ is locally asymptotically stable. As $\tau$ increases, $E_1^*$ may lose its stability only when some eigenvalues cross the imaginary axis to the right. In view of

$0$ is not an eigenvalue for any $\tau\geq0$. For simplicity, we drop the subscript $1$ in the following arguments. Substituting $\lambda = i\omega (\omega > 0)$ into (3.1) and separating the real and imaginary parts, we obtain

Squaring and adding both the above equations lead to

where

Then $i\omega \; (\omega > 0)$ is an imaginary root of (3.1) if and only if $G(\omega, \tau) = 0$. Since $c_1(\tau) > 0$ for all $0\leq\tau\le\tau_{max}$. we know that $G(\omega, \tau) = 0$ has exactly one positive root if and only if $c_0(\tau) < 0$, or equivalently, $a_0 < b_0$.

If $a_0 < b_0$, the implicit function theorem implies that there exists a unique $C^1$ function $\omega = \omega^*(\tau)$ such that $G(\omega^*(\tau), \tau) = 0$ for $0 < \tau\le\tau_{max}$. For $i\omega^*(\tau)$ to be a root of (3.1), $\omega^*(\tau)$ needs to satisfy the system

Set

For $\tau\in I$, let $\theta(\tau)$ be the unique solution of $\sin\theta = g_1$ and $\cos\theta = g_2$ in $(0, 2\pi]$; that is,

Following Beretta and Kuang ^[17], we define

One can check that $\pm i\omega^*(\tau^*)$ are a pair of purely imaginary eigenvalues of (3.1) if and only if $\tau^*$ is a zero of function $S_n(\tau)$ for some $n\in\mathbb{N}$. From [17,Theorem 2.2], we have

Note that ${\partial G\over\partial\omega}(\omega^*(\tau^*), \tau^*) > 0$, thus we have the following result concerning the transversality condition.

Lemma 3.4. Assume that $S_n(\tau)$ has a positive root $\tau^*\in I$ for some $n\in\mathbb{N}$, then a pair of simple purely imaginary roots $\pm i\omega^*(\tau^*)$ of (3.1) exist at $\tau = \tau^*$, and

Moreover, this pair of simple purely imaginary roots $\pm i\omega^*(\tau^*)$ cross the imaginary axis from left to right at $\tau = \tau^*$ if $S'_n(\tau^*) > 0$, and from right to left if $S'_n(\tau^*) < 0$.

It is easily seen that $S_n(\tau) > S_{n+1}(\tau)$ for all $\tau\in I, n\in\mathbb{N}$. When $\tau = 0$, the asymptotic stability of $E_1^*$ implies that $S_0(0) < 0$. Denote

If $\hat\tau < \tau_{max}$, then $c_0(\hat\tau) = 0$, and thus $\omega(\tau)\rightarrow0$ as $\tau\rightarrow\hat\tau^-.$ This, together with (3.3), yields $\lim\limits_{\tau\rightarrow\hat\tau^-}\sin\theta(\tau) = 0$ and $\lim\limits_{\tau\rightarrow\hat\tau^-}\cos\theta(\tau) = -1$. Therefore, $\lim\limits_{\tau\rightarrow\hat\tau^-}\theta(\tau) = \pi, \ \lim\limits_{\tau\rightarrow\hat\tau^-}S_n(\tau) = -\infty.$

If $\sup\limits_{\tau\in I} S_0(\tau) < 0$, $S_n(\tau)$ has no zeros in $I$ for all $n\in\mathbb{N}$, which excludes the existence of purely imaginary eigenvalues and thus implies that $E_1^*$ is locally asymptotically stable for all $\tau\in[0, \tau_{max}]$.

If $\sup\limits_{\tau\in I} S_0(\tau) = 0$, $S_0(\tau)$ has a double zero in $I$, denoted by $\tau_*$, and $S_0'(\tau_*) = 0$. This, together with Lemma 3.4, implies that the transversality condition at $\tau_*$ is not satisfied and all eigenvalues remain to the left of the pure imaginary axis. Thus, $E_1^*$ is still locally asymptotically stable for $\tau\in[0, \tau_{max}]$.

If $\sup\limits_{\tau\in I} S_0(\tau) > 0$, it then follows from $S_0(0) < 0$ and $\lim\limits_{\tau\rightarrow\hat\tau^-}S_n(\tau) = -\infty$ that $S_0(\tau)$ has at least two zeros in $I$. For simplicity, we assume that

(H4) $\sup\limits_{\tau\in I} S_0(\tau) > 0$ and all zeros of $S_n(\tau)$ with odd multiplicity are simple.

Under assumption ($\textbf{H}_{4}$), each $S_n(\tau)$ has even number of simple zeros. Now, we collect all simple zeros of $S_n(\tau)$ with $n\ge0$ and list them in increasing order: $0 < \tau_0 < \tau_1 < \cdots < \tau_{2K-1} < \hat\tau$ $(K\in\mathbb{N}^+)$. For each integer $0\leq i\leq K-1$, we apply Lemma 3.4 to obtain $S'_m(\tau_i) > 0$ and $S'_m(\tau_{2K-i-1}) < 0$ for some $m\in\mathbb{N}$. Therefore, the pair of simple conjugate purely imaginary eigenvalues $\pm i\omega(\tau_i)$ crosses the imaginary axis from left to right, and the pair of simple conjugate purely imaginary eigenvalues $\pm i\omega(\tau_{2K-i-1})$ crosses the imaginary axis from right to left. Thus, system (1.3) undergoes a Hopf bifurcation at $E_1^*$ when $\tau = \tau_j \ (0\le j\le 2K-1)$. Moreover, $E_1^*$ is asymptotically stable for $\tau\in[0, \tau_0)\cup(\tau_{2K-1}, \tau_{max}]$, and unstable for $\tau\in(\tau_0, \tau_{2K-1})$.

For each $k = 0, \cdots, 2K-1$, there exists $n_k$ such that $\tau_k$ is a simple zero of $S_{n_k}(\tau)$. Let $T_k$ be the period of periodic solution bifurcated at $\tau_k$. Applying the Hopf bifurcation theorem for delay differential equations ^[11,18], we have

This, together with $\theta(\tau_k)\in(0, 2\pi]$, implies that $T_k > \tau_k$ if $n_k = 0$, and ${\tau_k\over n_k+1}\le T_k < {\tau_k\over n_k}$ if $n_k > 0$. To summarize, we have the following results on stability of $E_1^*$ and Hopf bifurcation.

Theorem 3.5. Assume that ($\textbf{H}_{1}$)-($\textbf{H}_{3}$) hold, $y_1^* < \tilde y$ and $\hat\tau < \tau_{max}$. Let $I, S_n(\tau), \tau_{max}, \hat\tau$ be defined in (3.4), (3.5), ($\textbf{H}_{3}$) and (3.6), respectively.

(i) If either $I = \emptyset$ or $\sup\limits_{\tau\in I} S_0(\tau)\le 0$, then $E_1^*$ is locally asymptotically stable for all $\tau\in[0, \tau_{max}]$.

(ii) If ($\textbf{H}_{4}$) holds, then there exist exactly $2K$ local Hopf bifurcation values, namely, $0 < \tau_0 < \tau_1 < \cdots < \tau_{2K-1} < \hat\tau$ such that model (3.1) undergoes a Hopf bifurcation at $E_1^*$ when $\tau = \tau_j$ for $0\le j\le 2K-1$. $E_1^*$ is locally asymptotically stable for $\tau\in[0, \tau_0)\cup(\tau_{2K-1}, \tau_{max}]$, and unstable for $\tau\in(\tau_0, \tau_{2K-1})$. Moreover, the period $T_k$ of periodic solution bifurcated at $\tau_k$ satisfies $T_k > \tau_k$ if $\tau_k$ is a simple zero of $S_0(\tau)$, and ${\tau_k\over n_k+1}\le T_k < {\tau_k\over n_k}$ if $\tau_k$ is a simple zero of $S_{n_k}(\tau)$ with $n_k > 0$.

Theorem 3.6. Assume that ($\textbf{H}_{1}$)-($\textbf{H}_{3}$) are satisfied. If $y_1^* < y_2^* < \tilde{y}$ holds, then $E_2^*$ is unstable for all $\tau\geq 0$. Moreover, if $a_{0, 2}\le b_{0, 2}$, then the characteristic equation (3.1) with $i = 2$ has no purely imaginary eigenvalues.

Proof. Note that in (3.1) with $i = 2$, $\Delta_2(0) = p y_{2}^{*}z_{2}^{*}(h'(y_{2}^{*})e^{-\delta\tau}-q) < 0$. For any $\tau\geq 0$, we know that $\Delta_2(\infty) > 0$, thus as long as $E_2^*$ exists, the associated characteristic equation must have at least one positive real root and hence $E_2^*$ is always unstable for $\tau\geq 0$.

Using a similar argument in the study of root distribution for (3.1) at $E_1^*$, we note that $c_{1, 2}(\tau) > 0, c_{0, 2}(\tau) = a_{0, 2}^2-b_{0, 2}^2$ and $a_{0, 2}+b_{0, 2} = p y_{2}^{*}z_{2}^{*}(h'(y_{2}^{*})e^{-\delta\tau}-q) < 0$. Thus, $G(\omega, \tau)$ in (3.2) has no positive zeros if $a_{0, 2}\le b_{0, 2}$. Therefore, (3.1) with $i = 2$ has no purely imaginary eigenvalues if $a_{0, 2}\le b_{0, 2}$.

Since $E_2^*$ is unstable (and thus biologically irrelevant) whenever it exists, we name $E_1^*$, instead of $E_2^*$, the immune control equilibrium (ICE).

3.3.   Onset and termination of Hopf branches

Note that Theorem 3.5 states that if $y_1^* < \tilde y$ and ($\textbf{H}_{4}$) holds, then periodic solutions can bifurcate from $E_1^*$ when $\tau$ is near the local Hopf bifurcation values $\tau_k, k = 0, 1, 2, \ldots, 2K-1.$ For integer $n\ge0$, we denote

According to Theorem 3.5, $J_n$ has even number of bifurcation values, and the period of any periodic solution bifurcated near a bifurcation point in $J^0$ is bounded from both below and above, while the periodic solution bifurcated near a bifurcation point in $J_0$ has a lower bound for its period. As we shall see later in the numerical simulation, it seems impossible to find an upper bound for such period. In what follows, we will restrict our investigation on the set $J^0$ and assume that $J^0\neq\emptyset$. Especially, we will discuss the global continuation of periodic solutions bifurcated from the point $(E_1^*, \tau_k)$ with $\tau_k\in J^0$ as the bifurcation parameter $\tau$ varies. We shall use the global Hopf bifurcation theorem for delay differential equations ^[19] and show that model (1.3) admits periodic solutions globally for all $\tau\in(\underline\tau, \bar\tau)$, where $\underline\tau = \min J^0$ and $\bar\tau = \max J^0$.

Let $x(t) = (y(\tau t), z(\tau t))^T$, model (1.3) can be rewritten as a general functional differential equation:

where $x_{t}(\theta) = x(t+\theta)$ for $\theta\in[-1, 0]$, and $x_{t}\in \mathcal{X}: = C([-1, 0], \mathbb{R}_{+}^2)$, and

with $x_t = (x_{1t}, x_{2t})\in\mathcal{X}$. Identifying the subspace of $\mathcal{X}$ consisting of all constant functions with $\mathbb{R}_{+}^2$, we obtain a restricted function

Obviously, $\widetilde{F}$ is twice continuously differentiable, i.e., the condition $(\textbf{A1})$ in ^[19] holds. We denote the set of stationary solutions of system (3.8) by ${\bf{E}}(F) = \{(\widetilde{x}, \widetilde{\tau}, \widetilde{T})\in\mathbb{R}_{+}^2\times I\times\mathbb{R}_{+}: \; \widetilde{F}(\widetilde{x}, \widetilde{\tau}, \widetilde{T}) = 0 \}.$ It follows from Proposition 2.2 that (ⅰ) if $y_1^* < \tilde{y} < y_2^*$ holds,

(ⅱ) if $y_1^* < y_2^* < \tilde{y}$ holds,

For any stationary solution $(\widetilde{x}, \tau, T)$, the characteristic function is

Note that if either $y_1^* < \tilde{y} < y_2^*$ or $y_1^* < y_2^* < \tilde{y}$ holds, $0$ is not a eigenvalue of any stationary solution of (3.8) and hence the condition (A2) in ^[19] holds. It can be checked easily from (3.9) that the smoothness condition (A3) in ^[19] is satisfied.

Theorem 3.5 implies that if $y_1^* < \tilde y$ and ($\textbf{H}_{4}$) holds, then for each integer $0\le k\le 2K-1$ the stationary solution $(E_1^*, \tau_k, 2\pi/(\omega_k\tau_k))$ is an isolated center of (3.8), where $\omega_k = \omega(\tau_k)$ is the unique positive zero of $G(\omega, \tau)$ in (3.2), and there is only one purely imaginary characteristic value of the form $im (2\pi/\widetilde{T})$ with $m = 1$ and $\widetilde{T} = 2\pi/(\omega_k\tau_k)$. Moreover, if follows from Lemma 3.4 that the crossing number $\zeta_{1}(E_1^*, \tau_k, 2\pi/(\omega_k\tau_k))$ at each of these centers is

Thus the condition (A4) in ^[19] holds. We then define a closed subset $\Sigma(F)$ of $\mathcal{X}\times I \times \mathbb{R}_+$ by

and for each integer $0\le k\le 2K-1$, let ${\bf{C}}(E_1^*, \tau_k, 2\pi/(\omega_k\tau_k))$ denote the connected component of $(E_1^*, \tau_k, 2\pi/(\omega_k\tau_k))$ in $\Sigma(F)$. By Theorem 3.5, ${\bf{C}}(E_1^*, \tau_k, 2\pi/(\omega_k\tau_k))$ is a nonempty subset of $\Sigma(F)$.

To find the interval of $\tau$ in which periodic solutions exist, we shall further investigate the properties of periodic solutions of (3.8).

Lemma 3.7. All nonconstant and nonnegative periodic solutions of (3.8) are uniformly bounded. Actually, we have $0 < y(t), z(t)\le M$ for any $t\in\mathbb{R}$, where $M = \max\{\bar y, h'(0)M_f/(p\gamma)\}$.

Proof. Since $y(t)$ is nonconstant and nonnegative, there exists $t_0 > 0$ such that $y(t_0) > 0$. Integrating the equation for $y'(t)$ gives

Hence, $y(t) > 0$ for all $t > t_0$. Since $y$ is periodic, we have $y(t) > 0$ for all $t > 0$. Similarly, if $z(t_0) > 0$ for some $t_0 > 0$, we integrate the equation for $z'(t)$ to obtain

for all $t > t_0$. This together with periodicity of $z$ implies that $z(t) > 0$ for all $t > 0$.

By Proposition 2.1, we have $\limsup\limits_{t\rightarrow\infty}y(t)\le\bar y$. We claim $y(t)\le M$ for all $t\ge 0$. Otherwise, if $y(t_1) > M$ for some $t_1\ge0$, then

where $T$ is the period of the nonconstant and nonnegative periodic solutions. This contradicts the fact that $y(t)$ is eventually bounded above by $M$ as $t\rightarrow\infty$. Thus, $M$ is the upper bound of $y(t)$. Similarly, from Proposition 2.1, we have $\limsup\limits_{t\rightarrow\infty}z(t)\le h'(0)M_f/(p\gamma)\le M$. This implies that $z(t)$ is uniformly bounded by $M$.

Lemma 3.8. System (3.8) has no nontrivial periodic solution of period 1.

Proof. Note that any nontrivial $1$-periodic solution $x(t)$ of system (3.8) is also a nontrivial periodic solution of model (2.4), which does not admit any nontrivial periodic solution via Lemma 2.3.

We are now in the position to state our result concerning the properties of the global Hopf branches.

Theorem 3.9. Assume that ($\textbf{H}_{1}$)-($\textbf{H}_{4}$), either $y_1^* < \tilde{y} < y_2^*$ or $y_1^* < y_2^* < \tilde{y}$ holds, $\hat\tau < \tau_{max}, J^0\neq\emptyset$ and $a_{0, 2}\le b_{0, 2}$. Then for system (3.8), we have the following results:

(i) The connected component ${\bf{C}}(E_1^*, \tau_k, 2\pi/(\omega_k\tau_k))$ is bounded for all $\tau_k\in J^0$.

(ii) Let $n\ge1$. If $J_n$ has $2k_n$ bifurcation values, ordered as $\tau_{n, 1} < \cdots < \tau_{n, 2k_n}$, then for each $i = 1, \cdots, k_n$, there exists $j > k_n$ such that $\tau_{n, i}$ and $\tau_{n, j}$ are connected by a global Hopf branch. Similarly, for each $j = k_n+1, \cdots, 2k_n$, there exists $i\le k_n$ such that $\tau_{n, i}$ and $\tau_{n, j}$ are connected by a global Hopf branch. Especially, if $k_n = 1$, then the two bifurcation values of $J_n$ are connected.

(iii) For each $\tau\in(\min J_n, \max J_n)$ with integer $n\ge1$, system (3.8) has at least one periodic solution with period in $(1/(n+1), 1/n)$.

Proof. Lemma 3.7 implies that the projection of ${\bf{C}}(E_1^*, \tau_k, 2\pi/(\omega_k\tau_k))$ onto $\mathcal{X}$ is bounded. Note that system (3.8) has no periodic solutions of period $1$, and thus no periodic solutions of period $1/n_{k}$ or $1/(n_{k}+1)$ for any positive integer $n_{k}$. It follows from Lemma 3.8 that the period $T$ of a periodic solution on the connected component ${\bf{C}}(E_1^*, \tau_k, 2\pi/(\omega_k\tau_k))$ satisfies

with $\tau_{k}\in J_{n_k}$ and $n_k\ge1$. Hence, the projection of ${\bf{C}}(E_1^*, \tau_k, 2\pi/(\omega_k\tau_k))$ onto the $T$-space is bounded for $\tau_k\in J^0$. Note that $\tau\in I$ and $I$ is a bounded interval. Therefore, ${\bf{C}}(E_1^*, \tau_k, 2\pi/(\omega_k\tau_k))$ is bounded in $\mathcal{X}\times I\times \mathbb{R}_+$ for any $\tau_k\in J^0$. This proves (ⅰ).

From the proof of Theorems 3.1 and 3.6, we know that the stationary solutions $(E_0, \tau, T)$, $(E_b, \tau, T)$ and $(E_2^*, \tau, T)$ cannot be a center for any $\tau$ and $T$ if either $y_1^* < \tilde{y} < y_2^*$ or $y_1^* < y_2^* < \tilde{y}$ holds and $a_{0, 2}\le b_{0, 2}$. It then follows from the global bifurcation theorem [19,Theorem 3.3] that $\mathcal{E}: = {\bf{C}}(E_1^*, \tau_k, 2\pi/(\omega_k\tau_k))\cap {\bf{E}}(F)$ is finite and

If $J_n$ has $2k_n$ bifurcation values, ordered as $\tau_{n, 1} < \cdots < \tau_{n, 2k_n}$, then by (3.10), $\zeta_1(\widetilde{x}, \tau_{n, i}, T) = -1$ for each $i = 1, \cdots, k_n$, and $\zeta_1(\widetilde{x}, \tau_{n, j}, T) = 1$ for each $j = k_n+1, \cdots, 2k_n$. This together with the above equation implies (ⅱ).

Note that $\tau_{n, 1} = \min J_n$ and $\tau_{n, 2k_n} = \max J_n$. The $\tau$ projection of the global bifurcation branch bifurcated from $\tau_{n, 1}$ includes $(\tau_{n, 1}, \tau_{n, k_n+1})$, while the $\tau$ projection of the global bifurcation branch bifurcated from $\tau_{n, 2k_n}$ includes $(\tau_{n, k_n}, \tau_{n, 2k_n})$. Thus, for any $\tau\in(\tau_{n, 1}, \tau_{n, 2k})$, system (3.8) has at least one periodic solution. According to the proof of (ⅰ), the period of this periodic solution lies in $(1/(n+1), 1/n)$. Thus, (ⅲ) is proved.

4.   Numerical exploration

In this section, we present some numerical simulations to demonstrate our theoretical results. We will use the Matlab package DDE-BIFTOOL developed by Engelborghs et al. ^[20,21] to sketch the global Hopf branches ${\bf{C}}(E_1^*, \tau_k, 2\pi/(\omega_k\tau_k))$ for $0\le k\le 2K-1$. Following the study of ODE model in ^[6], we choose

and set the parameter values as follows:

Here, the units for viral loads $y$ and immune cell density $z$ are virus per $\text{mm}^3$ and cells per $\text{mm}^3$, respectively. The delay $\tau$ has the unit in days. It is easy to calculate

From Proposition 2.2, there are exactly two equilibria $E_{0} = (0, 0)$ and $E_{b} = (1.5, 0)$ when $\tau > \tau_{max}$, and, if $\tau\leq\tau_{max}$, there exists a third equilibrium $E_1^* = (y_1^*, z_1^*)$, where $y_1^*, z_1^*$ depend on $\tau$. This is the unique positive equilibrium when $\tau\in[0, \tau_e]$ and there exists another positive equilibrium $E_2^* = (y_2^*, z_2^*)$ when $\tau\in(\tau_e, \tau_{max})$; see Figure 1. Here, $\tau_e$ is the critical value when $y_2^* = \tilde y$; that is,

At $\tau = \tau_{max}$, the two positive equilibria coincide: $y_1^* = y_2^*$.

Theorem 3.1 implies that $E_0$ is a saddle point, and $E_b$ is locally asymptotically stable for $\tau\ge\tau_e$, and unstable for $\tau < \tau_e$, Moreover, by Theorem 3.3, $E_b$ is globally asymptotically stable in the region $X_1$ provided that $\tau\geq\tau_b$. It can be verified that the conditions of case (ⅱ) in Theorem 3.5 are satisfied only if $0\leq\tau < \hat\tau$. By Theorem 3.5, we know that there are exactly $4$ local Hopf bifurcation values with $K = 2$, namely,

as shown in Figure 2. Correspondingly,

Clearly, $E_1^*$ is stable for $\tau\in[0, \tau_0)\cup(\tau_3, \tau_{max}]$ and unstable for $\tau\in(\tau_0, \tau_3)$, which implies that $\tau_0$ and $\tau_3$ are stability switches. Remark that $\tau_3 < \tau_e$.

In Figure 3, we plot the four global Hopf branches ${\bf{C}}(E_1^*, \tau_k, 2\pi/(\omega_k\tau_k))$ with $0\le k\le 3$. It is observed that the branches ${\bf{C}}(E_1^*, \tau_1, 2\pi/(\omega_1\tau_1))$ and ${\bf{C}}(E_1^*, \tau_2, 2\pi/(\omega_2\tau_2))$ are bounded and connected, which agrees with Theorem 3.9. The periodic solutions on these global Hopf branches are plotted on the phase plane of $y$ and $z$; see Figures 4–6. As we know, the stability of the periodic solution is completely determined by the associated principal Floquet multiplier: if the principal Floquet multiplier is larger than $1$, then the corresponding periodic solution is unstable, otherwise, the bifurcated periodic solution is stable ^[11].

By using DDE-BIFTOOL, we can further calculate the associated principal Floquet multipliers; see Figure 7. Note that the periodic solution on the first branch ${\bf{C}}(E_1^*, \tau_0, 2\pi/(\omega_0\tau_0))$ is stable for small $\tau$, and becomes unstable as $\tau$ increases; the periodic solution on the second branch ${\bf{C}}(E_1^*, \tau_1, 2\pi/(\omega_1\tau_1))$–or equivalently, the third branch ${\bf{C}}(E_1^*, \tau_2, 2\pi/(\omega_2\tau_2))$–is always unstable; the periodic solution on the fourth branch ${\bf{C}}(E_1^*, \tau_3, 2\pi/(\omega_3\tau_3))$ is stable for $\tau$ near $\tau_3$, and becomes unstable as $\tau$ decreases. Moreover, for any $\tau\in(\tau_1, \tau_2)$, model (1.3) has at least one periodic solution on the second (or third) global Hopf branch with period in the interval $(\tau/2, \tau)$, as shown in Figure 8.

It is interesting to note from Figure 3 that the first and the fourth Hopf branches do not connect. This phenomenon is quite different from that in the existing literature. We have proved that both the periodic solutions and delays on any global Hopf branch are uniformly bounded. According to the global Hopf bifurcation theory ^[19], either the first and the fourth Hopf branches connect, or these two branches have unbounded periods. From Figure 8, we observe that the periods of periodic oscillations on the fourth branch ${\bf{C}}(E_1^*, \tau_3, 2\pi/(\omega_3\tau_3))$ seem to be unbounded. To further explore this interesting phenomenon, we numerically sketch in Figure 9 a bifurcation diagram for model (1.3). It is suggested that model (1.3) can have chaotic solutions and thus allow aperiodic oscillations when $\tau$ lies in the interval where the first and fourth Hopf branches disappear. Figure 9 also confirms that $\tau_0$ and $\tau_3$ are stability switches: (ⅰ) when $\tau$ crosses $\tau_0$ from left to right, the positive equilibrium $E_1^*$ becomes unstable, while a stable periodic solution bifurcates; and (ⅱ) when $\tau$ crosses $\tau_3$ from left to right, $E_1^*$ regains its stability while there exists a stable periodic solution when $\tau$ is close to $\tau_3$ from the left. It should be mentioned that the bifurcation diagram has a jump at $\tau_{max}$, which is due to the fact that there exists no positive equilibrium when $\tau > \tau_{max}$ but the positive equilibrium at $\tau = \tau_{max}$ is nonzero.

5.   Summary and discussion

Considering the fact the immune response process involves a sequence of events such as antigenic activation, selection, and proliferation of the immune cells and the mortality of immune cells during the time lag, we derive an immunosuppressive infection model from a stage structured population model. We find that the time lag $\tau$ plays a key role in the infection and immune response dynamics. As $\tau$ increases, the model dynamics shifts through four possible outcomes: (ⅰ) virus persists without immune response when $\tau > \tau_{max}$; (ⅱ) the outcome is initial condition dependent and delay dependent, either the immunity can be completely destroyed, or the immune mediated control can be established when $\tau\in(\tau_e, \tau_{max}]$; (ⅲ) immune response controls the virus when $\tau\in[0, \tau_0)\cup(\tau_{2K-1}, \tau_e)$; and (ⅳ) sustained oscillation appears or chaotic dynamics occurs when $\tau\in(\tau_0, \tau_{2K-1})$. Our theoretical results coincide with the phenomenon of oscillatory viral loads and immune cells observed in clinical data ^[9].

From numerical simulation, we also observe an unbounded global Hopf bifurcation branch with bounded periodic solutions, bounded delays, but unbounded periods. To the best of our knowledge, this phenomenon is new, because, in the existing literature, either the global Hopf branch is bounded or it is unbounded with unbounded delays. A detailed theoretical study of this interesting numerical phenomenon should enrich the global Hopf bifurcation theory for delay differential equations, and we leave it as a future work.

Acknowledgments

We are very grateful to the anonymous referees for their valuable comments and suggestions, which help to improve the presentation of this manuscript. H. Shu was partially supported by the National Natural Science Foundation of China (No.11971285, No.11601392), and the Fundamental Research Funds for the Central Universities.

Conflict of interest

The authors declare there is no conflict of interest.