Stability and Hopf bifurcation of HIV-1 model with Holling II infection rate and immune delay

Abstract

This paper aims to analyse stability and Hopf bifurcation of the HIV-1 model with immune delay under the functional response of the Holling II type. The global stability analysis has been considered by Lyapunov–LaSalle theorem. And stability and the sufficient condition for the existence of Hopf Bifurcation of the infected equilibrium of the HIV-1 model with immune response are also studied. Some numerical simulations verify the above results. Finally, we propose a novel three dimension system to the future study.

Keywords:

• Stability
• Hopf bifurcation
• HIV-1 model
• immune delay

AMS Mathematics Subject Classifications:

• 34K20
• 34K60
• 92B05
• 92D30

1. Introduction

In recent years, people pay more attention to damages of the immune system caused by HIV virus. According to recent studies, they found that latent infected cells could transform themselves into health cells by autoimmune response before that viral genome is integrated into cellular genome (e.g. see [15]). Therefore, some scholars began to study HIV-1 models with latent infected cells and the corresponding dynamic properties (e.g. see [1–14, 16–33]). The literature (e.g. see [1]) considered the following model: (1) $\{\begin{array}{l}{\displaystyle \frac{\mathrm{d}x(t)}{\mathrm{d}t}}=s-d_{1}x(t)-\beta x(t)v(t)+\delta w(t),\\ {\displaystyle \frac{\mathrm{d}w(t)}{\mathrm{d}t}}=\beta x(t)v(t)-(\delta +d_2+q)w(t),\\ {\displaystyle \frac{\mathrm{d}y(t)}{\mathrm{d}t}}=qw(t)-d_{3}y(t),\\ {\displaystyle \frac{\mathrm{d}v(t)}{\mathrm{d}t}}=\sigma y(t)-\gamma v(t),\end{array}$(1) where $x(t),w(t),y(t),v(t)$ denote the concentration of the uninfected $CD{4}^{+}$T cells, latent infected cells, infected cells and virus at time t, respectively. And $s(s>0)$ is the recruitment rate of uninfected T cells, the $\beta xv$ is the bilinear incidence of the healthy cells caused by HIV virus, and δ represents the proportion that latent cells restore to healthy cells before integrating into viral genome. Moreover, $d_1,d_2,d_3,\gamma$ respectively denote the death rate of the uninfected T cells, the latent infected cells, infected cells and the virus. Finally, q is the rate of at which the latent infected cells change into infected cells and σ is the rate at which cells release the virus.

It is noticed that the disease incidence rate is bilinear in Model (1). However, studies have shown that when the number of the target cells is large enough (see [3]), bilinear incidence may not be a valid assumption, namely the virus and host cell is a nonlinear relationship. Therefore, we consider the Holling II $\frac{\beta x(t)y(t)}{1+y(t)}$ instead of bilinear incidence rate $\beta x(t)v(t)$, which is more in line with the actual situation. As we know that the body's immune is an important factor, in the process of inhibition and destroy infected cells. The body's immune system could delay production when the body accepts to produce in the process of lymphocyte antigen stimulation (see [7, 20]). So we establish the following model: (2) $\{\begin{array}{l}{\displaystyle \frac{\mathrm{d}x(t)}{\mathrm{d}t}=s-d_{1}x(t)-\frac{\beta x(t)y(t)}{1+y(t)}+\delta w(t),}\\ {\displaystyle \frac{\mathrm{d}w(t)}{\mathrm{d}t}=\frac{\beta x(t)y(t)}{1+y(t)}-(\delta +d_2+q)w(t),}\\ {\displaystyle \frac{\mathrm{d}y(t)}{\mathrm{d}t}}=qw(t)-d_{3}y(t)-hy(t)z(t),\\ {\displaystyle \frac{\mathrm{d}z(t)}{\mathrm{d}t}}=ky(t-\tau )z(t-\tau )-d_{4}z(t),\end{array}$(2) where $z(t)$ denotes the concentration of the Immune cells at time t, τ denotes the immune delay, $d_4$ denotes the death rate of the immune cell, $ky(t-\tau )z(t-\tau )$ denotes the birth rate of immune cell.

In this paper, we mainly study the global stability and Hopf bifurcation of HIV-1 model (2) with immune delay under the functional response of the Holling II type.

The remainders of this paper are as follows. In Section 2, we consider the positivity and boundaries of solutions and equilibria of model (2). In Section 3, we mainly study the global stability of the viral free equilibrium and infected equilibrium by Lyapunov–LaSalle theorem. In Section 4, we mainly study the global stability and the existence of Hopf Bifurcation of the infected equilibrium of the HIV-1 model with immune response. In Section 5, some numerical simulations are performed to illustrate the main results. The sixth part gives some conclusions and prospects. Here, we propose a novel three dimension system to the future study.

2. Positivity and boundaries of solutions and equilibria given

Considering the biological significance of the model, we assume that the initial conditions of system (2(2) $\{\begin{array}{l}{\displaystyle \frac{\mathrm{d}x(t)}{\mathrm{d}t}=s-d_{1}x(t)-\frac{\beta x(t)y(t)}{1+y(t)}+\delta w(t),}\\ {\displaystyle \frac{\mathrm{d}w(t)}{\mathrm{d}t}=\frac{\beta x(t)y(t)}{1+y(t)}-(\delta +d_2+q)w(t),}\\ {\displaystyle \frac{\mathrm{d}y(t)}{\mathrm{d}t}}=qw(t)-d_{3}y(t)-hy(t)z(t),\\ {\displaystyle \frac{\mathrm{d}z(t)}{\mathrm{d}t}}=ky(t-\tau )z(t-\tau )-d_{4}z(t),\end{array}$(2) ) are as follows: (3) $\{\begin{array}{l}x(\theta )={\phi }_1(\theta ),w(\theta )={\phi }_2(\theta ),y(\theta )={\phi }_3(\theta ),z(\theta )={\phi }_4(\theta ),\\ {\phi }_1\ge 0,{\phi }_i>0,i=2,3,4,\theta \in [-\tau ,0],\end{array}$(3) where $\phi =({\phi }_1,{\phi }_2,{\phi }_3,{\phi }_4{)}^T\in C([-\tau ,0],R_{+}^4)$. It expresses a continuous function from [$-\tau$,0] to $R_{+}^4$ and with supremum norm in Banach space. ($R_{+}^4$=($x_1$,$x_2$,$x_3$,$x_4$): $x_i\ge 0$, i = 1, 2, 3, 4 .)

Define the infection of the basic reproductive number $R_0$ and the basic immune response reproductive number $R_1$. Through calculation we can get $R_0=\frac{\beta qs}{d_1{d}_3(\delta +d_2+q)},R_1=\frac{kq\beta s}{d_3[\beta d_4(d_2+q)+d_1(k+d_4)(\delta +d_2+q)]}.$It is also easy to know system (2(2) $\{\begin{array}{l}{\displaystyle \frac{\mathrm{d}x(t)}{\mathrm{d}t}=s-d_{1}x(t)-\frac{\beta x(t)y(t)}{1+y(t)}+\delta w(t),}\\ {\displaystyle \frac{\mathrm{d}w(t)}{\mathrm{d}t}=\frac{\beta x(t)y(t)}{1+y(t)}-(\delta +d_2+q)w(t),}\\ {\displaystyle \frac{\mathrm{d}y(t)}{\mathrm{d}t}}=qw(t)-d_{3}y(t)-hy(t)z(t),\\ {\displaystyle \frac{\mathrm{d}z(t)}{\mathrm{d}t}}=ky(t-\tau )z(t-\tau )-d_{4}z(t),\end{array}$(2) ) has following three equilibrium points:

1. If $R_0<1$, there exists an uninfection equilibrium $E_0=(\frac{s}{d_1},0,0,0)$.

2. If $R_0>1$ and $R_1<1$, there exists an infected equilibrium without immune response $E_1(x_1,w_1,y_1,z_1)$, where $\begin{array}{rl}{x}_1& =\frac{sq(\delta +d_2+q)+(d_2+q)d_3(\delta +d_2+q)}{\beta q(d_2+q)+d_{1}q(\delta +d_2+q)},\\ w_1& =\frac{\beta sq-d_1{d}_3(\delta +d_2+q)}{\beta q(d_2+q)+d_{1}q(\delta +d_2+q)},\\ y_1& =\frac{\beta sq-d_1{d}_3(\delta +d_2+q)}{\beta d_3(d_2+q)+d_1{d}_3(\delta +d_2+q)},\\ z_1& =0.\end{array}$

3. If $R_1>1$, there exists an infected equilibrium $E_2(x_2,w_2,y_2,z_2)$ with immune response, where $\begin{array}{rl}{x}_2& =\frac{s(k+d_4)(\delta +d_2+q)}{\beta d_4(d_2+q)+d_1(k+d_4)(\delta +d_2+q)},\\ w_2& =\frac{s\beta d_4}{\beta d_4(d_2+q)+d_1(k+d_4)(\delta +d_2+q)},\\ y_2& =\frac{d_4}{k},\\ z_2& =\frac{\beta sqk-d_3[d_1(k+d_4)(\delta +d_2+q)+\beta d_4(d_2+q)]}{[\beta d_4(d_2+q)+d_1(k+d_4)(\delta +d_2+q)]h}.\end{array}$

Lemma 2.1

Suppose that $x(t),w(t),y(t),z(t)$ is the solution of system (2(2) $\{\begin{array}{l}{\displaystyle \frac{\mathrm{d}x(t)}{\mathrm{d}t}=s-d_{1}x(t)-\frac{\beta x(t)y(t)}{1+y(t)}+\delta w(t),}\\ {\displaystyle \frac{\mathrm{d}w(t)}{\mathrm{d}t}=\frac{\beta x(t)y(t)}{1+y(t)}-(\delta +d_2+q)w(t),}\\ {\displaystyle \frac{\mathrm{d}y(t)}{\mathrm{d}t}}=qw(t)-d_{3}y(t)-hy(t)z(t),\\ {\displaystyle \frac{\mathrm{d}z(t)}{\mathrm{d}t}}=ky(t-\tau )z(t-\tau )-d_{4}z(t),\end{array}$(2) ) satisfying initial conditions (3(3) $\{\begin{array}{l}x(\theta )={\phi }_1(\theta ),w(\theta )={\phi }_2(\theta ),y(\theta )={\phi }_3(\theta ),z(\theta )={\phi }_4(\theta ),\\ {\phi }_1\ge 0,{\phi }_i>0,i=2,3,4,\theta \in [-\tau ,0],\end{array}$(3) ). Then $x(t)>0,w(t)>0,y(t)>0,z(t)>0$ for all $t\ge 0$.

Proof.

Assume that $t_1$ is the first point satisfy $t_1$=min{$t>0:x(t)\times w(t)\times y(t)\times z(t)=0$}.

(1) If $x(t_1)=0$, from the first equation of system (2(2) $\{\begin{array}{l}{\displaystyle \frac{\mathrm{d}x(t)}{\mathrm{d}t}=s-d_{1}x(t)-\frac{\beta x(t)y(t)}{1+y(t)}+\delta w(t),}\\ {\displaystyle \frac{\mathrm{d}w(t)}{\mathrm{d}t}=\frac{\beta x(t)y(t)}{1+y(t)}-(\delta +d_2+q)w(t),}\\ {\displaystyle \frac{\mathrm{d}y(t)}{\mathrm{d}t}}=qw(t)-d_{3}y(t)-hy(t)z(t),\\ {\displaystyle \frac{\mathrm{d}z(t)}{\mathrm{d}t}}=ky(t-\tau )z(t-\tau )-d_{4}z(t),\end{array}$(2) ) we can know $\dot{x}(t_1)=s-d_{1}x(t_1)-\frac{\beta x(t_1)y(t_1)}{1+y(t_1)}+\delta w(t_1)=s+\delta w(t_1),$because $t_1$ is the first time meet $x(t)\times w(t)\times y(t)\times z(t)=0$, so $w(t_1)\ge 0,y(t_1)\ge 0,z(t_1)\ge 0$. It is easy to know $\dot{x}(t_1)>0$. So for any sufficiently small ${\epsilon }_1$, when $t\in (t_1-{\epsilon }_1,t_1)$, we have $x(t)<0$, But on the other hand, $x(t)>0$, $t\in [0,t_1)$ the assumption $x(t_1)=0$ does not hold. Hence $x(t)>0$, $t\in (0,\mathrm{\infty })$.

(2) If $w(t_1)=0$, from the second equation of system (2(2) $\{\begin{array}{l}{\displaystyle \frac{\mathrm{d}x(t)}{\mathrm{d}t}=s-d_{1}x(t)-\frac{\beta x(t)y(t)}{1+y(t)}+\delta w(t),}\\ {\displaystyle \frac{\mathrm{d}w(t)}{\mathrm{d}t}=\frac{\beta x(t)y(t)}{1+y(t)}-(\delta +d_2+q)w(t),}\\ {\displaystyle \frac{\mathrm{d}y(t)}{\mathrm{d}t}}=qw(t)-d_{3}y(t)-hy(t)z(t),\\ {\displaystyle \frac{\mathrm{d}z(t)}{\mathrm{d}t}}=ky(t-\tau )z(t-\tau )-d_{4}z(t),\end{array}$(2) ) we can know $w(t_1)={\mathrm{e}}^{-(\delta +d_2+q)t_1}w(0)+{\int }_0^{t_1}\frac{\beta x(\eta )y(\eta )}{1+y(\eta )}{\mathrm{e}}^{-(\delta +d_2+q)(t_1-\eta )},\mathrm{d}\eta >0,$This is contradictory with $w(t_1)=0$, so we can't find any $t_1$ to meet $w(t_1)=0$. Therefore $w(t)>0$, $t\in (0,+\mathrm{\infty })$. By a recursive demonstration and initial conditions, we can easily get $y(t)>0,z(t)>0,t\in (0,+\mathrm{\infty })$. The proof is completed.

Lemma 2.2

Suppose that $x(t),w(t),y(t),z(t)$ are the solutions of system (2(2) $\{\begin{array}{l}{\displaystyle \frac{\mathrm{d}x(t)}{\mathrm{d}t}=s-d_{1}x(t)-\frac{\beta x(t)y(t)}{1+y(t)}+\delta w(t),}\\ {\displaystyle \frac{\mathrm{d}w(t)}{\mathrm{d}t}=\frac{\beta x(t)y(t)}{1+y(t)}-(\delta +d_2+q)w(t),}\\ {\displaystyle \frac{\mathrm{d}y(t)}{\mathrm{d}t}}=qw(t)-d_{3}y(t)-hy(t)z(t),\\ {\displaystyle \frac{\mathrm{d}z(t)}{\mathrm{d}t}}=ky(t-\tau )z(t-\tau )-d_{4}z(t),\end{array}$(2) ), each of them is bounded.

Proof.

We define $F(t)=x(t)+w(t)+y(t)+\frac{h}{k}z(t+\tau ),m=min\{d_1,d_2,d_3,d_4\}.$For boundedness of the solution, calculating the derivative of $F(t)$, we get $\begin{array}{rl}{F}^{\mathrm{\prime }}(t)& =x^{\mathrm{\prime }}(t)+w^{\mathrm{\prime }}(t)+y^{\mathrm{\prime }}(t)+\frac{h}{k}{z}^{\mathrm{\prime }}(t+\tau )\\ & =s-[d_{1}x(t)+d_{2}w(t)+d_{3}y(t)+\frac{d_{4}h}{k}z(t+\tau )]\\ & <s-m[x(t)+w(t)+y(t)+\frac{h}{k}z(t+\tau )]\\ & =s-mF(t),\end{array}$where $F(t)<\epsilon +\frac{s}{m}$ (positive number ε can be arbitrarily small). This implies that $F(t)$ is bounded by the comparison theorem, and so are $x(t),w(t),y(t)$ and $z(t)$. The proof is completed.

3. Stability analysis of equilibrium point $E_0$ and $E_1$

In this section, we mainly consider the stability of the viral free equilibrium $E_0$ by employing Lyapunov function.

Theorem 3.1

If $R_0<1$, the uninfected equilibrium $E_0$ of system (2(2) $\{\begin{array}{l}{\displaystyle \frac{\mathrm{d}x(t)}{\mathrm{d}t}=s-d_{1}x(t)-\frac{\beta x(t)y(t)}{1+y(t)}+\delta w(t),}\\ {\displaystyle \frac{\mathrm{d}w(t)}{\mathrm{d}t}=\frac{\beta x(t)y(t)}{1+y(t)}-(\delta +d_2+q)w(t),}\\ {\displaystyle \frac{\mathrm{d}y(t)}{\mathrm{d}t}}=qw(t)-d_{3}y(t)-hy(t)z(t),\\ {\displaystyle \frac{\mathrm{d}z(t)}{\mathrm{d}t}}=ky(t-\tau )z(t-\tau )-d_{4}z(t),\end{array}$(2) ) is globally asymptotical stable for any time delay $\tau \ge 0$.

Proof.

Define the Lyapunov function $V_0$ as follows: $\begin{array}{rl}{V}_0=& x(t)-x_0-{\int }_{x_0}^{x(t)}\frac{x_0}{\theta }d\theta +w(t)+\frac{\delta }{2(d_1+d_2+q)x_0}[(x(t)-x_0)+w(t){]}^2\\ & +\frac{\delta +d_2+q}{q}y(t)+\frac{(\delta +d_2+q)h}{qk}z(t)+\frac{(\delta +d_2+q)h}{q}{\int }_{t-\tau }^{t}y(\theta )z(\theta )\mathrm{d}\theta .\end{array}$Through $x_0=\frac{s}{d_1}$, we can push that $s=d_1{x}_0$. Calculating the derivative of $V_0$ along the solution of system (2(2) $\{\begin{array}{l}{\displaystyle \frac{\mathrm{d}x(t)}{\mathrm{d}t}=s-d_{1}x(t)-\frac{\beta x(t)y(t)}{1+y(t)}+\delta w(t),}\\ {\displaystyle \frac{\mathrm{d}w(t)}{\mathrm{d}t}=\frac{\beta x(t)y(t)}{1+y(t)}-(\delta +d_2+q)w(t),}\\ {\displaystyle \frac{\mathrm{d}y(t)}{\mathrm{d}t}}=qw(t)-d_{3}y(t)-hy(t)z(t),\\ {\displaystyle \frac{\mathrm{d}z(t)}{\mathrm{d}t}}=ky(t-\tau )z(t-\tau )-d_{4}z(t),\end{array}$(2) ), we get $\begin{array}{rl}{V}_0^{\mathrm{\prime }}& =(d_1{x}_0-d_{1}x(t)-\frac{\beta x(t)y(t)}{1+y(t)}+\delta w(t))(1-\frac{x_0}{x(t)})+\frac{\beta x(t)y(t)}{1+y(t)}-(\delta +d_2+q)w(t)\\ & +\frac{\delta }{(d_1+d_2+q)x_0}(x(t)-x_0+w(t))[d_1{x}_0-d_{1}x(t)+\delta w(t)-(\delta +d_2+q)w(t)]\\ & +\frac{\delta +d_2+q}{q}(qw(t)-d_{3}y(t)-hy(t)z(t))\\ & +\frac{(\delta +d_2+q)h}{qk}[ky(t-\tau )z(t-\tau )-d_{4}z(t)]\\ & +\frac{(\delta +d_2+q)h}{qk}y(t)z(t)-\frac{(\delta +d_2+q)h}{qk}[ky(t-\tau )z(t-\tau )-d_{4}z(t)]\\ & =\frac{-d_1(x(t)-x_0)}{x(t)}+\delta w(t)(1-\frac{x_0}{x(t)})+\frac{\beta x_{0}y(t)}{1+y(t)}+\frac{-\delta d_1(x(t)-x_0{)}^2}{(d_1+d_2+q)x_0}\\ & -\delta w(t)(\frac{x(t)}{x_0}-1)-\frac{\delta (d_2+q)w(t{)}^2}{(d_1+d_2+q)x_0}-\frac{\delta +q+d_2}{q}{d}_{3}y(t)-\frac{(\delta +q+d_2)h}{qk}{d}_{4}z(t)\\ & =-(\frac{1}{x(t)}+\frac{\delta }{(d_1+d_2+q)x_0})d_1(x(t)-x_0{)}^2+\delta w(t)(2-\frac{x_0}{x(t)}-\frac{x(t)}{x_0})\\ & +\beta x_{0}y(t)(\frac{1}{1+y(t)}-1)+(R_0-1)\frac{\delta +q+d_2}{q}{d}_{3}y(t)\\ & -\frac{(\delta +q+d_2)h}{qk}{d}_{4}z(t)-\frac{\delta (d_2+q)w(t{)}^2}{(d_1+d_2+q)x_0}.\end{array}$Since $2-\frac{x_0}{x}-\frac{x}{x_0}\le 0$ and $R_0<1$ hold, we can get the above $V_0^{\mathrm{\prime }}\le 0$. In addition to that only if $(x(t),w(t),y(t),z(t))=(\frac{s}{d_1},0,0,0)$, we obtain $V_0^{\mathrm{\prime }}=0$. The uninfected equilibrium $E_0$ of system (2(2) $\{\begin{array}{l}{\displaystyle \frac{\mathrm{d}x(t)}{\mathrm{d}t}=s-d_{1}x(t)-\frac{\beta x(t)y(t)}{1+y(t)}+\delta w(t),}\\ {\displaystyle \frac{\mathrm{d}w(t)}{\mathrm{d}t}=\frac{\beta x(t)y(t)}{1+y(t)}-(\delta +d_2+q)w(t),}\\ {\displaystyle \frac{\mathrm{d}y(t)}{\mathrm{d}t}}=qw(t)-d_{3}y(t)-hy(t)z(t),\\ {\displaystyle \frac{\mathrm{d}z(t)}{\mathrm{d}t}}=ky(t-\tau )z(t-\tau )-d_{4}z(t),\end{array}$(2) ) is globally asymptotically stable according to the Lyapunov–LaSalle theorem in [2]. The proof is complete.

Theorem 3.2

If $R_0\in (1,1+\frac{sq+(d_2+q)d_3}{\delta d_3}]$ and $R_1<1$ hold, the infected equilibrium $E_1$ without immune response of system (2(2) $\{\begin{array}{l}{\displaystyle \frac{\mathrm{d}x(t)}{\mathrm{d}t}=s-d_{1}x(t)-\frac{\beta x(t)y(t)}{1+y(t)}+\delta w(t),}\\ {\displaystyle \frac{\mathrm{d}w(t)}{\mathrm{d}t}=\frac{\beta x(t)y(t)}{1+y(t)}-(\delta +d_2+q)w(t),}\\ {\displaystyle \frac{\mathrm{d}y(t)}{\mathrm{d}t}}=qw(t)-d_{3}y(t)-hy(t)z(t),\\ {\displaystyle \frac{\mathrm{d}z(t)}{\mathrm{d}t}}=ky(t-\tau )z(t-\tau )-d_{4}z(t),\end{array}$(2) ) is globally asymptotical stable for any time delay $\tau \ge 0$.

Proof.

Define the Lyapunov function $V_1$ as follows: $\begin{array}{rl}{V}_1& =x(t)-x_1-{\int }_{x_1}^{x(t)}\frac{x_1}{\theta }\mathrm{d}\theta +w(t)-w_1-w_1\ln \frac{w(t)}{w_1}\\ & +\frac{\delta }{2(d_1+d_2+q)x_1}[(x(t)-x_1)+(w(t)-w_1){]}^2\\ & +\frac{\delta +d_2+q}{q}(y(t)-y_1-y_1\ln \frac{y(t)}{y_1})\\ & +\frac{(\delta +d_2+q)h}{qk}z(t)+\frac{(\delta +d_2+q)h}{q}{\int }_{t-\tau }^{t}y(\theta )z(\theta )\mathrm{d}\theta ,\end{array}$where

$s=d_1{x}_1+\frac{\beta x_1{y}_1}{1+y_1}-\delta w_1$, $n_1=\frac{\beta x_1{y}_1}{1+y_1}=(\delta +d_2+q)w_1$, $1-\frac{x_1}{x}=\frac{-(x-x_1{)}^2}{x{x}_1}+\frac{x-x_1}{x_1}$.

Calculating the derivative of $V_1$ along the solution of the system (2(2) $\{\begin{array}{l}{\displaystyle \frac{\mathrm{d}x(t)}{\mathrm{d}t}=s-d_{1}x(t)-\frac{\beta x(t)y(t)}{1+y(t)}+\delta w(t),}\\ {\displaystyle \frac{\mathrm{d}w(t)}{\mathrm{d}t}=\frac{\beta x(t)y(t)}{1+y(t)}-(\delta +d_2+q)w(t),}\\ {\displaystyle \frac{\mathrm{d}y(t)}{\mathrm{d}t}}=qw(t)-d_{3}y(t)-hy(t)z(t),\\ {\displaystyle \frac{\mathrm{d}z(t)}{\mathrm{d}t}}=ky(t-\tau )z(t-\tau )-d_{4}z(t),\end{array}$(2) ), we obtain: $\begin{array}{rl}{V}_1^{\mathrm{\prime }}& =(d_1{x}_1+\frac{\beta x_1{y}_1}{1+y_1}-\delta w_1-d_{1}x(t)-\frac{\beta x(t)y(t)}{1+y(t)}+\delta w(t))(1-\frac{x_1}{x(t)})\\ & +[\frac{\beta x(t)y(t)}{1+y(t)}-(\delta +d_2+q)w(t)](1-\frac{w_1}{w(t)})\\ & +\frac{\delta }{(d_1+d_2+q)x_1}(x(t)-x_1+w(t)-w_1)[d_1{x}_1-\frac{\beta x_1{y}_1}{1+y_1}-\delta w_1-d_{1}x(t)\\ & +\delta w(t)-(\delta +d_2+q)w(t)]+\frac{\delta +d_2+q}{q}[qw(t)-d_{3}y(t)-hy(t)z(t)](1-\frac{y_1}{y(t)})\\ & +\frac{(\delta +d_2+q)h}{qk}[ky(t-\tau )z(t-\tau )-d_{4}z(t)]\\ & +\frac{(\delta +d_2+q)h}{qk}y(t)z(t)-\frac{(\delta +d_2+q)h}{qk}[ky(t-\tau )z(t-\tau )-d_{4}z(t)]\\ & =\frac{-d_1(x(t)-x_1)}{x(t)}+n_1[1-\frac{x_1}{x(t)}+\frac{y(t)(1+y_1)}{y_1(1+y(t))}]\\ & +\delta (w(t)-w_1)(\frac{-(x(t)-x_1{)}^2}{x(t)x_1}+\frac{x(t)-x_1}{x(t)})-\frac{(\delta +q+d_2)h}{qk}{d}_{4}z(t)\\ & -\frac{\delta d_1(x(t)-x_1{)}^2}{(d_1+d_2+q)x_1}-\frac{\delta (d_2+q)(w(t)-w_1{)}^2}{(d_1+d_2+q)x_1}-\delta (w(t)-w_1)\frac{x(t)-x_1}{x(t)}\\ & +n_1(1-\frac{x(t)y(t)(1+y_1)w_1}{x_1{y}_1(1+y(t)w(t))})+n_1(-\frac{d_{3}y(t)}{q{w}_1}-\frac{w(t)y_1}{w_{1}y(t)}+1+\frac{h{y}_{1}z(t)}{q{w}_1})\\ & =-[(d_1{x}_1-\delta w_1)+\delta w(t)+\frac{\delta d_{1}x(t)}{d_1+d_2+q}]\frac{(x(t)-x_1{)}^2}{x(t)x_1}-\frac{\delta (d_2+q)(w(t)-w_1{)}^2}{(d_1+d_2+q)x_1}\\ & +n_1(4-\frac{x_1}{x(t)}-\frac{x(t)y(t)(1+y_1)w_1}{x_1{y}_1(1+y(t)w(t))}-\frac{w(t)y_1}{w_{1}y(t)}-\frac{1+y(t)}{1+y_1})\\ & +n_1\frac{-(y(t)-y_1{)}^2}{y_1(1+y(t))(1+y_1)}\\ & +\frac{\beta x_1{y}_{1}hz(t)}{q{w}_1}\cdot \frac{d_1(k+d_4)(\delta +d_2+q)+d_4\beta (d_2+q)}{k\beta (d_2+q)+k{d}_1(\delta +d_2+q)}(R_1-1),\end{array}$where $d_1{x}_1-\delta w_1\ge 0$ can be formulated as $R_0\le 1+\frac{sq+(d_2+q)d_3}{\delta d_3}$. Since the arithmetic mean is greater than or equal to the geometric mean, it follows that $4-\frac{x_1}{x(t)}-\frac{x(t)y(t)(1+y_1)w_1}{x_1{y}_1(1+y(t)w(t))}-\frac{w(t)y_1}{w_{1}y(t)}-\frac{1+y(t)}{1+y_1}\le 0.$If $R_0\in (1,1+\frac{sq+(d_2+q)d_3}{\delta d_3})$ and $R_1<1$, we can get the above $V_1^{\mathrm{\prime }}\le 0$. In addition, if and only if $(x(t),w(t),y(t),z(t))=(x_1,w_1,y_1,0)$, we obtain $V_1^{\mathrm{\prime }}=0$. According to Lyapunov–LaSalle, we can know that the infected equilibrium $E_1$ of system (2(2) $\{\begin{array}{l}{\displaystyle \frac{\mathrm{d}x(t)}{\mathrm{d}t}=s-d_{1}x(t)-\frac{\beta x(t)y(t)}{1+y(t)}+\delta w(t),}\\ {\displaystyle \frac{\mathrm{d}w(t)}{\mathrm{d}t}=\frac{\beta x(t)y(t)}{1+y(t)}-(\delta +d_2+q)w(t),}\\ {\displaystyle \frac{\mathrm{d}y(t)}{\mathrm{d}t}}=qw(t)-d_{3}y(t)-hy(t)z(t),\\ {\displaystyle \frac{\mathrm{d}z(t)}{\mathrm{d}t}}=ky(t-\tau )z(t-\tau )-d_{4}z(t),\end{array}$(2) ) without immune is globally asymptotically stable. The proof is complete.

4. Stability analysis and the existence of Hopf bifurcation of equilibrium point $E_2$

In this section, we mainly discuss the stability and the existence of Hopf bifurcation of the infected equilibrium $E_2$ of system (2(2) $\{\begin{array}{l}{\displaystyle \frac{\mathrm{d}x(t)}{\mathrm{d}t}=s-d_{1}x(t)-\frac{\beta x(t)y(t)}{1+y(t)}+\delta w(t),}\\ {\displaystyle \frac{\mathrm{d}w(t)}{\mathrm{d}t}=\frac{\beta x(t)y(t)}{1+y(t)}-(\delta +d_2+q)w(t),}\\ {\displaystyle \frac{\mathrm{d}y(t)}{\mathrm{d}t}}=qw(t)-d_{3}y(t)-hy(t)z(t),\\ {\displaystyle \frac{\mathrm{d}z(t)}{\mathrm{d}t}}=ky(t-\tau )z(t-\tau )-d_{4}z(t),\end{array}$(2) ) with immune response.

Theorem 4.1

If $R_1>1$ and (H1) hold, the infected equilibrium $E_2$ of system (2(2) $\{\begin{array}{l}{\displaystyle \frac{\mathrm{d}x(t)}{\mathrm{d}t}=s-d_{1}x(t)-\frac{\beta x(t)y(t)}{1+y(t)}+\delta w(t),}\\ {\displaystyle \frac{\mathrm{d}w(t)}{\mathrm{d}t}=\frac{\beta x(t)y(t)}{1+y(t)}-(\delta +d_2+q)w(t),}\\ {\displaystyle \frac{\mathrm{d}y(t)}{\mathrm{d}t}}=qw(t)-d_{3}y(t)-hy(t)z(t),\\ {\displaystyle \frac{\mathrm{d}z(t)}{\mathrm{d}t}}=ky(t-\tau )z(t-\tau )-d_{4}z(t),\end{array}$(2) ) with immune response is globally stable when $\tau =0$. (H1) $\delta d_4\beta \le d_1(k+d_4)(\delta +d_2+q)$(H1)

Proof.

We define the Lyapunov function $V_2$ as follows: $\begin{array}{rl}{V}_2& =x(t)-x_2-{\int }_{x_2}^{x(t)}\frac{x_2}{\theta }\mathrm{d}\theta +w(t)-w_2-w_2\ln \frac{w(t)}{w_2}\\ & +\frac{\delta }{2(d_1+d_2+q)x_2}[(x(t)-x_2)+(w(t)-w_2){]}^2\\ & +\frac{\delta +d_2+q}{q}(y(t)-y_2-y_2\ln \frac{y(t)}{y_2})\\ & +\frac{(\delta +d_2+q)h}{qk}(z(t)-z_2-z_2\ln \frac{z(t)}{z_2}).\end{array}$Calculating the derivative of $V_2$ along the solution of the system (2(2) $\{\begin{array}{l}{\displaystyle \frac{\mathrm{d}x(t)}{\mathrm{d}t}=s-d_{1}x(t)-\frac{\beta x(t)y(t)}{1+y(t)}+\delta w(t),}\\ {\displaystyle \frac{\mathrm{d}w(t)}{\mathrm{d}t}=\frac{\beta x(t)y(t)}{1+y(t)}-(\delta +d_2+q)w(t),}\\ {\displaystyle \frac{\mathrm{d}y(t)}{\mathrm{d}t}}=qw(t)-d_{3}y(t)-hy(t)z(t),\\ {\displaystyle \frac{\mathrm{d}z(t)}{\mathrm{d}t}}=ky(t-\tau )z(t-\tau )-d_{4}z(t),\end{array}$(2) ), we obtain $\begin{array}{rl}{V}_2^{\mathrm{\prime }}& =(s-d_{1}x(t)-\frac{\beta x(t)y(t)}{1+y(t)}+\delta w(t))(1-\frac{x_2}{x(t)})\\ & +[\frac{\beta x(t)y(t)}{1+y(t)}-(\delta +d_2+q)w(t)](1-\frac{w_2}{w(t)})\\ & +\frac{\delta }{(d_1+d_2+q)x_2}[(x(t)-x_2)\\ & +(w(t)-w_2)][s-d_{1}x(t)+\delta w(t)-(\delta +d_2+q)w(t)]\\ & +\frac{\delta +d_2+q}{q}(qw(t)-d_{3}y(t)-hy(t)z(t))(1-\frac{y_2}{y(t)})\\ & +\frac{(\delta +d_2+q)h}{qk}(ky(t)z(t)-d_{4}z(t))(1-\frac{z_2}{z(t)})\end{array}$where $\begin{array}{rl}s& =d_1{x}_2+\frac{\beta x_2{y}_2}{1+y_2}-\delta w_2,d_4=k{y}_2,h{y}_2{z}_2=q{w}_2-d_3{y}_2,\\ n_2& =\frac{\beta x_2{y}_1}{1+y_2}=(\delta +d_2+q)w_2,1-\frac{x_2}{x}=\frac{-(x-x_2{)}^2}{x{x}_2}+\frac{x-x_2}{x_2}.\\ V_2^{\mathrm{\prime }}& =-[(d_1{x}_2-\delta w_2)+\delta w(t)+\frac{d_1\delta x(t)}{d_1+d_2+q}]\frac{(x(t)-x_2{)}^2}{x(t)x_2}\\ & -\frac{\delta (d_2+q)}{(d_1+d_2+q)x(t)}(w(t)-w_2{)}^2\\ & +n_2[4-\frac{x_2}{x(t)}-\frac{(1+y_2)x(t)y(t)w_2}{(1+y(t))x_2{y}_{2}w(t)}-\frac{w(t)y_2}{w_{2}y(t)}]-n_2\frac{(y(t)-y_2{)}^2}{y_2(1+y(t))(1+y_2)},\end{array}$where $d_1{x}_2-\delta w_2\ge 0$ can be formulated as (H1), we can get the above $V_2^{\mathrm{\prime }}\le 0$, In addition, if and only if $(x(t),w(t),y(t),z(t))=(x_2,w_2,y_2,z_2)$ we obtain $V_2^{\mathrm{\prime }}=0$. The infected equilibrium $E_2$ of system (2(2) $\{\begin{array}{l}{\displaystyle \frac{\mathrm{d}x(t)}{\mathrm{d}t}=s-d_{1}x(t)-\frac{\beta x(t)y(t)}{1+y(t)}+\delta w(t),}\\ {\displaystyle \frac{\mathrm{d}w(t)}{\mathrm{d}t}=\frac{\beta x(t)y(t)}{1+y(t)}-(\delta +d_2+q)w(t),}\\ {\displaystyle \frac{\mathrm{d}y(t)}{\mathrm{d}t}}=qw(t)-d_{3}y(t)-hy(t)z(t),\\ {\displaystyle \frac{\mathrm{d}z(t)}{\mathrm{d}t}}=ky(t-\tau )z(t-\tau )-d_{4}z(t),\end{array}$(2) ) with immune is globally asymptotically stable from Lyapunov–LaSalle.

Next, when $\tau >0$, we linearize system (2(2) $\{\begin{array}{l}{\displaystyle \frac{\mathrm{d}x(t)}{\mathrm{d}t}=s-d_{1}x(t)-\frac{\beta x(t)y(t)}{1+y(t)}+\delta w(t),}\\ {\displaystyle \frac{\mathrm{d}w(t)}{\mathrm{d}t}=\frac{\beta x(t)y(t)}{1+y(t)}-(\delta +d_2+q)w(t),}\\ {\displaystyle \frac{\mathrm{d}y(t)}{\mathrm{d}t}}=qw(t)-d_{3}y(t)-hy(t)z(t),\\ {\displaystyle \frac{\mathrm{d}z(t)}{\mathrm{d}t}}=ky(t-\tau )z(t-\tau )-d_{4}z(t),\end{array}$(2) ) at $E_2$ to obtain (4) $\{\begin{array}{l}{\displaystyle \frac{\mathrm{d}x}{\mathrm{d}t}=(-d_1-\frac{\beta y_2}{1+y_2})x(t)+\delta w(t)-\frac{\beta x_2}{(1+y_2{)}^2}y(t),}\\ {\displaystyle \frac{\mathrm{d}w}{\mathrm{d}t}=\frac{\beta y_2}{1+y_2}x(t)-(\delta +d_2+q)w(t)+\frac{\beta x_2}{(1+y_2{)}^2}y(t),}\\ {\displaystyle \frac{\mathrm{d}y}{\mathrm{d}t}=qw(t)-(d_3+h{z}_2)y(t)-h{y}_{2}z(t),}\\ {\displaystyle \frac{\mathrm{d}z}{\mathrm{d}t}=k{z}_{2}y(t-\tau )+k{y}_{2}z(t-\tau )-d_{4}z(t).}\end{array}$(4) The associated characteristic equation of system (4(4) $\{\begin{array}{l}{\displaystyle \frac{\mathrm{d}x}{\mathrm{d}t}=(-d_1-\frac{\beta y_2}{1+y_2})x(t)+\delta w(t)-\frac{\beta x_2}{(1+y_2{)}^2}y(t),}\\ {\displaystyle \frac{\mathrm{d}w}{\mathrm{d}t}=\frac{\beta y_2}{1+y_2}x(t)-(\delta +d_2+q)w(t)+\frac{\beta x_2}{(1+y_2{)}^2}y(t),}\\ {\displaystyle \frac{\mathrm{d}y}{\mathrm{d}t}=qw(t)-(d_3+h{z}_2)y(t)-h{y}_{2}z(t),}\\ {\displaystyle \frac{\mathrm{d}z}{\mathrm{d}t}=k{z}_{2}y(t-\tau )+k{y}_{2}z(t-\tau )-d_{4}z(t).}\end{array}$(4) ) at $E_2$ becomes (5) $H(\lambda ;\tau )={\lambda }^4+b_1{\lambda }^3+b_2{\lambda }^2+b_3\lambda +b_4+(c_1{\lambda }^3+c_2{\lambda }^3+c_3\lambda +c_4){\mathrm{e}}^{-\lambda \tau }=0,$(5) where $\begin{array}{rl}{b}_1& =A+B+C+d_4,\\ b_2& =AB+BC+(d_2+q)C+\delta d_1+d_4(A+B+C)-qD,\\ b_3& =ABC+\delta d_{1}B-\delta BC-d_{1}qD+d_4(AB+BC+AC+\delta d_1-\delta C-qD),\\ b_4& =d_4(ABC+\delta d_{1}B-\delta BC-d_{1}qD),\\ c_1& =-k{y}_2,\\ c_2& =hk{y}_2{z}_2-(A+B+C)k{y}_2,\\ c_3& =Ahk{y}_2{z}_2+Chk{y}_2{z}_2-(AB+BC+AC+\delta d_1-\delta C-qD)k{y}_2,\\ c_4& =hk{y}_2{z}_2(AC+\delta d_1-\delta C)-(ABC+\delta d_{1}B-\delta BC-d_{1}qD)k{y}_2,\\ A& =\delta +d_2+q,B=d_3+h{z}_2,C=d_1+\frac{\beta y_2}{1+y_2},D=\beta X_2(1+Y_2{)}^2.\end{array}$We suppose (5(5) $H(\lambda ;\tau )={\lambda }^4+b_1{\lambda }^3+b_2{\lambda }^2+b_3\lambda +b_4+(c_1{\lambda }^3+c_2{\lambda }^3+c_3\lambda +c_4){\mathrm{e}}^{-\lambda \tau }=0,$(5) ) has a purely imaginary root $\lambda =i\omega$, then we obtain $\begin{array}{rl}& {\omega }^4-i{b}_1{\omega }^3-b_2{\omega }^2+i{b}_3\omega +b_4\\ & +[\cos (\omega \tau )-i\sin (\omega \tau )](i{c}_1{\omega }^3-c_2{\omega }^2+i{c}_3\omega +c_3)=0.\end{array}$Separating the real parts and imaginary parts of the above equation, we can get (6) $\{\begin{array}{l}{\omega }^4-b_2{\omega }^2+b_4=-\cos (\omega \tau )(-c_2{\omega }^2+c_4)+\sin (\omega \tau )(c_1{\omega }^3+c_3\omega ),\\ -b_1{\omega }^3+b_3\omega =cos(\omega \tau )(c_1{\omega }^3-c_3\omega )+\sin (\omega \tau )(-c_2{\omega }^2+c_4).\end{array}$(6) Then we have (7) ${\omega }^8+l_1{\omega }^6+l_2{\omega }^4+l_3{\omega }^2+l_4=0,$(7) where $\begin{array}{rl}{l}_1& =b_1^2-2b_2-c_1^2,\\ l_2& =b_2^2+2b_4-2b_1{b}_3-c_2^2+2c_1{c}_3,\\ l_3& =b_3^2-2b_2{b}_4-c_3^2+2c_2{c}_4,\\ l_4& =h{d}_4^2{z}_2(d_{2}C+qC+\delta d_1)[(d_3+B)(d_{2}C+qC+\delta d_1)-2d_{1}qD].\end{array}$Denote (8) $G(\varpi )={\varpi }^4+l_1{\varpi }^3+l_2{\varpi }^2+l_3\varpi +l_4.$(8) If Equation (5(5) $H(\lambda ;\tau )={\lambda }^4+b_1{\lambda }^3+b_2{\lambda }^2+b_3\lambda +b_4+(c_1{\lambda }^3+c_2{\lambda }^3+c_3\lambda +c_4){\mathrm{e}}^{-\lambda \tau }=0,$(5) ) has a purely imaginary root $i\omega$, equation (9) $G(\varpi )={\varpi }^4+l_1{\varpi }^3+l_2{\varpi }^2+l_3\varpi +l_4=0.$(9) will have a positive real root ${\omega }^2$.

If $l_4<0$, we can obtain the following inequality: (H2) $(d_3+B)(d_{2}C+qD+\delta d_1)<2d_{1}qD.$(H2) The above formulation implies that Equation (9(9) $G(\varpi )={\varpi }^4+l_1{\varpi }^3+l_2{\varpi }^2+l_3\varpi +l_4=0.$(9) ) has one positive real root at least.

Suppose that Equation (9(9) $G(\varpi )={\varpi }^4+l_1{\varpi }^3+l_2{\varpi }^2+l_3\varpi +l_4=0.$(9) ) has $n(1\le n\le 4)$ positive real roots, then Equation (7(7) ${\omega }^8+l_1{\omega }^6+l_2{\omega }^4+l_3{\omega }^2+l_4=0,$(7) ) has n positive real roots ${\omega }_1=\sqrt{{\varpi }_1},{\omega }_2=\sqrt{{\varpi }_2}\cdots {\omega }_n=\sqrt{{\varpi }_n}(1\le n\le 4)$. Through Equation (6(6) $\{\begin{array}{l}{\omega }^4-b_2{\omega }^2+b_4=-\cos (\omega \tau )(-c_2{\omega }^2+c_4)+\sin (\omega \tau )(c_1{\omega }^3+c_3\omega ),\\ -b_1{\omega }^3+b_3\omega =cos(\omega \tau )(c_1{\omega }^3-c_3\omega )+\sin (\omega \tau )(-c_2{\omega }^2+c_4).\end{array}$(6) ), we get $\{\begin{array}{l}\sin (\omega \tau )=F_n={\displaystyle \frac{(c_1{\omega }^3-c_3\omega )({\omega }^4-b_2{\omega }^2+b_4)+(-c_1{\omega }^2+c_4)(-b_1{\omega }^3+b_3\omega )}{(c_1{\omega }^3-c_3\omega {)}^2+(-c_2{\omega }^2+c_4{)}^2}},\\ \cos (\omega \tau )=J_n={\displaystyle \frac{(c_1{\omega }^3-c_3\omega )(-b_1{\omega }^3+b_3\omega )-(-c_2{\omega }^2+c_4)({\omega }^4-b_2{\omega }^2+b_4)}{(c_1{\omega }^3-c_3\omega {)}^2+(-c_2{\omega }^2+c_4{)}^2}}.\end{array}$Then, we have (10) ${\tau }_n^{(j)}=\frac{1}{{\omega }_n}\arccos (J_n)+\frac{2\pi j}{{\omega }_n}(1\le n\le 4,j=0,1,2,3,\dots ).$(10) It is easy to show that $\pm i{\omega }_n$ is a pair of purely imaginary root of Equation (5(5) $H(\lambda ;\tau )={\lambda }^4+b_1{\lambda }^3+b_2{\lambda }^2+b_3\lambda +b_4+(c_1{\lambda }^3+c_2{\lambda }^3+c_3\lambda +c_4){\mathrm{e}}^{-\lambda \tau }=0,$(5) ), for every integer j and n, let ${\lambda }_n^{(j)}(\tau )={\alpha }_n^{(j)}(\tau )+i{\omega }_n^{(j)}(\tau )$ be the roots of (5(5) $H(\lambda ;\tau )={\lambda }^4+b_1{\lambda }^3+b_2{\lambda }^2+b_3\lambda +b_4+(c_1{\lambda }^3+c_2{\lambda }^3+c_3\lambda +c_4){\mathrm{e}}^{-\lambda \tau }=0,$(5) ) near $\tau ={\tau }_n^{(j)}$ satisfying ${\alpha }_n^{(j)}=0,{\omega }_n^{(j)}(\tau )={\omega }_n$. Then, we have the following theorem.

Theorem 4.2

The $\frac{\mathrm{d}Re(\lambda )}{\mathrm{d}\tau }{\mid }_{\tau ={\tau }_n^{(j)}}$ and $G^{\mathrm{\prime }}({\omega }_n^2)$ have the same sign.

Proof.

Put ${\lambda }_n^{(j)}(\tau )$ into Equation (5(5) $H(\lambda ;\tau )={\lambda }^4+b_1{\lambda }^3+b_2{\lambda }^2+b_3\lambda +b_4+(c_1{\lambda }^3+c_2{\lambda }^3+c_3\lambda +c_4){\mathrm{e}}^{-\lambda \tau }=0,$(5) ) we get (11) $p(\lambda )+f(\lambda ){\mathrm{e}}^{-\lambda \tau }=0,$(11) where $\begin{array}{rl}p(\lambda )& ={\lambda }^4+b_1{\lambda }^3+b_2{\lambda }^2+b_3\lambda +b_4,\\ f(\lambda )& =c_1{\lambda }^3+c_2{\lambda }^2+c_3\lambda +c_4.\end{array}$Differentiating (11(11) $p(\lambda )+f(\lambda ){\mathrm{e}}^{-\lambda \tau }=0,$(11) ) with respect to τ we obtain that $p^{\mathrm{\prime }}(\lambda )\frac{\mathrm{d}\lambda }{\mathrm{d}\tau }+f^{\mathrm{\prime }}(\lambda )\frac{\mathrm{d}\lambda }{\mathrm{d}\tau }{\mathrm{e}}^{-\lambda \tau }-(\lambda +\tau \frac{\mathrm{d}\lambda }{\mathrm{d}\tau })f(\lambda ){\mathrm{e}}^{-\lambda \tau }=0.$Hence, we get ${\left[\frac{\mathrm{d}\lambda }{\mathrm{d}\tau }\right]}^{-1}=\frac{f^{\mathrm{\prime }}(\lambda )+p^{\mathrm{\prime }}(\lambda ){\mathrm{e}}^{\lambda \tau }}{\lambda f(\lambda )}-\frac{\tau }{\lambda }.$But $p(i{\omega }_n)+f(i{\omega }_n){\mathrm{e}}^{-i{\omega }_n{\tau }_n^{(j)}}=0$, we have $\begin{array}{rl}Re{\left[\frac{d\lambda }{d\tau }{|}_{\tau ={\tau }_n}\right]}^{-1}& =Re\left[\frac{f^{\mathrm{\prime }}(i{\omega }_n)+p^{\mathrm{\prime }}(i{\omega }_n){\mathrm{e}}^{i{\omega }_n{\tau }_n^j}}{i{\omega }_{n}f(i{\omega }_n)}\right]\\ & =Re\left[\frac{-f^{\mathrm{\prime }}(i{\omega }_n)}{{\omega }_{n}f(i{\omega }_n)}i\right]+Re\left[\frac{p^{\mathrm{\prime }}(i{\omega }_n)}{{\omega }_{n}p(i{\omega }_n)}i\right]\\ & =Im[\frac{f^{\mathrm{\prime }}(i{\omega }_n)}{{\omega }_{n}f(i{\omega }_n)}-\frac{p^{\mathrm{\prime }}(i{\omega }_n)}{{\omega }_{n}p(i{\omega }_n)}].\end{array}$On the other hand, we define $\phi (\omega )=|p(i\omega ){|}^2-|f(i\omega ){|}^2$. By calculating we can know that $\phi (\omega )=G(\omega {)}^2$. Calculating the derivative of $|p(i\omega ){|}^2$ with respect to ω, we obtain $\begin{array}{rl}\frac{\mathrm{d}}{\mathrm{d}\omega }(|p(i\omega ){|}^2)& =\frac{\mathrm{d}}{\mathrm{d}\omega }\{[Rep(i\omega {)}^2]-[Imp(i\omega ){]}^2\}\\ & =2Rep(i\omega )\cdot Re[p^{\mathrm{\prime }}(i\omega )i]+2Imp(i\omega )\cdot Im[p^{\mathrm{\prime }}(i\omega )i]\\ & =2Re[p(i\omega )p^{\mathrm{\prime }}(i\omega )i]\\ & =-2Im[\overline{p(i\omega )}{p}^{\mathrm{\prime }}(i\omega )].\end{array}$Then $\begin{array}{rl}\frac{1}{2\omega }\frac{\mathrm{d}\phi }{\mathrm{d}\omega }& =\frac{1}{2\omega }\frac{\mathrm{d}}{\mathrm{d}\omega }(|p(i\omega ){|}^2-|f(i\omega ){|}^2)\\ & =\frac{1}{\omega }\cdot Im[\overline{f(i\omega )}{f}^{\mathrm{\prime }}(i\omega )-\overline{p(i\omega )}{p}^{\mathrm{\prime }}(i\omega )]\\ & =Im[|f(i\omega ){|}^2\frac{f^{\mathrm{\prime }}(i\omega )}{\omega f(i\omega )}-|p(i\omega ){|}^2\frac{p^{\mathrm{\prime }}(i\omega )}{\omega p(i\omega )}].\end{array}$Since $|f(i{\omega }_n){|}^2=|p(i{\omega }_n){|}^2$, we get $\left(\frac{1}{2\omega }\frac{\mathrm{d}\phi }{\mathrm{d}\omega }\right){|}_{\omega ={\omega }_n}=|f(i{\omega }_n){|}^{2}Im[\frac{f^{\mathrm{\prime }}(i{\omega }_n)}{{\omega }_{n}f(i{\omega }_n)}-\frac{p^{\mathrm{\prime }}(i{\omega }_n)}{{\omega }_{n}p(i{\omega }_n)}]=|f(i{\omega }_n){|}^{2}Re{\left[\frac{\mathrm{d}\lambda }{\mathrm{d}\tau }{|}_{\tau ={\tau }_n}\right]}^{-1}.$Then $sign[G^{\mathrm{\prime }}({\omega }_n^2)]=sign\left[\left(\frac{1}{2\omega }\frac{\mathrm{d}\phi }{\mathrm{d}\omega }\right){|}_{\omega ={\omega }_n}\right]=signRe{\left[\frac{\mathrm{d}\lambda }{\mathrm{d}\tau }{|}_{\tau ={\tau }_n}\right]}^{-1}.$Because $sign\left[\frac{\mathrm{d}Re(\lambda )}{\mathrm{d}\tau }{|}_{\tau ={\tau }_n}\right]=signRe\left[\frac{\mathrm{d}\lambda }{\mathrm{d}\tau }{|}_{\tau ={\tau }_n}\right]=signRe{\left[\frac{\mathrm{d}\lambda }{\mathrm{d}\tau }{|}_{\tau ={\tau }_n}\right]}^{-1},$Therefore, we have $sign\left[\frac{\mathrm{d}Re(\lambda )}{\mathrm{d}\tau }{|}_{\tau ={\tau }_n}\right]=sign[G^{\mathrm{\prime }}({\omega }_n^2)].$It is obvious that if $G^{\mathrm{\prime }}({\omega }_n^2)\ne 0$, then $\frac{dRe(\lambda )}{d\tau }{|}_{\tau ={\tau }_n}\ne 0$. So, according to the above analysis and Hopf Bifurcation theorems given in the literature [6], we have the following conclusions.

Theorem 4.3

Assume that $R_1>1$, (H1) and (H2) hold. Then

1. the infected equilibrium with immune response $E_2$ $\tau \in [0,{\tau }_0)$, where ${\tau }_0={\tau }_{n0}^{(0)}=min\{{\tau }_n^j|1\le n\le 4,j=0,1,2,3,\dots \}$, ${\omega }_0={\omega }_{n0}$. (${\tau }_n^j$ is defined by (10(10) ${\tau }_n^{(j)}=\frac{1}{{\omega }_n}\arccos (J_n)+\frac{2\pi j}{{\omega }_n}(1\le n\le 4,j=0,1,2,3,\dots ).$(10) ))

2. if $G^{\mathrm{\prime }}({\omega }_{N0}^2)\ne 0$, there is a Hopf bifurcation for the system (2(2) $\{\begin{array}{l}{\displaystyle \frac{\mathrm{d}x(t)}{\mathrm{d}t}=s-d_{1}x(t)-\frac{\beta x(t)y(t)}{1+y(t)}+\delta w(t),}\\ {\displaystyle \frac{\mathrm{d}w(t)}{\mathrm{d}t}=\frac{\beta x(t)y(t)}{1+y(t)}-(\delta +d_2+q)w(t),}\\ {\displaystyle \frac{\mathrm{d}y(t)}{\mathrm{d}t}}=qw(t)-d_{3}y(t)-hy(t)z(t),\\ {\displaystyle \frac{\mathrm{d}z(t)}{\mathrm{d}t}}=ky(t-\tau )z(t-\tau )-d_{4}z(t),\end{array}$(2) ) near as τ is increased past ${\tau }_0$.

If (H1), (H2), (H3), (H4) hold, then (i) the positive equilibrium $E^{\ast }$ of system (1.7) is locally asymptotically stable for $0\le \tau <{\tau }_0$; (ii) $E^{\ast }$ is unstable for $\tau >{\tau }_0$; (iii) system (1.7) undergoes a Hopf bifurcation at $E^{\ast }$ for $\tau ={\tau }_0$.

5. Numerical simulations

In order to illustrate feasibility of the results of Theorem 4.3, we use the software Matlab to perform numerical simulations. Considering the following special system (12(12) $\{\begin{array}{l}{\displaystyle \frac{\mathrm{d}x(t)}{\mathrm{d}t}=7-0.03x(t)-0.011\frac{x(t)y(t)}{1+y(t)+0.001w(t)},}\\ {\displaystyle \frac{\mathrm{d}w(t)}{\mathrm{d}t}=0.011\frac{x(t)y(t)}{1+y(t)}-(0.001+0.21+0.89)w(t),}\\ {\displaystyle \frac{\mathrm{d}y(t)}{\mathrm{d}t}=0.89w(t)-0.26y(t)-0.02y(t)z(t)),}\\ {\displaystyle \frac{\mathrm{d}z(t)}{\mathrm{d}t}=0.62y(t-\tau )z(t-\tau )-0.25z(t).}\end{array}$(12) ) of system (2(2) $\{\begin{array}{l}{\displaystyle \frac{\mathrm{d}x(t)}{\mathrm{d}t}=s-d_{1}x(t)-\frac{\beta x(t)y(t)}{1+y(t)}+\delta w(t),}\\ {\displaystyle \frac{\mathrm{d}w(t)}{\mathrm{d}t}=\frac{\beta x(t)y(t)}{1+y(t)}-(\delta +d_2+q)w(t),}\\ {\displaystyle \frac{\mathrm{d}y(t)}{\mathrm{d}t}}=qw(t)-d_{3}y(t)-hy(t)z(t),\\ {\displaystyle \frac{\mathrm{d}z(t)}{\mathrm{d}t}}=ky(t-\tau )z(t-\tau )-d_{4}z(t),\end{array}$(2) ) : (12) $\{\begin{array}{l}{\displaystyle \frac{\mathrm{d}x(t)}{\mathrm{d}t}=7-0.03x(t)-0.011\frac{x(t)y(t)}{1+y(t)+0.001w(t)},}\\ {\displaystyle \frac{\mathrm{d}w(t)}{\mathrm{d}t}=0.011\frac{x(t)y(t)}{1+y(t)}-(0.001+0.21+0.89)w(t),}\\ {\displaystyle \frac{\mathrm{d}y(t)}{\mathrm{d}t}=0.89w(t)-0.26y(t)-0.02y(t)z(t)),}\\ {\displaystyle \frac{\mathrm{d}z(t)}{\mathrm{d}t}=0.62y(t-\tau )z(t-\tau )-0.25z(t).}\end{array}$(12) For the parameters from (12(12) $\{\begin{array}{l}{\displaystyle \frac{\mathrm{d}x(t)}{\mathrm{d}t}=7-0.03x(t)-0.011\frac{x(t)y(t)}{1+y(t)+0.001w(t)},}\\ {\displaystyle \frac{\mathrm{d}w(t)}{\mathrm{d}t}=0.011\frac{x(t)y(t)}{1+y(t)}-(0.001+0.21+0.89)w(t),}\\ {\displaystyle \frac{\mathrm{d}y(t)}{\mathrm{d}t}=0.89w(t)-0.26y(t)-0.02y(t)z(t)),}\\ {\displaystyle \frac{\mathrm{d}z(t)}{\mathrm{d}t}=0.62y(t-\tau )z(t-\tau )-0.25z(t).}\end{array}$(12) ), we can calculate $d_1(k+d_4)(\delta +d_2+q)-\delta d_4\beta =0.0287>0$ and ${\omega }_0=0.2313$ is a simple root of Equation (9(9) $G(\varpi )={\varpi }^4+l_1{\varpi }^3+l_2{\varpi }^2+l_3\varpi +l_4=0.$(9) ), so (H1), (H2) hold. We can also calculate that ${\tau }_0=4.1244,R_1=5.1452>1$ by using the software Matlab. If $\tau =2<{\tau }_0$, we can get Figure 1; If $\tau =6>{\tau }_0$, we can get Figure 2. From Figures 1 and 2, we can know that Theorem 4.3 holds.

Figure 1. Numerical simulations show that the equilibrium $E_2$ of system (12) is locally asymptotically stable when $\tau =2<{\tau }_0$ holds.

Figure 2. Numerical simulations show that the equilibrium $E_2$ of system (12) is unstable when $\tau =6>{\tau }_0$ holds.

6. Conclusion and prospects

In this paper, we established a mathematical model for HIV-1 with the immune delay and Holling II infection rate. In this model, we define the infection of the basic reproductive number $R_0$ and the basic immune response reproductive number $R_1$ by calculating and we identify the three equilibrium of the model. Then, it will show that uninfected equilibrium $E_0$ of system (2(2) $\{\begin{array}{l}{\displaystyle \frac{\mathrm{d}x(t)}{\mathrm{d}t}=s-d_{1}x(t)-\frac{\beta x(t)y(t)}{1+y(t)}+\delta w(t),}\\ {\displaystyle \frac{\mathrm{d}w(t)}{\mathrm{d}t}=\frac{\beta x(t)y(t)}{1+y(t)}-(\delta +d_2+q)w(t),}\\ {\displaystyle \frac{\mathrm{d}y(t)}{\mathrm{d}t}}=qw(t)-d_{3}y(t)-hy(t)z(t),\\ {\displaystyle \frac{\mathrm{d}z(t)}{\mathrm{d}t}}=ky(t-\tau )z(t-\tau )-d_{4}z(t),\end{array}$(2) ) is globally asymptotically stable for any time delay $\tau \ge 0$ by using the Lyapunov–LaSalle theorem when the infection of the basic reproductive number $R_0$; the infected equilibrium without immune response of system (2(2) $\{\begin{array}{l}{\displaystyle \frac{\mathrm{d}x(t)}{\mathrm{d}t}=s-d_{1}x(t)-\frac{\beta x(t)y(t)}{1+y(t)}+\delta w(t),}\\ {\displaystyle \frac{\mathrm{d}w(t)}{\mathrm{d}t}=\frac{\beta x(t)y(t)}{1+y(t)}-(\delta +d_2+q)w(t),}\\ {\displaystyle \frac{\mathrm{d}y(t)}{\mathrm{d}t}}=qw(t)-d_{3}y(t)-hy(t)z(t),\\ {\displaystyle \frac{\mathrm{d}z(t)}{\mathrm{d}t}}=ky(t-\tau )z(t-\tau )-d_{4}z(t),\end{array}$(2) ) is globally asymptotical stable for any time delay $\tau \ge 0$, when $R_0\in (1,1+\frac{sq+(d_2+q)d_3}{\delta d_3}]$ and $R_1<0$ hold; the infected equilibrium with immune response $E_2$ of system (2(2) $\{\begin{array}{l}{\displaystyle \frac{\mathrm{d}x(t)}{\mathrm{d}t}=s-d_{1}x(t)-\frac{\beta x(t)y(t)}{1+y(t)}+\delta w(t),}\\ {\displaystyle \frac{\mathrm{d}w(t)}{\mathrm{d}t}=\frac{\beta x(t)y(t)}{1+y(t)}-(\delta +d_2+q)w(t),}\\ {\displaystyle \frac{\mathrm{d}y(t)}{\mathrm{d}t}}=qw(t)-d_{3}y(t)-hy(t)z(t),\\ {\displaystyle \frac{\mathrm{d}z(t)}{\mathrm{d}t}}=ky(t-\tau )z(t-\tau )-d_{4}z(t),\end{array}$(2) ) is globally stable for $\tau =0$, when $R_1>1$ and $(H1)$ hold. And we also give sufficient condition of the Hopf bifurcation in equilibrium $E_2$ of system (2(2) $\{\begin{array}{l}{\displaystyle \frac{\mathrm{d}x(t)}{\mathrm{d}t}=s-d_{1}x(t)-\frac{\beta x(t)y(t)}{1+y(t)}+\delta w(t),}\\ {\displaystyle \frac{\mathrm{d}w(t)}{\mathrm{d}t}=\frac{\beta x(t)y(t)}{1+y(t)}-(\delta +d_2+q)w(t),}\\ {\displaystyle \frac{\mathrm{d}y(t)}{\mathrm{d}t}}=qw(t)-d_{3}y(t)-hy(t)z(t),\\ {\displaystyle \frac{\mathrm{d}z(t)}{\mathrm{d}t}}=ky(t-\tau )z(t-\tau )-d_{4}z(t),\end{array}$(2) ). Our analysis provides the effective reference value of the prevention and treatment of AIDS.

By using a novel reducing dimension modelling idea proposed in [29], see also [30, 31], we may reduce (2(2) $\{\begin{array}{l}{\displaystyle \frac{\mathrm{d}x(t)}{\mathrm{d}t}=s-d_{1}x(t)-\frac{\beta x(t)y(t)}{1+y(t)}+\delta w(t),}\\ {\displaystyle \frac{\mathrm{d}w(t)}{\mathrm{d}t}=\frac{\beta x(t)y(t)}{1+y(t)}-(\delta +d_2+q)w(t),}\\ {\displaystyle \frac{\mathrm{d}y(t)}{\mathrm{d}t}}=qw(t)-d_{3}y(t)-hy(t)z(t),\\ {\displaystyle \frac{\mathrm{d}z(t)}{\mathrm{d}t}}=ky(t-\tau )z(t-\tau )-d_{4}z(t),\end{array}$(2) ) into a three dimension system if we regard $x(t)$ or $z(t)$ as a known function rather than an independent variable satisfying an independent dynamical equation, which will be a challenging topic to the future study.

Disclosure statement

No potential conflict of interest was reported by the author(s).

Additional information

Funding

This work is supported partly by Hunan Provincial Natural Science Foundation of China (No. 2020JJ4516), Fund of Education Department of Hunan Province (No. 17A181), Hunan Key Laboratory of Mathematical Modeling and Scientific Computing (No. 2019).