Double Hopf Bifurcation of a Hematopoietic System with State Feedback Control

The dynamics of a system composed of hematopoietic stem cells and its relationship with neutrophils is ubiquitous due to periodic oscillating behavior induce cyclical neutropenia. Underlying the methodology of state feedback control with two time delays, double Hopf bifurcation occurs as varying time delay to reach its threshold value. By applying center manifold theory, the analytical analysis of system exposed the different dynamical feature in the classified regimes near double Hopf point. The novel dynamics as periodical solution and quasi-periodical attractor coexistence phenomena are explored and verified by numerical simulation.


Introduction
Mathematically and Biologically, the discussion of the complex dynamics of hematopoietic cell model is ubiquitous and interesting due to its inherent highly nonlinear attribute.The hematopoietic stem cells(HSCs) originates from the bone marrow to proliferate and differentiate into all other types of blood cells.Underlying biomedical discpline, people introduced population dynamics to describe the relationship between cells compartments and the application of delay differential equations is the necessary useful tool on hand.
The definition of a HSC has two properties: self-renewal and multipotent differentiation to cells with different function capabilities.People set forth all kinds of mathematcical models with the good understanding of how introduce compartment cells of proliferating phase into comparment cells of nonproliferating phase.Therefore, HSCs system is build on the basis that both nonproliferating phase and proliferating phase can produce, by successive divisions, all types of blood cells(white cells, red blood cells and platelets).Mackeys model proposed at the end of the seventies describe a decoupled quiescent phase of HSCS to study the dynamics of HSCs and its related diseases (Mackey, 1978).Since then it has been improved and analyzed by many authors, including Mackey and coauthors (Pujo−Menjouet & Mackey, 2004;Pujo−Menjouet et al., 2005;Fowler & Mackey, 2002).In their models, delay is an salient factor in introducing HSCs population models.Proliferating cells can die with apoptosis and divide after a given time to produce two daughter cells, then immedieatelly enter a new nonproliferating phase.
Due to the dysfunction in regulatory control process of blood cell production, some haematological disease emerge and cyclical neutropenia(CN) is the most extensively studied with low level oscillating phenomena happening.Medical observation have exhibited that the neutrophils fall down into abnormally level in a period about 19 to 21 days, and longer period happen in some patients with 40 to 80 days in human being.With the exception of period oscillating from 11 to 15 days, the similar oscillation character of cyclical neutropenia have been observed in the contrast experiment result of grey collie.With the assumption that neutrophils experience maturation time to release into circulation through the body, people have put forth the cyclical neutropenia dynamical system which is based on the negative feedack regulation mechanism (Fowler & Mackey, 2002;Bernard et al., 2003;Santillan et al., 2000;Ma et al.,2010).Haurie suggests that oscillation mechanism in cyclical neutropenia is due to destabilization of HSCs regulatory would explain the fact that other cell lineages oscillate with the same period as cyclical neutropenia (Haurie et al., 1998;Haurie et al., 2000;Hearn et al.,1998).Hopf bifurcation phenomenon has been observed in cyclical neutropenia model which lead to periodical oscillation of neutrophil numbers and the stable periodic solution branches out from equilibrium (Haurie et al., 1998).Biological meaningfully, Lei have developed CN model by introducing granulocyte-colony stimulating factors (G-CSF) into CN administration with the assumption that G-CSF may decrease the amplification coefficient in neutrophil line for the treatment of cyclical neutropenia (Lei & Mackey, 2014;Zhuge et al., 2012).Recently, the dependence of neutrophil response on the period of simulated chemotherapy and the secondary response to G-CSF adminstaration is also reported physiologically and mathematically (Haurie et al., 1998).
Lei developed cyclical neutropenia model including G-CSF administration which can be described as follows: Wherein, HSCs exist in its proliferating phase with nonlinear rate β, and enter into CN compartment with nonlinear rate k N .Commonly, both β and k N are given by the Holling function which listed In model (1.1), HSCs can be self-renewal and differentiate to produce all kinds of blood cells including white blood cells, platelets and erythrocytes with rate k δ .γ s denotes apoptosis rate and τ s denotes proliferation time.The delay in circulating neutrophil is τ N which can be decomposed as different time period experienced, that is, the proliferation time τ NP (days) of its progenitor cells and the maturation time τ N M (days) of cyclical neutropenia.η NP expresses the proliferation rate and γ 0 is the death rate in its maturation phase.Therefore, the amplification coefficient in circulating neutropenia is In this paper, we introduce state feedback control mechanism to describe system (1.1) as Hopf bifurcation and the instability of system (1.2) may bring out the complex dynamical behavior as varying parameter continuously.However, it is tremendous rough work to work out the threshold of Hopf bifurcation as varying time delay.Recently, DDE-Biftool is developed to effectively solve the root attribution of the related characteristic equation and compute the stability of the equlibrium solution of system (1.2).Hopf bifurcation occurs as a pair of imaginary roots of the characteristic equation cross the imaginary axis from left half plane to right half plane transversally and gives rise to oscillating periodic solutions with small amplitudes.In system (1.1),Hopf bifurcation point is tracked by regarding amplification coefficient as bifurcation parameters and bistable regime of steady states and periodic solution is found by using DDE-biftool software (Engelborghs et al., 2001;Green et al., 2002).
In this paper, we explore double Hopf bifurcation point of system (1.2) by tracking the secondary Hopf bifurcation lines which form the resonant region near "death island" (Zhang & Xu, 2013;Zhang & Xu, 2011;Ma et al., 2009).Different features as asymptotically stability in "death island" and periodic or chaotic attractor after double Hopf bifurcation exhibit the richness of system dynamics inherently.In section 2, we discuss the stability of the related linearized system and double Hopf bifurcation point which is highlighted by "death island" is explored by using DDE-Biftool.In section 3, by using center manifold method and normal form theory, dynamics of the linearized system near the double Hopf bifurcation point is analyzed.In section 4, dynamics of nonlinear system near the double Hopf point is analyzed.Numerical simulation further verifies the classified dynamical results, and coexistence of periodical solutions and torus are detected near double Hopf bifurcation point.
The linearized equation of system (1.2) is written as where The related characteristic equation of linearized system (2.1) is

