Bifurcation and Stability Analysis of HIV Infectious Model with Two Time Delays

S. Q. Ma

Department of Mathematics, China Agricultural University, Beijing, China.

DOI: 10.4236/ijmnta.2021.102004                

Abstract

The HIV problem is studied by version of delay mathematical models which consider the apoptosis of uninfected CD4⁺ T cells which cultured with infected T cells in big volume. The opportunistic infection and the apoptosis of uninfected CD4⁺ T cells are caused directly or indirectly by a toxic substance produced from HIV genes. Ubiquitously, the nonlinear incidence rate brings forth the increasing number of infected CD4⁺ T cells with introduction of small time delay, and in addition, there also exists a natural time delay factor during the process of virus replication. With state feedback control of time delay, the bifurcating periodical oscillating phenomena is induced via Hopf bifurcation. Mathematically, with the geometrical criterion applied in the stability analysis of delay model, the critical threshold of Hopf bifurcation in multiple delay differential equations which satisfy the transversal condition is derived. By applying reduction dimensional method combined with the center manifold theory, the stability of the bifurcating periodical solution is analyzed by the perturbation near Hopf point.

Keywords

Hopf Bifurcation, Multiple Time Delay, HIV Model

1. Introduction

The serious infectious problems of the human immunodeficiency virus (HIV) which can attack human immune system to infect human body healthy cells have drawn more attentions in the fields of dynamical investigation. Through its slowly replicating retrovirus process to cause the acquired immunodeficiency syndrom (AIDs) of human health problem, the body becomes gradually very susceptible to opportunistic infections while CD4⁺ T cells fall below a critical level or oscillate for a long period with the lower level [1] [2] [3]. The most devastating thing is that people have made big efforts in mathematical models to give a deep understanding of the inherent dynamical mechanism of HIV viral infection [4] [5]. Especially, the classic mathematical model of HIV infection dynamics which is originally modeled by Anderson and May and developed subsequently by Nowak and Bangham is well accepted as the basic model(ODEs) to be investigated both by theorists and experimentalists [6] [7] [8].

As is well known, HIV gene expression products can produce toxic substances which indirectly or directly lead to apoptosis of uninfected CD4⁺ T cells [1] [2] [8]. It is verified in a laboratory experiment that the apoptosis may occur without virus replication while both uninfected cells and infected cells are cultured together [9] [10]. It is suggested that viral proteins associate with uninfected CD4⁺ T cells can induce an apoptotic signal which induces the death of uninfected CD4⁺ T cells [11] [12]. People also have made consideration of the interaction relationship acts between infected cells and uninfected T cells via Holling type III functional response, Beddington-Deangelis functional response, and bilinear infection rate or more general nonlinear infection rate as suggested in papers [7] [13] [14] [15].

Accordingly, based on the biological meanings, we develop the following HIV infection model with a viral infection and replication kinetics

$\begin{array}{l}{x}^{\prime }\left(t\right)=s-dx(t)-c\frac{x\left(t\right)}{\theta +x\left(t\right)}y\left(t\right),\\ y^{\prime }\left(t\right)=\alpha c\frac{x\left(t-{\tau }_1\right)}{\theta +x\left(t-{\tau }_1\right)}y\left(t-{\tau }_1\right)+n{\text{e}}^{-\mu {\tau }_2}v\left(t-{\tau }_2\right)-\left(d+k\right)y\left(t\right),\\ v^{\prime }\left(t\right)=ky\left(t\right)-mv\left(t\right)\end{array}$ (1.1)

where $x\left(t\right)\mathrm{,}y\left(t\right)$ denote uninfected and infected CD4⁺T cells respectively, whilst $v\left(t\right)$ represents the concentration of virus. The explainition of system (1.1) biologically meaningful lies that virus replication is fractionally positive to the increasing number of infected T cells. In addition, the bilinear incidence rate of uninfected T cells and infected cells induces the increasing number of infected T cells with time period ${\tau }_1$. The virus replication exaggerates the number of infected T cells through virus replicating time ${\tau }_2$ though the apoptosis of virus is reduced with time ${\tau }_2$. In comparison with other HIV infection model, system (1.1) is delay differential equations(DDEs) which contains double time delay ${\tau }_1$ and ${\tau }_2$. With multiple time delay, the dynamics of system is complicated since instability produces period oscillation phenomena. With the attempt to track the period varying and bifurcation phenomena of the limit cycle as varying parameters continuously, the DDE-Biftool is a technology tool of art in the dynamical investigation field [16] [17].

The extending geometrical criterion to work out the eigenvalue problem of DDEs with multiple time delay is recently given in papers [18] [19] [20]. Hopf lines are tracked by varying time delay in parameter space and the transversal condition is also tidily outlined. The state of instability switching phenomenon happening is associated with the characteristic roots with zero real parts appearing in complex plane and limit cycle always arise near the threshold value. With time delay, the quanititive property of the evolution dynamical behavior of system (1.1) is investigated with the estimated value of loss rate k varying.

Periodical oscillation phenomena in every respect of T cells population arise as system equilibrium loses its stability. It is verified that T cells population oscillating in a time period with delayed feedback control of state variable. It is explained that difference estimation of T healthy cells of past history and present time may induce periodically evolutional behavior of system states. In another respect, in some extent, the difference estimation is regarded as disease symptom proof. Therefore, we introduce the delayed feedback control in system (1.1) as the following

$\begin{array}{l}{x}^{\prime }\left(t\right)=s-dx\left(t\right)-c\frac{x\left(t\right)}{\theta +x\left(t\right)}y\left(t\right)+r\left(x\left(t-{\tau }_1\right)-x\left(t\right)\right),\\ y^{\prime }\left(t\right)=\alpha c\frac{x\left(t-{\tau }_1\right)}{\theta +x\left(t-{\tau }_1\right)}y\left(t-{\tau }_1\right)+n{\text{e}}^{-\mu {\tau }_2}v\left(t-{\tau }_2\right)-\left(d+k\right)y\left(t\right),\\ v^{\prime }\left(t\right)=ky\left(t\right)-mv\left(t\right)\end{array}$ (1.2)

Based on the fundamental theory of DDEs, the dynamics of disease model (1.2) are studied as varying time delays. We also apply the geometrical criterion to deduce the quantitative property of system (1.2) to derive the instability condition with k regraded as free parameter. The periodical solution is bifurcating from the threshold value of Hopf bifurcation hence Hopf bifurcation lines are tracked as varying free parameters continuously on the parameter plane. The associated normal form is also computed via using Schmidt dimension reduction scheme combined with center manifold theory [21] [22].

The whole paper is organized as the following listed. In Section 2, the stability property of disease infectious equilibrium is analyzed by applying the geometrical criterion to derive Hopf bifurcation with regarded as time delay varying. In Section 3, the normal form is computed by applying center manifold theory and combined with dimension reduction method near Hopf point. Finally, the periodical solutions bifurcated from Hopf point and the continuous calculation of periodical solution is carried out as varying time delay continuously.

2. Study of System Stability

System (1.2) has two equilibrium solutions which are denoted as $E_1=\left(\frac{s}{d},0,0\right)$ and $E_2=\left(x^{*},y^{*},v^{*}\right)$ with the expression

$\begin{array}{l}{x}^{*}=\frac{\theta \left(-kne-\mu {\tau }_2+dm+km\right)}{\alpha cm+kn{\text{e}}^{-\mu {\tau }_2}-dm-km},\\ y^{*}=\frac{{\text{e}}^{\mu {\tau }_2}\alpha m\left(d\ast x-s\right)}{d{\text{e}}^{\mu {\tau }_2}m+k{\text{e}}^{\mu {\tau }_2}m-nk},\\ v^{*}=\frac{k}{m}{y}^{*},\end{array}$

Do axis transformation by setting

$x=x-x^{*},y=y-y^{*},v=v-v^{*}$

then the linearized equation of Equation (1.2) is listed as

$\begin{array}{l}{x}^{\prime }=\left(-d-c\frac{y^{*}}{\theta +x^{*}}+\frac{c{x}^{*}{y}^{*}}{{\left(\theta +x^{*}\right)}^2}\right)x-c\frac{x^{*}}{\theta +x^{*}}y+r\left(x\left(t-{\tau }_1\right)-x\left(t\right)\right),\\ y^{\prime }=\frac{\alpha c\theta y^{*}}{{\left(\theta +x^{*}\right)}^2}x\left(t-{\tau }_1\right)+\frac{\alpha c\theta }{\theta +x^{*}}y\left(t-{\tau }_1\right)-\left(d+k\right)y+n{\text{e}}^{-\mu {\tau }_2}v\left(t-{\tau }_2\right),\\ v^{\prime }=ky-mv,\end{array}$ (2.1)

The characteristic coefficient matrix of Equation (2.1) is

$A=\left(\begin{array}{ccc}\frac{-c\theta y^{*}-d{\theta }^2-2d\theta x^{*}-d{x}^{*2}}{{\left(\theta +x^{*}\right)}^2}+r{\text{e}}^{-\lambda {\tau }_1}-r-\lambda & -\frac{c{x}^{*}}{\theta +x^{*}}& 0\\ \frac{\alpha c\theta y^{*}}{{\left(\theta +x^{*}\right)}^2}{\text{e}}^{-\lambda {\tau }_1}& \frac{\alpha c{x}^{*}}{\theta +x^{*}}{\text{e}}^{-\lambda {\tau }_1}-d-k-\lambda & n{e}^{(-\lambda -\mu ){\tau }_2}\\ 0& k& -m-\lambda \end{array}\right)$ (2.2)

Compute the determination of the matrix A, the characteristic equation of the equilibrium solution $E_2$ has root ${\lambda }_1=-d$ and another branch is

$\begin{array}{c}\Delta \left(\lambda ,k,{\tau }_1,{\tau }_2\right)=\left(-d\theta -s\right){\lambda }^2+\left({\text{e}}^{-\lambda {\tau }_1}\alpha cs-d^2\theta -dk\theta -dm\theta -ds-ks-ms\right)\lambda \\ +{\text{e}}^{-\lambda {\tau }_1}\alpha cms+{\text{e}}^{-\left(\lambda +\mu \right){\tau }_2}dkn\theta +{\text{e}}^{-\left(\lambda +\mu \right){\tau }_2}kns\\ -d^{2}m\theta -dkm\theta -dms-kms\\ =0\end{array}$ (2.3)

