Generalized Hopf Bifurcation in a Delay Model of Neutrophil Cells Model

Suqi Ma¹, S. J. Hogan²

¹Department of Mathematics, China Agricultural University, Beijing, China.
²Department of Mathematics, University of Bristol, Bristol, UK.

DOI: 10.4236/ijmnta.2024.132002                

Abstract

The DDE-Biftool software is applied to solve the dynamical stability and bifurcation problem of the neutrophil cells model. Based on Hopf point finding with the stability property of the equilibrium solution loss, the continuation of the bifurcating periodical solution starting from Hopf point is exploited. The generalized Hopf point is tracked by seeking for the critical value of free parameter of the switching phenomena of the open loop, which describes the lineup of bifurcating periodical solutions from Hopf point. The normal form near the generalized Hopf point is computed by Lyapunov-Schimdt reduction scheme combined with the center manifold analytical technique. The near dynamics is classified by geometrically different topological phase portraits.

Keywords

Generalized Hopf Bifurcation, DDE-Biftool Software, Norm Form

1. Introduction

Mathematical modelling is powerful designed in understanding hematological stem cells (HSCs) regulation mechanism to explore its complex behaviors underlying dynamical bifurcation mechanism ([1]-[3]). Normally, the circulating white cells, red cells and platelets are maintained at a homeostatic level, hence no obvious evidence of any oscillation behavior occurring. However, the neutrophil dynamics is further discussed based on the circulating cells dynamics with the discipline of the cell lineage. When neutrophil numbers fall to sufficiently low levels, the periodic oscillation phenomena of neutrophil counts are observed ([4]-[6]). The often hematopoietic diseases are modeled with several delays history in one period, which describes neutrophil dynamics by DDEs with nonlinearity ([4] [7] [8]).

Even brilliantly with some known models, the administration of the periodic diseases and the production of blood cells number is partially understood yet. Biologically and mathematically, the outline of the circulating cells feedback mechanism is schemed, which originated from HSCs number, as shown in Figure 1. In mammals, hematopoietic stem cells can proliferate and differentiate into one of the three major cells, and end with the release of the mature blood cells into the circulation ([9]-[11]). The computation models of HSCs dynamics obey the Hill rules with the introduction rate manifesting time delay feedback in proliferating phase, which is described as ([9]):

$Q^{\prime }\left(t\right)=-\beta \left(Q\right)Q+2{\text{e}}^{-{\gamma }_0{\tau }_s}\beta \left(Q\left(t-{\tau }_s\right)\right)Q\left(t-{\tau }_s\right)$ (1.1)

Here in Q denotes HSCs number, and giving Hill function ([10] [12] [13])

$\beta \left(Q\right)=\frac{{\theta }^n}{Q^n+{\theta }^n}$ (*)

With all sorts of bifurcating solutions, and even chaos, system (1.1) has shed light on the pathophysiology dynamically development of the related hematopoietic diseases. In this paper, we also get a brand new function to replace of Hill function,

$\beta \left(Q\right)={\beta }_0\frac{1-p^n}{Q^n\left(1+p^n+2p^{n/2}\cos \left(Q\right)\right)},$ (F)

with $p={\text{e}}^{-\theta }$ .

Figure 1. The outline of the two department cells model of neutropenia dynamics.

The control of the circulating blood cells is mediated by a delayed negative feedback mechanism, such as granulocyte colony stimulating factor (G-CSF) for the white blood cells ([14]-[16]). As is known, the apoptosis is the programmed death of cells happening during the proliferation and differentiation process. In paper ([5] [17]-[19]), the feedback control of G-CSF has supposed to experience two time segments ${\tau }_{NP}$ and ${\tau }_{NM}$ , and the amplifying coefficients in the peripheral cells cycle are set forth, which is formulated as

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

With the guide of the outline schemed in Figure 1, the studied neutrophil model is put forth as the follows

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

System (1.3) is of highly nonlinearity of DDEs with history of multi-delays, wherein with positive integer $m=1,n=2$ . Hill function of $K\left(N\right)$ is often given as ([7] [20]).

$K\left(N\right)=f_0\frac{{\theta }_1^m}{{\theta }_1^m+N^m}$ (1.4)

We set $c\left({\tau }_s\right)=1-{\displaystyle {\int }_0^{{\tau }_s}}{g}_{\gamma }^j\left(s\right)\text{d}s$ , which is increasing with time delay ${\tau }_s$ . In general, the function $g_{\gamma }^j\left(s\right)$ denotes Gamma-distribution with formula

$g_{\gamma }^j\left(s\right)=\frac{s^{j-1}{\gamma }_0^j}{{\text{e}}^{{\gamma }_{0}s}\Gamma \left(j\right)}$ (1.5)

where in ${\gamma }_0$ is positive constant and $j\in N_0$ , $\Gamma \left(j\right)=\left(j-1\right)\mathrm{<mo>!</mo>}$ . It is easily seen that if $j=1$ , $c\left({\tau }_s\right)={\text{e}}^{-{\gamma }_0{\tau }_s}$ , and herein after we choose $j=2$ to get

$c\left({\tau }_s\right)={\text{e}}^{-{\gamma }_0{\tau }_s}+{\gamma }_0{\tau }_s{\text{e}}^{-{\gamma }_0{\tau }_s}$ (1.6)

and ${\lim }_{{\tau }_s\to 0}c\left({\tau }_s\right)=1$ .

The stability property of system loss as Hopf bifurcation is observed, and double Hopf bifurcation occurs as Hopf lines intersect itself with varying free parameters. We discuss double Hopf bifurcation further in another paper. The often analytical technique of Hopf bifurcation is completed by the imaginary roots computation of the related characteristic equation of the linearized system ([21] [22]). And the bifurcation direction of the periodical solutions arising from Hopf point is calculated by the normal form on the center manifold ([23] [24]). With the Lyapunov-Schmidt reduction method referred, the near dynamics of the double Hopf point is classified geometrically in topological phase portraits by the derived normal form. Generalized Hopf bifurcation happens at the switch point of super-Hopf bifurcation and sub-Hopf bifurcation. And in general, the calculated first Lyapunov coefficient in normal form is zero ([25]-[27]). Hence after the expanded normal form to five degree is demanded and the limit point cycle bifurcation line is computed by normal form coefficients.

DDE-Biftool is an artificial hand-tool in DDEs bifurcation computation and has the big compatibility with every sort of delay systems ([28] [29]). As another preserved method, DDE-Biftool discovers the generalized Hopf point by a special observation of loops of limit cycles. With free parameter ${\gamma }_0$ varying, the bifurcating periodical solutions arise from Hopf point. However, near GH point, the super-critical Hopf point and the sub-critical Hopf point are close to each other, and surely candidates for a lineup of the limit cycles continued from one Hopf point to another. The featured open loops of lineup of periodical solutions continuation are described. At GH point, the collide of stable limit cycles and unstable limit cycles also happens, which extends a limit point cycle bifurcation line from GH point. Whether the limit point cycle bifurcation line is in coincidence with the above referred bifurcation line calculated by normal form coefficients?