Double Hopf Bifurcation
As analyzed in section 2, Double Hopf bifurcation happened at pairs of critical values (k c0 , τ N0 ).For convenience, set where , which is a Banach space composed as all continuous and differential function by mapping from [−τ N , 0] → R 2 .For every ϕ ∈ C , one can define the linear operator as Furthermore, we define with and Suppose the linear operator at double Hopf bifurcation point (k c0 , τ N0 ) is and The nonlinear operator R is defined as and ) .
System (3.1) can be rewritten as its operator form as (Hale, 2003) with u = (x, y) T , and R 2 ) define the adjoint operator H * (0) of operator H(0)£and adjoint operator H * (k ϵ , τ ϵ ) of H(k ϵ , τ ϵ ) ,respectively as By the above discussion, there are two pairs of imaginary roots at the double Hopf bifurcation point (k ϵ , τ ϵ ) , and other characteristic roots with negative real parts.We denote the collection set Λ = {iω 1 , −iω 1 , iω 2 , −iω 2 } and decompose the phase space C to be the direct sum of its subspaces as C = P Λ ⊕ Q Λ , with P Λ being the characteristic space of Λ and Q Λ being the corresponding complement subspace.System (3.1) is an ordinary differential system on the Banach space C of functions from [−τ N , 0] to R 2 , which is bounded and continuous on (−τ N , 0) with a possible jump discontinuity at 0 (Buono et al., 2003).We choose Φ(θ) = (q 1 (θ), q1 (θ), q 2 (θ), q2 (θ)), as the bases of P Λ and P * Λ , where for −τ N ≤ θ < 0, are respectively the eigenvectors of operator H(0) with respect to the characteristic root iω 1 , iω 2 , and the eigenvectors corresponding to adjoint operator A * (0) with respect to −iω 1 and −iω 2 are respectively as (3.17) for 0 ≤ s ≤ τ N ,where 1 is satisfied and further we can prove that Let z = (z 1 , z1 , z 2 , z2 ) T , for every u t ∈ C, and with the reduction form on its center manifold as where The matrix B generates the torus group T 2 = S 1 × S 1 whose action on C 2 is given by Then the T 2 -equivalent normal form, which is truncated to the quadratic order, becomes ( ż1 ż2 , we express the associated amplitude equation generated from system (3.19) as

Center Manifold Reduction
Furthermore, consider the nonlinear system ut = H(0)u t + Ru t (4.1) It is directly computed that the reduction form derived from Eq.(4.1) is with I identity and Π is the projection mapping from C([−τ, 0], R 2 ) to the subspace P Λ , and X 0 denotes We write v t into its Taylor series with coefficient vector w i jkl only depend on θ.We now expand f into the Taylor series and From system (4.1),we obtain The first equation in Eq.(4.6) can be rewritten as with the coefficients Taking account of formula of operator A, we have Integrating Eq.(4.9) under the initial condition given be the second equation in Eq.( 4.6), we can compute the coefficient vectors w i jkl (θ) with i Using the standard technique as invertible parameter-dependent complex coordinates transformation which is introduced in (Buono et al., 2003), the normal form at (k c0 , τ N0 ) of Eq.(4.1) is described as (4.10) In polar coordination (ρ 1 , θ 1 , ρ 2 , θ 2 ), the amplitude equation generated from Eq.(4.10) is Combining the results of Eqs(3.19) and Eqs(4.10), the formal normal form arising from system (4.11) is Using polar coordinate transformation z 1 = ρ 1 e iθ 1 , z 2 = ρ 2 e iθ 2 , system (4.12) is transformed into Using Maple software, we compute the coefficients g i jlk with i + j + k + l = 2, 3 ,and further obtain the following formulas (Zhang, S. & Xu, J., 2011).
At double Hopf point DH with Q * = 0.3082205820115, N * = 388.4842602444987,ω 1 = 0.254524564411299,ω 2 = 0.633751246398309, we compute the coefficients g 1,2 i jlk , (i + j + k + l = 2, 3).Therefore, by formula(4.14),chooseϵ = 0.1, then we have We also have By time continuous transformation with direction preserved, system (4.13) is equivalent to System (4.17) may have equilibrium solutions (ρ * 1 , 0), (0, ρ * 2 ) and (ρ * 1 , ρ * 2 ) and stability of equilibrium solution can be derived via computation of its corresponding Jacobin determinant.Therefore, based on system(4.16),we deduce the following results (I)Neutral saddle bifurcation of equilibrium (ρ = 0 and correspondingly, the associated torus of system (4.17)getting into torus instability state.As shown in Figure2, system (4.17) has the equilibrium solution (0, ρ * 2 ) with parameter pairs(k c , τ N ) chosen in the regime above boundary line L1 in green color, and the regime bounded by line in red color and line in black color exhibits torus solution.We discuss system dynamics in detail with the partition of plane into different regimes.The stable periodic solution observed in regime I, V I is shown in Figure3(a)-(f).The periodic solution exhibited in regime I further bifurcates into quasi-periodical solution after crossing neutral saddle line HL1 shown in cyan color.The scenarios of the bifurcating quasi-periodical solutions observed in regimes from (II) to V is plotted as shown in Figure4 .Due to neutral line HL2 in purple color , system may have multi attractors coexistence in regime IV.We also obtain the observed coexisted attractors in regime II and V alike, as shown in Figure5.

Discussion
A blood cell model which composed of two compartment cells known as hematopoietic stem cells and its relationship with neutrophil was investigated.The dysfunction in regulatory control process of blood cell production induce the observed hematological disease and cyclical neutropenia model including G-CSF adminstration was developed.As system state dependent delay reflects different time stage during feedback regulation process, we introduced state feedback control with time delay originally arising in system which is subtle disturbance around stationary equilibrium solution.We use DDE-Biftool software to explore double Hopf bifurcation of system via tracking second Hopf bifurcaton lines as varying parameters continuously and dynamics near double Hopf point was discussed.A classified method of dynamics near double Hopf point was developed by using center manifold theory.The analytical norm form near double Hopf point is computed based on fundamental theory of functional differential equation and the observed numerical simulation results were in coincidence with dynamics of the classified method which further verified the correctness of theoretical analytical results.

Figure 1 .
Figure 1.Double Hopf bifurcation occurs in system(1.2). (a)Hopf bifurcation lines drawn on (k 0 , τ N ) parameter plane£(b)The so-called "death island" is found near double Hopf bifurcation point and denoted by the shaded regime

Figure 2 .Figure 3 .
Figure 2. The dynamical classify results near the double Hopf point DH in parameter space (k 0 , τ N ) plane we have the expression u t = Φz + v t and the linearized system corresponds to(3.11)isut= H(0)u t + H(k ϵ , τ ϵ )u t ,(3.18)