Bifurcation Analysis of a Neutrophil Periodic Oscillation Model with State Feedback Control

Suqi Ma

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

DOI: 10.4236/ijmnta.2023.121001                

Abstract

The mathematical model of stem cells is discussed with its motivation to describe the tissue relationship by technically introducing a two compartments model. The clear link between the proliferation phase of stem cells and the circulating neutrophil phase is set forth after delay feedback control of the state variable of stem cells. Hopf bifurcation is discussed with varying free parameters and time delays. Based on the center manifold theory, the normal form near the critical point is computed and the stability of bifurcating periodical solution is rigorously discussed. With the aids of the artificial tool on-hand which implies how much tedious work doing by DDE-Biftool software, the bifurcating periodic solution after Hopf point is continued by varying time delay.

Keywords

Neutrophil Phase, Time Delay, Hopf Bifurcation, DDE-Biftool, Fold Periodical Bifurcation, Period-Doubling Bifurcation

1. Introduction

With concentration on the mathematical model of stem cells in understanding its tissue organization, the discussion of the cell’s differentiation and proliferation discipline has insight people to get the interesting view sights in dynamics of stem cells. The investigation work of Mackey and his collaborators has shown some examples of hematopoietic system dynamics within the limitation of the topic of DDEs (delay differential equations), which is infinite-dimensional and its present dynamical behavior also depends on the past history [1] [2] [3] [4] . Refer to [2] [5] [6] [7] , some mathematical models with dynamical description which is originates from stem cells compartment to periodical oscillation neutrophil disease, and combined with the production and regulation of blood cells, are set forth by authors to be governed by DDEs versions.

The clear link between the proliferation phase of stem cells and the circulating neutrophil compartment is described technically as shown in Figure 1. With the introduction rate $\beta$ from resting phase, stem cells can enter into their proliferating phase, then after time delay ${\tau }_s$ , either go to programmed death due to its apoptosis or reentrily go into the resting phase since experiencing the maturation phase. After a cell division, the neutrophil number which gets through the proliferation phase and maturation phase is released into circulation through the body. Examples are as often seen in reference of Mackey’s work, which investigates a tissue system composed of the stem cells compartment and neutrophil compartment, wherein the governing discipline and the relationship between two compartments are described by DDEs.

Alike Lei but a few other people’s points of view [8] [9] [10] , the mechanism applying for periodic oscillating neutrophil numbers is very complex, however, the granulocyte-colony stimulation factor in its amplification coefficient is introduced. For the treatment of cyclical neutropenia experiencing time ${\tau }_{NM}$ , the amplification coefficient is also time delay-dependent. Hence, the time delay in neutrophil coupling is composed of two time segments, that is, the time duration respectively of proliferation phase and maturation phase of a neutrophil precursor. Lei emphasises that the level of stem cells decreases periodically in the treatment

Figure 1. The cell lineage of the stem cells with its clear link to the neutrophils.

of neutrophil oscillation phenomena, and the corresponding amplification coefficient can be written as

$A={\text{e}}^{{\eta }_{NP}{\tau }_{NP}-{\gamma }_0{\tau }_{NM}}$ (1.1)

herein, the proliferation rate is denoted as ${\eta }_{NP}$ , and the death rate during maturation phase is ${\gamma }_0$ , and ${\tau }_N={\tau }_{NP}+{\tau }_{NM}$ . Therefore, DDEs model to govern the discipline of neutrophil cell lines is set forth [6] [11]

$\begin{array}{l}\frac{\text{d}Q}{\text{d}t}=-\left(\beta \left(Q\right)+k_N\left(N\right)+k_{\delta }\right)Q+2{\text{e}}^{-r_s{\tau }_s}\beta \left(Q_{{\tau }_s}\right)Q_{{\tau }_s}\\ \frac{\text{d}N}{\text{d}t}=-{\gamma }_{N}N+A{k}_N\left(N_{{\tau }_N}\right)Q_{{\tau }_N}\end{array}$ (1.2)

The hematopoietic stem cells can enter into the proliferation phase with the rate $\beta$ (days⁻¹), and differentiate into the committed neutrophil compartment at the rate $k_N$ (days⁻¹), which is listed as

$\begin{array}{l}{k}_N\left(N\right)=\frac{f_0{\theta }_1^m}{{\theta }_1^m+N^m},\beta \left(Q\right)=\frac{k_0{\theta }_2^n}{{\theta }_2^n+Q^n}\\ k_N\left(N\left(t-{\tau }_N\right)\right)=\frac{f_0{\theta }_1^m}{{\theta }_1^m+N{\left(t-{\tau }_N\right)}^m},\beta \left(Q_{{\tau }_s}\right)=\frac{k_0{\theta }_2^n}{{\theta }_2^n+Q{\left(t-{\tau }_s\right)}^n}\end{array}$

The estimation of parameter is valued as the following ${\gamma }_N=2.4$ , ${\gamma }_0=0.01$ , $k_0=2.5$ , $f_0=0.4$ , ${\theta }_2=0.3$ , ${\theta }_1=36$ , $k_{\delta }=0.01$ , $r_s=0.2$ , ${\eta }_{NP}=2.5379$ , ${\tau }_N=11$ , ${\tau }_s=2.8238$ . Some known results as the fold bifurcation phenomena of limit cycle, since the collision of stable and unstable limit cycles happened with Takens-Bagnov sigularity of the equilibrium solution, have discussed in paper [11] .

$\begin{array}{l}\frac{\text{d}Q}{\text{d}t}=-\left(\beta \left(Q\right)+k_N\left(N\right)+k_{\delta }\right)Q+2{\text{e}}^{-r_s{\tau }_s}\beta \left(Q_{{\tau }_s}\right)Q_{{\tau }_s}\\ +K_1\left(x\left(t-{\tau }_s\right)-x\left(t\right)\right)+K_2\left(x\left(t-{\tau }_N\right)-x\left(t\right)\right)\\ \frac{\text{d}N}{\text{d}t}=-{\gamma }_{N}N+A{k}_N\left(N_{{\tau }_N}\right)Q_{{\tau }_N}\end{array}$ (1.3)

For further discussing the periodical oscillating mechanism, we set forth the state feedback control strategy which changes the system stability according to Hopf bifurcation criterion, hence, the small amplitude periodical solutions appear. Herein, we introduce state feedback control via different ways with delay signals, therefore, $K_1\mathrm{,}{K}_2$ are weak feedback coefficients.

Whilst the real parts of the corresponding characteristic roots become zero, DDEs system may lose the stability property of the positive equilibrium solution. In the general case, Hopf bifurcation occurs which is born with the new periodical solution under the assumption that Hopf bifurcation satisfies the transversal condition [12] [13] . We mainly focus on Hopf bifurcation of hematopoietic system with respect to time delay ${\tau }_s$ . As free parameter ${\tau }_s$ varies continuously, the bifurcating periodical solution changes its stability character which happened with one real Flouquet multiplier or a pair of conjugate complex multipliers getting through the unit circle [14] [15] [16] [17] [18] . If and only if the Flouquet multiplier attains at −1, the period-doubling bifurcation of periodical solution arises at the threshold value. The collision phenomena of the limit cycles also happen as Flouquet multiplier is positive 1. The hysteresis phenomena of the limit cycle are observed as the subcritical Hopf bifurcation occurs at some threshold value. The interesting dynamical bifurcation is also discussed with feedback control coefficient $K_2$ . Simultaneously, the routes of period-doubling bifurcation to chaos are discovered.