The whole paper is organized as follows. In Section 2, the continuation of equilibrium solution is done with varying free parameter ${\gamma }_0$ , and the bifurcating periodical solution arising Hopf point is finished. In Section 3, the preserved method of discover generalized Hopf point is programmed, which display both the open loops of lineup solutions nearby and the closed lineup solutions of periodical oscillation at the critical values. In Section 4, the programmed lineup solutions of limit cycle are also observed by varying parameter $k_{\delta }$ and ${\tau }_N$ , which are either open loops or closed loops continued by DDE-Biftool; In Section 5, the normal form of GH are computed by dimensional reduction method.

2. Hopf Bifurcation Analysis of System

Suppose $E^{\mathrm{<mo>*</mo>}}=\left(Q^{\mathrm{<mo>*</mo>}}\mathrm{,}{N}^{\mathrm{<mo>*</mo>}}\right)$ is the equilibrium solution of system (1.3), then it is satisfied the right side of the equations. Applying DDE-Biftool with sys_ rhs command to setup system (1.3), the continuation solution of the equilibrium is often carried out. The stability property of equilibrium solutions is analyzed by br_ stable command. As shown in Figure 2(a), applying DDE-Biftool software, we continue the equilibrium solution and compute the stability property as varying ${\gamma }_0$ . Hopf points are found due to the loss of stability of the equilibrium solution as varying free parameter. We use ${\gamma }_{\mathrm{0<mo>*</mo>}}$ to denote the threshold of Hopf bifurcation and Hopf points are listed

(a) (b)

Figure 2. The continuation of the equilibrium solution and bifurcating periodical solutions arise from Hopf points v.s. ${\gamma }_0$ , the other parameters are fixed as $k_{\delta }=0.143$ , ${\tau }_N=0.577$ . (a) The equilibrium solutions and Hopf points; (b) The continuation of periodical solutions with maximal magnitude and minimal magnitude versus ${\gamma }_0$ .

The bifurcating periodical solution arise from Hopf point, as shown in Figure 2(b), in which using solid line to denote the maximal amplitude and dash line to represent the minimal amplitude of periodical solution. The bifurcating periodical solution is devastating to find nearby possible generalized Hopf(GH) bifurcation.

3. Generalized Hopf Point with ${\gamma }_0$

(a) (b)

Figure 3. The detection of GH point and the continuation of periodical solutions. (a) The imaginary roots calculation of GH point; (b) The amplitude picture with respect to ${\gamma }_0$ is continued.

Figure 4. (a) The continuation of periodical solution near GH point in X-Z view; (b) The picture near GH point in X-Y-Z view; (c) The continuation of periodical solution of GH point in X-Z view; (d) The picture of GH point in X-Y-Z view.

We give the preserved method to detect Generalized Hopf point by programmed code in DDE-Biftool. That is, the periodical solution from left and right Hopf point on equilibrium solutions line collide. Therefore, with the amination to detect GH point, we pick up two periodical solutions to do continuation of periodical solutions as varying ${\gamma }_0$ , as shown in Figure 6. Originated from the following two periodical solutions psol1 and psol2 respectively,

The continue of periodical solutions is completed from Hopf point to Hopf point with respect to ${\gamma }_0$ . As shown in Figure 3(a), the nearby GH are found with

As seen in Figure 3(b), the periodical solution continued from GH is completed. The amplitude of periodical solutions versus ${\gamma }_0$ is shown which hints both subcritical Hopf bifurcation and supercritical bifurcation occurring. However, the calculation of bifurcating direction of periodical solution of GH point should be further discussed subsequently.

The line up of the bifurcating periodical solutions is looped by the continuation method in DDE-Biftool. The red color and blue color separate the periodical solutions before and after the limit point cycle bifurcation. As seen in Figure 4(a) and Figure 4(b), the open loop of line up of solution is obtained with the continuation from the originated periodical solution “psol1”; however, the closed loop of line up of solution which started from the solution “psol2” gets into sight of Figure 4(c) and Figure 4(d). Hence GH point is found as listed above. Other parameter are fixed as ${\gamma }_N=2.5$ , ${\eta }_{NP}=2.542$ , $f_0=0.8$ , ${\theta }_1=0.36$ , $r_s=0.0631$ , ${\tau }_{NP}=5$ , ${\beta }_0=2.3$ , $p=0.43$ .

4. GH Point with Parameter ${\tau }_N$ and $k_{\delta }$

It is claimed that GH point is found if a super-critical Hopf point collide with a sub-critical Hopf point. It is concluded that before and after threshold value of GH point, the lineup of limit cycle by continuation method manifests different routes design. Either lineup of periodical solutions from Hopf point continuation to Hopf point, or only lineup of periodical solutions from limit point cycle to limit point cycle without Hopf point found. For example, we choose the following two different periodical solutions, which respectively continued to be the referred different lineup solutions of system. With $k_{\delta }=0.0124$ , the continuation of limit cycle is done with varying delay ${\tau }_N$ , the closed loop of lineup of periodical solutions are observed. We increase the value of $k_{\delta }$ until attach at GH point, about $k_{\delta }=0.013431205338737$ , the loop of lineup of periodical solutions shows connecting cusp point at GH point. The value of $k_{\delta }$ is increased further, the loop of lineup of periodical solutions is done, which is open lineup loop. As shown in Figure 5, three open loop of lineup solutions are observed when varying ${\tau }_N$ with $k_{\delta }=\mathrm{0.0124,0.0133,0.0134}$ ; three closed loop of lineup solutions are also exhibited with $k_{\delta }=\mathrm{0.0135,0.0137,0.0141}$ .

At GH point listed as following, we setup the new solutions with profile method (example is illustrated at DDE-Biftool), the continuation of lineup of limit cycles is shown in Figure 6, as varying $k_{\delta }$ .

Figure 5. The loops set of the continuation of periodical solutions as increasing free parameter $k_{\delta }$ . The lineup solutions are obtained by varying time delay ${\tau }_N$ .

Figure 6. The loops set of the continuation of periodical solutions as increasing free parameter $k_{\delta }$ . The lineup solutions are obtained by varying time delay ${\tau }_N$ .

5. The Norm Form of GH Point

As discussed in the above sections, the generalized Hopf bifurcation occurs at $\left(k_{\delta 0}\mathrm{,}{\tau }_{N0}\right)$ with other parameter fixed. On Hopf line, the generalized Hopf point separates the supercritical Hopf points from the subcritical Hopf points. Since it happens with a limit point cycle bifurcation line arise from bifurcation point, two limit cycles coexist due to GH bifurcation. In Normal form, we can compute the first lyapunov exponent equals zero. The bifurcating periodical solutions arising from GH point is determined by the second Lyapunov exponent, which is stable if the exponent is less than zero, and otherwise versus. Applying Lyapunov-Schimdt reduction method, the norm form is often computed by the projection of solution operator on the phase space into the center manifold, and the bifurcating direction is determined.

