Stability and Bifurcation Analysis of a Type of Hematopoietic Stem Cell Model

Suqi Ma

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

DOI: 10.4236/ijmnta.2021.101002                

Abstract

The observed dynamical property illustrates that state feedback control may stabilize invariant attractor to stable state in a simple version of hematopoietic stem cell model. The stability character of the positive steady state is analyzed by the computation of the rightmost characteristic roots in complex plane. Hopf bifurcation points are tracked as the roots curve crossing imaginary axis from the left half plane to the right half plane continuously. The bifurcation direction and stability of the bifurcating periodical solution are discussed by norm form computation combined with the center manifold theory. Furthermore, the numerical simulation verifies that instead of chaos, system is stabilized to period-1, 2, 3, 4 and period-7 periodical solutions in some delay windows, and the continuous of periodical solutions is also numerical simulated with varying free parameters continuously.

Keywords

State Feedback Control, Time Delay, Hopf Bifurcation

1. Introduction

Delay differential equations have been broadly focused in every fields of scientific investigation work since time delay is a natural factor in the reality life [1] [2] [3]. For example, cell’s maturation time in biology, the reflection time for driving in road traffic, etc. In engineering control fields, people try to design delay control apparatus to let the original system stabilized to stable periodical orbits or produce new bifurcation behavior [4] [5]. The more attention work is the state feedback control in engineering control studying, the periodical oscillation behavior is induced via either time delay state feedback or difference feedback control [6] [7].

As is well known, the new periodical oscillation phenomena occur at Hopf point since system lost its stability. As for the eigenvalue problem of DDEs, Hopf bifurcation occurs as the imaginary roots cross the imaginary axis from the left half plane to the right half plane. With single time delay effects, people try to analyze Hopf bifurcation of linear DDEs by computing imaginary roots by the algebra method [8] [9]. Wang and Hu have shown the high analyzing technique by applying Sturm criterion with Maple language computation [10] [11]. Kuang computes the eigenvalue problem of a type DDEs by locating imaginary roots of Hopf bifurcation with consideration of delay-dependant physical parameter [12]. Together with Cookie’s work [13], Kuang’s work invokes people’s big interest in studying biological model with delay-dependent nonlinear coefficient. In fact, it is ubiquitous to introduce delay-dependent physical parameter in biological models since nonlinear birth rate to specify the loss rate in specie’s growth stage [14] [15] [16].

The stability switching always brings forth the periodical oscillation phenomena and complex dynamical behavior. A standard and much studied work of DDE is the Mackey-Glass equation which is proposed to model the production of white blood cells and given by

$\stackrel{\dot{}}{x}=\frac{ax\left(t-\tau \right)}{1+x{\left(t-\tau \right)}^c}-bx$ (1)

For which an invariant attractor is observed with parameters such as $a=0.2,b=0.1,c=10$ and $\tau =16.5$. With state feedback control of time delay, system (1) is described as

$\stackrel{\dot{}}{x}=\frac{ax\left(t-\tau \right)}{1+x{\left(t-\tau \right)}^c}-bx+K\left(x\left(t-\tau \right)-x\left(t\right)\right)$ (2)

Herein, x denotes the concentration of white blood cells, a, c are Hill coefficients, and b represents death rate.

Since in one respect, with state feedback control to perturb system’s dynamics to produce new bifurcation behavior; and in another respect, mathematically, to reflect the function of perturbation of state difference biologically since migration phenomena in model analysis.

The different dynamical character of hematopoietic stem cell model is displayed as shown in Figure 1(a) and Figure 1(b). The observed invariant attractor is becoming asymptotically stable steady state as exerting feedback control with strength $K=-0.2$. As varying the free parameter, the steady state may loss its stability into instability state to induce the oscillation phenomena arising in system. To further analyze the bifurcation mechanism of the stability of the equilibrium solution, the characteristic root with zero real part appearing in its characteristic equation is analyzed analytically. Hence, we endeavour to calculate the rightmost characteristic root wisdomly from geometrical point view which to decide the asymptotically stability of the equilibrium solution.

With the consideration of delay-dependent physical parameter in DDEs, we solve the eigenvalue problem of Hopf bifurcation by a new method which is

(a)(b)

Figure 1. The dynamics of hematopoietic stem cell model (1) and (1) with chosen time delay $\tau =16.5$, (a) The observed invariant attractor of system (1); (b) System is stabilized to the steady state under delay feedback control with strength $K=-0.2$.

illustrated as geometrical criterion in paper [17] [18]. Numerically, DDE-Biftool can compute the rightmost roots of the linear characteristic equation. We devote to compute the rightmost characteristic roots of the linear DDE of system (1) and (1) to give out the asymptotically stable condition, and make earnest endeavors to draw curve of the rightmost characteristic roots continuously as varying parameter K. Hopf bifurcation of Equation (1) and (1) are also analyzed to show the bifurcation of periodical oscillation solutions with small amplitude. Furthermore, based on the fundamental theory of functional differential equations [19] [20] [21], the bifurcating periodical solutions is computed via norm form analytical technique combined center manifold theory. The numerical simulation displays that the observed perioid-1, 2, 3, 4 and perioid-7 solution and chaos solution in different delay window underlying state feedback control.

The whole paper is organized as the listed. In Section 2, the distribution of characteristic roots in a band is calculated via geometrical analyze technique and the rightmost characteristic roots determines the stability of the positive steady state, and Hopf bifurcation is analyzed in Section 3. Based on the fundamental theory of DDEs, the dimension reduction system of system (1) and (1) is computed and analyzed combined with the center manifold theory. The bifurcation direction and the stability of the bifurcating periodical solution are derived via formal norm analytical technique. The numerical simulation has shown the continuous of oscillation solutions as varying free parameters.