As shown in Figure 2, the coexistence phenomena of periodical solutions are observed in system (1.3) undertaking the strategy of the state feedback control. The coexistence limit cycle can be produced by the fold bifurcation or period-doubling bifurcation of the limit cycle. Generally, the hysteresis phenomena of limit cycle are also observed since fold bifurcation with limit cycle collision phenomena. In either case of bifurcation scenario, the limit cycle with a doubly period arises near the period-doubling bifurcation point. DDE-Biftool software program solves the window visualization of limit cycle continuation as varying free parameters. As shown in Figure 2(a), three different period-2 solutions are observed with chosen feedback coefficient $K_1=-0.1,K_2=-0.005$ . We also observe period-2 and period-4 solutions coexistence phenomena with chosen feedback coefficient $K_1=-0.1,K_2=-0.015$ , as shown in Figure 2(b). The period-doubling bifurcation occurs in system (1.3) while varying feedback coefficient which rich the system bifurcation behaviors.

The whole paper is organized as the listed. In Section 2, the routes to chaos are discussed via period-doubling bifurcation and the Poincare section simulation verifies the coincident results as varying $K_2$ feedback coefficient continuously. In Section 3, Hopf bifurcation is discussed and the transversal condition is derived. In Section 4, the stability of bifurcating periodical solution from Hopf bifurcation is further discussed by dimension reduction technique combined with the center manifold theory [19] [20] [21] . And the simulation results verify the

Figure 2. The coexitence phenomena of stable and unstable periodical solutions finding in system (1.3) with different feedback control coeffients (a) $K_1=-0.1,K_2=-0.005$ ; (b) $K_1=-0.1,K_2=-0.015$ .

analysis results. The limit cycle bifurcation is discussed further with the artificial DDE-Biftool computation technique. Finally, a conclusion is given.

2. Routes to Chaos

As Hopf bifurcation happens, the new bifurcating periodical solution occurs which may bring forth the dynamical property of system (1.3) qualitatively changed. As Flouquet multiplier of the periodical solutions has attained at −1, the period-doubling bifurcation of periodical solution arises which brings forth a new periodical solution having its doubly period. The most interesting phenomenon is the scenario cascades that periodical bifurcation produces new solution with enlarged period continuously until chaos appears, which is called the routes of period-doubling bifurcation to chaos [14] . With chosen $K_1=-0.1$ , whilst varying $K_2$ continuously, we get period-2, period-4, period-8 solution, etc. of system (1.3), as shown in Figure 3. The period window is observed by the numerical analysis with the Poincare section $N\left(t-{\tau }_N\right)=5000$ , as shown in Figure 4.

3. Hopf Bifurcation

As a pair of imaginary roots cross the imaginary axis from the left half plane to the right half plane, Hopf bifurcation occurs with the non-degenerate transversal condition. To address the result derive in paper [13] , we write out the linearized version of system (1.2) at the positive equilibrium solution $E^{\mathrm{<mo>\; *\; </mo>}}=\left(Q^{\mathrm{<mo>\; *\; </mo>}}\mathrm{,}{N}^{\mathrm{<mo>\; *\; </mo>}}\right)$ as the following listed

$\begin{array}{l}{x}^{\prime }=a_{11}x+a_{12}y+b_{1}x\left(t-{\tau }_s\right)+K_1\left(x\left(t-{\tau }_s\right)-x\left(t\right)\right)+K_2\left(x\left(t-{\tau }_N\right)-x\left(t\right)\right),\\ y^{\prime }=-{\gamma }_{N}y+b_{2}x\left(t-{\tau }_N\right)+b_{3}y\left(t-{\tau }_N\right)\end{array}$ (2.1)

with

$\begin{array}{l}{a}_{11}=-{\beta }^{\prime }\left(Q^{*}\right)Q^{*}-\beta \left(Q^{*}\right)-k_N\left(N^{*}\right)-k_{\delta }-K_1-K_2\\ a_{12}=-k_{N^{\prime }}\left(N^{*}\right)Q^{*},\\ b_1=2{\text{e}}^{-r_s{\tau }_s}\left(\beta \left(Q^{*}\right)+{\beta }^{\prime }\left(Q^{*}\right)Q^{*}\right)+K_1,\\ b_2=Ak\left(N^{*}\right),b_3=A{k}^{\prime }\left(N^{*}\right)Q^{*}\end{array}$

The characteristic equation of the linear equation (2.1) is written as

$\begin{array}{c}\Delta \left(\lambda ,r_s,{\tau }_s,{\tau }_N\right)=K_2{b}_3{\text{e}}^{-2\lambda {\tau }_N}+{\text{e}}^{-\lambda \left({\tau }_s+{\tau }_N\right)}{b}_1{b}_3-{\text{e}}^{-\lambda {\tau }_s}\left({\gamma }_N{b}_1+b_1\lambda \right)\\ +{\text{e}}^{-\lambda {\tau }_N}\left(a_{11}{b}_3-a_{12}{b}_2-b_3\lambda -K_2{\gamma }_N-K_2\lambda \right)\\ -a_{11}{\gamma }_N-a_{11}\lambda +{\gamma }_N\lambda +{\lambda }^2\\ =0\end{array}$ (2.2)

We aims to compute Hopf bifurcation point by considering parameter on $\left(r_s\mathrm{,}{\tau }_s\right)$ -plane, and assume the characteristic equation (2.2) has the imaginary root $\lambda =i\omega \left(\omega >0\right)$ . Define three new angle variables $\theta =\omega {\tau }_0,\varphi =\omega {\tau }_s,\psi =\omega {\tau }_N$ since the introduction of new time delay ${\tau }_0$ , which is an auxiliary delay though with the assumption ${\tau }_0=1$ as discussed subsequently.

Figure 3. (a) P-1 solution with $K_2=-0.005$ ; (b) P-2 solution with $K_2=-0.005$ ; (c) P-4 solution with $K_2=-0.015$ ; (d) P-8 solution with $K_2=-0.026$ ; (e) strange attractors with $K_2=-0.028$ ; (f) chaos with $K_2=-0.03$ .

Figure 4. The continuation of equilibrium solutions and bifurcating periodical solutions as varying free parameter ${\tau }_s$ continuously. Other parameters are fixed with (a) $K_2=-0.005$ ; (b) $K_2=-0.015$ .

Hence after, by Euler formula, one gets

$\begin{array}{l}{\text{e}}^{-\lambda \left({\tau }_N+{\tau }_s\right)}=\cos \left(\varphi +\psi \right)-i\sin \left(\varphi +\psi \right),\\ {\text{e}}^{-\lambda {\tau }_s}=\cos \left(\varphi \right)-i\sin \left(\varphi \right),\\ {\text{e}}^{-\lambda {\tau }_N}=\cos \left(\psi \right)-i\sin \left(\psi \right),\\ {\text{e}}^{-2\lambda {\tau }_N}=\cos \left(2\psi \right)-i\sin \left(2\psi \right)\end{array}$ (2.3)

Substitute them into the characteristic Equation (2.2) and separate the real part from the imaginary part, we get