It is seen that the stability of $E_1$ is determined at time delay ${\tau }_1={\tau }_2=0$ by the above coefficient matrix with

$Cond{i}_1=\alpha cs-d^2\theta -dk\theta -dm\theta -ds-ks-ms$

That is, Equilibrium solution $E_1$ is asymptotically stable if $Cond{i}_1<0$ ; or otherwise unstable as $Cond{i}_1>0$. Hence, we have the conclusion that $E_1$ is also asymptotically stable if $Cond{i}_1<0$ for any ${\tau }_1=0,{\tau }_2>0$ or ${\tau }_1>0,{\tau }_2=0$.

The stability of equilibrium solution $E_2$ is successively studied as varying time delay ${\tau }_1$ which is the period of the occurrence of disease. Compute the determination of the matrix A, the characteristic equation of the equilibrium solution $E_2$ is written as

$\begin{array}{c}\Delta \left(\lambda ,k,{\tau }_1,{\tau }_2\right)=P\left(\lambda ,k,{\tau }_1,{\tau }_2\right)+Q\left(\lambda ,k,{\tau }_1,{\tau }_2\right){\text{e}}^{-\lambda {\tau }_1}+R\left(\lambda ,k,{\tau }_1,{\tau }_2\right){\text{e}}^{-\left(\lambda +\mu \right){\tau }_2}\\ +S\left(\lambda ,k,{\tau }_1,{\tau }_2\right){\text{e}}^{-\lambda {\tau }_1-\left(\lambda +\mu \right){\tau }_2}+W\left(\lambda ,k,{\tau }_1,{\tau }_2\right){\text{e}}^{-2\lambda {\tau }_1}\end{array}$ (2.4)

with

$\begin{array}{c}P\left(\lambda ,k,{\tau }_1,{\tau }_2\right)=\left(-{\theta }^2-2\theta x^{*}-x^{*2}\right){\lambda }^3+(-c\theta y^{*}-2d{\theta }^2-4d\theta x^{*}-2d{x}^{*2}-k{\theta }^2\\ -2k\theta x^{*}-k{x}^{*2}-m{\theta }^2-2m\theta x^{*}-m{x}^{*2}-r{\theta }^2-2r\theta x^{*}-r{x}^{*2}){\lambda }^2\\ +(-cd\theta y^{*}-ck\theta y^{*}-cm\theta y^{*}-d^2{\theta }^2-2d^2\theta x^{*}-d^2{x}^{*2}-dk{\theta }^2\end{array}$

$\begin{array}{c}-2dk\theta x^{*}-dk{x}^{*2}-2dm{\theta }^2-4dm\theta x^{*}-2dm{x}^{*2}-dr{\theta }^2-2dr\theta x^{*}\\ -dr{x}^{*2}-km{\theta }^2-2km\theta x^{*}-km{x}^{*2}-kr{\theta }^2-2kr\theta x^{*}-kr{x}^{*2}\\ -mr{\theta }^2-2mr\theta x^{*}-mr{x}^{*2})\lambda -cdm\theta y^{*}-ckm\theta y^{*}-d^{2}m{\theta }^2\\ -2d^{2}m\theta x^{*}-d^{2}m{x}^{*2}-dkm{\theta }^2-2dkm\theta x^{*}-dkm{x}^{*2}-dmr{\theta }^2\\ -2dmr\theta x^{*}-dmr{x}^{*2}-kmr{\theta }^2-2kmr\theta x^{*}-kmr{x}^{*2},\end{array}$

$\begin{array}{c}Q\left(\lambda ,k,{\tau }_1,{\tau }_2\right)=\alpha cd\lambda \theta x^{*}+\alpha cd\lambda x^{*2}+\alpha cdm\theta x^{*}+\alpha cdm{x}^{*2}+\alpha c{\lambda }^2\theta x^{*}\\ +\alpha c{\lambda }^2{x}^{*2}+\alpha c\lambda m\theta x^{*}+\alpha c\lambda m{x}^{*2}+\alpha c\lambda r\theta x^{*}+\alpha c\lambda r{x}^{*2}\\ +\alpha cmr\theta x^{*}+\alpha cmr{x}^{*2}+d\lambda r{\theta }^2+2d\lambda r\theta x^{*}+d\lambda r{x}^{*2}+dmr{\theta }^2\\ +2dmr\theta x^{*}+dmr{x}^{*2}+k\lambda r{\theta }^2+2k\lambda r\theta x^{*}+k\lambda r{x}^{*2}+kmr{\theta }^2\\ +2kmr\theta x^{*}+kmr{x}^{*2}0c+{\lambda }^{2}r{\theta }^2+2{\lambda }^{2}r\theta x^{*}+{\lambda }^{2}r{x}^{*2}+\lambda mr{\theta }^2\\ +2\lambda mr\theta x^{*}+\lambda mr{x}^{*2},\end{array}$

$\begin{array}{c}R\left(\lambda ,k,{\tau }_1,{\tau }_2\right)=ckn\theta y^{*}+dkn{\theta }^2+2dkn\theta x^{*}+dkn{x}^{*2}+k\lambda n{\theta }^2+2k\lambda n\theta x^{*}\\ +k\lambda n{x}^{*2}+knr{\theta }^2+2knr\theta x^{*}+knr{x}^{*2},\end{array}$

$\begin{array}{c}S\left(\lambda ,k,{\tau }_1,{\tau }_2\right)=ckn\theta y^{*}+2dkn\theta x^{*}+2k\lambda n\theta x^{*}+2knr\theta x^{*}-knr{\theta }^2\\ -2knr\theta x^{*}-knr{x}^{*2},\end{array}$

$W\left(\lambda ,k,{\tau }_1,{\tau }_2\right)=-\alpha c\lambda r\theta x^{*}-\alpha c\lambda r{x}^{*2}-\alpha cmr\theta x^{*}-\alpha cmr{x}^{*2},$

Set $Y={\text{e}}^{-\mu {\tau }_2}$, with the assumption $\lambda =I\omega$, then separate the real part from the imaginary part in Equation (2.2), one has

$\begin{array}{l}\left(R_I\cos \left(\theta \right)+R_R\sin \left(\theta \right)+S_I\right)Y\sin \left(S\right)+\left(\cos \left(\theta \right)R_R-\sin \left(\theta \right)R_I+S_R\right)Y\cos \left(S\right)\\ +P_R\cos \left(\theta \right)+W_R\cos \left(\theta \right)-P_I\sin \left(\theta \right)+W_I\sin \left(\theta \right)+Q_R=0,\\ \left(-\cos \left(\theta \right)R_R+\sin \left(\theta \right)R_I-S_R\right)Y\sin \left(S\right)+\left(R_I\cos \left(\theta \right)+R_R\sin \left(\theta \right)+S_I\right)Y\cos \left(S\right)\\ +P_I\cos \left(\theta \right)+W_I\cos \left(\theta \right)+P_R\sin \left(\theta \right)-W_R\sin \left(\theta \right)+Q_I=0\end{array}$ (2.5)

and introducing two angle variables $S=\omega {\tau }_2$ and $\theta =\omega {\tau }_1$, by solving $Y\cos \left(S\right)\mathrm{,}Y\sin \left(S\right)$ from Equation (2.4), one gets Hopf surfaces

$\begin{array}{l}Y\cos \left(S\right)={\Phi }_1\left({\tau }_1,{\tau }_2,\mu ,\theta \right),\\ Y\sin \left(S\right)={\Phi }_2\left({\tau }_1,{\tau }_2,\mu ,\theta \right),\end{array}$ (2.6)

with

$\begin{array}{l}{\Phi }_1\left({\tau }_1,{\tau }_2,k,S\right)\\ =\frac{\left(P_I{S}_R-P_R{S}_I-Q_I{R}_R+Q_R{R}_I+S_I{W}_R-S_R{W}_I\right)\sin \left(\theta \right)}{{\left(R_I\cos \left(\theta \right)+R_R\sin \left(\theta \right)+S_I\right)}^2+{\left(\cos \left(\theta \right)R_R-\sin \left(\theta \right)R_I+S_R\right)}^2}\\ +\frac{\left(-P_I{S}_I-P_R{S}_R-Q_I{R}_I-Q_R{R}_R-S_I{W}_I-S_R{W}_R\right)\cos \left(\theta \right)}{{\left(R_I\cos \left(\theta \right)+R_R\sin \left(\theta \right)+S_I\right)}^2+{\left(\cos \left(\theta \right)R_R-\sin \left(\theta \right)R_I+S_R\right)}^2}\\ +\frac{\left(R_I{W}_R-R_R{W}_I\right)\sin \left(2\theta \right)+\left(-R_I{W}_I-R_R{W}_R\right)\cos \left(2\theta \right)}{{\left(R_I\cos \left(\theta \right)+R_R\sin \left(\theta \right)+S_I\right)}^2+{\left(\cos \left(\theta \right)R_R-\sin \left(\theta \right)R_I+S_R\right)}^2}\\ +\frac{-P_I{R}_I-P_R{R}_R-Q_I{S}_I-Q_R{S}_R}{{\left(R_I\cos \left(\theta \right)+R_R\sin \left(\theta \right)+S_I\right)}^2+{\left(\cos \left(\theta \right)R_R-\sin \left(\theta \right)R_I+S_R\right)}^2},\end{array}$

$\begin{array}{l}{\Phi }_2\left({\tau }_1,{\tau }_2,k,S\right)\\ =\frac{\left(P_I{S}_I+P_R{S}_R-Q_I{R}_I-Q_R{R}_R-S_I{W}_I-S_R{W}_R\right)\sin \left(\theta \right)}{{\left(R_I\cos \left(\theta \right)+R_R\sin \left(\theta \right)+S_I\right)}^2+{\left(\cos \left(\theta \right)R_R-\sin \left(\theta \right)R_I+S_R\right)}^2}\\ +\frac{\left(P_I{S}_R-P_R{S}_I+Q_I{R}_R-Q_R{R}_I-S_I{W}_R+S_R{W}_I\right)\cos \left(\theta \right)}{{\left(R_I\cos \left(\theta \right)+R_R\sin \left(\theta \right)+S_I\right)}^2+{\left(\cos \left(\theta \right)R_R-\sin \left(\theta \right)R_I+S_R\right)}^2}\\ +\frac{\left(-R_I{W}_I-R_R{W}_R\right)\sin \left(2\theta \right)+\left(-R_I{W}_R+R_R{W}_I\right)\cos \left(2\theta \right)}{{\left(R_I\cos \left(\theta \right)+R_R\sin \left(\theta \right)+S_I\right)}^2+{\left(\cos \left(\theta \right)R_R-\sin \left(\theta \right)R_I+S_R\right)}^2}\\ +\frac{P_I{R}_R-P_R{R}_I+Q_I{S}_R-Q_R{S}_I}{{\left(R_I\cos \left(\theta \right)+R_R\sin \left(\theta \right)+S_I\right)}^2+{\left(\cos \left(\theta \right)R_R-\sin \left(\theta \right)R_I+S_R\right)}^2},\end{array}$