2. The Computation of the Rightmost Characteristic Root

The characteristic equation of the linear DDEs of Equation (1) and Equation (1) can be written as

$\Delta \left(\lambda \mathrm{,}\sigma \mathrm{,}\tau \right)=P\left(\lambda \mathrm{,}\sigma \right)+Q\left(\lambda \mathrm{,}\sigma \right){\text{e}}^{-\lambda \tau }$ (3)

For the fixed value of $\sigma ={\sigma }^{*},\tau ={\tau }^{*}$, suppose $\lambda =\alpha +i\beta$, substitute it into the characteristic Equation (2) then seperate the imaginary part from the real part to get

$\begin{array}{l}{P}_R\left(\alpha ,\beta ,{\sigma }^{*}\right)+Q_R\left(\alpha ,\beta ,{\sigma }^{*}\right){\text{e}}^{-\alpha {\tau }^{*}}\cos \left(\beta {\tau }^{*}\right)+Q_I\left(\alpha ,\beta ,{\sigma }^{*}\right){\text{e}}^{-\alpha {\tau }^{*}}\sin \left(\beta {\tau }^{*}\right)=0,\\ P_I\left(\alpha ,\beta ,{\sigma }^{*}\right)+Q_I\left(\alpha ,\beta ,{\sigma }^{*}\right){\text{e}}^{-\alpha {\tau }^{*}}\cos \left(\beta {\tau }^{*}\right)-Q_R\left(\alpha ,\beta ,{\sigma }^{*}\right){\text{e}}^{-\alpha {\tau }^{*}}\sin \left(\beta {\tau }^{*}\right)=0\end{array}$ (4)

By Equation (2) to get

$\begin{array}{l}{\text{e}}^{-\alpha {\tau }^{*}}\cos \left(\beta {\tau }^{*}\right)=-\frac{P_I{Q}_I+P_R{Q}_R}{Q_I^2+Q_R^2},{\text{e}}^{-\alpha {\tau }^{*}}\sin \left(\beta {\tau }^{*}\right)=\frac{P_I{Q}_R-P_R{Q}_I}{Q_I^2+Q_R^2}\hfill \end{array}$ (5)

Set $S=\beta {\tau }^{*}$, the equivalent equation of Equations (2) is obtained as

$\begin{array}{l}{\text{e}}^{-2\alpha {\tau }^{*}}={\left(\frac{P_I{Q}_I+P_R{Q}_R}{Q_I^2+Q_R^2}\right)}^2+{\left(\frac{P_I{Q}_R-P_R{Q}_I}{Q_I^2+Q_R^2}\right)}^2,\\ G\left(\alpha ,S\right)=-\frac{P_I{Q}_R-P_R{Q}_I}{P_I{Q}_I+P_R{Q}_R}-\tan \left(S\right)\end{array}$ (6)

Hence we derive the following lemma 2.1,

Lemma 2.1: Suppose ${G^{\prime }}_S\left(\alpha \mathrm{,}S\right)\ne 0$, then there’s implicit function $\alpha =\alpha \left(S\right)$ determined by $G\left(\alpha \mathrm{,}S\right)\equiv 0$.

Define the line

$L:Y_1={\text{e}}^{-2\alpha \left(S\right){\tau }^{*}}$ (7)

and the line

$R:Y_2={\left(\frac{P_I{Q}_I+P_R{Q}_R}{Q_I^2+Q_R^2}\right)}^2+{\left(\frac{P_I{Q}_R-P_R{Q}_I}{Q_I^2+Q_R^2}\right)}^2$ (8)

then the intersection point $S^{\mathrm{<mo>\; *\; </mo>}}$ determines the corresponding characteristic roots $\lambda =\alpha \left(S^{\mathrm{<mo>\; *\; </mo>}}\right)+i\beta \left(S^{\mathrm{<mo>\; *\; </mo>}}\right)$.

Setting $x\left(t-\tau \right)=x\left(t\right)$, the positive steady state of Equation (1) and Equation

(1) is calculated as $x^{*}=\exp \left(\frac{1}{C}\ln \left(\frac{a-b}{b}\right)\right)$. Doing the axis transformation $x=x-x^{*}$, $x\left(t-\tau \right)=x\left(t-\tau \right)-x^{*}$, one obtains the linear DDEs

$x^{\prime }\left(t\right)=Ax\left(t\right)+Bx\left(t-\tau \right)$ (9)

with

$A=-b+K,B=-\frac{\left(-1+\left(c-1\right){\left({\left(\left(a-b\right)/b\right)}^{\frac{1}{c}}\right)}^c\right)a}{{\left(1+{\left({\left(\left(a-b\right)/b\right)}^{\frac{1}{c}}\right)}^c\right)}^2}-K$ (10)

Hence one has

$P\left(\lambda ,b\right)=A-\lambda ,Q\left(\lambda ,b\right)=B$ (11)

Equation (2) is rewritten as

$\begin{array}{l}{\text{e}}^{-2\alpha \tau }={\left(\frac{A-\alpha }{B}\right)}^2+{\left(\frac{-S}{\tau B}\right)}^2,\\ G\left(\alpha ,S\right)=\frac{S}{\tau \left(A-\alpha \right)}-\tan \left(S\right)\end{array}$ (12)

By the second equation in Equations (2) one gets