$\begin{array}{l}{b}_1{b}_3\cos \left(\varphi +\psi \right)-{\gamma }_N{b}_1\cos \left(\varphi \right)-b_1\omega \sin \left(\varphi \right)\\ =-\left(a_{11}{b}_3-a_{12}{b}_2-K_2{\gamma }_N\right)\cos \left(\psi \right)-K_2{b}_3\cos \left(2\psi \right)\\ +\left(K_2\omega +b_3\omega \right)\sin \left(\psi \right)+a_{11}{\gamma }_N+{\omega }^2\\ b_1{b}_3\sin \left(\varphi +\psi \right)+b_1\omega \cos \left(\varphi \right)-{\gamma }_N{b}_1\sin \left(\varphi \right)\\ =-\left(a_{11}{b}_3-a_{12}{b}_2-K_2{\gamma }_N\right)\sin \left(\psi \right)-K_2{b}_3\sin \left(2\psi \right)\\ -a_{11}\omega +{\gamma }_N\omega -\left(K_2\omega +b_3\omega \right)\cos \left(\psi \right)\end{array}$ (2.4)

Solving $\sin \left(\varphi \right)\mathrm{,}\cos \left(\varphi \right)$ from Equation (2.4) to get

$\begin{array}{l}\cos \left(\varphi \right)=-\frac{\left(2K_2{b}_3\omega \cos \left(\psi \right)-\omega \left(-2a_{11}{b}_3+a_{12}{b}_2\right)\right)\sin \left(\psi \right)-2K_2\cos {\left(\psi \right)}^2{\gamma }_N{b}_3}{b_1\left(b_3^2+2b_3\omega \sin \left(\psi \right)-2\cos \left(\psi \right){\gamma }_N{b}_3+{\gamma }_N^2+{\omega }^2\right)}\\ -\frac{\left(K_2{\gamma }_N^2+K_2{b}_3^2+K_2{\omega }^2-2a_{11}{\gamma }_N{b}_3+a_{12}{\gamma }_N{b}_2\right)\cos \left(\psi \right)+a_{11}{b}_3^2-a_{12}{b}_2{b}_3+a_{11}\left({\gamma }_N^2+{\omega }^2\right)}{b_1\left(b_3^2+2b_3\omega \sin \left(\psi \right)-2\cos \left(\psi \right){\gamma }_N{b}_3+{\gamma }_N^2+{\omega }^2\right)}\\ \sin \left(\varphi \right)=-\frac{\left(-2K_2\cos \left(\psi \right){\gamma }_N{b}_3+\left(K_2+2b_3\right){\omega }^2+a_{12}{\gamma }_N{b}_2+K_2{b}_3^2+K_2{\gamma }_N^2\right)\sin \left(\psi \right)}{b_1\left(b_3^2+2b_3\omega \sin \left(\psi \right)-2\cos \left(\psi \right){\gamma }_N{b}_3+{\gamma }_N^2+{\omega }^2\right)}\\ -\frac{-2K_2{b}_3\cos {\left(\psi \right)}^2\omega +\left(a_{12}{b}_2-2{\gamma }_N{b}_3\right)\omega \cos \left(\psi \right)+\omega \left(2K_2{b}_3+{\gamma }_N^2+b_3^2+{\omega }^2\right)}{b_1\left(b_3^2+2b_3\omega \sin \left(\psi \right)-2\cos \left(\psi \right){\gamma }_N{b}_3+{\gamma }_N^2+{\omega }^2\right)}\end{array}$ (2.5)

By the triangle relationship, people set up the following equalities

$\begin{array}{l}F\left(r_s,{\tau }_s,\theta ,\varphi ,\psi \right)=\cos {\left(\varphi \right)}^2+\sin {\left(\varphi \right)}^2-1\equiv 0,\\ G\left(r_s,{\tau }_s,\theta ,\varphi ,\psi \right)=\varphi -\arctan \left(\frac{\sin \left(\varphi \right)}{\cos \left(\varphi \right)}\right)-2n\pi \equiv 0\end{array}$ (2.6)

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

Set the initial value of ${\tau }_0$ , for example ${\tau }_0=1$ , hence after, angle variable $\theta$ is chosen as an independent varying parameter. That is, as $\theta$ varying over its feasible value sets, one draw Hopf bifurcation curve by tracking the pathway with the given parameter value $r_s=r_s\left(\theta \right),{\tau }_s={\tau }_s\left(\theta \right)$ .

Differentiate both the side of Equation (2.5) with respect to $\theta$ , one gets

$\begin{array}{l}-\sin \left(\varphi \right)\left(\frac{1}{{\tau }_0}{\tau }_s+\frac{\theta }{{\tau }_0}{{\tau }^{\prime }}_s\left(\theta \right)\right)=\frac{\partial \cos \left(\varphi \right)}{\partial \theta }+\frac{\partial \cos \left(\varphi \right)}{\partial r_s}{r^{\prime }}_s\left(\theta \right),\\ \cos \left(\varphi \right)\left(\frac{1}{{\tau }_0}{\tau }_s+\frac{\theta }{{\tau }_0}{{\tau }^{\prime }}_s\left(\theta \right)\right)=\frac{\partial \sin \left(\varphi \right)}{\partial \theta }+\frac{\partial \sin \left(\varphi \right)}{\partial r_s}{r^{\prime }}_s\left(\theta \right)\end{array}$ (2.7)

Define the matrixes by

$\begin{array}{l}{D}_0=\left(\begin{array}{cc}\frac{\partial \cos \left(\varphi \right)}{\partial r_s}& \frac{\theta }{{\tau }_0}\sin \left(\varphi \right)\\ \frac{\partial \sin \left(\varphi \right)}{\partial r_s}& -\frac{\theta }{{\tau }_0}\cos \left(\varphi \right)\end{array}\right),D_1=\left(\begin{array}{cc}-\sin \left(\varphi \right)\frac{{\tau }_s}{{\tau }_0}-\frac{\partial \cos \left(\varphi \right)}{\partial \theta }& \frac{\theta }{{\tau }_0}\sin \left(\varphi \right)\\ \cos \left(\varphi \right)\frac{{\tau }_s}{{\tau }_0}-\frac{\partial \sin \left(\varphi \right)}{\partial \theta }& -\frac{\theta }{{\tau }_0}\cos \left(\varphi \right)\end{array}\right),\\ D_2=\left(\begin{array}{cc}\frac{\partial \cos \left(\varphi \right)}{\partial r_s}& -\sin \left(\varphi \right)\frac{{\tau }_s}{{\tau }_0}-\frac{\partial \cos \left(\varphi \right)}{\partial \theta }\\ \frac{\partial \sin \left(\varphi \right)}{\partial r_s}& \cos \left(\varphi \right)\frac{{\tau }_s}{{\tau }_0}-\frac{\partial \sin \left(\varphi \right)}{\partial \theta }\end{array}\right),\end{array}$ (2.8)

then we have

${r^{\prime }}_s\left(\theta \right)=\frac{D_1}{D_0},{{\tau }^{\prime }}_s\left(\theta \right)=\frac{D_2}{D_0}$ (2.9)

We compute the transversal condition of Hopf bifurcation subsequently. Without loss of generality, differentiate Equation (2.2) with respect to ${\tau }_s$ to get

$\begin{array}{l}\frac{\partial \Delta (\cdot )}{\partial \lambda }\frac{\text{d}\lambda }{\text{d}{\tau }_s}+\frac{\partial \Delta (\cdot )}{\partial {\tau }_s}=0\hfill \end{array}$ (2.10)

Hence after, we have