We set $x=Q-Q^{*},y=N-N^{*}$ , then applying the perturbation parameters method, the trunction system of system (1.3) with Taylor expansion to fifth order is written as

$U\left(t\right)=L{U}_t+F{U}_t$ (5.1)

With $U={\left(x\mathrm{,}y\right)}^{\text{T}}$ , and

$L{U}_t={\displaystyle {\int }_{-{\tau }_N}^0\text{d}\eta \left(\theta \right)\varphi \left(\theta \right)}$ (5.2)

With $\eta \mathrm{<mo>:</mo>}\left[-{\tau }_N\mathrm{,0}\right]\to R^2$ is a bounded variation function vector. We also have

$F{U}_t={\displaystyle \sum _{i+j+k\ge 2}}{B}_{ijk}\left(U_t{\left(0\right)}^i{U}_t{\left(-{\tau }_s\right)}^j{U}_t{\left(-{\tau }_s\right)}^k\right)$

where in the above expression represents the multi-linear form style as often referred.

System (5.1) is defined on the phase space of a Banach space $C=\left[-{\tau }_N\mathrm{,0}\right]\to R^2$ with super norm $\Vert \varphi \Vert ={\sup }_{-{\tau }_N\le \theta \le 0}\left|\varphi \left(\theta \right)\right|$ , with $U_t\left(\theta \right)=U\left(t+\theta \right)$ and $U_t\in C$ .

As often calculated in Section 3, the linearized system of Equation (1.3) has eigen roots with $\pm i\omega$ and eigenvector $q{\text{e}}^{i\omega \theta }\mathrm{,}\overline{q}{\text{e}}^{-i\omega \theta }$ for $-{\tau }_N\le \theta \le 0$ . Based on the fundamental theory of DDEs ([20]), the solutions of the linear operator are a strong continuous semigroup which has infinitesimal generator