$\alpha =\frac{\tan \left(S\right)\tau A-S}{\tan \left(S\right)\tau }$ (13)

Therefore by Lemma 2.1, one defines $L:Y_1=\exp \left(-2\frac{\tan \left(S\right)\tau A-S}{\tan \left(S\right)}\right)$ and $R:Y_2=\frac{S^2}{{\tau }^2{B}^2\sin {\left(S\right)}^2}$. We also draw the following conclusion:

Lemma 2.2: Suppose the intersection point of line L and line R is $\left(S^{\mathrm{<mo>\; *\; </mo>}}\mathrm{,}{A}^{\mathrm{<mo>\; *\; </mo>}}\mathrm{,}{B}^{\mathrm{<mo>\; *\; </mo>}}\right)$, then the only positive equilibrium solution is asymptotically stable if $Y_1=Y_1\left(S^{*},A^{*},B^{*}\right)<1$ and versus if $Y_1=Y_1\left(S^{*},A^{*},B^{*}\right)>1$. In addition, Hopf bifurcation is induced at $\left(S^{\mathrm{<mo>\; *\; </mo>}}\mathrm{,}{A}^{\mathrm{<mo>\; *\; </mo>}}\mathrm{,}{B}^{\mathrm{<mo>\; *\; </mo>}}\right)$ if $Y_1=Y_1\left(S^{*},a^{*},B^{*}\right)=1$.

The switching of stability of the only equilibrium solution is shown in Figure 2(a) and Figure 2(b). The green line describes the continuous of asymptotically stable equilibrium solution with varying b and the red line represents the unstable equilibrium. Chosen different feedback strength $K<0$, it is observed that the equilibrium solution is robust stable with smaller value of K. Figure 2(b) also shows the stability property of equilibrium solution as varying time delay $\tau$.

(a)(b)

Figure 2. The location of the imaginary roots in complex plane, with chosen time delay $\tau =16.5$. (a) Without output state feedback control since $K=0$ ; (b) The appearance of output state feedback control since $K=-0.2$.

The solution is asymptotically stable for small time delay then change to be unstable as increasing strength K. We draw the stable regime with green as shown in Figure 3 and red regime represents the unstable regime. It can be seen Hopf bifurcation occurs as stability property of the equilibrium solution is changed with varying parameter a and b and the amplifying picture denotes the stability

(a)(b)

Figure 3. The continuous of equilibrium solution with different strength of $K=0,K=-0.02$ and $K=-0.08$. (a) The stability property of the equilibrium solution as varying b; (b) The stability property of the equilibrium solution as varying time delay.

property of equilibrium solution with the specified value of $K=-0.1$ as varying b continuously. Hopf bifurcation line is also drawn with blue line and the discussion of Hopf bifurcation is given in Section 3.

3. Hopf Bifurcation

As shown in Figure 2 and Figure 3, Hopf bifurcation occurs respectively with varying parameter b and feedback strength K, and time delay $\tau$. The formula of Hopf bifurcation is deduced by setting $\alpha =0$. Set $\lambda =i\omega$, substitute it into Equation (2) then separate the real part from the imaginary part to get

$\begin{array}{l}{P}_R\left(\omega \mathrm{,}\sigma \right)+Q_R\left(\omega \mathrm{,}\sigma \right)\cos \left(\omega \tau \right)=\mathrm{0,}\\ P_I\left(\omega \mathrm{,}\sigma \right)-Q_R\left(\omega \mathrm{,}\sigma \right)\sin \left(\omega \tau \right)=0\end{array}$ (14)

Hence Hopf bifurcation occurs if and only if

$P_R^2+P_I^2-Q_R^2=0$ (15)

Set $S=\omega \tau$, by Equation (3), one gets

$\begin{array}{l}{\left(\frac{A}{B}\right)}^2+{\left(\frac{-S}{\tau B}\right)}^2-1=\mathrm{0,}\\ \frac{S}{\tau \left(A\right)}-\tan \left(S\right)=0\end{array}$ (16)

and further to obtain Hopf line in $\left(b\mathrm{,}K\right)$ -plane. Hopf line is shown in Figure 4 with blue color. By Equation (16), one obtains that

$\begin{array}{l}\tau =\frac{\left(1\pm \sqrt{\tan {\left(S\right)}^2+1}\right)S}{\tan \left(S\right)\left(B_1-b\right)}\\ K=\frac{\tan \left(S\right)\tau b+S}{\tan \left(S\right)\tau }\end{array}$ (17)

Following we will compute the transversal condition for Hopf bifurcation. The characteristic Equation (2) can be rewritten as

Figure 4. The equilibrium plane as varying parameter b and feedback strength K. The green regime denotes the asymptotical stability of the equilibrium solution whilst the red regime is unstable one.

$\Delta \left(\lambda \mathrm{,}K\mathrm{,}\tau \right)=P\left(\lambda \mathrm{,}K\right)+Q\left(\lambda \mathrm{,}K\right){\text{e}}^{-\lambda \tau }$ (18)

Differentiate Equation (3) with respect to time delay $\tau$, one has

${P^{\prime }}_{\lambda }\frac{\text{d}\lambda }{\text{d}\tau }+\left({Q^{\prime }}_{\lambda }-\tau Q\right){\text{e}}^{-\lambda \tau }\frac{\text{d}\lambda }{\text{d}\tau }-\lambda Q{\text{e}}^{-\lambda \tau }=0$ (19)

That is,