By the relationship of triage function, one gets the following equality

$F\left(\omega \mathrm{,}{\tau }_1\mathrm{,}{\tau }_2\mathrm{,}k\right)={\Phi }_1{\left({\tau }_1\mathrm{,}{\tau }_2\mathrm{,}k\mathrm{,}\theta \right)}^2+{\Phi }_2{\left({\tau }_1\mathrm{,}{\tau }_2\mathrm{,}k\mathrm{,}\theta \right)}^2-Y^2\equiv 0$ (2.7)

We also assume

$S=S_0+2l_1\pi ,\theta ={\theta }_0+2l_2\pi$

for $l_1=0,1,2,\cdots ,l_2=0,1,2,\cdots$, and define new mapping

$k=k\left(\theta \right)=k\left({\theta }_0+2l_1\pi \right)$ (2.8)

For given $k=k_{*}$, we seek for the finite values of $S_{\mathrm{<mo>\; *\; </mo>}}$ which determined by the inverse function of $k\left(\theta \right)=k_{\mathrm{<mo>\; *\; </mo>}}$ with more than one branch. However, what’s necessary is to find the admissible value of the range that parameter $\mu$ chosen. With this attempt, by differentiating Y with respect to $\theta$, one gets

$\frac{Y\text{d}Y}{\text{d}\theta }={\Phi }_1{{\Phi }^{\prime }}_{1\theta }+{\Phi }_2{{\Phi }^{\prime }}_{2\theta }$ (2.9)

with the assumption $\frac{\text{d}Y}{\text{d}\theta }=0$, we get the value of ${\theta }^{\mathrm{<mo>\; *\; </mo>}}$ of the bottom of the curve $Y\left(\theta \right)$, for some $l_1$ and $\theta ={\theta }_0+2l_1\pi$. In another respect, ${\mu }^{*}=-\frac{1}{2{\tau }_2}\ln \left(Y^2\left({\theta }^{*}\right)\right)$, and we assume $0<\mu <{\mu }^{*}$ to assure the solvability condition for Equation (2.6).

Subsequently, for some $l_2>0$, for given $k_{\mathrm{<mo>\; *\; </mo>}}$, we derive the solution ${\theta }_{\mathrm{<mo>\; *\; </mo>}}$ of mapping (2.8) since the interSection of curve $Y^2\left(\theta \right)$ with the curve $Y={\text{e}}^{-2\mu {\tau }_2}$. Therefore, one gets the threshold value of time delay ${\tau }_1$ determined by Equation (2.5) with the formula