$\frac{1}{\frac{\text{d}\lambda }{\text{d}{\tau }_s}}=-\frac{\frac{\partial \Delta (\cdot )}{\partial \lambda }}{\frac{\partial \Delta (\cdot )}{\partial {\tau }_s}}$ (2.11)

We also differentiate Equation (2.2) with respect to independent parameter $\theta$ to get

$\frac{\partial \Delta (\cdot )}{\partial \lambda }\frac{\text{d}\lambda }{\text{d}\theta }+\frac{\partial \Delta (\cdot )}{\partial {\tau }_s}\frac{\text{d}{\tau }_s}{\text{d}\theta }+\frac{\partial \Delta (\cdot )}{\partial r_s}\frac{\text{d}{r}_s}{\text{d}\theta }=0$ (2.12)

By Equation (2.12), one has

$\frac{1}{\frac{\text{d}\lambda }{\text{d}{\tau }_s}}=\frac{\frac{\partial \Delta (\cdot )}{\partial {\tau }_s}\frac{\text{d}{\tau }_s}{\text{d}\theta }+\frac{\partial \Delta (\cdot )}{\partial r_s}\frac{\text{d}{r}_s}{\text{d}\theta }}{\frac{\partial \Delta (\cdot )}{\partial {\tau }_s}\frac{\text{d}\lambda }{\text{d}\theta }}=-i{\tau }_0\left(\frac{\text{d}{\tau }_s}{\text{d}\theta }+\frac{\frac{\partial \Delta (\cdot )}{\partial r_s}\frac{\text{d}{r}_s}{\text{d}\theta }}{\frac{\partial \Delta (\cdot )}{\partial {\tau }_s}}\right)$ (2.13)

Therefore, the transversal condition is determined by the sign of the real part of the left side of Equation (2.13), which is written as

$\begin{array}{c}sign\Re \left(\frac{\text{d}\lambda }{\text{d}{\tau }_s}\right)=sign\Im \left(\frac{\frac{\partial \Delta (\cdot )}{\partial r_s}}{\frac{\partial \Delta (\cdot )}{\partial {\tau }_s}}\right)sign{r^{\prime }}_s\left(\theta \right)\\ =sign\Im \left(\frac{\partial \Delta (\cdot )}{\partial r_s}\overline{\frac{\partial \Delta (\cdot )}{\partial {\tau }_s}}\right)sign{r^{\prime }}_s\left(\theta \right)\end{array}$ (2.14)

Therefore, we have the following result.

Theorem 2.1. There is one pair of imaginary roots cross over the imaginary axis form left half plane to the right half plane if the condition $\Im \left(\frac{\partial \Delta (\cdot )}{\partial r_s}\overline{\frac{\partial \Delta (\cdot )}{\partial {\tau }_s}}\right){r^{\prime }}_s\left(\theta \right)>0$ is satisfied, or from the right half plane to the left half plane if the condition $\Im \left(\frac{\partial \Delta (\cdot )}{\partial r_s}\overline{\frac{\partial \Delta (\cdot )}{\partial {\tau }_s}}\right){r^{\prime }}_s\left(\theta \right)<0$ is satisfied.

Note: DDE-biftool however provide us a facilitate hand-tool to solve Hopf bifurcation curve with computation and drawing work. As usual as often, when we use DDE-Biftool software, we can look the advantage of our algorithm that the transversal condition is given explicitly, and Hopf bifurcation curve is shown in Figure 6(a).

4. Continuation of Bifurcating Solution as Varying ${\tau }_s$

We continue the equilibrium solution as varying delay ${\tau }_s$ . As shown in Figure 5(a) and Figure 5(b), Hopf bifurcation occurs at the star point since the stability switching push the character roots cross the imaginary axis then forward to the right half plane or backward to the left half plane. The bifurcating periodical solutions are continued with flouquet multiplier computed continuously. Thankful for the solving method of DDE-Biftool manuscript of Flouquet multiplier since it is right to find the key value on the unit circle. As two examples, choosing $K_2=-0.005$ and $K_2=-0.015$ , the bifurcation diagram is plotted in Figure 5(a) and Figure 5(b), with green color denotes the stable states of equilibrium solution or periodical solution, while red color and pink color explains the unstable states. We met the difficulties during the process of computing period-doubling bifurcation solution sometimes. As shown in Figure 5, the period-2 solution bifurcates as bifurcating limit cycle (denoted as P-1) arising from Hopf points switches its stability. As shown in Figure 5(a), three branches of P-1 solution are plotted and one branch of P-2 solution bifurcating via period-doubling bifurcation. We also get P-1, P-2 solutions in Figure 5(b), then P-4 solution bifurcates further again with flouquet multiplier gets close to −1 as varying control coefficient $K_2$ . How about two branches of P-4 solution, it is computed by DDE-Biftool software. The fold limit cycle bifurcation occurs as the stable and unstable limit cycle collision. We record it by fold limit cycle bifurcation line with $K_2=0.005$ . Try our best, with locating the real flouquet multiplier of maximal absolute values equating positive 1 in addition a common multiplier usually equals 1, the limit cycle bifurcation line is drawn on parameter plane $r_s-{\tau }_s$ plane. As shown in Figure 6(a), Hopf bifurcation curve is plotted with the transversal condition is given in Section 2. Furthermore, the period-doubling bifurcation line is plotted in Figure 6(b) and the fold limt cycle bifurcation line is drawn in Figure 6(c)

Figure 5. The continuation of equilibrium solutions and bifurcating periodical solutions as varying free parameter ${\tau }_s$ continuously. (a) $K_2=-0.005$ ; (b) $K_2=-0.015$ .

Figure 6. The bifurcation curve of equilibrium solution and periodical solution on $r_s-{\tau }_s$ plane with chosen control coefficent $K_2=0.005$ . (a) Hopf bifurcation curve of the equilibrium solution; (b) the period-doubling bifurcation curve of periodical solution; (c) the fold limt cycle bifurcation curve of periodical solution.

Figure 7. The time portraits of P-2 solution. (a) With chosen $K_2=0.005$ ; (b) The continuation of P-2 solution with varying $K_2$ continuously.

rigorously. Bautin bifurcation occurs as the intersection of fold limit cycle bifurcation with Hopf line. Without feedback control, that is $K_{1,2}=0$ , Bautin bifurcation point is analyzed with codimension 2 singularity in [12] . We also demonstrate the stability of bifurcating periodical solution follows subsequently in Subsection 4.1. As shown in Figure 7(a) and Figure 7(b), the time portraits of a P-2 solution is plotted and also the continuation of portraits are plotted in Figure 7(b) with varying $K_2$ continuously.

Normal Form Computation

In this section, we apply Schmidt-Lyapunov reduction scheme combined with center manifold technique to compute the normal form near Hopf point [13] [18] [21] . Assume that Hopf bifurcation occurs at the critical point $\left(r_s^{\mathrm{<mo>\; *\; </mo>}}\mathrm{,}{\tau }_s^{\mathrm{<mo>\; *\; </mo>}}\right)$ , with the equilibrium solution $\left(Q^{\mathrm{<mo>\; *\; </mo>}}\mathrm{,}{N}^{\mathrm{<mo>\; *\; </mo>}}\right)$ . Originated from system (1.3), by doing axis transformation $x=Q-Q^{*},y=N-N^{*}$ , one gets the 3rd order trunction of its Taylor expansion as