$\delta =signRe\frac{\text{d}\lambda }{\text{d}\tau }=signRe\frac{1}{\frac{\text{d}\lambda }{\text{d}\tau }}=signRe\frac{{P^{\prime }}_{\lambda }+\left({Q^{\prime }}_{\lambda }-\tau Q\right){\text{e}}^{-\lambda \tau }}{\lambda Q{\text{e}}^{-\lambda \tau }}$ (20)

On another respect, differentiate Equation (3) with respect to S to get

$\begin{array}{l}{P^{\prime }}_{\lambda }i{\omega }^{\prime }\left(S\right)+{P^{\prime }}_K{K}^{\prime }\left(S\right)+\left({Q^{\prime }}_{\lambda }i{\omega }^{\prime }\left(S\right)+{Q^{\prime }}_K{K}^{\prime }\left(S\right)\right){\text{e}}^{-\lambda \tau }\\ +Q{\text{e}}^{-\lambda \tau }\left(-i{\omega }^{\prime }\left(S\right)\tau -i\omega {\tau }^{\prime }\left(S\right)\right)=0\end{array}$ (21)

By Equation (3) and Equation (3) one can compute that

$\begin{array}{c}\delta \left(S\right)=signRe\frac{1}{-i{\omega }^{\prime }\left(S\right)}\frac{{{\Delta }^{\prime }}_K{K}^{\prime }\left(S\right)+Q{\text{e}}^{-iS}\left(-i\omega \right){\tau }^{\prime }\left(S\right)}{i\omega Q{\text{e}}^{-iS}}\\ =signRe\frac{{{\Delta }^{\prime }}_K{K}^{\prime }\left(S\right)\bar{Q}{\text{e}}^{iS}}{{\omega }^{\prime }\left(S\right)\omega {\left|Q\right|}^2}\end{array}$ (22)

with ${\omega }^{\prime }\left(S\right)=\frac{1-{\tau }^{\prime }\left(S\right)\omega }{\tau }$.

4. Bifurcation of Periodical Solutions

As we have discussed in Section 3, Hopf bifurcation occurs with one pair of imaginary roots $i\omega$ with zero real part of the characteristic equation. Suppose $\left(K^{\mathrm{<mo>\; *\; </mo>}}\mathrm{,}{\tau }^{\mathrm{<mo>\; *\; </mo>}}\right)$ is Hopf point, then we compute the bifurcation direction and stability of periodic solution by perturbation method combined with analytical technique on the center manifold. Set $K=K^{*}+\epsilon {\tau }_{\epsilon },\tau ={\tau }^{*}+\epsilon {\tau }_{\epsilon }$, then Equation (1) can be written as its 3rd truncated expansion which listed as

$\begin{array}{c}{x}^{\prime }\left(t\right)=Ax+Bx\left(t-{\tau }^{\mathrm{<mo>\; *\; </mo>}}\right)+\epsilon K_{\epsilon }x\left(t\right)-\epsilon K_{\epsilon }x\left(t-{\tau }^{\mathrm{<mo>\; *\; </mo>}}\right)\\ +B\left(x\left(t-{\tau }^{\mathrm{<mo>\; *\; </mo>}}-\epsilon {\tau }_{\epsilon }\right)-x\left(t-{\tau }^{\mathrm{<mo>\; *\; </mo>}}\right)\right)+s_{1}x{\left(t-{\tau }^{\mathrm{<mo>\; *\; </mo>}}\right)}^2+s_{2}x{\left(t-{\tau }^{\mathrm{<mo>\; *\; </mo>}}\right)}^3\mathrm{,}\end{array}$ (23)

with $s_1=\frac{ca\left(\left(-c-1\right)\left(\frac{a-b}{b}\right)+\left(c-1\right){\left(\frac{a-b}{b}\right)}^2\right){\left(\frac{a-b}{b}\right)}^{-\frac{1}{c}}}{{\left(1+\frac{a-b}{b}\right)}^3}$ and

$s_2=\frac{ca{\left(\frac{a-b}{b}\right)}^{-\frac{2}{c}}\left(4{\left(\frac{a-b}{b}\right)}^2{c}^2-{\left(\frac{a-b}{b}\right)}^3{c}^2-\left(\frac{a-b}{b}\right)c^2+2{\left(\frac{a-b}{b}\right)}^2+{\left(\frac{a-b}{b}\right)}^3+\left(\frac{a-b}{b}\right)\right)}{{\left(1+\frac{a-b}{b}\right)}^4},$ (24)

with the phase space $C=\left\{\varphi \left(\theta \right)\mathrm{<mo>\; |\; </mo>}\varphi \in C\left(\left[-\tau \mathrm{,0}\right]\to R\right)\right\}$, $\varphi$ is left continuous at $\theta =0$, the supremum norm is defined as $\Vert \varphi \Vert ={\sup }_{-\tau \le \theta \le 0}\varphi \left(\theta \right)$, then based on the fundamental theory of funcational differential equations, there exists matrix function to write the linearized equation of Equation (4) as

$L\left(\epsilon \right)\varphi ={\displaystyle {\int }_{-{\tau }^{\mathrm{<mo>\; *\; </mo>}}}^0}\text{d}\eta \left(\theta \right)\varphi \left(\theta \right)+\epsilon {\displaystyle {\int }_{-{\tau }^{\mathrm{<mo>\; *\; </mo>}}}^0}\text{d}{\eta }_1\left(\theta \right)\varphi \left(\theta \right)+{\displaystyle {\int }_{-\tau }^0}\text{d}{\eta }_2\left(\theta \right)\varphi \left(\theta \right)$ (25)

with

