Within its scope of dynamical discussion and biomedical discipline, stem cells are of great interest to give insight in understanding research work in its tissue organization. As a mathematical work, people devotes to use different methods to investigate both its spatial and temporal dynamics based on the understanding of stem cell proliferation and differentiation. Mathematical model which originates from hematopoietic stem cell compartment to periodic neutrophil blood diseases     , also the production and regulation of blood cells are governed by delay differential equations. Account for the variable success of granulocyte-colony stimulating factor, reduce of the amplification coefficient in neutrophil line for the treatment of cyclical neutropenia in reducing oscillation is proposed  . Recently, physiologically and mathematically, the dependence of neutrophil response on the period of simulated chemotherapy and the secondary response to G-CSF administration is reported     .
To support its function on tissue coupling, a clear link between the stem cell compartment and the differentiated mature cell lineages are motivated by delayed items to develop into fully hematopoietic system. The differentiation of hematopoietic stem cells take place in
phase. Quiescent phase HSCs (hematopoietic stem cells) can enter into its proliferative phase with the assumption of undergoing mitosis during time
. After a cell division, the total duration of both proliferative phase and maturation phase of neutropenia are experienced to release into circulation through the body. In general, considering tissue system composed of stem cells and cyclical neutrophils population, the mathematical model is formed by two equations with time delays which beyond illustration with schematic picture tools as shown in Figure 1.
The simple version related to DDEs (delay differential equations) to govern the system dynamics with above description have been reported as the following    :
where holling function
Figure 1. A cartoon representation of hematopoietic system response to stem cell population and neutrophil population.
HSCs (hematopoietic stem cells) enter into the proliferative phase with rate
(days−1), and differentiate into the committed neutrophil compartment at a rate
(days−1). The physiological meaning of other parameters and the corresponding values are explained as the following:
1) The hematopoiesis process consists of mechanism of triggering differentiation and maturation. HSCs produce differentiated cells via cell division to form into all kinds of cell types(white cells, red blood cells and platelets) and capable of self-renewal.
denotes the apoptosis rate (days−1) and the duration of the proliferation phase is taken to be
2) Different to the committed neutrophil line, cells in HSC proliferative phase are assumed to enter into the combined megaryocyte/erythrocyte lines at a rate
3) The time delay considered in neutrophil coupling is
, which denotes the total duration of the proliferation and maturation phases of a neutrophil precursor, that is,
. The proliferation rate is assumed as
whilst the death rate during maturation phase is
. Therefore, the amplification coefficient in neutrophil combined phases is listed as formula
Notice that the rates
and maturation time
are dependent on the effects of chemotherapy and G-CSF administration. The estimated parameter values are listed in Table 1.
From point view of dynamics, stable equilibrium solution can enter into instability regime to bifurcate into oscillation mode periodically via super-critical Hopf/sub-critical Hopf bifurcation  -  . Commonly, people use method of linearization model near the steady states to get the exponential polynomial characteristic equation to analyze its asymptotic stability via roots attribution to such polynomials   . Mathematically, the difficulty resides in the presence of two independent delays and the fact that some coefficients in the model equations depend upon these delays   . Therefore, there is very
Table 1. The values of parameters used in system (1).
few studies on this topic theoretically. With the fast development of computation methods in all kinds of mathematical questions, authors detect the different dynamic regime corresponding to stable steady states/periodic oscillating modes/bistable periodic oscillating modes and phase-locked resonance phenomena, etc.      .
The observation of bistability regime of steady state and periodic solution is found in Model (1). To understand the dynamics of the hematopoietic system, we take regard the amplification coefficient in stem cells as bifurcation parameter. Hopf bifurcation is tracked at its unique critical value. By means of the perturbation procedure such as the method of multiple scales, a codimension-2 bifurcation of periodic motion is theoretically discussed. Finally, varying maturation time delay simultaneously, dynamical bifurcations respond to periodic solutions is discovered and verified near condimension-2 bifurcation point by applying numerical simulation methods.
2. Linear Stability Analysis
Assume system (1) has an equilibrium solution
which is assumed to be positive. Notice coefficients as
are expressed by Hilling function which has obvious different analytical characters when varying
are assumed. By simple calculation, it is proved the formula
Do transformation of phase variables
, Equation (1) can be rewritten and truncated into its three order form as
We introduce variable vector z by assuming
, and parameter vector
. Equation (1) is written as its vector form
The hematopoietic model (1) has two delays in its dynamical description. As a common point view, we use the related linear system to analyze its asymptotic stability.
The corresponding characteristic equation is derived from Equation (11) as
At critical values of parameters, the onset of instability of linear system (11) may lead to oscillation solutions of nonlinear system. Hopf bifurcation can be tracked via the change of direction of characteristic root distribution which determined by Equation (12). Herein we solve possible values of Hopf bifurcation by software DDE-Biftool. The imaginary roots in complex plane are plotted in Figure 2(a) which shows Hopf bifurcation occurs since a pair of pure imaginary roots with zero real part appears in characteristic equation. If denote critical value of time delay as
, the direction of imaginary roots transverse in complex plane is shown in Figure 2(b). As time delay
, the head imaginary root(which is complex) transverse from right
plane to left plane and
is satisfied. Therefore, periodic
solutions may appear near the Hopf point.
Figure 2. The characteristic roots of Equation (12) are plotted in complex plane: (a) While
; (b) Varying time delay, head imaginary roots are observed to transverse through the vertical axis from right plane to left plane.
Figure 3 depicts a branch diagram of periodic solution emanating from the Hopf point as a function of
, respectively. Variation of solutions along the branch is characterized by their maximal and minimal values over the period for each computed point on the branch. As
grows from its Hopf point value, an instable periodic solution emanates which designates Hopf point as a subcritical Hopf point which brings bi-stability of system and hysteresis.
3. Multiple Scales Analysis
We consider the related vector form (8) in this section to discuss periodic solutions of system (1). Variable vector z depends on control vector parameter
and has subcritical Hopf bifurcation point at critical value
. In the following, the multiple scale method is used to investigate system dynamical behavior around subcritical Hopf point.
Without loss of generality, set characteristic matrix
then the corresponding critical eigenvectors are governed by
with the assumption
Subsequently, we introduce two time scales
in Equation(8) to rewrite it into the following formula
Figure 3. Branch of periodic solutions is seen to emanate from its subcritical Hopf point to form hysteresis which give rise to a stable and an unstable periodic solutions simultaneously: (a) Subcritical Hopf bifurcation occurs at critical value
; (b) Subcritical Hopf bifurcation occurs at critical value
According to multiple scale analysis, the mono-parameter family of solutions of Equation (16) may expand into its two scale form
. We set
, the above delay terms can both expand into their Taylor’s series as
We substitute Equation (17) and Equation (19) into the vector differential equation (16), and equating separately coefficients of
to obtain the perturbation equation
Therefore, we have
The sequence of perturbation equations can be obtained by equating separately coefficients of
Further, we suppose Equation (8) has the following general solution
is assumed to be the complex conjugate terms of the preceding terms in equation and frequency
satisfies characteristic Equation (12) by root
into Equation (21) and Equation (22), we have
respectively represent the resonance terms and non-resonance terms. It is well known that resonance terms can produce the secular terms and it satisfy the solvability condition
A particular solution of Equation (27) is given by
Therefore, solution is calculated as
The above process is repeated in the following analysis. Considering Equation (23) and Equation (24), we have
There is also a solvability condition for Equation (37) as
and one obtain the complex valued norm form as
Finally, doing the nonlinear transform
The first Lyapunov coefficient of the normal form is calculated    as
4. Bautin Bifurcation
Bautin bifurcation occurs if
, as shown in Figure 4, near Bautin point, Super-critical Hopf bifurcation occurs at point 2 and subcritical Hopf bifurcation occurs at Point 1 respectively. As shown in Figure 5, in region (1),
Figure 4. Bautin bifurcation of periodic solutions determined by its first Lyapunov coefficient.
Figure 5. System dynamics near Bautin point.
periodic solution with big amplitude appears due to Hopf bifurcation which switches the stability character of the positive steady state. In region (2), the coexistence of both stable and unstable periodic solutions is inspired due to the Bautin bifurcation which further form the bistable dynamic regime and hysteresis. The long time evolution dynamical behavior jumps to oscillating periodic limit cycle and the steady state is stable as shown in region (3). Choose parameters in system (1) with values as shown in Table 1, We calculate lyapunov coefficient which is referred above to be vicinity at Bautin point with codimension 2. The suitable values of parameter
are chosen and the corresponding eigenvector are given as listed in Table 2.
Table 2. Parameters, eigenvalues and eigenvectors at Bautin point.
In this work, we discussed the subcritical Hopf bifurcation of a hematopoiesis system which models two compartments of cell lines. The bistable regimes and hysteresis which are well known produced by subcritical Hopf bifurcation were detected near codimesion-2 bifurcation point (Bautin point).