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