$\begin{array}{l}\text{d}\eta \left(\theta \right)=\left[A\delta \left(\theta \right)+B\delta \left(\theta +{\tau }^{*}\right)\right],\text{d}{\eta }_1\left(\theta \right)=\left[K_{\epsilon }\delta \left(\theta \right)-K_{\epsilon }\delta \left(\theta +{\tau }^{*}\right)\right]\\ \text{d}{\eta }_2\left(\theta \right)=\left[B\delta (\theta +{\tau }^{*}+\epsilon {\tau }_{\epsilon })-B\delta \left(\theta +{\tau }^{*}\right)\right]\end{array}$ (26)

Equation (3) can be written as its opearator differential form

$x\left(t\right)=L\left(\epsilon \right)x_t+F\left(x_t\right)$ (27)

with nonlinear part

$F\left(\varphi \right)=\left(s_1\varphi {\left(-{\tau }^{\mathrm{<mo>\; *\; </mo>}}\right)}^2+s_2\varphi {\left(-{\tau }^{\mathrm{<mo>\; *\; </mo>}}\right)}^3\right)$ (28)

The solution operator of Equation (4) is a strong continuous semigroup with infinitesimal generator

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

The adjoint operator in the conjugate space $C^{\mathrm{<mo>\; *\; </mo>}}\left(\left[\mathrm{0,}\tau \right]\mathrm{,}R\right)$ is also defined as

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

For any $\varphi \in C\mathrm{,}\psi \in C^{\mathrm{<mo>\; *\; </mo>}}$, define the bilinear form as

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

Define the collection set $\Lambda =\left\{i\omega \mathrm{,}-i\omega \right\}$, and the corresponding eigenspace $P_{\Lambda }=\left\{cq\left(\theta \right)\mathrm{,}c\in ℂ\right\}$ given that

$\begin{array}{l}\mathcal{A}q\left(\theta \right)=i\omega q\left(\theta \right)\mathrm{,}-{\tau }^{\mathrm{<mo>\; *\; </mo>}}\le \theta \le \mathrm{0,}\\ {\mathcal{A}}^{\mathrm{<mo>\; *\; </mo>}}p\left(s\right)=-i\omega p\left(s\right)\mathrm{,}0\le s\le {\tau }^{\mathrm{<mo>\; *\; </mo>}}\end{array}$ (32)

and $\langle p,q\rangle =1,\langle \bar{p},q\rangle =0$. It is easily calculated that

$q\left(\theta \right)={\text{e}}^{i\omega \theta }\mathrm{,}-{\tau }^{\mathrm{<mo>\; *\; </mo>}}\le \theta \le 0$ (33)

and

$p\left(s\right)=N{\text{e}}^{i\omega s}\mathrm{,0}\le s\le {\tau }^{\mathrm{<mo>\; *\; </mo>}}$ (34)

with $N=\frac{1}{1+B{\tau }^{\mathrm{<mo>\; *\; </mo>}}{\text{e}}^{i\omega {\tau }^{\mathrm{<mo>\; *\; </mo>}}}}$. Suppose Q is the complementary subspace of $P_{\Lambda }$ and

the phase space is decomposed into $C=P_{\Lambda }\oplus Q$. With a possible discontinuous jump at $\theta =0$, we also define the map