$\begin{array}{l}{x}^{\prime }=a_{1}x+a_{2}x\left(t-{\tau }_s\right)+K_{2}x\left(t-{\tau }_N\right)+c_{1}y+g_0{x}^2+g_1{x}^3+g_{2}x{\left(t-{\tau }_s\right)}^2\\ +g_{3}x{\left(t-{\tau }_s\right)}^3+g_{4}xy+g_5{y}^2+g_{6}x{y}^2+g_7{y}^3,\\ y^{\prime }=-{\gamma }_{N}y+b_{1}x\left(t-{\tau }_N\right)+b_{2}y\left(t-{\tau }_N\right)+h_{0}y{\left(t-{\tau }_N\right)}^2+h_{1}y{\left(t-{\tau }_N\right)}^3\\ +h_{2}y\left(t-{\tau }_N\right)x\left(t-{\tau }_N\right)+h_{3}y{\left(t-{\tau }_N\right)}^{2}x\left(t-{\tau }_N\right),\end{array}$ (4.1)

with

$\begin{array}{l}{a}_1=-k_{\delta }-\frac{k_0{\theta }_2^2}{{\theta }_2^2+Q^{*2}}-\frac{f_0{\theta }_1}{{\theta }_1+N^{*}}+\frac{2k_0{\theta }_2^2{Q}^{*2}}{{\left({\theta }_2^2+Q^{*2}\right)}^2}-K_1-K_2,\\ a_2=2{\text{e}}^{-r_s{\tau }_s}\frac{k_0{\theta }_2^2}{{\theta }_2^2+Q^{*2}}-4{\text{e}}^{-r_s{\tau }_s}\frac{k_0{\theta }_2^2{Q}^{*2}}{{\left({\theta }_2^2+Q^{*2}\right)}^2}+K_1,b_1=A\frac{f_0{\theta }_1}{{\theta }_1+N^{*}},\end{array}$

$\begin{array}{l}{b}_2=-A\frac{f_0{\theta }_1{Q}^{*}}{{\left({\theta }_1+N^{*}\right)}^2},c_1=\frac{f_0{\theta }_1{Q}^{*}}{{\left({\theta }_1+N^{*}\right)}^2},\\ g_0=-\frac{k_0{\theta }_2^2{Q}^{*}\left(-3{\theta }_2^2+Q^{*2}\right)}{{\left({\theta }_2^2+Q^{*2}\right)}^3},g_1=\frac{k_0{\theta }_2^2\left({\theta }_2^4-6{\theta }_2^2{Q}^{*2}+Q^{*4}\right)}{{\left({\theta }_2^2+Q^{*2}\right)}^4},\\ g_2=\frac{2{\text{e}}^{-r_s{\tau }_s}{k}_0{\theta }_2^2{Q}^{*}\left(-3{\theta }_2^2+Q^{*2}\right)}{{\left({\theta }_2^2+Q^{*2}\right)}^3},g_3=-\frac{2{\text{e}}^{-r_s{\tau }_s}{k}_0{\theta }_2^2\left({\theta }_2^4-6{\theta }_2^2{Q}^{*2}+Q^{*4}\right)}{{\left({\theta }_2^2+Q^{*2}\right)}^4},\end{array}$

$\begin{array}{l}{g}_4=\frac{f_0{\theta }_1}{{\left({\theta }_1+N^{*}\right)}^2},g_5=-\frac{f_0{\theta }_1{Q}^{*}}{{\left({\theta }_1+N^{*}\right)}^3},\\ g_6=-\frac{f_0{\theta }_1}{{\left({\theta }_1+N^{*}\right)}^3},g_7=\frac{f_0{\theta }_1{Q}^{*}}{{\left({\theta }_1+N^{*}\right)}^4},\\ h_0=\frac{A{f}_0{\theta }_1{Q}^{*}}{{\left({\theta }_1+N^{*}\right)}^3},h_1=-\frac{A{f}_0{\theta }_1{Q}^{*}}{{\left({\theta }_1+N^{*}\right)}^4},\\ h_2=-\frac{A{f}_0{\theta }_1}{{\left({\theta }_1+N^{*}\right)}^2},h_3=\frac{A{f}_0{\theta }_1}{{\left({\theta }_1+N^{*}\right)}^3}.\end{array}$

Suppose $u_t\left(\theta \right)=u\left(t+\theta \right)$ with $u={\left(x\mathrm{,}y\right)}^{\text{T}}$ , system (4.1) is a DDEs with its phase space defined on the Banach space $C\left(\left[-{\tau }_N\mathrm{,0}\right]\mathrm{,}{R}^2\right)$ with super norm $\Vert u\Vert ={\max }_{-{\tau }_N\le \theta \le 0}\left|u\left(t+\theta \right)\right|$ . Based on the fundamental theory of DDEs, there exists the bounded variation matrix $\eta \left(\theta \right)$ to represent the linearized version of Equation (4.1) as its Rieze representation

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

with

$\text{d}\eta \left(\theta \right)=\left[\begin{array}{cc}{a}_1\delta \left(\theta \right)+a_2\delta \left(\theta +{\tau }_s\right)+K_2\delta \left(\theta +{\tau }_N\right)& c_1\delta \left(\theta \right)\\ b_1\delta \left(\theta +{\tau }_N\right)& -{\gamma }_N\delta \left(\theta \right)+b_2\delta \left(\theta +{\tau }_N\right)\end{array}\right]$ (4.3)

For $\varphi \in C$ , the solution operator of the linear Equation (4.2) is a strong continuous semigroup with its infinitesimal generator