$\mathcal{A}\varphi =\{\begin{array}{ll}\frac{\text{d}\varphi }{\text{d}\theta }\mathrm{,}\hfill & -{\tau }_N\le \theta <\mathrm{0,}\hfill \\ ℒ\varphi \left(\theta \right)\mathrm{,}\hfill & \theta =0\hfill \end{array}$ (5.3)

For $\varphi \in C\left(\left[-{\tau }_N\mathrm{,0}\right]\to R^2\right)$ .

For $\psi \in C^{\mathrm{<mo>*</mo>}}\left(\left[\mathrm{0,}{\tau }_N\right]\to R^2\right)$ , the adjoint operator is also defined as

${\mathcal{A}}^{\text{*}}\psi =\{\begin{array}{ll}-\frac{\text{d}\psi }{\text{d}s},\hfill & 0<s\le {\tau }_N,\hfill \\ -{ℒ}^{\text{*}}\psi \left(-s\right),\hfill & s=0.\hfill \end{array}$ (5.4)

For $\varphi \in C\mathrm{,}\psi \in C^{\mathrm{<mo>*</mo>}}$ , we define the bilinear function of inner product as

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

Therefore, since it is satisfied with

$\begin{array}{lll}\mathcal{A}q{\text{e}}^{i\omega \theta }=i\omega q{\text{e}}^{i\omega \theta },\hfill & \mathcal{A}\overline{q}{\text{e}}^{-i\omega \theta }=-i\omega \overline{q}{\text{e}}^{-i\omega \theta },\hfill & \text{for}-{\tau }_N\le \theta \le 0,\hfill \\ {\mathcal{A}}^{*}p{\text{e}}^{i\omega s}=-i\omega p{\text{e}}^{i\omega s},\hfill & {\mathcal{A}}^{*}\overline{p}{\text{e}}^{-i\omega s}=i\omega \overline{p}{\text{e}}^{-i\omega s},\hfill & \text{for}0\le s\le {\tau }_N.\hfill \end{array}$ (5.6)

The center space is defined with its bases as $P=span\left\{q{\text{e}}^{i\omega \theta }\mathrm{,}\overline{q}{\text{e}}^{-i\omega \theta }\right\}$. Suppose $Q\subset C$ being the complementary subspace, and we can decompose the phase space C into its direct summation $C=P\oplus Q$ . We define the base function $\Phi \left(\theta \right)=\left(q{\text{e}}^{i\omega \theta }\mathrm{,}\overline{q}{\text{e}}^{-i\omega \theta }\right)$ for $-{\tau }_N\le \theta \le 0$ , and also the base function $\Psi \left(s\right)$ in its conjugate space $C^{\mathrm{<mo>*</mo>}}$ , $\Psi \left(s\right)=\left(p{\text{e}}^{i\omega s},\overline{\psi }\left(s\right)=\overline{p}{\text{e}}^{-i\omega s}\right)$ for $0\le s\le {\tau }_N$ . Hence, we have

$\langle \Psi \mathrm{,}\Phi \rangle =1$ (5.7)

For $-{\tau }_N\le \theta \le 0$ , make the axis transformation

$U_t=\Phi \left(\theta \right)\left(\begin{array}{l}z\left(t\right)\hfill \\ \overline{z}\left(t\right)\hfill \end{array}\right)+W_t$ (5.8)

with $U_t={\left(x_t\mathrm{,}{y}_t\right)}^{\text{T}}$ .

The axis $z\left(t\right)$ is calculated with formula

$z\left(t\right)=\langle p{\text{e}}^{i\omega s},U_t\rangle ,\overline{z}\left(t\right)=\langle \overline{p}{\text{e}}^{-i\omega s},U_t\rangle$ (5.9)

and $W_t=U_t-z\left(t\right)q{\text{e}}^{i\omega \theta }-\overline{z}\left(t\right)\overline{q}{\text{e}}^{-i\omega \theta }$.

We define the projection operator on the center manifold $\Pi \mathrm{,}C\to P$ by

$\Pi U_t=\Phi \left(\theta \right)\langle \Psi \left(s\right),U_t\rangle$ (5.10)

and the complementary operator $I-\Pi \mathrm{,}C\to Q$ , wherein $Q\subset Ker\Pi$ . Then on the center manifold, one derives

$W_t=\left(I-\Pi \right)U_t$ (5.11)

Then by system (1.3), one gets that

$\left(\begin{array}{l}{z}^{\prime }\left(t\right)\hfill \\ {\overline{z}}^{\prime }\left(t\right)\hfill \end{array}\right)=B\left(\begin{array}{l}z\left(t\right)\hfill \\ \overline{z}\left(t\right)\hfill \end{array}\right)+{\overline{\Psi }}^{\text{T}}\left(0\right)F\left(z\left(t\right){\Phi }_1\left(\theta \right)+\overline{z}\left(t\right){\Phi }_2\left(\theta \right)+W_t\right)$ (5.12)

and $B=\left(\begin{array}{cc}i\omega & 0\\ 0& -i\omega \end{array}\right)$ . Set

(5.13)

for $m=\mathrm{2,3,}\cdots$ .

Combined with center manifold analytical technique, we write

$W_t=H_{20}\left(\theta \right)z^2+H_{11}\left(\theta \right)z\overline{z}+H_{02}\left(\theta \right){\overline{z}}^2+\cdots$ (5.14)

By definition (5.11), one has

${W^{\prime }}_t=\{\begin{array}{l}\mathcal{A}{W}_t\left(\theta \right)-\Phi \left(\theta \right){\Psi }^{\text{T}}\left(0\right)F\left(z\left(t\right){\Phi }_1\left(\theta \right)+\overline{z}\left(t\right){\Phi }_2\left(\theta \right)+W_t\right),-{\tau }_N\le \theta <0\hfill \\ \mathcal{A}{W}_t\left(\theta \right)+F-\Phi \left(\theta \right){\Psi }^{\text{T}}\left(0\right)F\left(z\left(t\right){\Phi }_1\left(\theta \right)+\overline{z}\left(t\right){\Phi }_2\left(\theta \right)+W_t\right),\theta =0\hfill \end{array}$ (5.15)

Therefore, we differentiate $W_t$ with respect to t to get

$\begin{array}{l}\mathcal{A}{H}_{20}\left(\theta \right)=2i\omega H_{20}\left(\theta \right)+{\Phi }_1\left(\theta \right)g_{200}+{\Phi }_2\left(\theta \right){\overline{g}}_{020},\\ \mathcal{A}{H}_{11}\left(\theta \right)={\Phi }_1\left(\theta \right)g_{110}+{\Phi }_2\left(\theta \right){\overline{g}}_{110},\\ \mathcal{A}{H}_{02}\left(\theta \right)=-2i\omega H_{02}\left(\theta \right)+{\Phi }_1\left(\theta \right)g_{020}+{\Phi }_2\left(\theta \right){\overline{g}}_{200}\end{array}$ (5.16)

With initial value condition determined by the case $\theta =0$ in Equation (5.15).

Integral Equation (5.16) with respect to $\theta$ to get the coefficients $H_{20}\left(\theta \right),H_{11}\left(\theta \right),H_{02}\left(\theta \right)$ in the expression of $W_t$ on the center manifold.

By the near identity transformation

$\tilde{z}=z+p\left(z\mathrm{,}\overline{z}\right)$

With

$p\left(z\mathrm{,}\overline{z}\right)=\frac{g_{20}}{-i\omega }{z}^2+\frac{g_{11}}{i\omega }z\overline{z}+\frac{g_{02}}{3i\omega }{\overline{z}}^2$

We get the first Lyapunov coefficient as

$\begin{array}{c}{l}_0=g_{210}+\frac{g_{200}{g}_{110}}{i\omega }-\frac{{\left|g_{110}\right|}^2}{i\omega }+\frac{2{\left|g_{020}\right|}^2}{3i\omega }+g_{101,0}\cdot H_{11}\left(0\right)\\ +g_{101,1}\cdot H_{11}\left(-{\tau }_s\right)+g_{101,2}\cdot H_{11}\left(-{\tau }_N\right)+g_{011,0}\cdot H_{20}\left(0\right)\\ +g_{011,1}\cdot H_{20}\left(-{\tau }_s\right)+g_{011,2}\cdot H_{20}\left(-{\tau }_N\right)\end{array}$ (5.17)

With

$g_{ij1,0}=g_{ij1,W_t\left(0\right)},g_{ij1,1}=g_{ij1,W_t\left(-{\tau }_s\right)},g_{ij1,2}=g_{ij1,W_t\left(-{\tau }_N\right)}$

For $l=1,i+j+l=m$ , with $m=\mathrm{2,3,}\cdots$ . and system (5.1) is equivalent to

(5.18)

where in with no doubt, we omit $\tilde{\cdot }$ in Equation (5.18). We have $g_{21}=l_0$ , and

$g_{30}=g_{300}+\frac{2}{9\omega }i{g}_{200}^2+H_{20}\left(0\right)\cdot g_{101,0}+H_{20}\left(-{\tau }_s\right)\cdot g_{101,1}+H_{20}\left(-{\tau }_N\right)\cdot g_{101,2},$

$\begin{array}{l}{g}_{12}=g_{120}+H_{11}\left(0\right)\cdot g_{011,0}+H_{11}\left(-{\tau }_s\right)\cdot g_{011,1}+H_{11}\left(-{\tau }_N\right)\cdot g_{011,2}+H_{11}\left(0\right)\cdot g_{011,0}\\ +H_{02}\left(0\right)\cdot g_{101,0}+H_{02}\left(-{\tau }_s\right)\cdot g_{101,1}+H_{02}\left(-{\tau }_N\right)\cdot g_{101,2}+\frac{i}{\omega }{g}_{110}^2-\frac{2i}{3\omega }{g}_{020}{g}_{200},\end{array}$

$g_{03}=g_{030}-\frac{i}{\omega }{g}_{020}{g}_{110}+H_{02}\left(0\right)\cdot g_{011,0}+H_{02}\left(-{\tau }_s\right)\cdot g_{011,1}+H_{02}\left(-{\tau }_N\right)\cdot g_{011,2},$

$\begin{array}{c}{g}_{40}=g_{400}+\frac{i}{\omega }{g}_{020}{g}_{21}+H_{30}\left(0\right)\cdot g_{101,0}+H_{30}\left(-{\tau }_s\right)\cdot g_{101,1}+H_{30}\left(-{\tau }_N\right)\cdot g_{101,2}\\ +H_{20}\left(0\right)\cdot g_{201,0}+H_{20}\left(-{\tau }_s\right)\cdot g_{201,1}+H_{20}\left(-{\tau }_N\right)\cdot g_{201,2},\end{array}$

$\begin{array}{c}{g}_{31}=g_{310}-\frac{i}{\omega }{g}_{110}{g}_{21}+H_{30}\left(0\right)\cdot g_{011,0}+H_{30}\left(-{\tau }_s\right)\cdot g_{011,1}+H_{30}\left(-{\tau }_N\right)\cdot g_{011,2}\\ +H_{21}\left(0\right)\cdot g_{101,0}+H_{21}\left(-{\tau }_s\right)\cdot g_{101,1}+H_{21}\left(-{\tau }_N\right)\cdot g_{101,2}+H_{11}\left(0\right)\cdot g_{201,0}\\ +H_{11}\left(-{\tau }_s\right)\cdot g_{201,1}+H_{11}\left(-{\tau }_N\right)\cdot g_{201,2}+H_{20}\left(0\right)\cdot g_{111,0}+H_{20}\left(-{\tau }_s\right)\cdot g_{111,1}\\ +H_{20}\left(-{\tau }_N\right)\cdot g_{111,2}+\frac{4i}{3\omega }{g}_{21}{g}_{200},\end{array}$

$\begin{array}{c}{g}_{22}=g_{220}+H_{02}\left(0\right)\cdot g_{201,0}+H_{02}\left(-{\tau }_s\right)\cdot g_{201,1}+H_{02}\left(-{\tau }_N\right)\cdot g_{201,2}+H_{11}\left(0\right)\cdot g_{111,0}\\ +H_{11}\left(-{\tau }_s\right)\cdot g_{111,1}+H_{11}\left(-{\tau }_N\right)\cdot g_{111,2}+H_{12}\left(0\right)\cdot g_{101,0}+H_{12}\left(-{\tau }_s\right)\cdot g_{101,1}\\ +H_{12}\left(-{\tau }_N\right)\cdot g_{101,2}+H_{20}\left(0\right)g_{021,0}+H_{20}\left(-{\tau }_s\right)g_{021,1}+H_{20}\left(-{\tau }_N\right)g_{021,2}\\ +H_{21}\left(0\right)\cdot g_{011,0}+H_{21}\left(-{\tau }_s\right)\cdot g_{011,1}+H_{21}\left(-{\tau }_N\right)\cdot g_{011,2}-\frac{i}{3\omega }{g}_{200}{g}_{21}\\ +\frac{3i}{\omega }{g}_{21}{g}_{110},\end{array}$

$\begin{array}{c}{g}_{13}=g_{130}+H_{02}\left(0\right)\cdot g_{111,0}+H_{02}\left(-{\tau }_N\right)\cdot g_{111,2}+H_{02}\left(-{\tau }_s\right)\cdot g_{111,1}+H_{12}\left(0\right)\cdot g_{011,0}\\ +H_{12}\left(-{\tau }_s\right)\cdot g_{011,1}+H_{12}\left(-{\tau }_N\right)\cdot g_{011,2}+H_{11}\left(0\right)\cdot g_{021,0}+H_{11}\left(-{\tau }_s\right)\cdot g_{021,1}\\ +H_{11}\left(-{\tau }_N\right)\cdot g_{021,2}+H_{03}\left(0\right)\cdot g_{101,0}+H_{03}\left(-{\tau }_s\right)\cdot g_{101,1}+H_{03}\left(-{\tau }_N\right)\cdot g_{101,2}\\ -\frac{2i}{\omega }{g}_{21}{g}_{020},\end{array}$

$\begin{array}{c}{g}_{04}=g_{040}+H_{03}\left(0\right)\cdot g_{011,0}+H_{03}\left(-{\tau }_s\right)\cdot g_{011,1}+H_{03}\left(-{\tau }_N\right)\cdot g_{011,2}+H_{02}\left(0\right)\cdot g_{021,0}\\ +H_{02}\left(-{\tau }_s\right)\cdot g_{021,1}+H_{02}\left(-{\tau }_N\right)\cdot g_{021,2}\end{array}$

$\begin{array}{c}{g}_{50}=g_{500}+H_{20}{\left(0\right)}^2{g}_{\mathrm{102,00}}+H_{20}{\left(-{\tau }_s\right)}^2{g}_{\mathrm{102,11}}+H_{20}{\left(-{\tau }_N\right)}^2{g}_{\mathrm{102,22}}\\ +2H_{20}\left(0\right)H_{20}\left(-{\tau }_s\right)g_{\mathrm{102,01}}+2H_{20}\left(0\right)H_{20}\left(-{\tau }_N\right)g_{\mathrm{102,02}}\\ +2H_{20}\left(-{\tau }_N\right)H_{20}\left(-{\tau }_s\right)g_{\mathrm{102,12}}+H_{20}\left(0\right)\cdot g_{\mathrm{301,0}}+H_{20}\left(-{\tau }_s\right)\cdot g_{\mathrm{301,1}}\\ +H_{20}\left(-{\tau }_N\right)\cdot g_{\mathrm{301,2}}+H_{30}\left(0\right)\cdot g_{\mathrm{201,0}}+H_{30}\left(-{\tau }_s\right)\cdot g_{\mathrm{201,1}}+H_{30}\left(-{\tau }_N\right)\cdot g_{\mathrm{201,2}}\\ +H_{40}\left(0\right)\cdot g_{\mathrm{101,0}}+H_{40}\left(-{\tau }_s\right)\cdot g_{\mathrm{101,1}}+H_{40}\left(-{\tau }_N\right)\cdot g_{\mathrm{101,2}}+\frac{4}{3{\omega }^2}{g}_{020}{g}_{21}{g}_{200}\mathrm{,}\end{array}$

$\begin{array}{c}{g}_{41}=g_{410}+2H_{11}\left(0\right)H_{20}\left(0\right)g_{\mathrm{102,00}}+2H_{11}\left(-{\tau }_s\right)H_{20}\left(-{\tau }_s\right)g_{\mathrm{102,11}}+2H_{11}\left(0\right)H_{20}\left(-{\tau }_s\right)g_{\mathrm{102,01}}\\ +2H_{11}\left(-{\tau }_N\right)H_{20}\left(-{\tau }_N\right)g_{\mathrm{102,22}}+2H_{11}\left(-{\tau }_s\right)H_{20}\left(0\right)g_{\mathrm{102,10}}+2H_{11}\left(0\right)H_{20}\left(-{\tau }_N\right)g_{\mathrm{102,02}}\\ +2H_{11}\left(-{\tau }_N\right)H_{20}\left(0\right)g_{\mathrm{102,20}}+2H_{11}\left(-{\tau }_N\right)H_{20}\left(-{\tau }_s\right)g_{\mathrm{102,21}}+2H_{11}\left(-{\tau }_s\right)H_{20}\left(-{\tau }_N\right)g_{\mathrm{102,12}}\\ +\frac{5}{9{\omega }^2}{g}_{21}{g}_{200}^2+H_{20}^2\left(0\right)g_{\mathrm{012,00}}+H_{20}^2\left(-{\tau }_s\right)g_{\mathrm{012,11}}+H_{20}^2\left(-{\tau }_N\right)g_{\mathrm{012,22}}\\ +2H_{20}\left(0\right)H_{20}\left(-{\tau }_s\right)g_{\mathrm{012,01}}+2H_{20}\left(0\right)H_{20}\left(-{\tau }_N\right)g_{\mathrm{012,02}}+2H_{20}\left(-{\tau }_N\right)H_{20}\left(-{\tau }_s\right)g_{\mathrm{012,21}}\\ +H_{30}\left(0\right)\cdot g_{\mathrm{111,0}}+H_{30}\left(-{\tau }_s\right)\cdot g_{\mathrm{111,1}}+H_{30}\left(-{\tau }_N\right)\cdot g_{\mathrm{111,2}}+H_{20}\left(0\right)\cdot g_{\mathrm{211,0}}\\ +H_{20}\left(-{\tau }_s\right)\cdot g_{\mathrm{211,1}}+H_{20}\left(-{\tau }_N\right)\cdot g_{\mathrm{211,2}}+H_{31}\left(0\right)\cdot g_{\mathrm{101,0}}+H_{31}\left(-{\tau }_s\right)\cdot g_{\mathrm{101,1}}\\ +H_{31}\left(-{\tau }_N\right)\cdot g_{\mathrm{101,2}}+H_{11}\left(0\right)\cdot g_{\mathrm{301,0}}+H_{11}\left(-{\tau }_s\right)\cdot g_{\mathrm{301,1}}+H_{11}\left(-{\tau }_N\right)\cdot g_{\mathrm{301,2}}\\ +H_{21}\left(0\right)\cdot g_{\mathrm{201,0}}+H_{21}\left(-{\tau }_s\right)\cdot g_{\mathrm{201,1}}+H_{21}\left(-{\tau }_N\right)\cdot g_{\mathrm{201,2}}+H_{40}\left(0\right)\cdot g_{\mathrm{011,0}}\\ +H_{40}\left(-{\tau }_s\right)\cdot g_{\mathrm{011,1}}+H_{40}\left(-{\tau }_N\right)\cdot g_{\mathrm{011,2}}+\frac{3}{{\omega }^2}{g}_{21}{g}_{110}{g}_{020}-\frac{4}{3{\omega }^2}{g}_{21}{g}_{110}{g}_{200}\mathrm{,}\end{array}$

$\begin{array}{c}{g}_{14}=g_{140}+H_{03}\left(0\right)\cdot g_{\mathrm{111,0}}+H_{03}\left(-{\tau }_s\right)\cdot g_{\mathrm{111,1}}+H_{03}\left(-{\tau }_N\right)\cdot g_{\mathrm{111,2}}+H_{02}\left(0\right)\cdot g_{\mathrm{121,0}}\\ +H_{02}\left(-{\tau }_s\right)\cdot g_{\mathrm{121,1}}+H_{02}\left(-{\tau }_N\right)\cdot g_{\mathrm{121,2}}+H_{04}\left(0\right)\cdot g_{\mathrm{101,0}}+H_{04}\left(-{\tau }_s\right)\cdot g_{\mathrm{101,1}}\\ +H_{04}\left(-{\tau }_N\right)\cdot g_{\mathrm{101,2}}+H_{13}\left(0\right)\cdot g_{\mathrm{011,0}}+H_{13}\left(-{\tau }_s\right)\cdot g_{\mathrm{011,1}}+H_{13}\left(-{\tau }_N\right)\cdot g_{\mathrm{011,2}}\\ +H_{11}\left(0\right)\cdot g_{\mathrm{031,0}}+H_{11}\left(-{\tau }_s\right)\cdot g_{\mathrm{031,1}}+H_{11}\left(-{\tau }_N\right)\cdot g_{\mathrm{031,2}}+H_{12}\left(0\right)\cdot g_{\mathrm{021,0}}\\ +H_{12}\left(-{\tau }_s\right)\cdot g_{\mathrm{021,1}}+H_{12}\left(-{\tau }_N\right)\cdot g_{\mathrm{021,2}}+\frac{2}{3{\omega }^2}{g}_{020}{g}_{21}{g}_{200}-\frac{4}{{\omega }^2}{g}_{21}{g}_{110}{g}_{020}\\ +H_{02}^2\left(0\right)g_{\mathrm{102,00}}+H_{02}^2\left(-{\tau }_s\right)g_{\mathrm{102,11}}+H_{02}^2\left(-{\tau }_N\right)g_{\mathrm{102,22}}+2H_{02}\left(0\right)H_{02}\left(-{\tau }_s\right)g_{\mathrm{102,01}}\\ +2H_{02}\left(0\right)H_{02}\left(-{\tau }_N\right)g_{\mathrm{102,02}}+2H_{02}\left(-{\tau }_s\right)H_{02}\left(-{\tau }_N\right)g_{\mathrm{102,12}}+2H_{02}\left(0\right)H_{11}\left(0\right)g_{\mathrm{012,00}}\\ +2H_{02}\left(-{\tau }_s\right)H_{11}\left(-{\tau }_s\right)g_{\mathrm{012,11}}+2H_{02}\left(-{\tau }_N\right)H_{11}\left(-{\tau }_N\right)g_{\mathrm{012,22}}+2H_{02}\left(0\right)H_{11}\left(-{\tau }_s\right)g_{\mathrm{012,01}}\\ +2H_{02}\left(-{\tau }_s\right)H_{11}\left(0\right)g_{\mathrm{012,10}}+2H_{02}\left(0\right)H_{11}\left(-{\tau }_N\right)g_{\mathrm{012,02}}+2H_{02}\left(-{\tau }_N\right)H_{11}\left(0\right)g_{\mathrm{012,20}}\\ +2H_{02}\left(-{\tau }_N\right)H_{11}\left(-{\tau }_s\right)g_{\mathrm{012,21}}+2H_{02}\left(-{\tau }_s\right)H_{11}\left(-{\tau }_N\right)g_{\mathrm{012,12}}\end{array}$

$\begin{array}{c}{g}_{05}=g_{050}+\frac{1}{{\omega }^2}{g}_{020}^2{g}_{21}+H_{04}\left(0\right)\cdot g_{\mathrm{011,0}}+H_{04}\left(-{\tau }_s\right)\cdot g_{\mathrm{011,1}}+H_{04}\left(-{\tau }_N\right)\cdot g_{\mathrm{011,2}}\\ +H_{02}\left(0\right)\cdot g_{\mathrm{031,0}}+H_{02}\left(-{\tau }_s\right)\cdot g_{\mathrm{031,1}}+H_{02}\left(-{\tau }_N\right)\cdot g_{\mathrm{031,2}}+H_{03}\left(0\right)\cdot g_{\mathrm{021,0}}\\ +H_{03}\left(-{\tau }_s\right)\cdot g_{\mathrm{021,1}}+H_{03}\left(-{\tau }_N\right)\cdot g_{\mathrm{021,2}}+H_{02}^2\left(0\right)g_{\mathrm{012,00}}+H_{02}^2\left(-{\tau }_s\right)g_{\mathrm{012,11}}\\ +H_{02}^2\left(-{\tau }_N\right)g_{\mathrm{012,22}}+2H_{02}\left(0\right)H_{02}\left(-{\tau }_s\right)g_{\mathrm{012,01}}+2H_{02}\left(0\right)H_{02}\left(-{\tau }_N\right)g_{\mathrm{012,02}}\\ +2H_{02}\left(-{\tau }_N\right)H_{02}\left(-{\tau }_s\right)g_{\mathrm{012,21}}\end{array}$

Equation (5.8) is an ODE with its coefficients have unknown coefficient $H_{ij}\mathrm{,}i+j\ge 3$ , we solve it by the same differential methods as referred in Equation (5.11) - Equation (5.16), that is

$\begin{array}{l}\mathcal{A}{H}_{30}\left(\theta \right)=3i\omega H_{30}\left(\theta \right)+{\Phi }_1\left(\theta \right)g_{30}+{\Phi }_2\left(\theta \right){\overline{g}}_{03},\hfill \\ \mathcal{A}{H}_{21}\left(\theta \right)=i\omega H_{21}\left(\theta \right)+{\Phi }_1\left(\theta \right)g_{21}+{\Phi }_2\left(\theta \right){\overline{g}}_{12},\hfill \\ \mathcal{A}{H}_{12}\left(\theta \right)=-i\omega H_{12}\left(\theta \right)+{\Phi }_1\left(\theta \right)g_{12}+{\Phi }_2\left(\theta \right){\overline{g}}_{21},\hfill \\ \mathcal{A}{H}_{03}\left(\theta \right)=-3i\omega H_{03}\left(\theta \right)+{\Phi }_1\left(\theta \right)g_{03}+{\Phi }_2\left(\theta \right){\overline{g}}_{30},\hfill \\ \mathcal{A}{H}_{40}\left(\theta \right)=4i\omega H_{30}\left(\theta \right)+{\Phi }_1\left(\theta \right)g_{40}+{\Phi }_2\left(\theta \right){\overline{g}}_{04},\hfill \\ \mathcal{A}{H}_{31}\left(\theta \right)=2i\omega H_{31}\left(\theta \right)+{\Phi }_1\left(\theta \right)g_{31}+{\Phi }_2\left(\theta \right){\overline{g}}_{13},\hfill \\ \mathcal{A}{H}_{22}\left(\theta \right)={\Phi }_1\left(\theta \right)g_{22}+{\Phi }_2\left(\theta \right){\overline{g}}_{22},\hfill \\ \mathcal{A}{H}_{13}\left(\theta \right)=-2i\omega H_{13}\left(\theta \right)+{\Phi }_1\left(\theta \right)g_{13}+{\Phi }_2\left(\theta \right){\overline{g}}_{31},\hfill \\ \mathcal{A}{H}_{04}\left(\theta \right)=-4i\omega H_{04}\left(\theta \right)+{\Phi }_1\left(\theta \right)g_{04}+{\Phi }_2\left(\theta \right){\overline{g}}_{40}\hfill \end{array}$ (5.19)

With initial value condition determined by the case $\theta =0$ in Equation (5.15).

Set

$p\left(z,\overline{z}\right)=\frac{1}{2\omega }i{g}_{30}{z}^3-\frac{i}{2\omega }{g}_{12}z{\overline{z}}^2-\frac{i}{4\omega }{g}_{03}{\overline{z}}^3$ (5.20)

And do near identity transformation

$\tilde{z}=z+p\left(z,\overline{z}\right)$ (5.21)

We get the second Lyapunov exponent

$l_2=g_{32}-\frac{3}{4\omega }i{\overline{g}}_{03}{g}_{03}+\frac{i}{\omega }{g}_{12}{g}_{30}-\frac{i}{\omega }{\overline{g}}_{12}{g}_{12}$ (5.22)

Therefore, the normal form of the reduced system is written as

$z^{\prime }\left(t\right)=i\omega z+S_{\epsilon }z+l_1{z}^2\overline{z}+l_2{z}^3{\overline{z}}^2$ (5.23)

Under the polar axis $\left(\rho \mathrm{,}\theta \right)$ , the normal form is written as

${\rho }^{\prime }\left(t\right)=\rho \left(\Re \left\{S_{\epsilon }\right\}+\Re \left\{l_1\right\}{\rho }^2+\Re \left\{l_2\right\}{\rho }^4\right)$ (5.24)

The limit cycle bifurcation is computed as

$\Re {\left\{l_1\right\}}^2-\epsilon \Re \left\{S_{\epsilon }\right\}\Re \left\{l_2\right\}=0$ (5.25)

The first Lyapunov exponent of GH point equals zero, and it is concluded that the nearby dynamics of GH is classified by at least two parameters perturbation.

Fixed ${\gamma }_0=0.283470299684549$ , we draw Hopf bifurcation line on parameters $\left(k_{\delta }\mathrm{,}{\tau }_N\right)$ plane, as shown in Figure 7(a). GH bifurcation is found at the bottom of the valley located in Hopf line, and the eigen roots of GH is pictured in Figure 1(a). GH is the separation point between supercritical Hopf points and subcritical Hopf points which locates on Hopf line on $\left(k_{\delta }\mathrm{,}{\tau }_N\right)$ plane. The computation of the normal form are exemplified by choosing Hopf points respectively. For Hopf point 1 with $\left({\tau }_N,k_{\delta }\right)=\left(8.187,0.0134316\right)$ , which is subcritical Hopf point with $\Re \left\{l_1\right\}>0$ ; For Hopf point 2 with $\left({\tau }_N,k_{\delta }\right)=\left(8.202,0.0134305\right)$ , which is subcritical Hopf point with $\Re \left\{l_1\right\}>0$ ; For Hopf point 3 with $\left({\tau }_N,k_{\delta }\right)=\left(8.242,0.0134307\right)$ , which is supercritical Hopf point with $\Re \left\{l_1\right\}>0$ ; For Hopf point 4 with $\left({\tau }_N,k_{\delta }\right)=\left(8.256,0.0134318\right)$ , which is supercritical Hopf point with $\Re \left\{l_1\right\}>0$ . The bifurcating periodical solutions from Hopf point 1, 2, 3, 4 are shown in Figure 7(b).

Specially, for GH point with $\left({\tau }_N,k_{\delta }\right)=\left(8.22,0.01343\right)$ , which is both supercritical and subcritical Hopf bifurcation with singularity of codimension 2. The near dynamics of phase portraits of GH point is topological classified, as shown in Figure 8(a). The continuation results of periodical solutions arising from GH point are pictured as varying ${\tau }_N$ . The continuation of bifurcating periodical solutions is completed by DDE-Biftool software, as shown in Figure 8(b).

(a) (b)

Figure 7. The generalized Hopf bifurcation point in parameters $\left(k_{\delta }\mathrm{,}{\tau }_N\right)$ plane. (a) Hopf line with the GH point, the imaginary roots is pictured in the middle view; (b) The bifurcating solutions from subcritical Hopf points 1, 2 and supercritical Hopf points 3, 4.

(a) (b)

Figure 8. GH point is both the supercritical Hopf point and subcritical Hopf point. (a) The phase portraits of near dynamics of GH point geometrically; (b) The maximal magnitude and minimal magnitude of bifurcating solutions continued when varying ${\tau }_N$ , which bifurcates from GH point.

6. Conclusion

With delay dependent coefficients in DDEs of neutrophil dynamics, the neutrophil system was numerical simulated with the continuation method proposed by DDE-Biftool software. The generalized Hopf point (GH) was exploited to separate the supercritical Hopf points from the subcritical Hopf points. The observed GH point was located in the bottom of the valley of Hopf line, which manifested the produced phase portraits of the classification of the near dynamics. The usual stable and unstable limit cycle collide to disappear beyond the equilibrium solution preserved its local stable property. The preserved program by applying DDE-Biftool also displayed the lineup solutions of bifurcating period limit cycles. Near GH point, the open lineup of periodical solutions was observed, which often bifurcates from one Hopf point and then gets to another Hopf point to form a bubble of continuation of limit cycles. However, once the free parameter crossing the critical threshold, the closed lineup of periodical solutions was exhibited which is a cycle again with no equilibrium solutions bifurcation any more. Nevertheless, the lineup of bifurcating limit cycles connects GH point, which features GH bifurcation with both supercritical Hopf and subcritical Hopf derived.

Availability of Data and Materials

All data generated or analyzed during this study are included in this published article.

Conflicts of Interest

The authors declare no conflicts of interest regarding the publication of this paper.

References

[1]	Kaushansky, K. (2016) Hematopoietic Stem Cells, Progenitors, and Cytokines. In: Kaushansky, K., et al., Eds., Williams Hematology, 9th Edition, McGraw-Hill.
[2]	Brunetti, M., Mackey, M.C. and Craig, M. (2021) Understanding Normal and Pathological Hematopoietic Stem Cell Biology Using Mathematical Modelling. Current Stem Cell Reports, 7, 109-120. https://doi.org/10.1007/s40778-021-00191-9
[3]	Mackey, M.C. (2001) Cell Kinetic Status of Haematopoietic Stem Cells. Cell Proliferation, 34, 71-83. https://doi.org/10.1046/j.1365-2184.2001.00195.x
[4]	Colijn, C., Foley, C. and Mackey, M.C. (2007) G-CSF Treatment of Canine Cyclical Neutropenia: A Comprehensive Mathematical Model. Experimental Hematology, 35, 898-907. https://doi.org/10.1016/j.exphem.2007.02.015
[5]	Foley, C., Bernard, S. and Mackey, M.C. (2006) Cost-Effective G-CSF Therapy Strategies for Cyclical Neutropenia: Mathematical Modelling Based Hypotheses. Journal of Theoretical Biology, 238, 754-763. https://doi.org/10.1016/j.jtbi.2005.06.021
[6]	Hearn, T., Haurie, C. and Mackey, M.C. (1998) Cyclical Neutropenia and the Peripheral Control of White Blood Cell Production. Journal of Theoretical Biology, 192, 167-181. https://doi.org/10.1006/jtbi.1997.0589
[7]	Haurie, C., Dale, D., Rudnicki, R., et al. (2000) Modeling Complex Neutrophil Dynamics in the Grey Collie. Journal of Theoretical Biology, 192, 167-181.
[8]	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. https://doi.org/10.1182/blood.v92.8.2629.420a35_2629_2640
[9]	De Souza, D.C. and Humphries, A.R. (2019) Dynamics of a Mathematical Hematopoietic Stem-Cell Population Model. SIAM Journal on Applied Dynamical Systems, 18, 808-852. https://doi.org/10.1137/18m1165086
[10]	Crauste, F. (2006) Global Asymptotic Stability and Hopf Bifurcation for a Blood Cell Production Model. Mathematical Biosciences and Engineering, 3, 325-346.
[11]	Pujo-Menjouet, L., Bernard, S. and Mackey, M.C. (2005) Long Period Oscillations in a G0 Model of Hematopoietic Stem Cells. SIAM Journal on Applied Dynamical Systems, 4, 312-332. https://doi.org/10.1137/030600473
[12]	Fortin, P. and Mackey, M.C. (1999) Periodic Chronic Myelogenous Leukaemia: 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
[13]	Fowler, A.C. and Mackey, M.C. (2002) Relaxation Oscillations in a Class of Delay Differential Equations. SIAM Journal on Applied Mathematics, 63, 299-323. https://doi.org/10.1137/s0036139901393512
[14]	Colijn, C. and Mackey, M.C. (2007) Bifurcation and Bistability in a Model of Hematopoietic Regulation. SIAM Journal on Applied Dynamical Systems, 6, 378-394. https://doi.org/10.1137/050640072
[15]	Colijn, C. and Mackey, M.C. (2007) Bifurcation and Bistability in a Model of Hematopoietic Regulation. SIAM Journal on Applied Dynamical Systems, 6, 378-394. https://doi.org/10.1137/050640072
[16]	Bernard, S., Bélair, J. and Mackey, M.C. (2003) Oscillations in Cyclical Neutropenia: New Evidence Based on Mathematical Modeling. Journal of Theoretical Biology, 223, 283-298. https://doi.org/10.1016/s0022-5193(03)00090-0
[17]	Adimya, M., Craustea, F. and Ruanb, S. (2006) Modelling Hematopoiesis Mediated by Growth Factors with Applications to Periodic Hematological Diseases. Bulletin of Mathematical Biology, 68, 2321-2351. https://doi.org/10.1007/s11538-006-9121-9
[18]	Talibi Alaoui, H. and Yafia, R. (2007) Stability and Hopf Bifurcation in an Approachable Haematopoietic Stem Cells Model. Mathematical Biosciences, 206, 176-184. https://doi.org/10.1016/j.mbs.2006.03.004
[19]	Andersen, L.K. and Mackey, M.C. (2001) Resonance in Periodic Chemotherapy: A Case Study of Acute Myelogenous Leukemia. Journal of Theoretical Biology, 209, 113-130. https://doi.org/10.1006/jtbi.2000.2255
[20]	Ma, S. (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
[21]	Hale, J.K. (1977) Theory of Functional Differential Equations. Springer.
[22]	Kolmanovskii, V. and Myshkis, A. (1999) Introduction to the Theory and Applications of Functional Differential Equations. Kluwer Academic Publishers.
[23]	Vainstein, V., Ginosar, Y., Shoham, M., Ranmar, D.O., Ianovski, A. and Agur, Z. (2005) The Complex Effect of Granulocyte Colony-Stimulating Factor on Human Granulopoiesis Analyzed by a New Physiologically-Based Mathematical Model. Journal of Theoretical Biology, 234, 311-327. https://doi.org/10.1016/j.jtbi.2004.11.026
[24]	Kuznetsov, Y.A. (2004) Elements of Applied Bifurcation Theory. 3rd Edition, Springer.
[25]	Hassard, B., Kazarinoff, N. and Wan, Y. (1981) Theory and Applications of Hopf Bifurcation. Cambridge Univ. Press.
[26]	Faria, T. and Magalhaes, L.T. (1995) Normal Forms for Retarded Functional Differential Equations with Parameters and Applications to Hopf Bifurcation. Journal of Differential Equations, 122, 181-200. https://doi.org/10.1006/jdeq.1995.1144
[27]	Ma, S.Q., Lu, Q. and Hogan, S. (2007) The Normal Form for Double Hopf Bifurcation for Stuart-Landau System with Nonlinear Delay Feedback and Delay-Dependent Parameters. Advances in Complex Systems.
[28]	Engelborghs, K., Luzyanina, T. and Roose, D. (2002) Numerical Bifurcation Analysis of Delay Differential Equations Using Ddebiftool. ACM Transactions on Mathematical Software, 28, 1-21. https://doi.org/10.1145/513001.513002
[29]	Engelborghs, K., Luzyanina, T. and Samaey, G. (2001) DDE-BIFTOOL v. 2.00: A MATLAB Package for Bifurcation Analysis of Delay Differential Equations. Tech. Rep., Department of Computer Science, K. U. Leuven.