${\tau }_{2*}=\{\begin{array}{ll}\frac{{\tau }_1}{{\theta }_{*}}\left[\arctan \left(\frac{{\Phi }_2}{{\Phi }_1}\right)+2l_1\pi \right],\hfill & \sin \left(\theta \right)>0,\cos \left(\theta \right)>0,\hfill \\ \frac{{\tau }_1}{{\theta }_{*}}\left[\arccos \left({\Phi }_1\right)+2l_1\pi \right],\hfill & \sin \left(\theta \right)>0,\cos \left(\theta \right)<0,\hfill \\ \frac{{\tau }_1}{{\theta }_{*}}\left[\arccos \left({\Phi }_1\right)+\pi +2l_1\pi \right],\hfill & \sin \left(\theta \right)<0,\cos \left(\theta \right)<0,\hfill \\ \frac{{\tau }_1}{{\theta }_{*}}\left[\arctan \left(\frac{{\Phi }_2}{{\Phi }_1}\right)+2\pi +2l_1\pi \right],\hfill & \sin \left(\theta \right)<0,\cos \left(\theta \right)>0\hfill \end{array}$ (2.10)

for $l_1=0,1,2,\cdots$.

To obtain the transversal condition at ${\tau }_2={\tau }_{2*}$, we apply the mathematical technique as shown by geometrical criterion given in paper [20] and rewrite Equation (2.2) as

$\begin{array}{l}P\left(\lambda ,k,{\tau }_1,{\tau }_2\right)+Q\left(\lambda ,k,{\tau }_1,{\tau }_2\right){\text{e}}^{-\lambda {\tau }_1}+R\left(\lambda ,k,{\tau }_1,{\tau }_2\right){\text{e}}^{-\left(\lambda +\mu \right){\tau }_2}\\ +V\left(\lambda ,k,{\tau }_1,{\tau }_2\right){\text{e}}^{-\lambda {\tau }_1-\left(\lambda +\mu \right){\tau }_2}+W\left(\lambda ,k,{\tau }_1,{\tau }_2\right){\text{e}}^{-2\lambda {\tau }_1}=0\end{array}$ (2.11)

Furthermore, differentiate Equation (2.11) with respect to ${\tau }_2$ to get

$\begin{array}{l}{P^{\prime }}_{\lambda }\frac{\text{d}\lambda }{\text{d}{\tau }_2}+{P^{\prime }}_{{\tau }_2}+\left({Q^{\prime }}_{\lambda }\frac{\text{d}\lambda }{\text{d}{\tau }_2}+{Q^{\prime }}_{{\tau }_2}-{\tau }_{1}Q\frac{\text{d}\lambda }{\text{d}{\tau }_2}\right){\text{e}}^{-\lambda {\tau }_1}\\ +\left(\left({R^{\prime }}_{\lambda }-{\tau }_{2}R\right)\frac{\text{d}\lambda }{\text{d}{\tau }_2}+\left({R^{\prime }}_{{\tau }_2}-\lambda R\right)\right){\text{e}}^{-\left(\lambda +\mu \right){\tau }_2}\\ +\left(\left({V^{\prime }}_{\lambda }-\left({\tau }_1+{\tau }_2\right)V\right)\frac{\text{d}\lambda }{\text{d}{\tau }_2}+\left({V^{\prime }}_{{\tau }_2}-\left(\lambda +\mu \right)V\right)\right){\text{e}}^{-\lambda {\tau }_1-\left(\lambda +\mu \right){\tau }_2}\\ +\left(\left({W^{\prime }}_{\lambda }-2{\tau }_{1}W\right)\frac{\text{d}\lambda }{\text{d}{\tau }_2}+{W^{\prime }}_{{\tau }_2}\right){\text{e}}^{-2\lambda {\tau }_1}=0\end{array}$ (2.12)

Solving $\frac{\text{d}\lambda }{\text{d}{\tau }_2}$ from Equation (2.12) to get

$\frac{\text{d}\lambda }{\text{d}{\tau }_2}=-\frac{F_1\left(\lambda \mathrm{,}k\mathrm{,}{\tau }_1\mathrm{,}{\tau }_2\right)}{F_2\left(\lambda \mathrm{,}k\mathrm{,}{\tau }_1\mathrm{,}{\tau }_2\right)}$ (2.13)

with

$\begin{array}{l}{F}_1\left(\lambda ,k,{\tau }_1,{\tau }_2\right)={P^{\prime }}_{{\tau }_2}+{Q^{\prime }}_{{\tau }_2}{\text{e}}^{-\lambda {\tau }_1}+\left({R^{\prime }}_{{\tau }_2}-\lambda R\right){\text{e}}^{-\left(\lambda +\mu \right){\tau }_2}\\ +\left({V^{\prime }}_{{\tau }_2}-\left(\lambda +\mu \right)V\right){\text{e}}^{-\lambda {\tau }_1-\left(\lambda +\mu \right){\tau }_2}+{W^{\prime }}_{{\tau }_2}{\text{e}}^{-2\lambda {\tau }_1},\\ F_2\left(\lambda ,k,{\tau }_1,{\tau }_2\right)={P^{\prime }}_{\lambda }+\left({Q^{\prime }}_{\lambda }-{\tau }_{1}Q\right){\text{e}}^{-\lambda {\tau }_1}+\left({R^{\prime }}_{\lambda }-{\tau }_{2}R\right){\text{e}}^{-\left(\lambda +\mu \right){\tau }_2}\\ +\left({V^{\prime }}_{\lambda }-\left({\tau }_1+{\tau }_2\right)V\right){\text{e}}^{-\lambda {\tau }_1-\left(\lambda +\mu \right){\tau }_2}+\left({W^{\prime }}_{\lambda }-2{\tau }_{1}W\right)\end{array}$

In another respect, Set $\lambda =i\omega$, differentiate $P\mathrm{,}Q\mathrm{,}R\mathrm{,}S\mathrm{,}W$ with respect to $\theta$ to get

$\begin{array}{l}{P^{\prime }}_{\theta }=i{P}^{\prime }\left(i\omega \right){\omega }^{\prime }\left(\theta \right)+{P^{\prime }}_{{\tau }_2}{{\tau }^{\prime }}_2\left(\theta \right)+{P^{\prime }}_k{k}^{\prime }\left(\theta \right),\\ {Q^{\prime }}_{\theta }=i{Q}^{\prime }\left(i\omega \right){\omega }^{\prime }\left(\theta \right)+{Q^{\prime }}_{{\tau }_2}{{\tau }^{\prime }}_2\left(\theta \right)+{Q^{\prime }}_k{k}^{\prime }\left(\theta \right),\\ {R^{\prime }}_{\theta }=i{R}^{\prime }\left(i\omega \right){\omega }^{\prime }\left(\theta \right)+{R^{\prime }}_{{\tau }_2}{{\tau }^{\prime }}_2\left(\theta \right)+{R^{\prime }}_k{k}^{\prime }\left(\theta \right),\\ {V^{\prime }}_{\theta }=i{V}^{\prime }\left(i\omega \right){\omega }^{\prime }\left(\theta \right)+{V^{\prime }}_{{\tau }_2}{{\tau }^{\prime }}_2\left(\theta \right)+{V^{\prime }}_k{k}^{\prime }\left(\theta \right),\\ {W^{\prime }}_{\theta }=i{W}^{\prime }\left(i\omega \right){\omega }^{\prime }\left(\theta \right)+{W^{\prime }}_{{\tau }_2}{{\tau }^{\prime }}_2\left(\theta \right)+{W^{\prime }}_k{k}^{\prime }\left(\theta \right)\end{array}$ (2.14)

and differentiate Equation (2.13) with respect to $\theta$ and k, respectively, to get

$\begin{array}{c}{{\Delta }^{\prime }}_{\theta }\left(\lambda ,k,{\tau }_1,{\tau }_2\right)={P^{\prime }}_{\theta }+\left({Q^{\prime }}_{\theta }-iQ\right){\text{e}}^{-i\theta }\\ +\left({R^{\prime }}_{\theta }-i{\omega }^{\prime }\left(\theta \right){\tau }_{2}R\right){\text{e}}^{-\left(i\omega +\mu \right){\tau }_2}\\ +\left({V^{\prime }}_{\theta }-iV-i{\omega }^{\prime }\left(\theta \right){\tau }_{2}V\right){\text{e}}^{-i\theta -\left(i\omega +\mu \right){\tau }_2}\\ +\left({W^{\prime }}_{\theta }-2iW\right){\text{e}}^{-2i\theta }\end{array}$ (2.15)

and

${{\Delta }^{\prime }}_k\left(\lambda ,k,{\tau }_1,{\tau }_2\right)={P^{\prime }}_k+{Q^{\prime }}_k{\text{e}}^{-i\theta }+{R^{\prime }}_k{\text{e}}^{-\left(i\omega +\mu \right){\tau }_2}+{V^{\prime }}_k{\text{e}}^{-i\theta -\left(i\omega +\mu \right){\tau }_2}+{W^{\prime }}_k{\text{e}}^{-2i\theta }$ (2.16)

Noticed that ${\omega }^{\prime }\left(\theta \right)=\frac{1}{{\tau }_1}$, then we have

$\frac{\text{d}\lambda }{\text{d}{\tau }_2}=\frac{-\frac{i}{{\tau }_1}{F}_1\left(i\omega ,k,{\tau }_1,{\tau }_2\right)}{{{\Delta }^{\prime }}_k{k}^{\prime }\left(\theta \right)+{{\tau }^{\prime }}_2\left(\theta \right)\left(F_1+i\omega R{e}^{-\left(i\omega +\mu \right){\tau }_2}+\left(i\omega +\mu \right)V{\text{e}}^{-i\theta -\left(i\omega +\mu \right){\tau }_2}\right)}$ (2.17)

Further, we have

$\begin{array}{l}\delta \left({\theta }_{\mathrm{<mo>\; *\; </mo>}}\right)=sign\left\{Re\left(\frac{\text{d}\lambda }{\text{d}{\tau }_2}\right)\right\}=sign\left\{Re\left(\frac{1}{\frac{\text{d}\lambda }{\text{d}{\tau }_2}}\right)\right\}\\ =sign\Re \{i\left({{\Delta }^{\prime }}_k{k}^{\prime }\left({\theta }_{\mathrm{<mo>\; *\; </mo>}}\right)+{{\tau }^{\prime }}_2\left({\theta }_{\mathrm{<mo>\; *\; </mo>}}\right)\left(F_1+i\omega R{\text{e}}^{-\left(i\omega +\mu \right){\tau }_2}+\left(i\omega +\mu \right)V{\text{e}}^{-i{\theta }_{\mathrm{<mo>\; *\; </mo>}}-\left(i\omega +\mu \right){\tau }_2}\right)\right)\\ \times \left(\Re \left(F_1\right)-i\Im \left(F_1\right)\right)\}\\ =sign\{\Im \left(F_1\right)\Re \left({{\Delta }^{\prime }}_k{k}^{\prime }\left({\theta }_{\mathrm{<mo>\; *\; </mo>}}\right)+{{\tau }^{\prime }}_2\left({\theta }_{\mathrm{<mo>\; *\; </mo>}}\right)\left(F_1+i\omega R{\text{e}}^{-\left(i\omega +\mu \right){\tau }_2}+\left(i\omega +\mu \right)V{\text{e}}^{-i\theta -\left(i\omega +\mu \right){\tau }_2}\right)\right)\\ -\Re \left(F_1\right)\Im \left({{\Delta }^{\prime }}_k{k}^{\prime }\left({\theta }_{\mathrm{<mo>\; *\; </mo>}}\right)+{{\tau }^{\prime }}_2\left({\theta }_{\mathrm{<mo>\; *\; </mo>}}\right)\left(F_1+i\omega R{\text{e}}^{-\left(i\omega +\mu \right){\tau }_2}+\left(i\omega +\mu \right)V{\text{e}}^{-i{\theta }_{\mathrm{<mo>\; *\; </mo>}}-\left(i\omega +\mu \right){\tau }_2}\right)\right)\}\end{array}$

$\begin{array}{l}=sign\{\Im \left(F_1\right)\Re \left({{\Delta }^{\prime }}_k\right)k^{\prime }\left({\theta }_{\mathrm{<mo>\; *\; </mo>}}\right)+{{\tau }^{\prime }}_2\left({\theta }_{\mathrm{<mo>\; *\; </mo>}}\right)(\Re \left(F_1\right)+{\text{e}}^{-\mu {\tau }_2}\omega (\Re \left(R\right)\sin \left(\frac{{\theta }_{\mathrm{<mo>\; *\; </mo>}}}{{\tau }_1}{\tau }_2\right)\\ -\Im \left(R\right)\cos \left(\frac{{\theta }_{\mathrm{<mo>\; *\; </mo>}}}{{\tau }_1}{\tau }_2\right)+\mu \left(\Re \left(V\right)\cos \left({\theta }_{\mathrm{<mo>\; *\; </mo>}}+\frac{{\theta }_{\mathrm{<mo>\; *\; </mo>}}}{{\tau }_1}{\tau }_2\right)+\Im \left(V\right)\sin \left({\theta }_{\mathrm{<mo>\; *\; </mo>}}+\frac{{\theta }_{\mathrm{<mo>\; *\; </mo>}}}{{\tau }_1}{\tau }_2\right)\right)\\ -\omega \left(\Im \left(V\right)\cos \left({\theta }_{\mathrm{<mo>\; *\; </mo>}}+\frac{{\theta }_{\mathrm{<mo>\; *\; </mo>}}}{{\tau }_1}{\tau }_2\right)+\Re \left(V\right)\sin \left({\theta }_{\mathrm{<mo>\; *\; </mo>}}+\frac{{\theta }_{\mathrm{<mo>\; *\; </mo>}}}{{\tau }_1}{\tau }_2\right)\right))\\ -\Re \left(F_1\right)(\Im \left({{\Delta }^{\prime }}_k\right)k^{\prime }\left({\theta }_{\mathrm{<mo>\; *\; </mo>}}\right)+{{\tau }^{\prime }}_2\left({\theta }_{\mathrm{<mo>\; *\; </mo>}}\right)\mathrm{<mo\; stretchy="false">\; (\; </mo>}\Im \left(F_1\right)+{\text{e}}^{-\mu {\tau }_2}\omega (\Im \left(R\right)\sin \left(\frac{{\theta }_{\mathrm{<mo>\; *\; </mo>}}}{{\tau }_1}{\tau }_2\right)\\ +\Re \left(R\right)\cos \left(\frac{{\theta }_{\mathrm{<mo>\; *\; </mo>}}}{{\tau }_1}{\tau }_2\right)+\mu \left(\Im \left(V\right)\cos \left({\theta }_{\mathrm{<mo>\; *\; </mo>}}+\frac{{\theta }_{\mathrm{<mo>\; *\; </mo>}}}{{\tau }_1}{\tau }_2\right)-\Re \left(V\right)\sin \left({\theta }_{\mathrm{<mo>\; *\; </mo>}}+\frac{{\theta }_{\mathrm{<mo>\; *\; </mo>}}}{{\tau }_1}{\tau }_2\right)\right)\\ +\omega \left(\Re \left(V\right)\cos \left({\theta }_{\mathrm{<mo>\; *\; </mo>}}+\frac{{\theta }_{\mathrm{<mo>\; *\; </mo>}}}{{\tau }_1}{\tau }_2\right)+\Im \left(V\right)\sin \left({\theta }_{\mathrm{<mo>\; *\; </mo>}}+\frac{{\theta }_{\mathrm{<mo>\; *\; </mo>}}}{{\tau }_1}{\tau }_2\right)\right))\}\end{array}$ (2.18)

Therefore we give the following result of the stability analysis of equilibrium solution $E_2$,

Theorem 2.1 $E_2$ losses its stability at the threshold value ${\tau }_1={\tau }_1^{*}$ with the fixed parameter $k=k_{*}$ with $\mu$ satisfies the solvability condition. The transversal condition $\delta \left({\tau }_{\mathrm{1\; <mo>\; *\; </mo>}}\right)\ne 0$ given by formula (2.19) determines the transverse direction of the characteristic roots associated with Hopf bifurcation. That is, there is a pair of imaginary roots that can transverse from the left half plane to the right half plane if $\delta \left({\tau }_{1*}\right)>0$, or otherwise transverse from the right half plane to the left half plane if $\delta \left({\tau }_{1*}\right)<0$.

3. Stability of Periodic Solution

In Section 2, stability of infectious disease equilibrium solution $E_2$ loss as parameter cross over Hopf line given that the transversal condition $\frac{\text{d}\Re \lambda }{\text{d}\tau }\ne 0$. The stability of bifurcating periodic solution is analyzed by perturbation technique near Hopf point $\left(k_c\mathrm{,}{\tau }_c\right)$. We suppose that Hopf bifurcation happens at the critical parameter pairs $\left(k_c\mathrm{,}{\tau }_{1c}\right)$. Do axis transformation $\bar{x}=x-x^{*},\bar{y}=y-y^{*},\bar{v}=v-v^{*}$ (for simplicity the overbar is ignored), to write Equation (1.1) as the following truncated system with Taylor expansion to three orders,

$\begin{array}{l}{x}^{\prime }\left(t\right)=\left(-d-c\frac{y^{*}}{\theta +x^{*}}+\frac{c{x}^{*}{y}^{*}}{{\left(\theta +x^{*}\right)}^2}\right)x-c\frac{x^{*}}{\theta +x^{*}}y+r\left(x\left(t-{\tau }_1\right)-x\left(t\right)\right)\\ +\frac{c\theta y^{*}}{{\left(\theta +x^{*}\right)}^3}{x}^2-\frac{c\theta y^{*}}{{\left(\theta +x^{*}\right)}^4}{x}^3+\frac{c\theta }{{\left(\theta +x^{*}\right)}^3}{x}^{2}y,\\ y^{\prime }\left(t\right)=\frac{\alpha c\theta y^{*}}{{\left(\theta +x^{*}\right)}^2}x\left(t-{\tau }_1\right)+\frac{\alpha c\theta }{\theta +x^{*}}y\left(t-{\tau }_1\right)-\left(d+k\right)y+n{\text{e}}^{-\mu {\tau }_2}v\left(t-{\tau }_2\right)\\ +\frac{\alpha c\theta y^{*}}{{\left(\theta +x^{*}\right)}^4}x{\left(t-{\tau }_1\right)}^3-\frac{\alpha c\theta }{{\left(\theta +x^{*}\right)}^3}{x}^2\left(t-{\tau }_1\right)y\left(t-{\tau }_1\right)\\ -\frac{\alpha c\theta y^{*}}{{\left(\theta +x^{*}\right)}^3}{x}^2\left(t-{\tau }_1\right)+\frac{\alpha c\theta }{{\left(\theta +x^{*}\right)}^2}x\left(t-{\tau }_1\right)y\left(t-{\tau }_1\right),\\ v^{\prime }\left(t\right)=ky\left(t\right)-mv\left(t\right)\end{array}$ (3.1)

Equation (3.1) is defined on Banach space $C\left(\left[-{\tau }_{\max },0\right],R^3\right)$ with the supreme norm defined as $\Vert \varphi \Vert ={\sup }_{\theta \in \left[-{\tau }_{\max },0\right]}\left|\varphi \left(\theta \right)\right|$, herein ${\tau }_{\max }=\max \left\{{\tau }_1,{\tau }_2\right\}$.

We define the phase space to be the extended phase space $C\left(\left[-{\tau }_{\max },0\right],R^3\right)$ with a possible jump discontinuity at $\theta =0$. Set ${\tau }_{1e}{v}_e={\tau }_1-{\tau }_{1c}$ and $k_e{v}_e=k-k_c$ and the linearized equation is written as

$\begin{array}{l}{x}^{\prime }\left(t\right)=a_{11}x+a_{12}y+r\left(x\left(t-{\tau }_{1c}\right)-x\right)+r\left(x\left(t-{\tau }_{1c}-{\tau }_{1e}{v}_e\right)-x\left(t-{\tau }_{1c}\right)\right),\\ y^{\prime }\left(t\right)=b_{1}x\left(t-{\tau }_{1c}\right)+b_{2}y\left(t-{\tau }_{1c}\right)-\left(d+k_c\right)y+n{\text{e}}^{-\mu {\tau }_2}v\left(t-{\tau }_2\right)-k_e{v}_{e}y\\ +b_1\left(x\left(t-{\tau }_{1c}-{\tau }_{1e}{v}_e\right)-x\left(t-{\tau }_{1c}\right)\right)+b_2\left(y\left(t-{\tau }_{1c}-{\tau }_{1e}{v}_e\right)-y\left(t-{\tau }_{1c}\right)\right),\\ v^{\prime }\left(t\right)=k_{c}y+k_e{v}_{e}y-mv\end{array}$ (3.2)

with

$\begin{array}{l}{a}_{11}=-d-c\frac{y^{*}}{\theta +x^{*}}+\frac{c{x}^{*}{y}^{*}}{{\left(\theta +x^{*}\right)}^2},a_{12}=-c\frac{x^{*}}{\theta +x^{*}},\\ b_1=\frac{\alpha c\theta y^{*}}{{\left(\theta +x^{*}\right)}^2},b_2=\frac{\alpha c\theta }{\theta +x^{*}}\end{array}$

Set $u\left(t\right)={\left(x\left(t\right),y\left(t\right),v\left(t\right)\right)}^{\text{T}}$ and $u_t=u\left(t+\theta \right)$ for $-{\tau }_{\max }\le \theta \le 0$, Equation (3.2) can be rewritten as

$u^{\prime }\left(t\right)=L{u}_t+L_e{u}_t$ (3.3)

By Rieze representation theorem, there exists a $3\times 3$ matrix function $\eta \left(\theta \right)$ of bounded variation and ${\eta }_e\left(\theta \right)$ to express

$L{u}_t={\displaystyle {\int }_{-{\tau }_{\max }}^0}\text{d}\eta \left(\theta \right)u_t\left(\theta \right)$ (3.4)

with

$\begin{array}{l}\text{d}\eta \left(\theta \right)=\left(\begin{array}{ccc}{a}_{11}\left(k_c\right)\delta \left(\theta \right)& a_{12}\left(k_c\right)\delta \left(\theta \right)& 0\\ 0& -\left(d+k_c\right)\delta \left(\theta \right)& n{\text{e}}^{-\mu {\tau }_2}\delta \left(\theta +{\tau }_2\right)\\ 0& k_c\delta \left(\theta \right)& -m\delta \left(\theta \right)\end{array}\right)\\ +\left(\begin{array}{ccc}r\delta \left(\theta +{\tau }_{1c}\right)-r\delta \left(\theta \right)& 0& 0\\ b_1\left(k_c\right)\delta \left(\theta +{\tau }_{1c}\right)& b_2\left(k_c\right)\delta \left(\theta +{\tau }_{1c}\right)& 0\\ 0& 0& 0\end{array}\right)\end{array}$

and

$L_e{u}_t={\displaystyle {\int }_{-{\tau }_{\max }}^0}\text{d}{\eta }_e\left(\theta \right)u_t\left(\theta \right)$ (3.5)

with

$\text{d}{\eta }_{e1}\left(\theta \right)=\left(\begin{array}{ccc}{a^{\prime }}_{11}\left(k_c\right)k_e\delta \left(\theta \right)& {a^{\prime }}_{12}\left(k_c\right)k_e\delta \left(\theta \right)& 0\\ -{b^{\prime }}_1\left(k_c\right)k_e\delta \left(\theta +{\tau }_{1c}\right)& -{b^{\prime }}_2\left(k_c\right)k_e\delta \left(\theta +{\tau }_{1c}\right)-k_e\delta \left(\theta \right)& 0\\ 0& k_e\delta \left(\theta \right)& 0\end{array}\right)$

and

$\text{d}{\eta }_{e2}\left(\theta \right)=\left(\begin{array}{ccc}r\left(\delta \left(\theta +{\tau }_{1c}+{\tau }_{1e}{v}_e\right)-\delta \left(\theta +{\tau }_{1c}\right)\right)& 0& 0\\ b_1\left(k_c\right)\left(\delta \left(\theta +{\tau }_{1c}+{\tau }_{1e}{v}_e\right)-\delta \left(\theta +{\tau }_{1c}\right)\right)& b_2\left(k_c\right)\left(\delta \left(\theta +{\tau }_{1c}+{\tau }_{1e}{v}_e\right)-\delta \left(\theta +{\tau }_{1c}\right)\right)& 0\\ 0& 0& 0\end{array}\right)$

Based on the fundamental theory of DDEs, define $\mathcal{A}$ to be the infinitesimal generator of the solution semigroup associated with linear operator $L\varphi (\cdot )$ such that

$\mathcal{A}\varphi =\{\begin{array}{l}\frac{\text{d}\varphi }{\text{d}\theta },-{\tau }_{\max }\le \theta <0,\hfill \\ {\displaystyle {\int }_{-{\tau }_{\max }}^0}\text{d}\eta \left(\theta \right)\varphi \left(\theta \right),\theta =0,\hfill \end{array}$ (3.6)

for $\varphi \in C\left(\left[-{\tau }_{\max },0\right],R^3\right)$. For $\psi \in C^{*}\left(\left[0,{\tau }_{\max }\right],R^3\right)$, the adjoint operator of $\mathcal{A}$ is defined as

${\mathcal{A}}^{*}\psi =\{\begin{array}{l}-\frac{\text{d}\psi }{\text{d}s},0<s\le {\tau }_{\max },\hfill \\ {\displaystyle {\int }_{-{\tau }_{\max }}^0}\text{d}\eta \left(-s\right)\psi \left(s\right),s=0,\hfill \end{array}$ (3.7)

Define the inner product

$\langle \psi ,\varphi \rangle ={\bar{\psi }}^{\text{T}}\left(0\right)\varphi \left(0\right)-{\displaystyle {\int }_{-{\tau }_{\max }}^0}{\displaystyle {\int }_0^{\theta }}{\bar{\psi }}^{\text{T}}\left(\xi +\theta \right)\text{d}\eta \left(\theta \right)\varphi \left(\xi \right)$ (3.8)

Suppose that the eigenvector q and $q^{\mathrm{<mo>\; *\; </mo>}}$ satisfy

$\mathcal{A}q\left(\theta \right)=i\omega q\left(\theta \right)\mathrm{,}{\mathcal{A}}^{\mathrm{<mo>\; *\; </mo>}}{q}^{\mathrm{<mo>\; *\; </mo>}}\left(s\right)=i\omega q^{\mathrm{<mo>\; *\; </mo>}}\left(s\right)$ (3.9)

with

$\langle q\mathrm{,}{q}^{\mathrm{<mo>\; *\; </mo>}}\rangle =\mathrm{1,}\langle \bar{q}\mathrm{,}{q}^{\mathrm{<mo>\; *\; </mo>}}\rangle =0$

For example, we choose

$q=\left(\begin{array}{c}\left(-m-i\omega \right)a_{12}\\ -\left(i\omega -r{\text{e}}^{-i\omega {\tau }_{1c}}-a_{11}+r\right)\left(i\omega +m\right)\\ k_c\left(a_{11}+r{\text{e}}^{-i\omega {\tau }_{1c}}-r-i\omega \right)\end{array}\right){\text{e}}^{i\omega \theta }$

$q^{*}=\frac{1}{M}\left(\begin{array}{c}-b_1{\text{e}}^{i\omega {\tau }_{1c}}\left(-m+i\omega \right)\\ \left(a_{11}+r{e}^{i\omega {\tau }_{1c}}-r+i\omega \right)\left(-m+i\omega \right)\\ \left(-n{\text{e}}^{-\mu {\tau }_2+i\omega {\tau }_2}\right)\left(a_{11}+r{\text{e}}^{i\omega {\tau }_{1c}}-r+i\omega \right)\end{array}\right){\text{e}}^{-i\omega s}$

with

$\begin{array}{l}M=-b_1{\text{e}}^{i\omega {\tau }_{1c}}\left(-m+i\omega \right)\left(a_{12}\left(m^2+{\omega }^2\right)\right)\\ -{\left(a_{11}+r{\text{e}}^{i\omega {\tau }_{1c}}-r+i\omega \right)}^2\left(-m+{\omega }^2\right)\left(-m-i\omega \right)\\ +\left(-n{\text{e}}^{-\mu {\tau }_2}{\text{e}}^{i\omega {\tau }_2}\left(a_{11}+r{\text{e}}^{i\omega {\tau }_{1c}}-r+i\omega \right)\right)k_c\left(a_{11}+r{\text{e}}^{i\omega {\tau }_{1c}}-r+i\omega \right)\left(-m-i\omega \right)\\ +(2(\left(1/2\right)b_2\left(-i\omega +m\right)r^2{\text{e}}^{i\omega {\tau }_{1c}}+((\left(1/2i\right){\omega }^3+\left(-\left(1/2\right)m-r+a_{11}\right){\omega }^2\end{array}$

$\begin{array}{l}-i\left(\left(1/2\right)r+m-\left(1/2\right)a_{11}\right)\left(r-a_{11}\right)\omega +\left(1/2\right)m{\left(r-a_{11}\right)}^2)b_2\\ +\left(1/2\right)b_1(-a_{12}{\omega }^2+\left(-im{a}_{12}-ir{a}_{12}+i{a}_{11}{a}_{12}\right)\omega +\left(r-a_{11}\right)a_{12}m)){\text{e}}^{-i\omega {\tau }_{1c}}\\ -b_{2}r\left(-{\omega }^2+\left(-im-ir+i{a}_{11}\right)\omega +m\left(r-a_{11}\right)\right)))\left(m^2+{\omega }^2\right){\tau }_{1c}\\ -{\tau }_2{\text{e}}^{{\tau }_2\left(-i\omega -\mu \right)}{\left(-i\omega -r{\text{e}}^{i\omega {\tau }_{1c}}-a_{11}+r\right)}^{2}n{k}_c\left(-m^2-{\omega }^2\right)\end{array}$

We also write Equation (3.1) into an operator differential equation

$u_{t^{\prime }}=\mathcal{A}{u}_t+\mathcal{R}{u}_t,-{\tau }_{\max }\le \theta \le 0$ (3.10)

with nonlinear terms

$\mathcal{R}{u}_t=\{\begin{array}{ll}0,\hfill & -{\tau }_{\max }\le \theta <0,\hfill \\ L_e{u}_t+F\left(u_t\right),\hfill & \theta =0\hfill \end{array}$ (3.11)

and

$\begin{array}{l}{L}_e{u}_t=\left(\begin{array}{c}{a^{\prime }}_{11}\left(k_c\right))k_e{v}_{e}x\left(t\right)+{a^{\prime }}_{12}\left(k_c\right)k_e{v}_{e}y\left(t\right)\\ -{b^{\prime }}_1\left(k_c\right)k_e{v}_{e}x\left(t-{\tau }_{1c}\right)-{b^{\prime }}_2\left(k_c\right)k_e{v}_{e}y\left(t-{\tau }_{1c}\right)-k_e{v}_{e}y\left(t\right)\\ k_e{v}_{e}y\left(t\right)\end{array}\right)\\ +\left(\begin{array}{c}rx\left(t-{\tau }_1\right)-rx\left(t-{\tau }_{1c}\right)\\ b_1\left(k_c\right)\left(x\left(t-{\tau }_1\right)-x\left(t-{\tau }_{1c}\right)\right)+b_2\left(k_c\right)\left(x\left(t-{\tau }_1\right)-x\left(t-{\tau }_{1c}\right)\right)\\ 0\end{array}\right)\mathrm{,}\end{array}$

$\begin{array}{c}F\left(u_t\right)=\left(\begin{array}{c}{c}_1{x}^2\\ c_2{x}^2\left(t-{\tau }_{1c}\right)+c_{3}x\left(t-{\tau }_{1c}\right)y\left(t-{\tau }_{1c}\right)\\ 0\end{array}\right)\\ +\left(\begin{array}{c}{c}_4{x}^3+c_5{x}^{2}y\\ c_{6}x{\left(t-{\tau }_{1c}\right)}^3+c_7{x}^2\left(t-{\tau }_{1c}\right)y\left(t-{\tau }_{1c}\right)\\ 0\end{array}\right)\end{array}$ (3.12)

$\begin{array}{l}{c}_1=\frac{c\theta y^{*}}{{\left(\theta +x^{*}\right)}^3},c_2=-\frac{\alpha c\theta y^{*}}{{\left(\theta +x^{*}\right)}^3},c_3=\frac{\alpha c\theta }{{\left(\theta +x^{*}\right)}^2},c_4=-\frac{c\theta y^{*}}{{\left(\theta +x^{*}\right)}^4},\\ c_5=\frac{c\theta }{{\left(\theta +x^{*}\right)}^3},c_6=\frac{\alpha c\theta y^{*}}{{\left(\theta +x^{*}\right)}^4},c_7=-\frac{\alpha c\theta }{{\left(\theta +x^{*}\right)}^3}\end{array}$

Since Hopf bifurcation occurs at the critical parameter pairs $\left(k_c\mathrm{,}{\tau }_{1c}\right)$ with the imaginary roots $\Lambda =\left\{i\omega \mathrm{,}-i\omega \right\}$, it is supposed $P_{\Lambda }$ be the corresponding eigenspace and Q is its complementary subspace. By decomposition, $C=P_{\Lambda }\oplus Q$, and we define $\pi \mathrm{<mo>\; :\; </mo>}C\to P_{\Lambda }$ be the projection operator and

$X_0=\{\begin{array}{ll}I,\hfill & \theta =0,\hfill \\ 0,\hfill & -{\tau }_{\max }\le \theta <0\hfill \end{array}$ (3.13)

Therefore, set $u_t=zq+\bar{z}\bar{q}+y$ with $y\in Q$, then by Equation (3.5) one gets

$\begin{array}{l}{z}^{\prime }=i\omega z+{\bar{q}}^{*}\left(0\right)R\left(zq+\bar{z}\bar{q}+y\right),\\ {\bar{z}}^{\prime }=-i\omega \bar{z}+q^{*}\left(0\right)R\left(zq+\bar{z}\bar{q}+y\right),\\ y^{\prime }=\mathcal{A}y+\left(I-\pi \right)X_{0}R\left(zq+\bar{z}\bar{q}+y\right)\end{array}$ (3.14)

henceforth, set $\Phi \left(\theta \right)=\left(q\left(\theta \right)\mathrm{,}\bar{q}\left(\theta \right)\right)$, Equation (3.14) is written as Taylor expansion to be truncated to 3 order with the expression

$\begin{array}{l}\left(\begin{array}{c}{z}^{\prime }\\ {\bar{z}}^{\prime }\end{array}\right)=Jz+f_2^{\left(1\right)}\left(z,\bar{z},0,v\right)+f_3^{\left(1\right)}\left(z,\bar{z},y,0\right),\\ y^{\prime }=\mathcal{A}y+\{\begin{array}{ll}-\Phi \left(\theta \right)f_2^{\left(1\right)}\left(z,\bar{z},0,0\right),\hfill & -{\tau }_{\max }\le \theta <0,\hfill \\ f_2^2\left(z,\bar{z},0,0\right)-\Phi \left(0\right)f_2^{\left(1\right)}\left(z,\bar{z},0,0\right)\hfill & \theta =0,\hfill \end{array}\end{array}$ (3.15)

with $J=\left(\begin{array}{cc}i\omega & 0\\ 0& -i\omega \end{array}\right)$, and

$F\left(zq+\bar{z}\bar{q}+y\right)=f_{2000}^{\left(2\right)}{z}^2+f_{1100}^{\left(2\right)}z\bar{z}+f_{0200}^{\left(2\right)}{\bar{z}}^2+\cdots$ (3.16)

Define the operator $M_2^{\mathrm{1,2}}$ on the space of homogeneous polynomial $H_2\left(z\mathrm{,}\bar{z}\right)$ by

$\begin{array}{l}{M}_2^1\left(\begin{array}{c}p\left(z,\bar{z}\right)\\ \bar{p}\left(z,\bar{z}\right)\end{array}\right)=J\left(\begin{array}{c}p\left(z,\bar{z}\right)\\ \bar{p}\left(z,\bar{z}\right)\end{array}\right)-\left(\begin{array}{cc}{p}_z& p_{\bar{z}}\\ {\bar{p}}_z& {\bar{p}}_{\bar{z}}\end{array}\right)J\left(\begin{array}{c}z\\ \bar{z}\end{array}\right),\\ M_2^{2}U\left(z,\bar{z}\right)=\mathcal{A}U\left(z,\bar{z}\right)-\left(U_z,U_{\bar{z}}\right)J\left(\begin{array}{c}z\\ \bar{z}\end{array}\right).\end{array}$ (3.17)

On the center manifold, the normal form of Equation (3.15) is expressed as

$\left(\begin{array}{c}{z}^{\prime }\\ {\bar{z}}^{\prime }\end{array}\right)=J\left(\begin{array}{c}z\\ \bar{z}\end{array}\right)+g_2^{\left(1\right)}\left(z,\bar{z},0,v\right)+g_3^{\left(1\right)}\left(z,\bar{z},0,0\right)$ (3.18)

with

$\begin{array}{l}{g}_2^{\left(1\right)}\left(z\mathrm{,}\bar{z}\mathrm{,0,}v\right)=Pro{j}_{Im{\left(M_2\right)}^c}{f}_2^{\left(1\right)}\left(z\mathrm{,}\bar{z}\mathrm{,0,}v\right)\mathrm{,}\\ g_3^{\left(1\right)}\left(z\mathrm{,}\bar{z}\mathrm{,0,0}\right)=Pro{j}_{Im{\left(M_3\right)}^c}{\tilde{f}}_3^{\left(1\right)}\left(z\mathrm{,}\bar{z}\mathrm{,0,0}\right)\end{array}$ (3.19)

By choosing

$\begin{array}{l}\left(\begin{array}{c}p\left(z\mathrm{,}\bar{z}\right)\\ \bar{p}\left(z\mathrm{,}\bar{z}\right)\end{array}\right)=M_1^{\left(1\right)-}{f}_2^{\left(1\right)}\left(z\mathrm{,}\bar{z}\mathrm{,0,}v\right)\mathrm{,}\\ U\left(z\mathrm{,}\bar{z}\right)=M_2^{\left(2\right)-}{f}_2^{\left(2\right)}\left(z\mathrm{,}\bar{z}\mathrm{,0,0}\right)\mathrm{,}\end{array}$

one gets

$\begin{array}{c}{\tilde{f}}_3^{\left(1\right)}\left(z\mathrm{,}\bar{z}\mathrm{,0,0}\right)=f_3^{\left(1\right)}\left(z\mathrm{,}\bar{z}\mathrm{,0,0}\right)-f_{2z}^{\left(1\right)}\left(z\mathrm{,}\bar{z}\mathrm{,0,0}\right)p\left(z\mathrm{,}\bar{z}\right)\\ -f_{2\bar{z}}^{\left(1\right)}\left(z\mathrm{,}\bar{z}\mathrm{,0,0}\right)\bar{p}\left(z\mathrm{,}\bar{z}\right)+f_{2y}\left(z\mathrm{,}\bar{z}\mathrm{,0,0}\right)U\left(z\mathrm{,}\bar{z}\right)\end{array}$ (3.20)

and the normal form (3.17) is derived as

$z^{\prime }=i\omega z+k_{10}z{v}_e+k_{21}{z}^2\bar{z}$ (3.21)

The bases of $Im\left(M_1^2\right)$ is expressed as

$\begin{array}{l}\left(\begin{array}{c}2i\omega \bar{z}{v}_e\\ 0\end{array}\right)\mathrm{,}\left(\begin{array}{c}0\\ -2i\omega z{v}_e\end{array}\right)\mathrm{,}\left(\begin{array}{c}-i\omega z^2\\ 0\end{array}\right)\mathrm{,}\left(\begin{array}{c}i\omega z\bar{z}\\ 0\end{array}\right)\mathrm{,}\\ \left(\begin{array}{c}3i\omega {\bar{z}}^2\\ 0\end{array}\right)\mathrm{,}\left(\begin{array}{c}0\\ -3i\omega z^2\end{array}\right)\mathrm{,}\left(\begin{array}{c}0\\ -i\omega z\bar{z}\end{array}\right)\mathrm{,}\left(\begin{array}{c}0\\ i\omega {\bar{z}}^2\end{array}\right)\end{array}$

The bases of $Im\left(M_1^3\right)$ is expressed as

$\left(\begin{array}{c}-2i\omega z^3\\ 0\end{array}\right),\left(\begin{array}{c}2i\omega z{\bar{z}}^2\\ 0\end{array}\right),\left(\begin{array}{c}4i\omega {\bar{z}}^3\\ 0\end{array}\right),\left(\begin{array}{c}0\\ -4i\omega z^3\end{array}\right),\left(\begin{array}{c}0\\ -2i\omega z^2\bar{z}\end{array}\right),\left(\begin{array}{c}0\\ 2i\omega {\bar{z}}^3\end{array}\right)$

Hence, we choose

$p\left(z\mathrm{,}\bar{z}\right)=\frac{f_{2000}^{\left(\mathrm{1,1}\right)}}{-i\omega }{z}^2+\frac{f_{1100}^{\left(\mathrm{1,1}\right)}}{i\omega }z\bar{z}+\frac{f_{0200}^{\left(\mathrm{1,1}\right)}}{3i\omega }{\bar{z}}^2\mathrm{,}$ (3.22)

Set

$U\left(z\mathrm{,}\bar{z}\right)=H_{20}{z}^2+H_{11}z\bar{z}+H_{02}{\bar{z}}^2$

By the above definition of $U\left(z\mathrm{,}\bar{z}\right)$ in Equation (3.7), one gets

$\begin{array}{l}\mathcal{A}{H}_{20}=2i\omega H_{20}+q\left(\theta \right)f_{2000}^{\left(1,1\right)}+\bar{q}\left(\theta \right){\bar{f}}_{0200}^{\left(1,1\right)},\\ \mathcal{A}{H}_{11}=q\left(\theta \right)f_{1100}^{\left(1,1\right)}+\bar{q}\left(\theta \right){\bar{f}}_{1100}^{\left(1,1\right)},\\ \mathcal{A}{H}_{02}=-2i\omega H_{02}+q\left(\theta \right)f_{0200}^{(1,1)}+\bar{q}\left(\theta \right){\bar{f}}_{2000}^{\left(1,1\right)}\end{array}$ (3.23)

which satisfies the initial condition

$\begin{array}{l}L{H}_{20}=2i\omega H_{20}\left(0\right)+q\left(0\right)f_{2000}^{\left(1,1\right)}+\bar{q}(0){\bar{f}}_{0200}^{\left(1,1\right)}-f_{2000}^{\left(2\right)},\\ L{H}_{11}=q\left(0\right)f_{1100}^{\left(1,1\right)}+\bar{q}\left(0\right){\bar{f}}_{1100}^{\left(1,1\right)}-f_{1100}^{\left(2\right)},\\ L{H}_{02}=-2i\omega H_{02}\left(0\right)+q\left(0\right)f_{0200}^{\left(1,1\right)}+\bar{q}\left(0\right){\bar{f}}_{2000}^{\left(1,1\right)}-f_{0200}^{\left(2\right)}\end{array}$ (3.24)

By the computation of integral of Equation (3.23), we derive all of the coeffcients $H_{20}\mathrm{,}{H}_{11}$ and $H_{02}$. Based on Equation (3.8), Equation (3.15), Equation (3.17), one gets

$\begin{array}{l}{k}_{10}=f_{1001}^{1,1},\\ k_{21}=\frac{i{f}_{2000}^{\left(1,1\right)}{f}_{1100}^{\left(1,1\right)}}{\omega }+f_{1010}^{\left(1\right)}{H}_{11}\left(0\right)+f_{0110}^{\left(1\right)}{H}_{20}\left(0\right)+h_{1010}^{\left(1\right)}{H}_{11}\left(-{\tau }_{1c}\right)\\ +h_{0110}^{\left(1\right)}{H}_{20}\left(-{\tau }_{1c}\right)+f_{2100}^{\left(1,1\right)}\end{array}$ (3.25)

By calculation, the coefficients in Equation (3.21) is derived as

$\begin{array}{c}{f}_{1001}^{1,1}={a^{\prime }}_{11}\left(k_c\right)k_e{q}^{\left(1\right)}\left(0\right){\bar{q}}^{*\left(1\right)}\left(0\right)+{a^{\prime }}_{12}\left(k_c\right)k_e{q}^{\left(2\right)}\left(0\right){\bar{q}}^{*\left(1\right)}\left(0\right)\\ -{\bar{q}}^{*\left(2\right)}\left(0\right){b^{\prime }}_1\left(k_c\right)k_e{q}^{\left(1\right)}\left(0\right){\text{e}}^{-i\omega {\tau }_{1c}}-{\bar{q}}^{*\left(2\right)}\left(0\right){b^{\prime }}_2\left(k_c\right)k_e{q}^{\left(2\right)}\left(0\right){\text{e}}^{-i\omega {\tau }_{1c}}\\ +{\bar{q}}^{*\left(2\right)}\left(0\right){b^{\prime }}_2\left(k_c\right)k_e{q}^{\left(2\right)}\left(0\right)+{\bar{q}}^{*\left(3\right)}\left(0\right)k_e{q}^{\left(2\right)}\left(0\right)\\ -i{\bar{q}}^{\left(1\right)}\left(0\right)r{q}^{\left(1\right)}\left(0\right)\omega {\tau }_e{\text{e}}^{-i\omega {\tau }_{1c}}-i{\bar{q}}^{*\left(2\right)}\left(0\right)b_1{q}^{\left(1\right)}\left(0\right)\omega {\tau }_e{\text{e}}^{-i\omega {\tau }_{1c}}\\ -i{\bar{q}}^{*\left(2\right)}\left(0\right)b_2{q}^{\left(2\right)}\left(0\right)\omega {\tau }_e{\text{e}}^{-i\omega {\tau }_{1c}}\end{array}$

and

$\begin{array}{c}{f}_{2000}^{1,1}={\bar{q}}^{*\left(1\right)}\left(0\right)c_1{q}^{\left(1\right)}{\left(0\right)}^2+{\bar{q}}^{*\left(2\right)}\left(0\right)c_2{q}^{\left(1\right)}{\left(0\right)}^2{\text{e}}^{-2i\omega {\tau }_{1c}}\\ +{\bar{q}}^{*\left(2\right)}\left(0\right)c_3{q}^{\left(1\right)}\left(0\right)q^{\left(2\right)}\left(0\right){\text{e}}^{-2i\omega {\tau }_{1c}},\end{array}$

$\begin{array}{c}{f}_{1100}^{1,1}=2{\bar{q}}^{\left(1\right)}\left(0\right)c_1{q}^{\left(1\right)}\left(0\right){\bar{q}}^{*\left(1\right)}\left(0\right)+2{\bar{q}}^{\left(1\right)}\left(0\right)c_2{q}^{\left(1\right)}\left(0\right){\bar{q}}^{*\left(2\right)}\left(0\right)\\ +{\bar{q}}^{\left(1\right)}\left(0\right)c_3{q}^{\left(2\right)}\left(0\right){\bar{q}}^{*\left(2\right)}\left(0\right)+{\bar{q}}^{\left(2\right)}\left(0\right)c_3{q}^{\left(1\right)}\left(0\right){\bar{q}}^{*\left(2\right)}\left(0\right),\end{array}$

$\begin{array}{c}{f}_{0200}^{1,1}={\bar{q}}^{*\left(1\right)}\left(0\right)c_1{\bar{q}}^{\left(1\right)}{\left(0\right)}^2+{\bar{q}}^{*\left(2\right)}\left(0\right)c_2{\bar{q}}^{\left(1\right)}{\left(0\right)}^2{\text{e}}^{2i\omega {\tau }_{1c}}\\ +{\bar{q}}^{*\left(2\right)}\left(0\right)c_3{\bar{q}}^{\left(1\right)}\left(0\right){\text{e}}^{2i\omega {\tau }_{1c}}{\bar{q}}^{\left(2\right)}\left(0\right),\end{array}$

$f_{1010}=\left(\begin{array}{c}2{\bar{q}}^{*\left(1\right)}\left(0\right)c_1{q}^{\left(1\right)}\left(0\right)\\ 0\\ 0\end{array}\right),f_{0110}=\left(\begin{array}{c}2{\bar{q}}^{\left(1\right)}\left(0\right)c_1{\bar{q}}^{*\left(1\right)}\left(0\right)\\ 0\\ 0\end{array}\right),$

$h_{1010}=\left(\begin{array}{c}2{\bar{q}}^{*\left(2\right)}\left(0\right)c_2{q}^{\left(1\right)}\left(0\right){\text{e}}^{-i\omega {\tau }_{1c}}+{\bar{q}}^{*\left(2\right)}\left(0\right)c_3{q}^{\left(2\right)}\left(0\right){\text{e}}^{-i\omega {\tau }_{1c}}\\ {\bar{q}}^{*\left(2\right)}\left(0\right)c_3{q}^{\left(1\right)}\left(0\right){\text{e}}^{-i\omega {\tau }_{1c}}\\ 0\end{array}\right),$

$h_{0110}=\left(\begin{array}{c}2{\bar{q}}^{*\left(2\right)}\left(0\right)c_2{\bar{q}}^{\left(1\right)}\left(0\right){\text{e}}^{i\omega {\tau }_{1c}}+{\bar{q}}^{*\left(2\right)}\left(0\right)c_3{\bar{q}}^{\left(2\right)}\left(0\right){\text{e}}^{i\omega {\tau }_{1c}}\\ {\bar{q}}^{*\left(2\right)}\left(0\right)c_3{\bar{q}}^{\left(1\right)}\left(0\right){\text{e}}^{i\omega {\tau }_{1c}}\\ 0\end{array}\right),$

$\begin{array}{c}{f}_{21000}^{1,1}=3{\bar{q}}^{\left(1\right)}\left(0\right)c_4{q}^{\left(1\right)}{\left(0\right)}^2{\bar{q}}^{*\left(1\right)}\left(0\right)+2{\bar{q}}^{\left(1\right)}(0)c_5{q}^{\left(1\right)}\left(0\right)q^{\left(2\right)}(0){\bar{q}}^{*\left(1\right)}\left(0\right)\\ +{\bar{q}}^{\left(2\right)}\left(0\right)c_5{q}^{\left(1\right)}{\left(0\right)}^2{\bar{q}}^{*\left(1\right)}\left(0\right)+3{\bar{q}}^{*\left(2\right)}\left(0\right)c_6{\bar{q}}^{\left(1\right)}\left(0\right)q^{\left(1\right)}{\left(0\right)}^2{\text{e}}^{-i\omega {\tau }_{1c}}\\ +2{\bar{q}}^{*\left(2\right)}\left(0\right)c_7{\bar{q}}^{\left(1\right)}\left(0\right)q^{\left(1\right)}\left(0\right){\text{e}}^{-i\omega {\tau }_{1c}}{q}^{\left(2\right)}\left(0\right)+{\bar{q}}^{*\left(2\right)}\left(0\right)c_7{q}^{\left(1\right)}{\left(0\right)}^2{\text{e}}^{-i\omega {\tau }_{1c}}{\bar{q}}^{\left(2\right)}(\; 0\; )\end{array}$

By the above analysis, we conclude that

Theorem 3.1 There is a periodical solution with small amplitude bifurcating from Hopf point $\left({\tau }_{1c}\mathrm{,}{k}_c\right)$ which is stable if the first Lyapunov exponent $l_1\left(0\right)=\Re \left(k_{21}\right)<0$, otherwise unstable if $l_1\left(0\right)=\Re \left(k_{21}\right)>0$ and the bifurcating direction is determined by the signature of $\mu =\frac{\Re \left(k_{21}\right)}{\Re \left\{{\frac{\text{d}\lambda }{\text{d}\tau }|}_{\tau ={\tau }_{1c}}\right\}}$, which is supercritical if $\mu <0$ or subcritical if $\mu >0$.

4. Conclusion

The stability property of the disease infectious equilibrium solution of a type of HIV mathematical model with delay feedback control was investigated by varying parameter pairs $\left(k\mathrm{,}{\tau }_1\right)$ on parameter space. The Hopf bifurcation lines were tracked via using geometrical criterion of DDEs with multiple time delays. By using Schmidt dimensional reduction method combined with center manifold analytical technique, the universal norm form was computed near Hopf point. As period time delay ${\tau }_1$ is prolonged, stability of the disease infectious equilibrium solution loss and the stable periodical oscillation solutions arise near Hopf point.

Conflicts of Interest

The authors declare no conflicts of interest.

References

[1]	Korobeinikov, A. (2004) Global Properties of Basic Virus Dynamics Models. Bulletin of Mathematical Biology, 66, 879-883. https://doi.org/10.1016/j.bulm.2004.02.001
[2]	Korobeinikov, A. (2007) Global Properties of Infectious Disease Models with Nonlinear Incidence. Bulletin of Mathematical Biology, 69, 1871-1886. https://doi.org/10.1007/s11538-007-9196-y
[3]	Li, M.Y. and Shu, H. (2010) Global Dynamics of an In-Host Viral Model with Intracellular Delay. Bulletin of Mathematical Biology, 72, 1492-1505. https://doi.org/10.1007/s11538-010-9503-x
[4]	Wang, K.F., Fan, A.J. and Torres, A. (2010) Global Properties of an Improved Hepatitis B Virus Model. Nonlinear Analysis: Real World Applications, 11, 3131-3138. https://doi.org/10.1016/j.nonrwa.2009.11.008
[5]	Perelson, A.S., Kirschner, D.E. and De Boer, R. (1993) Dynamics of HIV Infection of CD4+ T-Cells. Mathematical Biosciences, 114, 81-125. https://doi.org/10.1016/0025-5564(93)90043-A
[6]	Herz, A.V.M., Bonhoeffer, S., Anderson,R.M., May,R.M. and Nowak,M.A. (1996) Viral Dynamics in Vivo: Limitations on Estimates of Intracellular Delay and Virus Decay. Proceedings of the National Academy of Sciences of the United States of America, 93, 7247-7251. https://doi.org/10.1073/pnas.93.14.7247
[7]	Huang, G., Ma, W. and Takeuchi, Y. (2009) Global Properties for Virus Dynamics Model with Beddington-DeAngelis Functional Response. Applied Mathematics Letters, 22, 1690-1693. https://doi.org/10.1016/j.aml.2009.06.004
[8]	Srivastava, P.K. and Chandra, P. (2010) Modeling the Dynamics of HIV and CD4+ T Cells during Primary Infection. Nonlinear Analysis: Real World Applications, 11, 612-618. https://doi.org/10.1016/j.nonrwa.2008.10.037
[9]	Guo, S.B. and Ma, W.B. (2016) Global Behavior of Delay Differential Equations Model of HIV Infection with Apoptosis. Discrete and Continuous Dynamical Systems: Series B, 21, 103-119. https://doi.org/10.3934/dcdsb.2016.21.103
[10]	Song, X. and Neumann, A.U. (2007) Global Stability and Periodic Solution of the Viral Dynamics. Journal of Mathematical Analysis and Applications, 329, 281-297. https://doi.org/10.1016/j.jmaa.2006.06.064
[11]	Huang, G., Yokoi, H., Takeuchi, Y., Kajiwara, T. and Sasaki, T. (2011) Impact of Intracellular Delay, Immune Activation Delay and Nonlinear Incidence on Viral Dynamics. Japan Journal of Industrial and Applied Mathematics, 28, 383-411. https://doi.org/10.1007/s13160-011-0045-x
[12]	Hirsch, M.W., Smith, H.L. and Zhao, X.-Q. (2001) Chain Transitivity, Attractivity, and Strong Repellors for Semidynamical Systems. Journal of Dynamics and Differential Equations, 13, 107-131. https://doi.org/10.1023/A:1009044515567
[13]	De Boer, R.J. and Perelson, A.S. (1998) Target Cell Limited and Immune Control Models of HIV Infection: A Comparison. Journal of Theoretical Biology, 19, 201-214. https://doi.org/10.1006/jtbi.1997.0548
[14]	Kepler, T.B. and Perelson, A.S. (1993) Cyclic Re-Entry of Germinal Center B Cells and the Efficiency of Affinity Maturation. Immunology Today, 14, 412-415. https://doi.org/10.1016/0167-5699(93)90145-B
[15]	Kuang, Y. (1993) Delay Differential Equations with Applications in Population Dynamics. Academic Press, Inc., Boston.
[16]	Engelborghs, K., Luzynina, T. and Roose, D. (2002) Numerical Bifurcation Analysis of Delay Differential Equations Using DDE-BIFTOOL. ACM Transactions on Mathematical Software, 28, 1-21. https://doi.org/10.1145/513001.513002
[17]	Beretta, E. and Kuang, Y. (2002) Geometric Stability Switch Criteria in Delay Differential Systems with Delay Dependant Parameters. SIAM Journal on Mathematical Analysis, 33, 1144-1165. https://doi.org/10.1137/S0036141000376086
[18]	Ma, S.Q., Lu, Q.S. and Mei, S.L. (2005) Dynamics of a Logistic Population Model with Maturation Delay and Nonlinear Birth Rate. Discrete and Continuous Dynamical Systems: Series B, 5, 736-752. https://doi.org/10.3934/dcdsb.2005.5.735
[19]	Ma, S.Q. (2018) Application of Extended Geometrical Criterion to Population Model with Two Time Delays. Journal of Mathematical Researches, 10, 63-76. https://doi.org/10.5539/jmr.v10n3p63
[20]	Ma, S.Q., Feng, Z.S. and Lu, Q.S. (2008) A Two-Parameter Geometrical Criteria for Delay Differential Equations. Discrete and Continuous Dynamical Systems: Series B, 9, 397-413. https://doi.org/10.3934/dcdsb.2008.9.397
[21]	Hale, J.K. and Verduyn Lunel, S.M. (1993) Introduction to Functional Differential Equations. Springer-Verlag, New York. https://doi.org/10.1007/978-1-4612-4342-7_1
[22]	Song, Z.-G. and Xu, J. (2013) Stability Switches and Double Hopf Bifurcation in a Two Neural Network System with Multiple Delays. Cognitive Neurodynamics, 7, 505-521. https://doi.org/10.1007/s11571-013-9254-0