_{1}

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.

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 [

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 [

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

x ˙ = a x ( t − τ ) 1 + x ( t − τ ) c − b x (1)

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

x ˙ = a x ( t − τ ) 1 + x ( t − τ ) c − b x + K ( x ( t − τ ) − x ( t ) ) (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

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

illustrated as geometrical criterion in paper [

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.

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

Δ ( λ , σ , τ ) = P ( λ , σ ) + Q ( λ , σ ) e − λ τ (3)

For the fixed value of σ = σ * , τ = τ * , suppose λ = α + i β , substitute it into the characteristic Equation (2) then seperate the imaginary part from the real part to get

P R ( α , β , σ * ) + Q R ( α , β , σ * ) e − α τ * cos ( β τ * ) + Q I ( α , β , σ * ) e − α τ * sin ( β τ * ) = 0 , P I ( α , β , σ * ) + Q I ( α , β , σ * ) e − α τ * cos ( β τ * ) − Q R ( α , β , σ * ) e − α τ * sin ( β τ * ) = 0 (4)

By Equation (2) to get

e − α τ * cos ( β τ * ) = − P I Q I + P R Q R Q I 2 + Q R 2 , e − α τ * sin ( β τ * ) = P I Q R − P R Q I Q I 2 + Q R 2 (5)

Set S = β τ * , the equivalent equation of Equations (2) is obtained as

e − 2 α τ * = ( P I Q I + P R Q R Q I 2 + Q R 2 ) 2 + ( P I Q R − P R Q I Q I 2 + Q R 2 ) 2 , G ( α , S ) = − P I Q R − P R Q I P I Q I + P R Q R − tan ( S ) (6)

Hence we derive the following lemma 2.1,

Lemma 2.1: Suppose G ′ S ( α , S ) ≠ 0 , then there’s implicit function α = α ( S ) determined by G ( α , S ) ≡ 0 .

Define the line

L : Y 1 = e − 2 α ( S ) τ * (7)

and the line

R : Y 2 = ( P I Q I + P R Q R Q I 2 + Q R 2 ) 2 + ( P I Q R − P R Q I Q I 2 + Q R 2 ) 2 (8)

then the intersection point S * determines the corresponding characteristic roots λ = α ( S * ) + i β ( S * ) .

Setting x ( t − τ ) = x ( t ) , the positive steady state of Equation (1) and Equation

(1) is calculated as x * = exp ( 1 C ln ( a − b b ) ) . Doing the axis transformation x = x − x * , x ( t − τ ) = x ( t − τ ) − x * , one obtains the linear DDEs

x ′ ( t ) = A x ( t ) + B x ( t − τ ) (9)

with

A = − b + K , B = − ( − 1 + ( c − 1 ) ( ( ( a − b ) / b ) 1 c ) c ) a ( 1 + ( ( ( a − b ) / b ) 1 c ) c ) 2 − K (10)

Hence one has

P ( λ , b ) = A − λ , Q ( λ , b ) = B (11)

Equation (2) is rewritten as

e − 2 α τ = ( A − α B ) 2 + ( − S τ B ) 2 , G ( α , S ) = S τ ( A − α ) − tan ( S ) (12)

By the second equation in Equations (2) one gets

α = tan ( S ) τ A − S tan ( S ) τ (13)

Therefore by Lemma 2.1, one defines L : Y 1 = exp ( − 2 tan ( S ) τ A − S tan ( S ) ) and R : Y 2 = S 2 τ 2 B 2 sin ( S ) 2 . We also draw the following conclusion:

Lemma 2.2: Suppose the intersection point of line L and line R is ( S * , A * , B * ) , then the only positive equilibrium solution is asymptotically stable if Y 1 = Y 1 ( S * , A * , B * ) < 1 and versus if Y 1 = Y 1 ( S * , A * , B * ) > 1 . In addition, Hopf bifurcation is induced at ( S * , A * , B * ) if Y 1 = Y 1 ( S * , a * , B * ) = 1 .

The switching of stability of the only equilibrium solution is shown in

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

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.

As shown in

P R ( ω , σ ) + Q R ( ω , σ ) cos ( ω τ ) = 0, P I ( ω , σ ) − Q R ( ω , σ ) sin ( ω τ ) = 0 (14)

Hence Hopf bifurcation occurs if and only if

P R 2 + P I 2 − Q R 2 = 0 (15)

Set S = ω τ , by Equation (3), one gets

( A B ) 2 + ( − S τ B ) 2 − 1 = 0, S τ ( A ) − tan ( S ) = 0 (16)

and further to obtain Hopf line in ( b , K ) -plane. Hopf line is shown in

τ = ( 1 ± tan ( S ) 2 + 1 ) S tan ( S ) ( B 1 − b ) K = tan ( S ) τ b + S tan ( S ) τ (17)

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

Δ ( λ , K , τ ) = P ( λ , K ) + Q ( λ , K ) e − λ τ (18)

Differentiate Equation (3) with respect to time delay τ , one has

P ′ λ d λ d τ + ( Q ′ λ − τ Q ) e − λ τ d λ d τ − λ Q e − λ τ = 0 (19)

That is,

δ = s i g n R e d λ d τ = s i g n R e 1 d λ d τ = s i g n R e P ′ λ + ( Q ′ λ − τ Q ) e − λ τ λ Q e − λ τ (20)

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

P ′ λ i ω ′ ( S ) + P ′ K K ′ ( S ) + ( Q ′ λ i ω ′ ( S ) + Q ′ K K ′ ( S ) ) e − λ τ + Q e − λ τ ( − i ω ′ ( S ) τ − i ω τ ′ ( S ) ) = 0 (21)

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

δ ( S ) = s i g n R e 1 − i ω ′ ( S ) Δ ′ K K ′ ( S ) + Q e − i S ( − i ω ) τ ′ ( S ) i ω Q e − i S = s i g n R e Δ ′ K K ′ ( S ) Q ¯ e i S ω ′ ( S ) ω | Q | 2 (22)

with ω ′ ( S ) = 1 − τ ′ ( S ) ω τ .

As we have discussed in Section 3, Hopf bifurcation occurs with one pair of imaginary roots i ω with zero real part of the characteristic equation. Suppose ( K * , τ * ) 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 * + ε τ ε , τ = τ * + ε τ ε , then Equation (1) can be written as its 3rd truncated expansion which listed as

x ′ ( t ) = A x + B x ( t − τ * ) + ε K ε x ( t ) − ε K ε x ( t − τ * ) + B ( x ( t − τ * − ε τ ε ) − x ( t − τ * ) ) + s 1 x ( t − τ * ) 2 + s 2 x ( t − τ * ) 3 , (23)

with s 1 = c a ( ( − c − 1 ) ( a − b b ) + ( c − 1 ) ( a − b b ) 2 ) ( a − b b ) − 1 c ( 1 + a − b b ) 3 and

s 2 = c a ( a − b b ) − 2 c ( 4 ( a − b b ) 2 c 2 − ( a − b b ) 3 c 2 − ( a − b b ) c 2 + 2 ( a − b b ) 2 + ( a − b b ) 3 + ( a − b b ) ) ( 1 + a − b b ) 4 , (24)

with the phase space C = { ϕ ( θ ) | ϕ ∈ C ( [ − τ ,0 ] → R ) } , ϕ is left continuous at θ = 0 , the supremum norm is defined as ‖ ϕ ‖ = sup − τ ≤ θ ≤ 0 ϕ ( θ ) , then based on the fundamental theory of funcational differential equations, there exists matrix function to write the linearized equation of Equation (4) as

L ( ε ) ϕ = ∫ − τ * 0 d η ( θ ) ϕ ( θ ) + ε ∫ − τ * 0 d η 1 ( θ ) ϕ ( θ ) + ∫ − τ 0 d η 2 ( θ ) ϕ ( θ ) (25)

with

d η ( θ ) = [ A δ ( θ ) + B δ ( θ + τ * ) ] , d η 1 ( θ ) = [ K ε δ ( θ ) − K ε δ ( θ + τ * ) ] d η 2 ( θ ) = [ B δ ( θ + τ * + ε τ ε ) − B δ ( θ + τ * ) ] (26)

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

x ( t ) = L ( ε ) x t + F ( x t ) (27)

with nonlinear part

F ( ϕ ) = ( s 1 ϕ ( − τ * ) 2 + s 2 ϕ ( − τ * ) 3 ) (28)

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

A ϕ = { d ϕ d θ , − τ ≤ θ < 0, L ( 0 ) ϕ , θ = 0 (29)

The adjoint operator in the conjugate space C * ( [ 0, τ ] , R ) is also defined as

A * ψ = { − d ψ d s , 0 < s ≤ τ , − ∫ 0 τ * d η ( − s ) ψ ( s ) , s = 0 (30)

For any ϕ ∈ C , ψ ∈ C * , define the bilinear form as

〈 ϕ , ψ 〉 = ψ ¯ ( 0 ) ϕ ( 0 ) − ∫ − τ 0 ∫ 0 θ ψ ¯ ( ξ − θ ) d η ( θ ) ϕ ( ξ ) (31)

Define the collection set Λ = { i ω , − i ω } , and the corresponding eigenspace P Λ = { c q ( θ ) , c ∈ ℂ } given that

A q ( θ ) = i ω q ( θ ) , − τ * ≤ θ ≤ 0, A * p ( s ) = − i ω p ( s ) , 0 ≤ s ≤ τ * (32)

and 〈 p , q 〉 = 1 , 〈 p ¯ , q 〉 = 0 . It is easily calculated that

q ( θ ) = e i ω θ , − τ * ≤ θ ≤ 0 (33)

and

p ( s ) = N e i ω s ,0 ≤ s ≤ τ * (34)

with N = 1 1 + B τ * e i ω τ * . Suppose Q is the complementary subspace of P Λ and

the phase space is decomposed into C = P Λ ⊕ Q . With a possible discontinuous jump at θ = 0 , we also define the map

X 0 = { 0 , − τ max ≤ θ < 0 , I , θ = 0 (35)

Then for any x = φ + X 0 α ∈ C , define the projection operator Π : C ′ ( [ − τ ,0 ] → R ) → P Λ as Π ( φ + X 0 α ) = q ( θ ) [ 〈 p , φ 〉 + p ¯ ( 0 ) w ] . Therefore, for any x t ∈ C , it is written as x t = z q + z ¯ q ¯ + y t , substitute it into operator differential Equation (4) to obtain

z ′ = i ω z + p ¯ ( 0 ) ( L ( ε ) x t − L ( 0 ) x t + F ( x t ) ) = i ω z + p ¯ ( 0 ) ( ε ∫ − τ * 0 d η 1 ( θ ) x t ( θ ) + ∫ − τ 0 d η 2 ( θ ) x t ( θ ) + F ( x t ) ) = i ω z + N ¯ ε ( K ε ( z q ( 0 ) + z ¯ q ¯ ( 0 ) + y t ( 0 ) − z q ( − τ * ) − z ¯ q ¯ ( − τ * ) − y t ( − τ * ) ) + N ¯ ( B z q ( − τ * − ε τ ε ) + z ¯ q ¯ ( − τ * − ε τ ε ) + y t ( − τ * − ε τ ε )

+ N ¯ ( B z q ( − τ * ) + z ¯ q ¯ ( − τ * ) + y t ( − τ * ) ) + N ¯ ( − s 1 ) ( z q ( − τ * ) + z ¯ q ¯ ( − τ * ) + y t ( − τ * ) ) 2 + N ¯ s 2 ( z q ( − τ * ) + z ¯ q ¯ ( − τ * ) + y t ( − τ * ) ) 3 y ′ t = A y t + ( I − Π ) X 0 F ( z q ( − τ * ) + z ¯ q ¯ ( − τ * ) + y t ( − τ * ) (36)

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

z ′ = i ω z + N ¯ ε K ε ( z + z ¯ − z e − i ω τ * − z ¯ e i ω τ * ) − N ¯ B ε τ ε ( z e − i ω τ * + z ¯ e i ω τ * ) + N ¯ ( − s 1 ) ( z e − i ω τ * + z ¯ e i ω τ * + y t ( − τ * ) ) 2 + N ¯ s 2 ( z e − i ω τ * + z ¯ e i ω τ * + y t ( − τ * ) ) 3 y ′ t = A y t + { − 2 ℜ ( N ¯ ( − s 1 ) ( z e − i ω τ * + z ¯ e i ω τ * + y t ( − τ * ) ) 2 q ( θ ) ) , − τ * ≤ θ < 0, + ( − s 1 ) ( z e − i ω τ * + z ¯ e i ω τ * + y t ( − τ * ) ) 2 − 2 ℜ ( N ¯ ( − s 1 ) ( z e − i ω τ * + z ¯ e i ω τ * + y t ( − τ * ) ) 2 q ( θ ) ) , θ = 0, (37)

Set

y t = Y 20 ( θ ) z 2 + Y 11 ( θ ) z z ¯ + Y 02 ( θ ) z ¯ 2 (38)

By Equation (4) to get differentiation of all the coefficients Y i j ( θ ) , i + j = 2 , i , j ∈ N with respect to θ and the corresponding initial condition which is omitted. Hence we obtain

z ′ = i ω + ε c 1 z + ε ( g 20 z 2 + g 11 z z ¯ + g 02 z ¯ 2 ) + g 30 z 3 + g 21 z 2 z ¯ + g 12 z z ¯ 2 + g 03 z ¯ 3 + ⋯ , (39)

c 1 = N ¯ K ε ( 1 − e − i ω τ * ) − N ¯ B τ ε e − i ω τ * , g 20 = N ¯ ( − s 1 ) ( 2 Y 20 ( − τ * ) e − i ω τ * + e − 2 i ω τ * + Y 20 ( − τ * ) 2 )

g 11 = 2 N ¯ ( − s 1 ) ( 2 e i ω τ * Y 20 ( − τ * ) + 2 ) , g 02 = N ¯ ( − s 1 ) e 2 i ω τ * , g 21 = N ¯ ( − s 1 ) ( 2 Y 11 ( − τ * ) e − 2 i ω τ * + 2 Y 11 ( − τ * ) Y 20 ( − τ * ) ) + N ¯ s 2 e − i ω τ * . (40)

By set z = ρ e i θ , Equation (4) can be written into its polar form

ρ ′ = ε ρ ( R e ( c 1 ) + R e ( c 2 ) ρ 2 ) (41)

with

c 2 = 1 ω 2 R e ( i g 20 g 11 + ω g 21 ) (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 μ = R e ( c 2 ) d R e ( λ ) d τ , which is super-critical

Hopf bifurcation if μ < 0 or sub-critical Hopf bifurcation if μ > 0 ; The bifurcating periodical solution is stable if R e ( c 2 ) < 0 and unstable if R e ( c 2 ) > 0 .

For example, with fixed parameter a = 0.2 , c = 10 , two super-critical Hopf points are detected respectively at b = 0.148 , τ = 16.5 and b = 0.017 , τ = 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

point. By varying time delay and choosing Poincare section x ′ ( t ) = 0 , the senario of period 1, 2, 3, 4 and period 8 bifurcation of period solutions are pictured as shown in

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.

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

Ma, S.Q. (2021) Stability and Bifurcation Analysis of a Type of Hematopoietic Stem Cell Model International Journal of Modern Nonlinear Theory and Application, 10, 13-27. https://doi.org/10.4236/ijmnta.2021.101002