\documentclass[reqno]{amsart} \usepackage{hyperref} \usepackage{graphicx} \AtBeginDocument{{\noindent\small \emph{Electronic Journal of Differential Equations}, Vol. 2014 (2014), No. 27, pp. 1--13.\newline ISSN: 1072-6691. URL: http://ejde.math.txstate.edu or http://ejde.math.unt.edu \newline ftp ejde.math.txstate.edu} \thanks{\copyright 2013 Texas State University - San Marcos.} \vspace{9mm}} \begin{document} \title[\hfilneg EJDE-2014/27\hfil Hopf bifurcation] {Hopf bifurcation for tumor-immune competition systems with delay} \author[P. Bi, H. Xiao \hfil EJDE-2014/27\hfilneg] {Ping Bi, Heying Xiao} % in alphabetical order \address{Ping Bi \newline Department of Mathematics, Shanghai Key Laboratory of PMMP, East China Normal University, 500 Dongchuan RD, Shanghai 200241, China} \email{pbi@math.ecnu.edu.cn} \address{Heying Xiao \newline Department of Mathematics, Shanghai Key Laboratory of PMMP, East China Normal University, 500 Dongchuan RD, Shanghai 200241, China} \email{hawcajok701@163.com} \thanks{Submitted October 9, 2013. Published January 14, 2014.} \subjclass[2000]{34K18, 34K60, 92C37} \keywords{Tumor-immune system; Hopf bifurcation; delay} \begin{abstract} In this article, a immune response system with delay is considered, which consists of two-dimensional nonlinear differential equations. The main purpose of this paper is to explore the Hopf bifurcation of a immune response system with delay. The general formula of the direction, the estimation formula of period and stability of bifurcated periodic solution are also given. Especially, the conditions of the global existence of periodic solutions bifurcating from Hopf bifurcations are given. Numerical simulations are carried out to illustrate the the theoretical analysis and the obtained results. \end{abstract} \maketitle \numberwithin{equation}{section} \newtheorem{theorem}{Theorem}[section] \newtheorem{lemma}[theorem]{Lemma} \allowdisplaybreaks \section{Introduction} Recently, there has been much interest in mathematical model describing the interaction between tumor cells and effector cells of the immune system, see, for example, Bi et al \cite{ Bi}, d'Onofrio et al \cite{d'Onofrio2, d'Onofrio3, d'Onofrio4}, and Yafia et al \cite{RYafia7.1}. An ideal model can provide insights into the dynamics of interactions of the immune response with the tumor and may play a significant role in understanding of the corresponding cancer and developing effective drug therapy strategies against it. However, it is almost impossible to develop realistic models to describe such complex processes. In fact, mathematical models for the dynamics of the interaction of the immune components with tumor cells are very idealized. Thus it is feasible to propose simple models which are capable to display some of the essential immunological phenomena. Recently, delayed models of tumor and immune response interactions have been studied extensively, we refer to Bi et al \cite{Bi, xiao}, Galach \cite{MGalach6}, Mayer \cite{Mayer}, Yafia \cite{RYafia7.2,RYafia7.3,RYafia7.4} and the references cited therein. It would be interesting to consider the nonlinear dynamics of the delayed model. In 1994, Kuznetsov et al \cite{Kuznetsov} took into account the penetration of TCs by ECs, and proposed a model describing the response of effector cells (ECs) to the growth of tumor cells (TCs). They assumed that interactions between ECs and TCs in vitro can be described by the kinetic scheme shown in Figure \ref{fig1}, where $E, T, C, T^\ast, E^\ast$ are the local concentrations of ECs, TCs, EC-TC complexes, inactivated ECs, and lethally hit TCs, respectively. From the above kinetic scheme and the experimental observations, Kuznetsov et al \cite{Kuznetsov} claimed that the model can be reduced to two equations which describe the behavior of ECs and TCs only. Taking time scale as that in \cite{MGalach6} and \cite{xiao}, then Kuznetsov, Makalkin and Taylor's model can be written as \begin{equation}\label{80} \begin{gathered} \frac{\mathrm{d}x}{\mathrm{d}t}=\sigma +\zeta xy-\delta x, \\ \frac{\mathrm{d}y}{\mathrm{d}t}=\alpha y(1-\beta y)-xy, \end{gathered} \end{equation} where $x$ denotes the dimensionless density of ECs, $y$ stands for the dimensionless density of the population of TCs. All coefficients are positive, where $\zeta>0$ shows the stimulation coefficient of the immune system exceeds the neutralization coefficient of ECs in the process of the formation of EC-TC complexes, which is the most biologically meaningful. \begin{figure}[ht] \begin{center} \includegraphics[width=0.7\textwidth]{fig1} %Fig1.1.eps \end{center} \caption{Kinetic scheme describing interactions between ECs and TCs} \label{fig1} \end{figure} Yafia \cite{RYafia7.1} considered the model \eqref{80} and studied the linearizing stability of the equilibrium and the existence of the Hopf bifurcation. In \cite{MGalach6, RYafia7.2,RYafia7.3}, the authors obtained the same results as \cite{RYafia7.1} for the model \eqref{80} with delay $\tau$, as follows \begin{equation}\label{5} \begin{gathered} \frac{\mathrm{d}x}{\mathrm{d}t}=\sigma +\zeta x(t-\tau)y(t-\tau)-\delta x, \\ \frac{\mathrm{d}y}{\mathrm{d}t}=\alpha y(1-\beta y)-xy, \end{gathered} \end{equation} The rest of this paper is organized as follows. In section 2, we study the tumor model with delay \eqref{5}, and show the properties of the Hopf bifurcated periodic solutions of this system, including the direction of Hopf bifurcation, the period and stability of bifurcated periodic solutions. The numerical analysis and simulations are also given to illustrate the main results. Section 3 is devoted to the existence of the global Hopf bifurcation. A brief discussion is also given in this section. \section{Direction and stability of Hopf bifurcation} If $\alpha\delta>\sigma$, equations \eqref{80} and \eqref{5} have two possible nonnegative steady states $P_0(\frac{\sigma}{\delta}, 0)$ and $P_2(\frac{-\alpha(\beta\delta-\zeta)+\sqrt{\Delta}}{2\zeta}, \frac{\alpha(\beta\delta+\zeta)-\sqrt{\Delta}}{2\alpha\beta\zeta})$, where $\Delta=\alpha^2(\beta\delta-\zeta)^2+4\alpha\beta\zeta\sigma>0$. Galach\cite{MGalach6} and Yafia \cite{RYafia7.2} shows the stability of the equilibria and existence of the Hopf bifurcation for system \eqref{80} and \eqref{5}, the main results are summarized as follows. \begin{lemma}[\cite{MGalach6}] \label{lem2.1} If the point $P_2$ exists and has nonnegative coordinates, then it is stable. And there is no nonnegative periodic solution to equation \eqref{80}. \end{lemma} \begin{lemma}[\cite{RYafia7.2}] \label{lem2.2} Let $ \alpha\delta>\sigma$. (1) The equilibrium point $P_0$ is asymptotically stable for all $\tau>0$. (2) If $\alpha\beta>\zeta$, then there exists $\tau_l > 0$ such that $P_2$ is asymptotically stable for $\tau<\tau_l$ and unstable for $\tau>\tau_l$; (3) If $\alpha\beta>\zeta$, then there exists $\varepsilon_1>0$ such that, for each $0<\epsilon<\varepsilon_1$, system \eqref{5} has a family of periodic solutions $p_l(\varepsilon)$ with period $T_l=T_l(\varepsilon)$, for the parameter values $\tau=\tau(\varepsilon)$ satisfy $p_l(0)=P_2, T_l(0)=\frac{2\pi}{\omega_l}, \tau(0)=\tau_l$. Where \begin{gather*} \tau_l=\begin{cases} \frac{1}{\omega} \arccos\frac{q(\omega^2-r)-ps\omega^2} {s^2\omega^2+p^2}, & \text{if } s(r-\omega^2)-pq>0,\\ \frac{1}{\omega}\arccos\frac{q(\omega^2-r)-ps\omega^2} {s^2\omega^2+p^2}+\pi, & \text{if } s(r-\omega^2)-pq<0, \end{cases} \\ \omega^2=\frac{1}{2}(s^2-p^2+2r) +\frac{1}{2} \sqrt{(s^2-p^2+2r)^2-4(r^2-q^2)}. \end{gather*} \end{lemma} Let $\tau_j=\tau_l+\frac{2j\pi}{\omega_l},j=0,1,\dots$. Let $P_2(x_2, y_2)$ be the positive equilibrium, where $x_2=\frac{-\alpha(\beta\delta-\zeta)+\sqrt{\Delta}}{2\zeta}, y_2=\frac{\alpha(\beta\delta+\zeta)-\sqrt{\Delta}}{2\alpha\beta\zeta}$. From lemma \ref{lem2.1} we know that if $0<\frac{\zeta}{\beta}<\alpha, \alpha\delta>\sigma$, then \eqref{6.2} undergos Hopf bifurcation at the equilibrium $P_2(x_2, y_2)$, the corresponding purely imaginary roots are $\lambda=\pm i\tau_j\omega$, and the critical values are $\tau=\tau_j,j=0,1,\dots$. As pointed out by Hassard et al \cite{Hassard}, it is interesting to determine the direction, stability and period of these periodic solutions bifurcating from the steady state. In this section, we will study the stability and direction of the Hopf bifurcated periodic solution by using the center manifold reduction and normal form theory of retarded functional differential equations due to the ideals of Faria and Magalha\'es \cite{T.Faria8.1, T.Faria8.2}. Throughout this section, we assume that system \eqref{6.1} undergoes Hopf bifurcations at the equilibrium $P_2$ as the critical parameter $\tau=\tau_j$, $j=1,2,3\dots$ and the corresponding purely imaginary roots are $\pm i\omega$. Set $z_1(t)=x(t)-x_2, z_2(t)=y(t)-y_2$. Then system \eqref{5} can be written as \begin{equation}\label{6.1} \begin{gathered} z_1'(t)=\alpha_1z_1(t)+\alpha_2z_1(t-\tau)+\alpha_3z_2(t-\tau) +\zeta z_1(t-\tau)z_2(t-\tau), \\ z_2'(t)=\beta_1z_1(t)+\beta_2z_2(t) -\alpha\beta z_2^2(t)- z_1(t)z_2(t), \end{gathered} \end{equation} where $\alpha_1=-\delta<0$, $\alpha_2=\zeta y_2>0$, $\alpha_3=\zeta x_2>0$, $\beta_1=-y_2<0$, $\beta_2=\alpha-2\alpha\beta y_2-x_2=-\alpha\beta y_2<0$. Normalizing the delay $\tau$ in system \eqref{6.1} by the time-scaling $ t\to \frac{t}{\tau}$, then \eqref{6.1} is transformed into \begin{equation}\label{6.2} \begin{gathered} z_1'(t)=\tau(\alpha_1z_1(t)+\alpha_2z_1(t-1)+\alpha_3z_2(t-1) +\zeta z_1(t-1)z_2(t-1)), \\ z_2'(t)=\tau(\beta_1z_1(t)+\beta_2z_2(t) -\alpha\beta z_2^2(t) - z_1(t)z_2(t)). \end{gathered} \end{equation} Let $z(t)=(z_1(t), z_2(t))^T$. We transformed \eqref{6.2} into the FDE \begin{equation}\label{6.3} \dot{z}(t)=N(\tau)(z_t)+F(z_t, \tau), \end{equation} where $N(\varphi):C([-1,0],\mathbb{R}^2)\to\mathbb{R}^2,F(\varphi):C([-1,0], \mathbb{R}^2)\to\mathbb{R}^2, \varphi =\mathrm{col}(\varphi_1,\varphi_2) \in C([-1,0],\mathbb{R}^2)$ satisfy \begin{gather*} N(\tau)(\varphi)=\tau \begin{pmatrix} \alpha_1\varphi_1(0)+\alpha_2\varphi_1(-1) +\alpha_3\varphi_2(-1)\\ \beta_1\varphi_1(-1)+\beta_2\varphi_2(0) \end{pmatrix}, \\ F(\varphi, \tau) =\tau\begin{pmatrix} \zeta\varphi_1(-1)\varphi_2(-1)\\ -\alpha\beta\varphi^2_2(0)-\varphi_1(0)\varphi_2(0) \end{pmatrix}. \end{gather*} Setting the new parameter $\gamma=\tau-\tau_j$, then \eqref{6.3} can be written as \begin{equation}\label{6.4} \dot{z}(t)=N(\tau_j)(z_t)+\tilde{F}(z_t, \gamma), \end{equation} where $\tilde{F}(z_t, \gamma)=N(\gamma)(z_t)+F(z_t, \tau_j+\gamma)$. Let $\Lambda=\{i\omega, -i\omega\}$. Assuming $A$ is the infinitesimal generator of $\dot{z}(t)=N(\tau_k)(z_t)$ satisfying $A\Phi=\Phi B$ with \begin{equation}\label{35} B=\begin{pmatrix} i\omega & 0 \\ 0 & -i\omega \end{pmatrix}, \end{equation} and $A$ has a pair of conjugate purely imaginary roots $\pm i\omega$. Denote $P$ is the invariant space of $A$ associated with $\Lambda$, then $\dim P=2$. We can decompose $C:=C([-1,0],\mathbb{R}^2)$ as $ C=P\bigoplus Q$ by the formal adjoint theory for FDEs in \cite{J.K.Hale1}. Considering complex coordinates, we still denote $C([-1,0],\mathbb{C}^2)$ as $C$. Let $\Phi=(\Phi_1, \Phi_2)$ by the bases for $P$, where \begin{equation}\label{31} \Phi_1=e^{i\omega\theta}v, \quad \Phi_2=\bar{\Phi}_1, \quad \theta \in[-1,0], \end{equation} where $v= ( v_1, v_2)^T$ is a vector in $\mathbb{C}^2$ that satisfies $N(\tau_j)\Phi_1=i\omega v.$ Choose a basis $\Psi$ for the adjoint space $P^\ast$, where $\Psi=col(\Psi_1, \Psi_2)$, \begin{equation} \Psi_1=e^{-i\omega\tilde{\theta}}u^T, \quad \Psi_2=\bar{\Psi}_1, \quad u=\begin{pmatrix} u_1\\u_2\end{pmatrix}, \quad \tilde{\theta} \in[0,1]. \end{equation} Define $(\Psi,\Phi)=((\Psi_i,\Phi_j))_{i,j=1,2}$, \begin{equation} (\psi,\varphi)=\psi(0)\varphi(0)-\int_{-1}^{0}\int_{0}^{\theta} \psi(s-\theta)\mathrm{d}\eta(\theta)\varphi(s)\mathrm{d}s, \forall\varphi \in P,\psi \in P^\ast. \end{equation} Then $(\Psi,\Phi)$ can be transformed to the identify matrix $I_2$. Thus we can ensure \begin{equation} v=\begin{pmatrix} 1\\\frac{i\omega_j-(\alpha_1+\alpha_2e^{-i\omega_j}) \tau_j}{\tau_j\alpha_3e^{-i\omega_j}}\end{pmatrix}, \quad u=u_1\begin{pmatrix}1\\ \frac{i\omega_j-(\alpha_1+\alpha_2e^{-i\omega_j})\tau_j}{\tau_j\beta_1} \end{pmatrix}, \end{equation} with \[ \frac{1}{u_1}=1+v_2\frac{i\omega_j-(\alpha_1+\alpha_2e^{-i\omega_j}) \tau_j}{\tau_j\beta_1}+(\alpha_2+\alpha_3v_2)e^{-i\omega_j}. \] Define the enlarged phase space $BC$ as $$ BC:=\{\varphi:[-1,0]\to\mathbb{C}^2|\varphi \text{ is continuous on $[-1,0)$, $\lim_{\theta\to 0^-}\varphi(\theta)$ exists}\}. $$ The projection $\pi:BC\to P$ is defined as $\pi(\varphi+X_0 b)=\Phi[(\Psi,\varphi)+\Psi(0)b]$, for each $\varphi\in C$ and $b\in\mathbb{R}^2$, thus we have the decomposition $BC=P\bigoplus \ker \pi$. Let $z_t=\Phi x+y, x\in\mathbb{C}^2, y\in \ker (\pi)\cap C^1:=Q^1$, we can decompose \eqref{6.4} as \begin{equation}\label{6.6} \begin{gathered} \dot{x}=Bx+\Psi(0)\tilde{F}(\Phi x+y, \gamma),\\ \frac{\mathrm{d}y}{\mathrm{d}x}=A_{Q^1}y+(I-\pi)X_0\tilde{F}(\Phi x+y,\gamma), \end{gathered} \end{equation} where \begin{equation} X_0(\theta)=\begin{cases} I, &\theta=0;\\ 0, &-1\le\theta<0. \end{cases} \end{equation} We write the Taylor expansion \begin{gather*} \Psi(0)\tilde{F}(\Phi x+y, \gamma) =\frac{1}{2}f_2^1(x,y,\gamma)+\frac{1}{3!}f_3^1(x,y,\gamma)+\mathrm{h.o.t.},\\ (I-\pi)X_0\tilde{F}(\Phi x+y,\gamma)=\frac{1}{2}f_2^2(x,y,\gamma) +\frac{1}{3!}f_3^2(x,y,\gamma)+\mathrm{ h.o.t.}, \end{gather*} where $f_k^1,f_k^2$ are homogeneous polynomials in $x,y,\gamma$ of degree $k,k=2,3$, with coefficients in $\mathbb{C}^2$ and $\ker \pi$, respectively, and h.o.t. stands for higher order terms. The normal form method implies a normal form on the center manifold of the origin for \eqref{6.4} is \begin{equation}\label{6.12} \dot{x}=Bx+\frac{1}{2}g_2^1(x,0,\gamma)+\frac{1}{3!}g_3^1(x,0,\gamma) +\mathrm{h.o.t.}. \end{equation} where $g_2^1(x,0,\gamma),g_3^1(x,0,\gamma)$ are homogeneous polynomials in $x,\gamma$, and h.o.t. are the higher order terms. We follow the notation in \cite{T.Faria8.1}; that is, $$ V_j^{m+p}(X)=\left\{\Sigma_{|(q,l)|=j}C_{(q,l)}x^q\gamma^l:(q,l) \in \mathbb{N}_0^{m+p},C_{(q,l)}\in X\right\}, $$ with $x=(x_1,x_2,\dots,x_m),\gamma=(\gamma_1,\gamma_2,\dots,\gamma_p)$; \begin{equation} \begin{gathered} M_j(p,h)=(M_j^1p,M_j^2h),\\ M_j^1p(x,\gamma)=[B,p(\cdot,\gamma)](x)=D_xp(x,\gamma)Bx-Bp(x,\gamma),\\ M_j^2h(x,\gamma)=D_xh(x,\gamma)Bx-A_{Q^1}h(x,\gamma). \end{gathered} \end{equation} Since $V_j^3$ satisfies $V_j^3(\mathbb{C}^2)=\operatorname{Im}(M_j^1)\oplus\ker (M_j^1)$ and $$ \ker (M_j^1)=\operatorname{span}\big\{x^q\gamma^le_k:(q,\overline{\lambda}) =\lambda_j,\;\lambda_j\in \Lambda,\; j=1,2,q\in\mathbb{N}_0^2,\; l\in \mathbb{N}_0,\; |(q,l)|=j\big\}, $$ $e_1,e_2$ is the base of $\mathbb{C}^2$. Then we can obtain \begin{gather*} \ker (M_2^1)=\operatorname{span} \Big\{\begin{pmatrix}x_1\gamma\\ 0\end{pmatrix},\begin{pmatrix}0\\x_2\gamma\end{pmatrix}\Big\}, \\ \ker (M_3^1)=\operatorname{span}\Big\{\begin{pmatrix}x_1^2x_2\\0\end{pmatrix},\begin{pmatrix}x_1\gamma^2\\0\end{pmatrix}, \begin{pmatrix}0\\x_1x_2^2\end{pmatrix}, \begin{pmatrix}0\\x_2\gamma^2\end{pmatrix}\Big\}. \end{gather*} Noting \eqref{6.6}, one has $$ f_2^1(x,0,\gamma)=2\Psi(0)[N(\gamma)(\Phi x)+F(\Phi x,\tau_j)]; $$ i.e., \begin{equation}\label{6.7} f_2^1(x,0,\gamma)=2\begin{pmatrix}A_1x_1\gamma+A_2x_2\gamma+a_{20}x_1^2 +a_{11}x_1x_2+a_{02}x_2^2\\ \bar{A}_1x_2\gamma+\bar{A}_2x_1\gamma+\bar{a}_{02}x_1^2+\bar{a}_{11}x_1x_2 +\bar{a}_{20}x_2^2\end{pmatrix}, \end{equation} where \begin{gather*} A_1=\frac{i\omega_j}{\tau_j}u^Tv, A_2=\frac{-i\omega_j}{\tau_j}u^T\bar{v},\\ a_{20}=\tau_j[u_1e^{-2i\omega_j}\zeta v_1v_2+u_2(-\alpha\beta v_2^2-v_1v_2)],\\ a_{11}=\tau_j[\zeta u_1(v_1\bar{v}_2+\bar{v}_1v_2)+u_2(-2\alpha\beta v_2\bar{v}_2-(v_1\bar{v}_2+\bar{v}_1v_2)],\\ a_{02}=\tau_j[u_1e^{2i\omega_j}\zeta \bar{v}_1\bar{v}_2 +u_2(-\alpha\beta \bar{v}_2^2-\bar{v}_1\bar{v}_2)]. \end{gather*} Therefore, \begin{equation}\label{6.8} g_2^1(x,0,\gamma)=\operatorname{Proj}_{\ker (M_2^1)}f_2^1(x,0,\gamma) =\begin{pmatrix}2A_1x_1\gamma\\2\bar{A}_1x_2\gamma\end{pmatrix}. \end{equation} Since the terms $O(|x|\gamma^2)$ are irrelevant to determine the generic Hopf bifurcation, we assume $$ J=\operatorname{span}\Big\{\begin{pmatrix}x_1^2x_2\\0\end{pmatrix}, \begin{pmatrix}0\\x_1x_2^2\end{pmatrix}\Big\}, $$ then $$ g_3^1(x,0,\gamma)=\operatorname{Proj}_J\bar{f}_3^1(x,0,0)+o(|x|\gamma^2), $$ where $$ \bar{f}_3^1(x,0,0)=\frac{3}{2}[(D_xf_2^1)u_2^1-(D_xu_2^1) g_2^1]_{(x,0,0)}+\frac{3}{2}[(D_yf_2^1)u_2^2]_{(x,0,0)}. $$ Hence we will compute $g_3^1(x,0,\gamma)$ as follows. Firstly, noting \begin{gather*} f_2^1(x,0,0)=2\begin{pmatrix} a_{20}x_1^2+a_{11}x_1x_2+a_{02}x_2^2 \\\bar{a}_{02}x_1^2+\bar{a}_{11}x_1x_2+\bar{a}_{20}x_2^2 \end{pmatrix}, \\ u_2^1(x,0)=\frac{2}{i\omega_l}\begin{pmatrix} a_{20}x_1^2-a_{11}x_1x_2-\frac{1}{3}a_{02}x_2^2\\ \frac{1}{3}\bar{a}_{02}x_1^2+\bar{a}_{11}x_1x_2-\bar{a}_{20}x_2^2 \end{pmatrix}, \end{gather*} one has \begin{equation}\label{6.9} \begin{split} \quad \operatorname{Proj}_J[(D_xf_2^1)u_2^1]_{(x,0,0)} &=\frac{4}{i\omega_l}\begin{pmatrix}(-a_{20}a_{11}+\frac{2}{3}|a_{02}|^2 +|a_{11}|^2)x_1^2x_2\\ (-\frac{2}{3}|a_{02}|^2-|a_{11}|^2+\overline{a_{20}a_{11}})x_1x_2^2 \end{pmatrix}\\ &=4\begin{pmatrix}A_3x_1^2x_2\\ \bar{A}_3x_1x_2^2\end{pmatrix}. \end{split} \end{equation} Secondly, From \eqref{6.8}, we know $g_2^1(x,0,0)=0$, then $\operatorname{Proj}_J[(D_xu_2^1)g_2^1]_{(x,0,0)}=0$. Lastly, we will compute $\operatorname{Proj}_J[(D_yf_2^1)u_2^2]_{(x,0,0)}$ as follows. Let $$ h=u_2^2=h_{200}x_1^2+h_{020}x_2^2+h_{002}\gamma^2+ h_{110}x_1x_2+h_{101}x_1\gamma+h_{011}x_2\gamma. $$ Noting that $g_2^2=0, $ one has $$ M_2^2h(x,\gamma)=f_2^2=2(I-\pi)X_0\tilde{F}(\Phi x,\gamma) =2(I-\pi)X_0[N(\gamma)(\Phi x)+F(\Phi x,\tau_j)]. $$ On the other hand, we know \begin{equation} \begin{split} M_2^2h(x,\gamma) &=D_xh(x,\gamma)Bx-A_{Q^1}h(x,\gamma)\\ &=D_xh(x,\gamma)Bx-[\dot{h}(x,\gamma)+X_0(L(\tau_j)(h(x,\gamma))- \dot{h}(x,\gamma)(0))]. \end{split} \end{equation} If $\gamma=0$, then \begin{equation}\label{6.10} \begin{gathered} \dot{h}(x)-D_xh(x)Bx=2\Phi\Psi(0)F(\Phi x,\tau_j),\\ \dot{h}(x)(0)-L(\tau_j)(h(x))=2F(\Phi x,\tau_j). \end{gathered} \end{equation} Let \begin{gather*} W(\theta)=\Phi x+y=\Phi_1x_1+\Phi_2x_2+y(\theta)=e^{i\omega_l\theta}vx_1+ e^{-i\omega_l\theta}\bar{v}x_2+y(\theta), \\ \tilde{W}(\theta)=\Phi x=\Phi_1x_1+\Phi_2x_2=e^{i\omega_l\theta} vx_1+e^{-i\omega_l\theta}\bar{v}x_2. \end{gather*} From $$ f_2^1(x,y,0)=2\tau_j\begin{pmatrix}u^T\begin{pmatrix}\zeta W_1(-1)W_2(-1)\\ -\alpha\beta W_2^2(0)-W_1(-1)W_2(-1)\end{pmatrix}\\\bar{u}^T\begin{pmatrix}\zeta W_1(-1)W_2(-1)\\ -\alpha\beta W_2^2(0)-W_1(-1)W_2(-1)\end{pmatrix}\end{pmatrix}, $$ we obtain \begin{align*} &[(D_yf_2^1)h]_{(x,0,0)}\\ &=2\begin{pmatrix}\tau_j u^T\begin{pmatrix}\zeta \tilde{W}_2(-1)h^1(-1)+\zeta \tilde{W}_1(-1)h^2(-1)\\ -\tilde{W}_2(-1)h^1(-1)-\tilde{W}_1(-1)h^2(-1)-2\alpha\beta \tilde{W}_2(0)h^2(0)\end{pmatrix}\\ \tau_j\bar{u}^T\begin{pmatrix}\zeta \tilde{W}_2(-1)h^1(-1)+\zeta \tilde{W}_1(-1)h^2(-1)\\ -\tilde{W}_2(-1)h^1(-1)-\tilde{W}_1(-1)h^2(-1)-2\alpha\beta \tilde{W}_2(0)h^2(0)\end{pmatrix}\end{pmatrix} \end{align*} and \begin{equation}\label{6.11} \operatorname{Proj}_J[(D_yf_2^1)u_2^2]_{(x,0,0)} =2\begin{pmatrix}A_4x_1^2x_2\\\bar{A}_4x_1x_2^2\end{pmatrix}, \end{equation} where \begin{align*} A_4&=\tau_j\Big[u_1\zeta(e^{-i\omega_j}v_2h_{110}^1(-1) +e^{i\omega_j}\bar{v}_2h_{200}^1(-1) +e^{-i\omega_j}v_1h_{110}^2(-1)\\ &\quad +e^{i\omega_j}\bar{v}_1h_{200}^2(-1))\Big] +u_2\tau_j\Big[-e^{-i\omega_j}v_2h_{110}^1(-1) -e^{i\omega_j}\bar{v}_2h_{200}^1(-1)\\ &\quad -e^{-i\omega_k}v_1h_{110}^2(-1) -e^{i\omega_j}\bar{v}_1h_{200}^2(-1)\Big] -u_2\tau_j\big[2\alpha\beta(v_2h_{110}^2(0)+\bar{v}_2h_{200}^2(0))\big]. \end{align*} To compute $A_4$, we should get $h_{110}(\theta),h_{200}(\theta)$ firstly. From \eqref{6.10}, it follows \begin{equation}\label{zz} \begin{gathered} \dot{h}_{110}=2(\Phi_1,\Phi_2)\begin{pmatrix}a_{11}\\\bar{a}_{11}\end{pmatrix},\\ \dot{h}_{110}(0)-L(\tau_j)(h_{110})=\tau_j\begin{pmatrix}a_1\\b_1\end{pmatrix}, \end{gathered} \end{equation} and \begin{equation}\label{zz1} \begin{gathered} \dot{h}_{200}-2i\omega_j h_{200}=2(\Phi_1,\Phi_2)\begin{pmatrix}a_{20}\\\bar{a}_{02}\end{pmatrix},\\ \dot{h}_{200}(0)-L(\tau_j)(h_{200})=\tau_j\begin{pmatrix}a_2\\b_2\end{pmatrix}, \end{gathered} \end{equation} where $a_1=2[\zeta(v_1\bar{v}_2+\bar{v}_1v_2)]$, $b_1=2[-2\alpha\beta v_2\bar{v}_2-(v_1\bar{v}_2+\bar{v}_1v_2)]$, $a_2=2[\zeta v_1v_2e^{-2i\omega_k}]$, $b_2=2[-\alpha\beta v_2^2-e^{-2i\omega_j}v_1v_2]$. Solving the above equations \eqref{zz} and \eqref{zz1}, we obtain \begin{gather*} h_{110}={2[\frac{a_{11}}{i\omega_k }\Phi_1-\frac{\bar{a}_{11}}{i\omega_j}\Phi_2]+C_1},\\ h_{200}={2[\frac{a_{20}}{-i\omega_k }\Phi_1+\frac{\bar{a}_{02}}{-3i\omega_j}\Phi_2]+C_2e^{2i\omega_k \theta}}, \end{gather*} where \begin{gather*} C_1=\begin{pmatrix}C_1^1\\C_1^2\end{pmatrix}, \quad C_1^1=\frac{\begin{vmatrix} a_1&-\alpha_3\\ b_1&-(\beta_2+\beta_3) \end{vmatrix}} {\begin{vmatrix} -(\alpha_1+\alpha_2)&-\alpha_3\\ -\beta_1&-(\beta_2+\beta_3)\end{vmatrix}}, \\ C_1^2=\frac{\begin{vmatrix} -(\alpha_1+\alpha_2)&a_1\\ -\beta_1&b_1 \end{vmatrix}} {\begin{vmatrix} -(\alpha_1+\alpha_2)&-\alpha_3\\ -\beta_1&-(\beta_2+\beta_3)\end{vmatrix}},\quad C_2=\begin{pmatrix}C_2^1\\C_2^2\end{pmatrix}, \\ C_2^1=\frac{\begin{vmatrix} \tau_j a_2&-\tau_j\alpha_3e^{-2i\omega_j}\\ \tau_j b_2&2i\omega_k +\tau_j\beta_2+\tau_j\beta_3e^{-2i\omega_j} \end{vmatrix}} {\begin{vmatrix} 2i\omega_k -\tau_j\alpha_1-\tau_j\alpha_2e^{-2i\omega_j}&-\tau_j\alpha_3e^{-2i\omega_k }\\ -\tau_j\beta_1e^{-2i\omega_j}&2i\omega_k +\tau_j\beta_2+\tau_j\beta_3e^{-2i\omega_j}\end{vmatrix}}, \\ C_2^2=\frac{\begin{vmatrix} 2i\omega_k -\tau_j\alpha_1-\tau_j\alpha_2e^{-2i\omega_k }&\tau_j a_2\\ -\tau_j\beta_1e^{-2i\omega_j}&\tau_j b_2 \end{vmatrix}} {\begin{vmatrix} 2i\omega_k -\tau_j\alpha_1-\tau_j\alpha_2e^{-2i\omega_j}&-\tau_j\alpha_3e^{-2i\omega_j}\\ -\tau_j\beta_1e^{-2i\omega_j}&2i\omega_j+\tau_j\beta_2 +\tau_j\beta_3e^{-2i\omega_j}\end{vmatrix}}. \end{gather*} Hence, $$ g_3^1(x,0,0)=\begin{pmatrix}(6A_3+3A_4)x_1^2x_2\\(6\bar{A}_3 +3\bar{A}_4)x_1x_2^2\end{pmatrix}. $$ Thus, the normal form of the system \eqref{6.12} has the form \begin{equation}\label{normal} \dot{x}=Bx+\begin{pmatrix}A_1x_1\gamma\\\bar{A}_1x_2\gamma\end{pmatrix}+ \frac{1}{3!}\begin{pmatrix}(6A_3+3A_4)x_1^2x_2\\(6\bar{A}_3 +3\bar{A}_4)x_1x_2^2\end{pmatrix}+o(|x|^4+|x|\gamma^2). \end{equation} Let $x_1=\xi_1-i\xi_2$, $x_2=\xi_1+i\xi_2$, $\xi_1=\rho\cos\omega$, $\xi_2=\rho\sin\omega$. Then system \eqref{normal} can be written as \begin{equation}\label{normal1} \begin{gathered} \dot{\rho}=r_1\gamma\rho+r_2\rho^3+O(\gamma^2\rho+|(\rho,\gamma)|^4),\\ \dot{\omega}=-\omega_j-\operatorname{Im}(A_1)\gamma-\operatorname{Im} (A_3+\frac{1}{2}A_4)\rho^2+o(|(\rho^2,\gamma)|), \end{gathered} \end{equation} where $r_1=\operatorname{Re}A_1,r_2=\operatorname{Re}(A_3+\frac{1}{2}A_4)$. Summarizing, we have the following theorem. \begin{theorem} \label{thm2.3} The flow on the center manifold of the equilibrium $P_2$ at $\gamma=0$ is given by \eqref{normal1}. And also we can draw the following conclusion. \begin{itemize} \item[(1)] The Hopf bifurcation is supercritical if $r_1r_2<0$, and subcritical if $r_1r_2>0$; \item[(2)] The nontrivial periodic solution is stable if $r_2<0$, and unstable if $r_2>0$; \item[(3)] The period of the nontrivial solution is \[ P(\gamma)=\frac{2\pi}{\omega_j}(1-\frac{\operatorname{Im}(A_1)\gamma- \frac{r_1\gamma}{r_2}\operatorname{Im}(A_3+\frac{1}{2}A_4)}{\omega_j}) +O(\gamma^3) \] with $P(0)=2\pi/\omega_j$. \end{itemize} \end{theorem} To explain the result of Theorem \ref{thm2.3}, we provide the simulations of Hopf bifurcation at $P_2$. In this paper, we cite the parameters in \cite{Kuznetsov} to illustrate Theorem \ref{thm2.3}. Assume $\sigma=0.1181$, $\zeta=0.0031$, $\delta=0.3743$, $\alpha=1.636$, $\beta=0.002$, then the system \eqref{5} has a Tumor-free equilibrium $P_0(0.3155, 0)$, which is asymptotically stable and a positive equilibrium $P_2(1.33435, 92.1911)$, which is locally stable. We only simulate local properties of the positive equilibrium $P_2(1.33435, 92.1911)$ here in the following Figure \ref{fig2.1} and Figure \ref{fig2.2}, and the corresponding critical value is $\tau_0=1.8760$. \begin{figure}[ht] \begin{center} \includegraphics[width=0.45\textwidth]{fig2a} \quad % Fig2.11.eps \includegraphics[width=0.45\textwidth]{fig2b} %Fig2.12.eps \end{center} \caption{(left) Stable equilibrium $(1.3344, 92.1911)$ when $\tau=1<\tau_0$; (right) oscillation of the solution to \eqref{5} with respect to $t$ when $\tau=1<\tau_0$} \label{fig2.1} \end{figure} \begin{figure}[ht] \begin{center} \includegraphics[width=0.45\textwidth]{fig3a} \quad % Fig2.21.eps \includegraphics[width=0.45\textwidth]{fig3b} %Fig2.22.eps \end{center} \caption{(left) Bifurcated periodic solution at the positive equilibrium when $\tau=2.0>\tau_0$; (right) oscillation solution $x(t)$ and $y(t)$ of \eqref{5} with respect to $t$ when $\tau=2.0>\tau_0$} \label{fig2.2} \end{figure} \section{Existence of global Hopf bifurcation} In section 2, we discussed the direction and stability of the Hopf bifurcated periodic solutions of system \eqref{5} at the equilibrium $P_2(x_2, y_2)$ as $\tau=\tau_j$. However, this bifurcated periodic solution exists in a neighborhood of the positive equilibrium, so we want to know whether the non-constant periodic solution exist globally. In this section, we will study the existence of global bifurcated periodic solution using a global Hopf bifurcation result due to Wu\cite{wu16}. Throughout this section, we follow the notations in \cite{wu16} and rewrite system \eqref{5} as the functional differential equation \begin{equation}\label{2.1} \dot{z}=F(z_t, \tau), \end{equation} where $z=(x,y)^T,z_t(\theta)=z(t+\theta) \in C([-\tau,0],\mathbb{R}_+^2)$. Before stating the main results of this section, we will give a known lemma. Let \begin{equation} \begin{gathered} \mathbb{X}=C([-\tau,0],\mathbb{R}_+^2), \\ \Sigma=Cl\{(z,\tau,p)\in \mathbb{X}\times\mathbb{R}_+\times\mathbb{R}_+: z \text{$z$ is a nonconstant periodic solution of \eqref{2.1}}\},\\ N=\{(\bar{z},\tau, p): F(\bar{z},\tau,p)=0\}. \end{gathered} \end{equation} \begin{lemma}[\cite{wu16}] \label{lem3.1} Only one of the following two results holds \begin{itemize} \item[(1)] $l_{(P_2,\tau_j,2\pi/\omega_l)}$ is unbounded; \item[(2)] $l_{(P_2,\tau_j,2\pi/\omega_l)}$ is bounded and $\sum_{(\bar{z},\tau,p)\in l_{(P_2,\tau_j,2\pi/\omega_l)} \cap N}\gamma(\bar{z},\tau,p)=0$. \end{itemize} \end{lemma} Assume $l_{(P_2, \tau_j, 2\pi/\omega_l)}$ is the connected component of $(P_2, \tau_j, 2\pi/\omega_l)$ in $\Sigma$, from lemma \ref{lem2.2}, we know that $l_{(P_2, \tau_j, 2\pi/\omega_l)}$ is nonempty. Next we give two lemmas needed for the proof of the main result. \begin{lemma} \label{lem3.2} If $0<\zeta/\beta<\min \{\delta, \alpha\}, \alpha\delta>\sigma$, then all periodic solutions of \eqref{5} are uniformly bounded. \end{lemma} \begin{proof} From the second equation of system \eqref{5}, we have $$ y(t)=y(0)\exp\int_0^t( \alpha (1-\beta y(s))-x(s))ds, $$ noting the definition of $y(t)$ and $P_2$ is a positive equilibrium, one knows $x(0)\geq0$ and $y(0)\geq0$; that is, $y(t)\geq0$. It is easy to see $$ \frac{\mathrm{d}y}{\mathrm{d}t}=\alpha y(1-\beta y)-xy\leq\alpha y(1-\beta y), $$ it follows $y\leq\frac{1}{\beta}$, that is $0\leq y\leq\frac{1}{\beta}$. We will prove $x(t)\geq0$ as follows. Since $(x(t), y(t))$ is a periodic solution of \eqref{5}, then $x(t)$ must have the minimum $x(\eta_1),\ \eta_1>0$, that is $$ \sigma +\zeta x(\eta_1-\tau)y(\eta_1-\tau)-\delta x(\eta_1)=0, $$ hence $$ \delta x(\eta_1)\geq \zeta x(\eta_1-\tau)y(\eta_1-\tau). $$ If $x(\eta_1-\tau)\geq0$, then $x(t)\geq0$. If $x(\eta_1-\tau)<0$, then $x(\eta_1)\leq x(\eta_1-\tau)<0$, hence $$ \frac{\mathrm{d}x}{\mathrm{d}t}=\sigma +\zeta x(t-\tau)y(t-\tau)-\delta x > \frac{\zeta}{\beta} x(\eta_1)-\delta x(t). $$ By the constant variation formula and the comparison theorem, it implies that \begin{equation}\label{proof} x(t)> \frac{\zeta x(\eta_1)}{\beta\delta}(1-\exp^{-\delta t})+x(0) \exp^{-\delta t} , \quad t>0. \end{equation} noting $x(0)>0$ and $0<\frac{\zeta}{\beta}<\min \{\delta, \alpha\}$, it is easy to see that \eqref{proof} does not hold as $t=\eta_1$; this is a contradiction, thus $x(t)>0$. On the other hand, we have $$ \frac{\mathrm{d}x}{\mathrm{d}t}=\sigma +\zeta x(t-\tau)y(t-\tau) -\delta x\leq\sigma +\frac{\zeta}{\beta} x(t-\tau)-\delta x, $$ Then we prove the bound of the solution $x(t)$ in three cases: (1) If $x(t-\tau)\geq x(t)$ for $t>0$, since $(x(t), y(t))$ is a periodic solution of \eqref{5}, then $x(t)$ must have the maximum $x(\eta),\ \eta\in(0,\tau)$, that is $\sigma +\zeta x(\eta-\tau)y(\eta-\tau)-\delta x(\eta)=0$, hence $$ x(\eta)\leq \frac{\sigma}{\delta}+\frac{\zeta}{\beta\sigma}x(\eta-\tau)\leq \frac{\sigma}{\delta}+\frac{\zeta}{\beta\sigma}\|\varphi\|,\quad \eta-\tau\in(-\tau,0). $$ (2) If $x(t-\tau)0$, then $$ \frac{\mathrm{d}x}{\mathrm{d}t}\leq\sigma +\frac{\zeta}{\beta} x -\delta x, $$ from $0<\zeta/\beta<\delta$, one has $$ x(t)\leq\frac{\sigma}{\delta-\frac{\zeta} {\beta}}; $$ (3) If there exists $t>0$ such that $x(t-\tau)\geq x(t)$ and $t_1>0$ such that $x(t_1-\tau)\sigma$, then system \eqref{5} has no-nonconstant periodic solution with period $\tau$. \end{lemma} \begin{proof} Suppose \eqref{5} has no-constant periodic solution with period $\tau$, then the corresponding ordinary differential equations \begin{equation}\label{2.2} \begin{gathered} \frac{\mathrm{d}x}{\mathrm{d}t}=\sigma +\zeta xy-\delta x, \\ \frac{\mathrm{d}y}{\mathrm{d}t}=\alpha y(1-\beta y)-xy, \end{gathered} \end{equation} have a non-constant periodic solution. System \eqref{2.2} has a positive equilibrium $P_2(x_2, y_2)$ and a boundary equilibrium $P_0(\frac{\sigma}{\delta}, 0)$. From lemma \ref{lem2.1}, we know that $P_2$ is stable and \eqref{2.2} has no non-constant periodic solution. This is a contradiction. The proof is complete. \end{proof} \begin{theorem} \label{thm3.4} If $0<\zeta/\beta<\min \{\delta, \alpha\}, \alpha\delta>\sigma$, then \eqref{5} has at least $j-1$ periodic solutions for $\tau\geq\tau_j,\ j\geq 1$. \end{theorem} \begin{proof} The characteristic matrix of \eqref{2.1} at the equilibrium $\bar{z}=(\bar{z}_1,\bar{z}_2)^T\in \mathbb{R}^2$ is $$ \Delta(\bar{z},\tau,p)(\lambda)=\lambda Id -DF(\bar{z},\tau,p)(e^{\lambda\cdot}Id); $$ i.e. \begin{equation}\label{2.3} \Delta(\bar{z},\tau,p)(\lambda)= \begin{pmatrix} &\lambda-\zeta\bar{z}_2e^{-\lambda\tau}-\delta & -\zeta\bar{z}_1e^{-\lambda\tau}\\ &\bar{z}_2 &\lambda-\alpha+2\alpha\beta\bar{z}_2+\bar{z}_1 \end{pmatrix}, \end{equation} where $Id$ is the identity matrix, $DF(\bar{z},\tau,p)$ is the Fr\'echet derivative of $F$ with respect to $z_t$ evaluated at $(\bar{z},\tau,p)$. On the other hand, from the definition of \cite{wu16}, we know that $(\bar{z},\tau,p)$ is called a center of \eqref{2.1} if $(\bar{z},\tau,p)\in N$ and $\det (\Delta(\bar{z},\tau,p)(\frac{2\pi}{p}i))=0$ and the center is said to be isolated if it is the only center in the neighborhood of $(\bar{z},\tau,p)$. From \eqref{2.3} we know that $$ \det (\Delta(P_0,\tau,p)(\lambda)) =(\lambda-\delta)(\lambda-\alpha+\frac{\sigma}{\delta})=0, $$ obviously, this equation has no purely imaginary roots, thus \eqref{2.1} has no centers with the form $(P_0,\tau,p)$. From the proof of \cite[theorem 4.2]{xiao, RYafia7.2}, we know that there exist $\epsilon_1>0, \epsilon_2>0$ and a smooth curve $\lambda(\tau):[\tau_j-\epsilon_1,\tau_j+\epsilon_1]\to C$ such that $\det (\Delta(\lambda(\tau)))=0$ as $|\lambda(\tau)-i\omega_l|<\epsilon_2$ for all $\tau \in [\tau_j-\epsilon_1,\tau_j+\epsilon_1]$, furthermore $$ \lambda(\tau_j)=i\omega_l,\quad \frac{\operatorname{Re} \lambda(\tau)}{\mathrm{d}\tau}\big|_{\tau=\tau_j}>0; $$ thus, for all $\tau>0$, $(P_2, \tau_j, \frac{2\pi}{\omega_l})(j \in \mathrm{N})$ is the only center of \eqref{2.1}. Let $$ \Omega_{\epsilon,2\pi/\omega_l}=\Big\{(\eta,p) :0<\eta<\epsilon,|p-2\frac{\pi}{\omega_l}|<\epsilon\Big\}. $$ If $|\tau-\tau_j|\leq\epsilon_1$ and $(\eta,p)\in \partial \Omega_{\epsilon,2\pi/\omega_l}$, then $$ \det (\Delta(P_2,\tau,p)(\eta+\frac{2\pi}{p}i))=0\;\Leftrightarrow\; \eta=0, \tau=\tau_j,p=\frac{2\pi}{\omega_l}. $$ Define $$ H^\pm(P_2,\tau_j,2\pi/\omega_l)(\eta,p)=\det (\Delta(P_2,\tau\pm\epsilon_1,p)(\eta+\frac{2\pi}{p}i)). $$ Then the crossing number of isolated center is \begin{equation} \begin{split} &\gamma(P_2,\tau_j,2\pi/\omega_l)\\ &=\deg _B(H^-(P_2,\tau_j,\frac{2\pi} {\omega_l}),\Omega_{\epsilon,2\pi/\omega_l}) -\deg _B(H^+(P_2,\tau_j,2\pi/\omega_l), \Omega_{\epsilon,2\pi/\omega_l}) =-1; \end{split} \end{equation} thus $$ \sum_{(\bar{z},\tau,p)\in l_{(P_2,\tau_j,2\pi/\omega_l)}\cap N} \gamma(\bar{z},\tau,p)<0. $$ Noting lemma \ref{lem3.2}, we know that the connected component $l_{(P_2,\tau_j,2\pi/\omega_l)}$ passing through $(P_2,\tau_j,2\pi/\omega_l)$ is unbounded in $\Sigma$. From lemma \ref{lem2.1} we know that the projection of $l_{(P_2,\tau_j,2\pi/\omega_l)}$ onto $z$-space is bounded. Noting \cite{RYafia7.2}, one has \begin{equation}\label{2.4} \frac{2\pi}{\omega_l}<\tau_j. \end{equation} From the results of lemma \ref{lem3.3}, we know that \eqref{2.1} has no non-trivial periodic solutions for $\tau=0$. Hence the projection of $l_{(P_2,\tau_j,2\pi/\omega_l)}$ onto $\tau$-space is away from zero. If not, the projection of $l_{(P_2,\tau_j,2\pi/\omega_l)}$ onto $\tau$-space is included in the interval $(0,\tau^\ast), \tau^\ast>\tau_j$, with the help of \eqref{2.4}, we know $(\bar{z},\tau,p)\in l_{(P_2,\tau_j,2\pi/\omega_l)}$ as $p<\tau^\ast$. This implies the projection of $l_{(P_2,\tau_j,2\pi/\omega_l)}$ onto $p$-space is bounded. Thus if we want the connected component $l_{(P_2,\tau_j,2\pi/\omega_l)}$ is unbounded, the projection of $l_{(P_2,\tau_j,2\pi/\omega_l)}$ onto $\tau$-space must be unbounded; i.e., the projection of $l_{(P_2,\tau_j,2\pi/\omega_l)}$ onto $\tau$-space must conclude $[\tau_j,\infty)$, this implies that the system \eqref{5} has at least $j-1$ periodic solutions for all $\tau\geq\tau_j(j\geq1)$ . \end{proof} \subsection*{Discussion} We studied the nonlinear dynamics of a Kuznetsov Makalkin and Taylor's model with delay, which is analyzed in \cite{RYafia7.2}. In \cite{RYafia7.2}, Yafia only give the existence of the Hopf bifurcation. In this paper, we give the general formula of the direction, the estimation formula of period and stability of bifurcated periodic solution. Especially, using a global Hopf bifurcation due to Wu\cite{wu16}, the global existence of periodic solutions bifurcating from Hopf bifurcations is given. we show that the local Hopf bifurcation implies the global Hopf bifurcation after the second critical value of the delay. Numerical simulations are carried out to illustrate the the theoretical analysis and results. Our results on the existence and stability of the Hopf bifurcated periodic solutions of $P_2$ describe the equilibrium process. When a global stable periodic orbit exists, it can be understood that the tumor and the immune system can coexist for all the time although the cancer is not eliminated. The conditions for the parameters provide theories basis to control the development or progression of the tumors. The phenomena have been observed in some models d'Onofrio \cite{d'Onofrio2}, Kuznetsov et al \cite{Kuznetsov}, Bi et al \cite{xiao} . In particular, Bi et al. \cite{Bi} have shown that various bifurcations, including Hopf bifurcation, Bautin bifurcation and Hopf-Hopf bifurcation, can occur in such models. Finally, we point out that we have studied the dynamical behaviors of $P_2$. In fact, the dynamical behaviors of $P_0$ is more rich such as the existence of the Bogdanov-Takens bifurcation and steady-state bifurcation, which is studied in \cite{xiao}. The two equilibria may coexist, correspondingly, the system can exhibit more degenerate bifurcations including Hopf-Hopf and resonant higher codimension bifurcations et al. It would be interesting to consider these dynamics of the delayed model. \subsection*{Acknowledgments} This research was supported by National Natural Science Foundation of China (No. 11171110) and Shanghai Leading Academic Discipline Project (No. B407), 211 Project of Key Academic Discipline of East China Normal University. \begin{thebibliography}{00} \bibitem{Bi} P. Bi, S. Ruan; Bifurcations in Delay Differential Equations and Applications to Tumor and Immune System Interaction Models, \textit{ SIAM J. Appl. Dynamical Systems,} Vol. 12(4),(2013) 1847-1888. \bibitem{xiao} P. Bi, H. Xiao, Bifurcations of a Tumor-Immune Competition Systems with Delay, submitted. \bibitem{d'Onofrio2} A. d'Onofrio; A general framework for modeling tumour-immune system competition and immunotherapy: Mathematical analysis and biomedical inferences, \textit{Phys. D} \textbf{208} (2005), 220-235. \bibitem{d'Onofrio3} A. d'Onofrio; Tumour-immune system interaction: Modeling the tumour-stimulated proliferation of effectors and immunotherapy, \textit{Math. Models Methods Appl. Sci.} \textbf{16} (2006), 1375-1401. \bibitem{d'Onofrio4} A. d'Onofrio; Metamodeling tumour-immune system interaction, tumour evasion and immunotherapy, \textit{Math. Comput. Modelling} \textbf{47} (2008), 614-637. \bibitem{T.Faria8.1} T. Faria, L. Magalh\~aes; Normal forms for retarded functional differential equation and applications to Bogdanov-Takens singularity, \textit{J. Diff. Equ.} \textbf{122}(1995), 201--224. \bibitem{T.Faria8.2} T. Faria, L. Magalh\~aes, Normal forms for retarded functional differential equation with parameters and applications to Hopf bifurcation, \textit{J. Diff. Equ.} \textbf{122}(1995), 181--200. \bibitem{MGalach6} M. Galach; Dynamics of the tumor-immune system competition: the effect of time delay, \textit{Int. J. Appl. Math. Comput. Sci.} \textbf{13} 2003, 395--406. \bibitem{J.K.Hale1} J. K. Hale; Theory of functional differential equations, Springer-Verlag[M], New York, 1977. \bibitem{Hassard}B. Hassard, N. Kazarinoff, Y. Wan; Theory and Application of Hopf Bifurcation. Cambridge: Cambridge University Press, 1981. \bibitem{Kuznetsov} V. Kuznetsov, I. Makalkin, M. Taylor, A. Perelson; Nonlinear dynamics of immunogenic tumors: parameter estimationestimation and global bifurcation analysis, \textit{Bull. Math. Biol.} \textbf{56} (1994), 295-321. \bibitem{liu1} D. Liu, S. Ruan, D. Zhu; Stable periodic oscillations in a two-stage cancer model of tumor-immune interaction, \textit{Mathematical Biosciences and Engineering} \textbf{12} (2009)151--168. \bibitem{Mayer} H. Mayer, K. Zaenker, U. an der Heiden; A basic mathematical model of the immune response, \textit{Chaos},\textbf{5}(1995) 155--161. \bibitem{wu16} J. Wu.; Symmetric functional differential equations and neural networks with memory. \textit{ Trans. Amer. Math. Soc.,} \textbf{350}(1998) 4799--4838. \bibitem{RYafia7.1} R. Yafia; Hopf bifurcation analysis and numerical simulations in an ODE model of the immune system with positive immune response, \textit{Nonl. Anal.: Real World Appl.} \textbf{8}(2007) 1359--1369. \bibitem{RYafia7.2} R. Yafia; Hopf bifurcation in differential equations with delay for tumor-immune system competition model, \textit{SIAM J. Appl. Math,} \textbf{67}(2007) 1693--1703. \bibitem{RYafia7.3} R. Yafia; Hopf bifurcation in a delayed model for tumor-immune system competition with negative immune response, \textit{ Discrete Dynamics in Nature and Soc.} \textbf{11} (2006) 1--9. \bibitem{RYafia7.4} R. Yafia; Periodic solutions for small and large delays in a tumor-immune system model, \textit{Electron. J. of Diff. Equ.}, Conference \textbf{14} (2006) 241--248. \end{thebibliography} \end{document}