$X_0=\{\begin{array}{l}0,-{\tau }_{\max }\le \theta <0,\\ I,\theta =0\end{array}$ (35)

Then for any $x=\phi +X_0\alpha \in C$, define the projection operator $\Pi \mathrm{<mo>\; :\; </mo>}{C}^{\prime }\left(\left[-\tau \mathrm{,0}\right]\to R\right)\to P_{\Lambda }$ as $\Pi \left(\phi +X_0\alpha \right)=q\left(\theta \right)\left[\langle p\mathrm{,}\phi \rangle +\bar{p}\left(0\right)w\right]$. Therefore, for any $x_t\in C$, it is written as $x_t=zq+\bar{z}\bar{q}+y_t$, substitute it into operator differential Equation (4) to obtain

$\begin{array}{l}{z}^{\prime }=i\omega z+\bar{p}\left(0\right)\left(L\left(\epsilon \right)x_t-L\left(0\right)x_t+F\left(x_t\right)\right)\\ =i\omega z+\bar{p}\left(0\right)\left(\epsilon {\displaystyle {\int }_{-{\tau }^{\mathrm{<mo>\; *\; </mo>}}}^0}\text{d}{\eta }_1\left(\theta \right)x_t\left(\theta \right)+{\displaystyle {\int }_{-\tau }^0}\text{d}{\eta }_2\left(\theta \right)x_t\left(\theta \right)+F\left(x_t\right)\right)\\ =i\omega z+\bar{N}\epsilon (K_{\epsilon }\left(zq\left(0\right)+\bar{z}\bar{q}\left(0\right)+y_t\left(0\right)-zq\left(-{\tau }^{\mathrm{<mo>\; *\; </mo>}}\right)-\bar{z}\bar{q}\left(-{\tau }^{\mathrm{<mo>\; *\; </mo>}}\right)-y_t\left(-{\tau }^{\mathrm{<mo>\; *\; </mo>}}\right)\right)\\ +\bar{N}(Bzq\left(-{\tau }^{\mathrm{<mo>\; *\; </mo>}}-\epsilon {\tau }_{\epsilon }\right)+\bar{z}\bar{q}\left(-{\tau }^{\mathrm{<mo>\; *\; </mo>}}-\epsilon {\tau }_{\epsilon }\right)+y_t\left(-{\tau }^{\mathrm{<mo>\; *\; </mo>}}-\epsilon {\tau }_{\epsilon }\right)\end{array}$

$\begin{array}{l}+\bar{N}\left(Bzq\left(-{\tau }^{\mathrm{<mo>\; *\; </mo>}}\right)+\bar{z}\bar{q}\left(-{\tau }^{\mathrm{<mo>\; *\; </mo>}}\right)+y_t\left(-{\tau }^{\mathrm{<mo>\; *\; </mo>}}\right)\right)\\ +\bar{N}\left(-s_1\right){\left(zq\left(-{\tau }^{\mathrm{<mo>\; *\; </mo>}}\right)+\bar{z}\bar{q}\left(-{\tau }^{\mathrm{<mo>\; *\; </mo>}}\right)+y_t\left(-{\tau }^{\mathrm{<mo>\; *\; </mo>}}\right)\right)}^2\\ +\bar{N}{s}_2{\left(zq\left(-{\tau }^{\mathrm{<mo>\; *\; </mo>}}\right)+\bar{z}\bar{q}\left(-{\tau }^{\mathrm{<mo>\; *\; </mo>}}\right)+y_t\left(-{\tau }^{\mathrm{<mo>\; *\; </mo>}}\right)\right)}^3\\ {y^{\prime }}_t=\mathcal{A}{y}_t+\mathrm{<mo\; stretchy="false">\; (\; </mo>}I-\Pi \mathrm{<mo\; stretchy="false">\; )\; </mo>}{X}_{0}F\mathrm{<mo\; stretchy="false">\; (\; </mo>}zq\left(-{\tau }^{\mathrm{<mo>\; *\; </mo>}}\right)+\bar{z}\bar{q}\left(-{\tau }^{\mathrm{<mo>\; *\; </mo>}}\right)+y_t\left(-{\tau }^{\mathrm{<mo>\; *\; </mo>}}\right)\end{array}$ (36)

The normal form on the center manifold is represented by the following

$\begin{array}{l}{z}^{\prime }=i\omega z+\bar{N}\epsilon K_{\epsilon }\left(z+\bar{z}-z{\text{e}}^{-i\omega {\tau }^{*}}-\bar{z}{\text{e}}^{i\omega {\tau }^{*}}\right)-\bar{N}B\epsilon {\tau }_{\epsilon }\left(z{\text{e}}^{-i\omega {\tau }^{*}}+\bar{z}{\text{e}}^{i\omega {\tau }^{*}}\right)\\ +\bar{N}\left(-s_1\right){\left(z{\text{e}}^{-i\omega {\tau }^{*}}+\bar{z}{\text{e}}^{i\omega {\tau }^{*}}+y_t\left(-{\tau }^{*}\right)\right)}^2+\bar{N}{s}_2{\left(z{\text{e}}^{-i\omega {\tau }^{*}}+\bar{z}{\text{e}}^{i\omega {\tau }^{*}}+y_t\left(-{\tau }^{*}\right)\right)}^3\\ {y^{\prime }}_t=\mathcal{A}{y}_t+\{\begin{array}{l}-2\Re \left(\bar{N}\left(-s_1\right){\left(z{\text{e}}^{-i\omega {\tau }^{\mathrm{<mo>\; *\; </mo>}}}+\bar{z}{\text{e}}^{i\omega {\tau }^{\mathrm{<mo>\; *\; </mo>}}}+y_t\left(-{\tau }^{\mathrm{<mo>\; *\; </mo>}}\right)\right)}^{2}q\left(\theta \right)\right)\mathrm{,}-{\tau }^{\mathrm{<mo>\; *\; </mo>}}\le \theta <\mathrm{0,}\hfill \\ +\left(-s_1\right){\left(z{\text{e}}^{-i\omega {\tau }^{\mathrm{<mo>\; *\; </mo>}}}+\bar{z}{\text{e}}^{i\omega {\tau }^{\mathrm{<mo>\; *\; </mo>}}}+y_t\left(-{\tau }^{\mathrm{<mo>\; *\; </mo>}}\right)\right)}^2\hfill \\ -2\Re \left(\bar{N}\left(-s_1\right){\left(z{\text{e}}^{-i\omega {\tau }^{\mathrm{<mo>\; *\; </mo>}}}+\bar{z}{\text{e}}^{i\omega {\tau }^{\mathrm{<mo>\; *\; </mo>}}}+y_t\left(-{\tau }^{\mathrm{<mo>\; *\; </mo>}}\right)\right)}^{2}q\left(\theta \right)\right)\mathrm{,}\theta =\mathrm{0,}\hfill \end{array}\end{array}$ (37)

Set

$y_t=Y_{20}\left(\theta \right)z^2+Y_{11}\left(\theta \right)z\bar{z}+Y_{02}\left(\theta \right){\bar{z}}^2$ (38)

By Equation (4) to get differentiation of all the coefficients $Y_{ij}\left(\theta \right),i+j=2,i,j\in N$ with respect to $\theta$ and the corresponding initial condition which is omitted. Hence we obtain

$\begin{array}{c}{z}^{\prime }=i\omega +\epsilon c_{1}z+\epsilon \left(g_{20}{z}^2+g_{11}z\bar{z}+g_{02}{\bar{z}}^2\right)+g_{30}{z}^3\\ +g_{21}{z}^2\bar{z}+g_{12}z{\bar{z}}^2+g_{03}{\bar{z}}^3+\cdots \mathrm{,}\end{array}$ (39)

$\begin{array}{l}{c}_1=\bar{N}{K}_{\epsilon }\left(1-{\text{e}}^{-i\omega {\tau }^{\mathrm{<mo>\; *\; </mo>}}}\right)-\bar{N}B{\tau }_{\epsilon }{\text{e}}^{-i\omega {\tau }^{\mathrm{<mo>\; *\; </mo>}}}\mathrm{,}\\ g_{20}=\bar{N}\left(-s_1\right)\left(2Y_{20}\left(-{\tau }^{\mathrm{<mo>\; *\; </mo>}}\right){\text{e}}^{-i\omega {\tau }^{\mathrm{<mo>\; *\; </mo>}}}+{\text{e}}^{-2i\omega {\tau }^{\mathrm{<mo>\; *\; </mo>}}}+Y_{20}{\left(-{\tau }^{\mathrm{<mo>\; *\; </mo>}}\right)}^2\right)\end{array}$

$\begin{array}{l}{g}_{11}=2\bar{N}\left(-s_1\right)\left(2{\text{e}}^{i\omega {\tau }^{\mathrm{<mo>\; *\; </mo>}}}{Y}_{20}\left(-{\tau }^{\mathrm{<mo>\; *\; </mo>}}\right)+2\right)\mathrm{,}\\ g_{02}=\bar{N}\left(-s_1\right){\text{e}}^{2i\omega {\tau }^{\mathrm{<mo>\; *\; </mo>}}}\mathrm{,}\\ g_{21}=\bar{N}\left(-s_1\right)\left(2Y_{11}\left(-{\tau }^{\mathrm{<mo>\; *\; </mo>}}\right){\text{e}}^{-2i\omega {\tau }^{\mathrm{<mo>\; *\; </mo>}}}+2Y_{11}\left(-{\tau }^{\mathrm{<mo>\; *\; </mo>}}\right)Y_{20}\left(-{\tau }^{\mathrm{<mo>\; *\; </mo>}}\right)\right)+\bar{N}{s}_2{\text{e}}^{-i\omega {\tau }^{\mathrm{<mo>\; *\; </mo>}}}\mathrm{.}\end{array}$ (40)

By set $z=\rho {\text{e}}^{i\theta }$, Equation (4) can be written into its polar form

${\rho }^{\prime }=\epsilon \rho \left(Re\left(c_1\right)+Re\left(c_2\right){\rho }^2\right)$ (41)

with

$c_2=\frac{1}{{\omega }^2}Re\left(i{g}_{20}{g}_{11}+\omega g_{21}\right)$ (42)

Hence, we deduce the following conclusion,

Theorem 4.1. The bifurcation direction of the small amplitude periodical solution arise from Hopf point determined by $\mu =\frac{Re\left(c_2\right)}{\frac{\text{d}Re\left(\lambda \right)}{\text{d}\tau }}$, which is super-critical

Hopf bifurcation if $\mu <0$ or sub-critical Hopf bifurcation if $\mu >0$ ; The bifurcating periodical solution is stable if $Re\left(c_2\right)<0$ and unstable if $Re\left(c_2\right)>0$.

For example, with fixed parameter $a=0.2,c=10$, two super-critical Hopf points are detected respectively at $b=0.148,\tau =16.5$ and $b=0.017,\tau =16.5$. By formula (4), it is calculated that $c_2=-0.2184741763<0$ and

$c_2=-0.3382864105<0$. Hence, the bifurcating periodical solutions are stable and the continuous stable periodical solutions with varying parameter b is shown in Figure 5(a). The bifurcating stable periodical solution at $b=0.13$ is shown in Figure 5(b) with blue lines, and the bifurcating period-2 solution at $b=0.12$ is plotted by green lines. The picture of phase portraits in Figure 5(c) is produced by x versus $x\left(t-\tau \right)$ with sub-figures. And we deduce a possible period-doubling bifurcation arising at $b=0.125$. Comparing with the period of oscillating solution before bifurcation, a periodical solution with two times period is produced since floquent multiplier is arriving at −1 at the bifurcation

(a) (b) (b) (c) (e) (f)

Figure 5. Fixed $\tau =16.5,K=-0.02$, periodical oscillating solutions of Equation (1) and the continuous of periodical solutions of Equation (1) as varying free parameter b. (a) The vicinity of periodical solution at Hopf points; (b) The periodical solution at $b=0.08$ and the period-2 solution at $b=0.1$ with $\tau =16.5$. (c) The corresponding phase portraits of (b); (d) The period-4 solution at $b=0.125,\tau =21$ ; (e) Fixed $\tau =21$, the bifurcation of periodical solutions on Poincare section $\stackrel{\dot{}}{x}\left(t\right)=0$ ; (f) Fixed $\tau =23.5$, the bifurcation of periodical solutions on Poincare section $\stackrel{\dot{}}{x}\left(t\right)=0$.

point. By varying time delay and choosing Poincare section $x^{\prime }\left(t\right)=0$, the senario of period 1, 2, 3, 4 and period 8 bifurcation of period solutions are pictured as shown in Figure 5(d). The state feedback control in system (1) produce stable periodical orbits within some time delay windows, which may function hematopoietic stem cells system to visualize stable oscillation solution instead of chaos.

5. Conclusions

The dynamics of a hematopoietic stem cell model with delay state feedback control is discussed. Underlying super-critical Hopf bifurcation, system lost its stability to experience periodical oscillation behavior. The curve of the rightmost characteristic roots is continuously simulated with varying free parameters, hence Hopf point is found as the roots curve cross the imaginary axis from the left half plane to the right half plane. Furthermore, the stable and unstable regime of the steady state is partitioned by Hopf bifurcation curve. The continuous bifurcating periodical solution is carried out with varying free parameters. It is discovered that period-1, 2, 3, 4, and period-7 solution arises underlying adding-period bifurcation and period-doubling bifurcation of periodical solutions. The results visualized the stable periodical orbits instead of chaos under the state feedback control with time delay.

In this paper, we discussed the hematological system model underlying delay state feedback control. However, for simplicity, the feedback delay is uniform with the mature delay as in stem cells growth stage. We will further discuss the hematological system which contains two different time delays in later paper.

Conflicts of Interest

The authors declare no conflicts of interest.

References

[1]	Smith, H. (2011) An Introduction to Delay Differential Equations with Applications to the Life Sciences. Vol. 57, Springer, New York. https://doi.org/10.1007/978-1-4419-7646-8
[2]	Kuang, Y. (1992) Delay Differential Equations with Application in Population Dynamics. Springer, New York.
[3]	Davis, L.C. (2003) Modifications of the Optimal Velocity Traffic Model to Include Delay Due to Driver Reaction Time. Physica A: Statistical Mechanics and Its Applications, 319, 557-567. https://doi.org/10.1016/S0378-4371(02)01457-7
[4]	Beddington, J.R. and Ray, R.M. (1975) Time Delays Are Not Necessarily Destabilizing. Mathematical Biosciences, 27, 109-117. https://doi.org/10.1016/0025-5564(75)90028-0
[5]	Bengea, S.C., Li, X.Q. and DeCarlo, R.A. (2004) Combined Controller-Observer Design for Uncertain Time Delay Systems with Application to Engine Idle Speed Control. Journal of Dynamic Systems Measurement and Control, 126, 772-780. https://doi.org/10.1115/1.1849239
[6]	Franklin, G.F., Powell, J.D. and Emami-Naeini, A. (2005) Feedback Control of Dynamic Systems. Pearson Prentice Hall, Upper Saddle River.
[7]	Frost, M.G. (1982) Controllability, Observability and the Transfer Function Matrix for a Delay-Differential System. International Journal of Control, 35, 175-182. https://doi.org/10.1080/00207178208922610
[8]	Belair, J. and Campbell, S.A. (1994) Stability and Bifurcations of Equilibria in a Multiple-Delayed Differential Equation. SIAM Journal on Applied Mathematics, 54, 1402-1424. https://doi.org/10.1137/S0036139993248853
[9]	Huang, K.L. and Lu, Q.S. (1995) Some Theorems for a Class of Dynamical System with Delay and Their Applications. Acta Mathematic Application Sinica, 18, 422-428. (in Chinese).
[10]	Shi, M. and Wang, Z.H. (2011) An Effective Analytical Criterion for Stability Testing of Fractional-Delay Systems. Automatica, 47, 2001-2005. https://doi.org/10.1016/j.automatica.2011.05.018
[11]	Wang, Z.H. and Hu, H.Y. (1999) Delay-Independent Stability of Retarded Dynamic Systems of Multiple Degrees of Freedom. Journal of Sound and Vibration, 226, 57-81. https://doi.org/10.1006/jsvi.1999.2282
[12]	Beretta, E. and Kuang, Y. (2002) Geometric Stability Switch Criteria in Delay Differential Systems with Delay Dependent Parameters. SIAM Journal on Mathematical Analysis, 33, 1144-1165. https://doi.org/10.1137/S0036141000376086
[13]	Cooke, K.L. and van den Driesche, P. (1986) On Zeros of Some Transcendental Equaions. Funkcialaj Ekvacioj, 29, 77-90.
[14]	Ma, S.Q., Lu, Q.S. and Mei, S. L. (2005) Dynamics of a Logistic Population Model with Maturation Delay and Nonlinear Birth Rate. Discrete and Continuous Dynamical Systems-B, 5, 735-752. http://dx.doi.org/10.3934/dcdsb.2005.5.735
[15]	Aiello, W.G., Freedman, H.I. and Wu, J. (1992) A Model of Stage Structural Population Growth with Density Dependent Time Delay. SIAM Journal on Applied Mathematics, 52, 855-869. https://doi.org/10.1137/0152048
[16]	Beretta, E. and Kuang, Y. (2001) Modeling and Analysis of a Marine Bacteriophage Infection with Latency Period. Nonlinear Analysis: Real World Applications, 2, 35-74. https://doi.org/10.1016/S0362-546X(99)00285-0
[17]	Ma, S.Q., Feng, Z.S. and Lu, Q.S. (2008) A Two Parameter Criteria for Delay Differential Equations. Discrete & Continuous Dynamical Systems-B, 9, 397-413. http://dx.doi.org/10.3934/dcdsb.2008.9.397
[18]	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
[19]	Xu, J. and Chung, K.W. (2003) Effects of Time Delayed Position Feedback on a Van Der Pol-Duffing Oscillator. Physica D: Nonlinear Phenomena, 180, 17-39. https://doi.org/10.1016/S0167-2789(03)00049-6
[20]	Wang, Z.H., Hu, H.Y., Xu, Q. and Stepan, G. (2016) Effect of Delay Combinations on Stability and Hopf Bifurcation of an Oscillator with Acceleration-Derivative Feedback. International Journal of Nonlinear Mechanics, 94, 392-399. https://doi.org/10.1016/j.ijnonlinmec.2016.10.008
[21]	Hale, J.K. and Lunel, S.M.V. (1993) Introduction of Functional Differential Equations. Springer-Verlag, New York. https://doi.org/10.1007/978-1-4612-4342-7