$\mathcal{A}\varphi =\{\begin{array}{ll}\frac{\text{d}\varphi }{\text{d}\theta },\hfill & -{\tau }_N\le \theta <0\hfill \\ L\varphi ,\hfill & \theta =0\hfill \end{array}$ (4.4)

For $\psi \in C^{\mathrm{<mo>\; *\; </mo>}}$ with $C^{\mathrm{<mo>\; *\; </mo>}}\left(\left[\mathrm{0,}{\tau }_N\right]\mathrm{,}{R}^2\right)$ the conjugate space of C, the adjoint operator $A^{\mathrm{<mo>\; *\; </mo>}}$ of A is denoted as

${\mathcal{A}}^{*}\psi =\{\begin{array}{ll}-\frac{\text{d}\psi }{\text{d}s},\hfill & 0<s\le {\tau }_N\hfill \\ {\displaystyle {\int }_{-{\tau }_N}^0}\text{d}{\eta }^{\text{T}}\left(-s\right)\psi \left(s\right),\hfill & \theta =0\hfill \end{array}$ (4.5)

And for $\varphi \in C\mathrm{,}\psi \in C^{\mathrm{<mo>\; *\; </mo>}}$ , define the inner product as the following bilinear form

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

Furthermore, for $u_t\in C$ , we rewrite system (4.1) as its operator form of differential equation as follows

${u^{\prime }}_t=A{u}_t+R{u}_t$ (4.7)

with the nonlinear operator R defined as

$R{u}_t=\{\begin{array}{ll}0,\hfill & -{\tau }_N\le \theta <0,\hfill \\ F{u}_t,\hfill & \theta =0\hfill \end{array}$ (4.8)

with

$F{u}_t=\left(\begin{array}{c}{g}_0{x}^2+g_1{x}^3+g_{2}x{\left(t-{\tau }_s\right)}^2\\ +g_{3}x{\left(t-{\tau }_s\right)}^3+g_{4}xy+g_5{y}^2+g_{6}x{y}^2+g_7{y}^3,\\ b_{2}y\left(t-{\tau }_N\right)+h_{0}y{\left(t-{\tau }_N\right)}^2+h_{1}y{\left(t-{\tau }_N\right)}^3\\ +h_{2}y\left(t-{\tau }_N\right)x\left(t-{\tau }_N\right)+h_{3}y{\left(t-{\tau }_N\right)}^{2}x\left(t-{\tau }_N\right)\end{array}\right)$

With the assumption of Hopf bifurcation at the critical value, we define the finite characteristic roots set $\Lambda =\left\{i\omega \mathrm{,}-i\omega \right\}$ . The corresponding eigenvector function $\Phi$ for differential operator A and $\Psi$ for its adjoint operator $A^{\mathrm{<mo>\; *\; </mo>}}$ satisfies

$A\Phi =\left(\begin{array}{cc}i\omega & 0\\ 0& -i\omega \end{array}\right)\Phi ,A^{*}\Psi =\left(\begin{array}{cc}-i\omega & 0\\ 0& i\omega \end{array}\right)\Psi$ (4.9)

Denote the subspace P is expanded by the eigenvector function $\Phi$ , and the corresponding complementary subspace is denoted by Q, then decompose $C=P\oplus Q$ .

Consider the extended phase space $BC\left(\left[-{\tau }_N\mathrm{,0}\right]\mathrm{,}{R}^2\right)$ , which is left continuous with a possible jump at $\theta =0$ , for any $\phi \in BC$ , one has $\phi =\varphi +X_0\alpha$ with fundamental matrix

$X_0=\{\begin{array}{ll}0\hfill & -{\tau }_N\le \theta <0,\hfill \\ I_2\hfill & \theta =0\hfill \end{array}$ (4.10)

we define the projection mapping $\Pi \mathrm{<mo>\; :\; </mo>}BC\to P$ as

$\Pi \phi =\Phi \left(\theta \right)\langle \Psi ,\phi \rangle +\Phi \left(\theta \right)\Psi \left(0\right)X_0\xi$ (4.11)

then we have $Ker\Pi ⊊Q$ .

As stated in the paper, we depart the local center $W_{loc}^c$ as the sum of linear part belonging to center eigenspace P and nonlinear part belonging to complementary subspace Q taht is

$W_{loc}^c=\left\{\phi \mathrm{<mo>\; |\; </mo>}\phi =\mathrm{<mo>\; |\; </mo>}\Phi \left(\begin{array}{l}z\hfill \\ \bar{z}\hfill \end{array}\right)+h\left(\theta \mathrm{,}z\mathrm{,}\bar{z}\right)\right\}$ (4.12)

with $z\in C^2$ . It can be seen that $\left(\begin{array}{l}z\hfill \\ \bar{z}\hfill \end{array}\right)=\langle \Psi \mathrm{,}\phi \rangle$ , $h\left(\theta \mathrm{,}z\mathrm{,}\bar{z}\right)\in Ker\left(\Pi \right)$ , therefore, we get

$\langle \Psi \mathrm{,}h\left(\theta \mathrm{,}z\mathrm{,}\bar{z}\right)\rangle =0$ (4.13)

and

$\langle \Psi \mathrm{,}\left(\frac{\partial h}{\partial z}\left(z\left(t\right)\mathrm{,}\bar{z}\left(t\right)\right)\mathrm{,}\frac{\partial h}{\partial \bar{z}}\left(z\left(t\right)\mathrm{,}\bar{z}\left(t\right)\right)\right)\rangle =0$ (4.14)

The solution $u_t\left(\theta \right)=u\left(t+\theta \right)$ of operator differential Equation (4.7) on the center manifold can be written as

$u_t\left(\theta \right)=\Phi \left(\theta \right)\left(\begin{array}{l}z\left(t\right)\hfill \\ \bar{z}\left(t\right)\hfill \end{array}\right)+h\left(\theta \mathrm{,}z\left(t\right)\mathrm{,}\bar{z}\left(t\right)\right)$ (4.15)

and we have

$\begin{array}{l}\left(\Phi \left(\theta \right)+\left(\frac{\partial h}{\partial z}\left(\theta \mathrm{,}z\left(t\right)\mathrm{,}\bar{z}\left(t\right)\right)\mathrm{,}\frac{\partial h}{\partial \bar{z}}\left(\theta \mathrm{,}z\left(t\right)\mathrm{,}\bar{z}\left(t\right)\right)\right)\right)\stackrel{\dot{}}{z}\\ =\{\begin{array}{ll}\Phi \left(\theta \right)B\left(\begin{array}{l}z\left(t\right)\hfill \\ \bar{z}\left(t\right)\hfill \end{array}\right)+\frac{\partial h}{\partial \theta }\left(\theta \mathrm{,}z\left(t\right)\mathrm{,}\bar{z}\left(t\right)\right)\mathrm{,}\hfill & -\tau \le \theta <0\hfill \\ \Phi \left(0\right)B\left(\begin{array}{l}z\left(t\right)\hfill \\ \bar{z}\left(t\right)\hfill \end{array}\right)+Lh\left(\mathrm{0,}z\left(t\right)\mathrm{,}\bar{z}\left(t\right)\right)+F\left(\Phi \left(\theta \right)z\left(t\right)+h\left(\theta \mathrm{,}z\left(t\right)\mathrm{,}\bar{z}\left(t\right)\right)\right)\mathrm{,}\hfill & \theta =0\hfill \end{array}\end{array}$ (4.16)

Therefore, we have

$\Pi \left(\stackrel{\dot{}}{u}\left(t\right)-\mathcal{A}{u}_t-\mathcal{R}{u}_t\right)=\mathrm{0,}$ (4.17)

and

$\left(I-\Pi \right)\left(\stackrel{\dot{}}{u}\left(t\right)-\mathcal{A}{u}_t-\mathcal{R}{u}_t\right)=0$ (4.18)

By computation of Equation (4.17) and Equation (4.18), we have

$\left(\begin{array}{l}\stackrel{\dot{}}{z}\left(t\right)\hfill \\ \stackrel{\dot{}}{\bar{z}}\left(t\right)\hfill \end{array}\right)=B\left(\begin{array}{l}z\left(t\right)\hfill \\ \bar{z}\left(t\right)\hfill \end{array}\right)+\Psi \left(0\right)F\left(\Phi \left(\theta \right)\left(\begin{array}{l}z\left(t\right)\hfill \\ \bar{z}\left(t\right)\hfill \end{array}\right)+h\left(\theta ,z\left(t\right),\bar{z}\left(t\right)\right)\right)$ (4.19)

and

$\begin{array}{l}\frac{\partial h}{\partial z}\left(\theta ,z\left(t\right),\bar{z}\left(t\right)\right)\stackrel{\dot{}}{z}\left(t\right)=\mathcal{A}h\left(\theta ,z\left(t\right),\bar{z}\left(t\right)\right)\\ +\{\begin{array}{l}-\Phi \left(\theta \right)\Psi \left(0\right)F\left(\Phi \left(\theta \right)\left(\begin{array}{l}z\hfill \\ \bar{z}\hfill \end{array}\right)+h\left(\theta ,z,\bar{z}\right)\right),-\tau \le \theta <0\\ F\left(\Phi \left(0\right)\left(\begin{array}{l}z\hfill \\ \bar{z}\hfill \end{array}\right)+h\left(0,z,\bar{z}\right)\right)-\Phi \left(0\right)\Psi \left(0\right)F\left(\Phi \left(0\right)\left(\begin{array}{l}z\hfill \\ \bar{z}\hfill \end{array}\right)+h\left(0,z,\bar{z}\right)\right),\theta =0\end{array}\end{array}$ (4.20)

Set

$h\left(\theta \mathrm{,}z\left(t\right)\mathrm{,}\bar{z}\left(t\right)\right)=W_{20}\left(\theta \right)z^2+W_{11}\left(\theta \right)z\bar{z}+W_{02}\left(\theta \right){\bar{z}}^2+\cdots$ (4.21)

and at $\theta =0$ ,

$\begin{array}{c}F\left(\theta ,z\left(t\right),\bar{z}\left(t\right)\right)=f_{20,0}\left(0\right)z^2+f_{11,0}\left(0\right)z\bar{z}+f_{02,0}\left(0\right){\bar{z}}^2\\ +f_{20,1}\left(-{\tau }_s\right)z^2+f_{11,1}\left(-{\tau }_s\right)z\bar{z}+f_{02,1}\left(-{\tau }_s\right){\bar{z}}^2\\ +f_{20,2}\left(-{\tau }_N\right)z^2+f_{11,2}\left(-{\tau }_N\right)z\bar{z}+f_{02,2}\left(-{\tau }_N\right){\bar{z}}^2\\ +f_{21,0}\left(0\right)z^2\bar{z}+f_{21,1}\left(-{\tau }_s\right)z^2\bar{z}+f_{21,2}\left(-{\tau }_N\right)z^2\bar{z}+\cdots \end{array}$ (4.22)

Substitute it into Equation (4.20), then we get

$\begin{array}{l}\mathcal{A}{W}_{20}\left(\theta \right)=2i\omega W_{20}\left(\theta \right)+\Phi \left(\theta \right)\Psi \left(0\right)\left(f_{20,0}\left(0\right)+f_{20,1}\left(-{\tau }_s\right)+f_{20,2}\left(-{\tau }_N\right)\right)\\ \mathcal{A}{W}_{11}\left(\theta \right)=\Phi \left(\theta \right)\Psi \left(0\right)\left(f_{11,0}\left(0\right)+f_{11,1}\left(-{\tau }_s\right)+f_{11,2}\left(-{\tau }_N\right)\right)\\ \mathcal{A}{W}_{02}\left(\theta \right)=-2i\omega W_{02}\left(\theta \right)+\Phi \left(\theta \right)\Psi \left(0\right)\left(f_{02,0}\left(0\right)+f_{02,1}\left(-{\tau }_s\right)+f_{02,2}\left(-{\tau }_N\right)\right)\end{array}$ (4.23)

Equation (4.23) satisfies the initial value condition

$\begin{array}{l}L{W}_{20}\left(\theta \right)=2i\omega W_{20}\left(0\right)+\Phi \left(0\right)\Psi \left(0\right)\left(f_{20,0}\left(0\right)+f_{20,1}\left(-{\tau }_s\right)+f_{20,2}\left(-{\tau }_N\right)\right)\\ -\left(f_{20,0}\left(0\right)+f_{20,1}\left(-{\tau }_s\right)+f_{20,2}\left(-{\tau }_N\right)\right)\\ L{W}_{11}\left(\theta \right)=\Phi \left(0\right)\Psi \left(0\right)\left(f_{11,0}\left(0\right)+f_{11,1}\left(-{\tau }_s\right)+f_{11,2}\left(-{\tau }_N\right)\right)\\ -\left(f_{11,0}\left(0\right)+f_{11,1}\left(-{\tau }_s\right)+f_{11,2}\left(-{\tau }_N\right)\right)\\ L{W}_{02}\left(\theta \right)=-2i\omega W_{02}\left(0\right)+\Phi \left(0\right)\Psi \left(0\right)\left(f_{02,0}\left(0\right)+f_{02,1}\left(-{\tau }_s\right)+f_{02,2}\left(-{\tau }_N\right)\right)\\ -\left(f_{02,0}\left(0\right)+f_{02,1}\left(-{\tau }_s\right)+f_{02,2}\left(-{\tau }_N\right)\right)\end{array}$ (4.24)

Based on Equations (4.23) and (4.24), by integral Equation (4.23) with the given initial condition, then the coefficent matrix $W_{20}\left(\theta \right)\mathrm{,}{W}_{11}\left(\theta \right)\mathrm{,}{W}_{02}\left(\theta \right)$ can be easily computed. Hence after, the expression on the local center manifold is approximately derived. Substitute it into Equation (4.19), one also gets

$\begin{array}{c}\left(\begin{array}{l}\stackrel{\dot{}}{z}\left(t\right)\hfill \\ \stackrel{\dot{}}{\bar{z}}\left(t\right)\hfill \end{array}\right)=B\left(\begin{array}{l}z\left(t\right)\hfill \\ \bar{z}\left(t\right)\hfill \end{array}\right)+\Psi \left(0\right)\left(f_{20,0}\left(0\right)+f_{20,1}\left(-{\tau }_s\right)+f_{20,2}\left(-{\tau }_N\right)\right)\left(\begin{array}{l}{z}^2\left(t\right)\hfill \\ {\bar{z}}^2\left(t\right)\hfill \end{array}\right)\\ +\Psi \left(0\right)\left(f_{11,0}\left(0\right)+f_{11,1}\left(-{\tau }_s\right)+f_{11,2}\left(-{\tau }_N\right)\right)\left(\begin{array}{l}z\left(t\right)\bar{z}\left(t\right)\hfill \\ z\left(t\right)\bar{z}\left(t\right)\hfill \end{array}\right)\\ +\Psi \left(0\right)\left(f_{02,0}\left(0\right)+f_{02,1}\left(-{\tau }_s\right)+f_{02,2}\left(-{\tau }_N\right)\right)\left(\begin{array}{l}{\bar{z}}^2\left(t\right)\hfill \\ z^2\left(t\right)\hfill \end{array}\right)\\ +\Psi \left(0\right)\left(f_{21,0}\left(0\right)+f_{21,1}\left(-{\tau }_s\right)+f_{21,2}\left(-{\tau }_N\right)\right)\left(\begin{array}{l}{z}^2\left(t\right)\bar{z}\left(t\right)\hfill \\ z^2\left(t\right)\bar{z}\left(t\right)\hfill \end{array}\right)\end{array}$ (4.25)

Therefore, we have expression

$\stackrel{\dot{}}{z}\left(t\right)=i\omega z\left(t\right)+g_{20}{z}^2\left(t\right)+g_{11}z\left(t\right)\bar{z}\left(t\right)+g_{02}{\bar{z}}^2\left(t\right)+g_{21}{z}^2\left(t\right)\bar{z}\left(t\right)$ (4.26)

with

$\begin{array}{l}\left(\begin{array}{c}{g}_{20}\\ {\bar{g}}_{20}\end{array}\right)=\Psi \left(0\right)\left(f_{20,0}\left(0\right)+f_{20,1}\left(-{\tau }_s\right)+f_{20,2}\left(-{\tau }_N\right)\right),\\ \left(\begin{array}{c}{g}_{11}\\ {\bar{g}}_{11}\end{array}\right)=\Psi \left(0\right)\left(f_{11,0}\left(0\right)+f_{11,1}\left(-{\tau }_s\right)+f_{11,2}\left(-{\tau }_N\right)\right),\\ \left(\begin{array}{c}{g}_{02}\\ {\bar{g}}_{02}\end{array}\right)=\Psi \left(0\right)\left(f_{02,0}\left(0\right)+f_{02,1}\left(-{\tau }_s\right)+f_{02,2}\left(-{\tau }_N\right)\right),\\ \left(\begin{array}{c}{g}_{21}\\ {\bar{g}}_{21}\end{array}\right)=\Psi \left(0\right)\left(f_{21,0}\left(0\right)+f_{21,1}\left(-{\tau }_s\right)+f_{21,2}\left(-{\tau }_N\right)\right)\end{array}$

Furthermore, the normal form of system (4.1) is written as

$\stackrel{\dot{}}{z}\left(t\right)=i\omega z\left(t\right)+C\left(0\right)z^2\left(t\right)\bar{z}\left(t\right)$ (4.27)

with

$C\left(0\right)=g_{21}-\frac{g_{20}{g}_{11}}{i\omega }+\frac{{\left|g_{11}\right|}^2}{i\omega }+\frac{{\left|g_{02}\right|}^2}{3i\omega }$

5. Conclusion

The dynamics of stem cells were discussed in stem cells renewal procedure, with the discipline between proliferating stem cells and oscillating cyclic neutropenia, the coexistence phenomena of periodical solutions were exhibited. It was seen that chaos occurs via scenario cascades of period-doubling bifurcation phenomena with the control strategy via state feedback with time delay. Hopf bifurcation criterion was derived and a pair of imaginary roots cross the imaginary axis if and only if the transversal condition was not zero. Based on the center manifold theory, the normal form computation was done based on Schmidt-Lyapunov reduction technique. The continuous work of limit cycle was expanded with the benefits to adopt DDE-Biftool technique. The periodical oscillation phenomena of neutropenia were discussed which can provide some evidence further in the diagnosis treatment.

Conflicts of Interest

The author declares no conflicts of interest regarding the publication of this paper.

References

[1]	Adimy, M., Crauste, F. and Ruan, S. (2005) A Mathematical Study of the Hematopoiesis Process with Applications to Chronic Myelogenous Leukemia. SIAM Journal on Apllied Mathematics, 65, 1328-1352. https://doi.org/10.1137/040604698
[2]	Adimy, M., Crauste, F. and Ruan, S. (2005) Stability and Hopf Bifurcation in a Mathematical Model of Pluripotent Stem Cell Dynamics. Nonlinear Analysis: Real World Applications, 6, 651-670. https://doi.org/10.1016/j.nonrwa.2004.12.010
[3]	Daniel C. and Humphries, A.R. (2017) Dynamics of a Mathematical Hematopoietic Stem-Cell Population Model. ArXiv: 1712.08308. https://arxiv.org/abs/1712.08308
[4]	Fortin, P. and Mackey, M.C. (1999) Periodic Chronic Myelogenous Leukemia: Spectral Analysis of Blood Cell Counts and Aetiological Implications. British Journal of Haematology, 104, 336-345. https://doi.org/10.1046/j.1365-2141.1999.01168.x
[5]	Mackey, M.C. (1978) United Hypothesis for the Origin of Aplasric Anemia and Periodic Hematopoiesis. Blood, 51, 941-956.
[6]	Haurie, C., Dale, D.C. and Mackey, M.C. (1998) Cyclical Neutropenia and Other Periodic Hematological Disorders: A Review of Mechanisms and Mathematical Models. Blood, 92, 2629-2640.
[7]	Alaoui, T.H., Yafia, R. and Aziz-Alaoui, M.A. (2007) Hopf Bifurcation Direction in a Delayed Hematopoietic Stem Cells Model. Arab Journal of Mathematics and Mathematical Sciences, 1, 35-50.
[8]	Zhuge, C.J., Lei, J.Z. and Mackey, M.C. (2012) Neutrophil Dynamics in Response to Chemotherapy and G-SCF. Journal of Theory Biology, 293, 111-120. https://doi.org/10.1016/j.jtbi.2011.10.017
[9]	Brooks, G., Langlois, G.P., Lei, J.Z. and Mackey, M.C. (2012) Neutrophil Dynamics after Chemotherapy and G-CSF: The Role of Pharmacokinetics in Shaping the Response. Journal of Theory Biology, 315, 97-109. https://doi.org/10.1016/j.jtbi.2012.08.028
[10]	Lei, J.Z. and Mackey, M.C. (2014) Understanding and Treating Cytopenia through Mathematical Modeling. In: Corey, S., Kimmel, M. and Leonard, J., Eds., A Systems Biology Approach to Blood, Advances in Experimental Medicine and Biology, Vol. 844, Springer, New York, 279-302. https://doi.org/10.1007/978-1-4939-2095-2_14
[11]	Ma, S.Q. (2014) Stability and Bifurcation Analysis of a Type of Hematopoietic Stem Cell Model. International Journal of Modern Nonlinear Theory and Application, 10, 13-27. https://doi.org/10.4236/ijmnta.2021.101002
[12]	Ma, S.Q. (2017) Periodic Oscillation in Neutrophil Models with Time Delays. International Journal of Modern Nonlinear Theory and Application, 6, 119-133. https://doi.org/10.4236/ijmnta.2017.64011
[13]	Ma, S.Q. (2019) Hopf Bifurcation of a Type of Neuron Model with Multiple Time Delays. International Journal of Bifurcation and Chaos, 29, Article ID: 1950163. https://doi.org/10.1142/S0218127419501633
[14]	Ma, S.Q. (2022) Bifurcation Analysis of Periodic Oscillation in a Hematopoietic Stem Cells Model with Time Delay Control. Mathematical Problems in Engineering, 2022, Article ID: 7304280. https://doi.org/10.1155/2022/7304280
[15]	Sieber, J., Engelborghs, K., Luzyanina, T., Samaey, G. and Roose, D. (2016) DDE-BIFTOOL v. 3.1.1 Manual—Bifurcation Analysis of Delay Differential Equations. Dynamical Systems. ArXiv: 1406.7144. https://arxiv.org/abs/1406.7144
[16]	Engelborghs, K., Luzyanina, T. and Samae, G. (2001) DDE-BIFTOOL v. 2.00: A Matlab Package for Bifurcation Analysis of Delay Differential Equations. Technical Report TW330.
[17]	Verheyden, K., Luzyanina, T. and Roose, D. (2004) Location and Numerical Preservation of Characteristic Roots of Delay Differential Equations by LMS Methods. Technical Report TW-382.
[18]	Kuzenetsov, Y.A. (2004) Elements of Applied Bifurcation Theory Applied Mathematical Sciences. Vol. 112, Springer, New York.
[19]	Hale, J.K. and Lunel, S.M.V. (1993) Introduction of Functional Differential Equations. In: Bloch, A., Charles, L., Epstein, A.G. and Greengard, L., Eds., Applied Mathematical Sciences, Vol. 99, Springer, New York, 1-10. https://doi.org/10.1007/978-1-4612-4342-7
[20]	Hassard, B.D., Kazarinoff, N.D. and Wan, Y.H. (1981) Theory and Applications of Hopf Bifurcation. Combridge University Press, Cambridge.
[21]	Faria, T. and Magalhaes, M.L. (1995) Normal Form for Retarded Functional Differential Equations with Parameters and Applications to Hopf Bifurcation. Journal of Differential Equations, 122, 181-200. https://doi.org/10.1007/978-1-4612-4342-7