Global Asymptotic Stability and Hopf Bifurcation for a Blood Cell Production Model

We analyze the asymptotic stability of a nonlinear system of two differential equations with delay describing the dynamics of blood cell production. This process takes place in the bone marrow where stem cells differentiate throughout divisions in blood cells. Taking into account an explicit role of the total population of hematopoietic stem cells on the introduction of cells in cycle, we are lead to study a characteristic equation with delay-dependent coefficients. We determine a necessary and sufficient condition for the global stability of the first steady state of our model, which describes the population's dying out, and we obtain the existence of a Hopf bifurcation for the only nontrivial positive steady state, leading to the existence of periodic solutions. These latter are related to dynamical diseases affecting blood cells known for their cyclic nature.


Introduction
Blood cell production process is based upon the differentiation of so-called hematopoietic stem cells, located in the bone marrow. These undifferentiated and unobservable cells have unique capacities of differentiation (the ability to produce cells committed to one of the three blood cell types: red blood cells, white cells or platelets) and self-renewal (the ability to produce cells with the same properties).
Mathematical modelling of hematopoietic stem cells dynamics has been introduced at the end of the seventies by Mackey [21]. He proposed a system of two differential equations with delay where the time delay describes the cell cycle duration. In this model, hematopoietic stem cells are separated in proliferating and nonproliferating cells, these latter being introduced in the proliferating phase with a nonlinear rate depending only upon the nonproliferating cell population. The resulting system of delay differential equations is then uncoupled, with the nonproliferating cells equation containing the whole information about the dynamics of the hematopoietic stem cell population. The stability analysis of the model in [21] highlighted the existence of periodic solutions, through a Hopf bifurcation, describing in some cases diseases affecting blood cells, characterized by periodic oscillations [19].
The model of Mackey [21] has been studied by many authors, mainly since the beginning of the nineties. Mackey and Rey [23,24,25] numerically studied the behavior of a structured model based on the model in [21], stressing the existence of strange behaviors of the cell populations (like oscillations, or chaos). Mackey and Rudnicky [26,27] developed the description of blood cell dynamics through an age-maturity structured model, stressing the influence of hematopoietic stem cells on blood production. Their model has been further developed by Dyson et al. [13,14,15], Adimy and Pujo-Menjouet [7], Adimy and Crauste [2,3] and Adimy et al. [4]. Recently, Adimy et al. [5,6] studied the model proposed in [21] taking into account that cells in cycle divide according to a density function (usually gamma distributions play an important role in cell cycles durations), contrary to what has been assumed in the above-cited works, where the division has always been assumed to occur at the same time.
More recently, Pujo-Menjouet and Mackey [30] and Pujo-Menjouet et al. [29] gave a better insight into the model of Mackey [21], highlighting the role of each parameter of the model on the appearance of oscillations and, more particularly, of periodic solutions, when the model is applied to the study of chronic myelogenous leukemia [16].
Contrary to the assumption used in all of the above-cited works, we study, in this paper, the model introduced by Mackey [21] considering that the rate of introduction in the proliferating phase, which contains the nonlinearity of this model, depends upon the total population of hematopoietic stem cells, and not only upon the nonproliferating cell population. The introduction in cell cycle is partly known to be a consequence of an activation of hematopoietic stem cells due to molecules fixing on them. Hence, the entire population is in contact with these molecules and it is reasonable to think that the total number of hematopoietic stem cells plays a role in the introduction of nonproliferating cells in the proliferating phase.
The first consequence is that the model is not uncoupled, and the nonproliferating cell population equation does not contain the whole information about the dynamics of blood cell production, contrary to the model in [21,29,30]. Therefore, we are lead to the study of a modified system of delay differential equations (system (3)-(4)), where the delay describes the cell cycle duration, with a nonlinear part depending on one of the two populations.
Secondly, while studying the local asymptotic stability of the steady states of our model, we have to determine roots of a characteristic equation taking the form of a first degree exponential polynomial with delay-dependent coefficients. For such equations, Beretta and Kuang [9] developed a very useful and powerful technic, that we will apply to our model.
Our aim is to show, through the study of the steady states' stability, that our model, described in (3)-(4), exhibits similar properties than the model in [21] and that it can be used to model blood cells production dynamics with good results, in particularly when one is interested in the appearance of periodic solutions in blood cell dynamics models. We want to point out that the usually accepted assumption about the introduction rate may be limitative and that our model can display interesting dynamics, such as stability switches, that have never been noted before.
The present work is organized as follows. In the next section we present our model, stated in equations (3) and (4). We then determine the steady states of this model. In section 3, we linearize the system (3)-(4) about a steady state and we deduce the associated characteristic equation. In section 4, we establish necessary and sufficient conditions for the global asymptotic stability of the trivial steady state (which describes the extinction of the hematopoietic stem cell population). In section 5, we focus on the asymptotic stability of the unique nontrivial steady state. By studying the existence of pure imaginary roots of a first degree exponential polynomial with delay-dependent coefficients, we obtain the existence of a critical value of the time delay for which a Hopf bifurcation occurs at the positive steady state, leading to the appearance of periodic solutions. Using numerical illustrations, we show how these solutions can be related to periodic hematological diseases in section 6, and we note the existence of a stability switch. We conclude with a discussion.

A Nonlinear Model of Blood Cell Production
Let consider a population of hematopoietic stem cells, located in the bone marrow. These cells actually perform a succession of cell cycles, in order to differentiate in blood cells (white cells, red blood cells and platelets). According to early works, by Burns and Tannock [11] for example, we assume that cells in cycle are divided in two groups: proliferating and nonproliferating cells. The respective proliferating and nonproliferating cell populations are denoted by P and N .
All hematopoietic stem cells die with constant rates, namely γ > 0 for proliferating cells and δ > 0 for nonproliferating cells. These latter are introduced in the proliferating phase, in order to mature and divide, with a rate β. At the end of the proliferating phase, cells divide in two daughter cells which immediately enter the nonproliferating phase.
Then the populations P and N satisfy the following evolution equations (see Mackey [21] or Pujo-Menjouet and Mackey [30]), In each of the above equations, τ denotes the average duration of the proliferating phase. The term e −γτ then describes the survival rate of proliferating cells. The last terms in the right hand side of equations (1) and (2) account for cells that have performed a whole cell cycle and leave (enter, respectively) the proliferating phase (the nonproliferating phase, respectively). These cells are in fact nonproliferating cells introduced in the proliferating phase a time τ earlier. The factor 2 in equation (2) represents the division of each proliferating hematopoietic stem cell in two daughter cells. We assume that the rate of introduction β depends upon the total population of hematopoietic stem cells, that we denote by S. With our notations, S = P +N . This assumption stresses the fact that the nature of the trigger signal for introduction in the proliferating phase is the result of an action on the entire cell population. For example, it can be caused by molecules entering the bone marrow and fixing on hematopoietic stem cells, activating or inhibiting their proliferating capacity. This occurs in particularly for the production of red blood cells. Their regulation is mainly mediated by an hormone (a growth factor, in fact) called erythropoietin, produced by the kidneys under a stimulation by circulating blood cells (see Bélair et al. [8], Mahaffy et al. [28]).
The function β is supposed to be continuous and positive on [0, +∞), and strictly decreasing. This latter assumption describes the fact that the less hematopoietic stem cells in the bone marrow, the more cells introduced in the proliferative compartment [21,30]. Furthermore, we assume that lim S→∞ β(S) = 0.
Adding equations (1) and (2) we can then deduce an equation satisfied by the total population of hematopoietic stem cells S(t). We assume, for the sake of simplicity, that proliferating and nonproliferating cells die with the same rate, that is δ = γ. Then the populations N and S satisfy the following nonlinear system with time delay τ , corresponding to the cell cycle duration, From Hale and Verduyn lunel [18], for each continuous initial condition, system (3)-(4) has a unique continuous solution (S(t), N (t)), well-defined for t ≥ 0. Proof. First assume that there exists ξ > 0 such that N (ξ) = 0 and N (t) > 0 for t < ξ. Then, from (4) and since β is a positive function, Consequently, N (t) ≥ 0 for t > 0. If there exists ζ > 0 such that S(ζ) = 0 and S(t) > 0 for t < ζ, then the same reasoning, using (3), leads to dS dt (ζ) = e −δτ β(S(ζ − τ ))N (ζ − τ ) > 0, and we deduce that S(t) ≥ 0 for t > 0. Using a classical variation of constant formula, the solutions P (t) of (1) are given, for t ≥ 0, by Setting the change of variable σ = θ − τ , we obtain This condition is biologically relevant since 0 −τ e δθ β(S(θ))N (θ)dθ represents the population of cells that have been introduced in the proliferating phase at time θ ∈ [−τ, 0] and that have survived at time t = 0. Hence, from a biological point of view, the population in the proliferating phase at time t = 0 should be larger than this quantity.
In the stability analysis of system (3)-(4), the existence of stationary solutions, called steady states, is relevant since these particular solutions are potential limits of system (3)- (4).
A steady state of system (3)-(4) is a solution (S, N ) satisfying Let (S, N ) be a steady state of (3)-(4). Then We immediately notice that (0, 0) is a steady state of system (3)- (4), that we will denote, in the following, by E 0 . This steady state always exists. It describes the extinction of the hematopoietic stem cell population.
Assume that (S, N ) is a steady state of (3)-(4) with S, N = 0. Then, from (5) and (6), A necessary condition to obtain a nontrivial steady state is then Under this condition, since the function β is decreasing, positive, and tends to zero at infinity, there exists S > 0 satisfying if and only if (2e −δτ − 1)β(0) > δ.
One can easily check that condition (8) is equivalent to In this case, S > 0 solution of (7) is unique, and One can check that N ≤ S. These results are summed up in the next proposition.
In the next section, we linearize the system (3)-(4) about one of its steady states in order to analyze its local asymptotic stability.

Linearization and Characteristic Equation
We are interested in the asymptotic stability of the steady states of system (3)- (4). To that aim, we linearize system (3)-(4) about one of its steady state and we determine the associated characteristic equation. We assume that β is continuously differentiable on [0, +∞).
Let (S, N ) be a steady state of system (3) The linearization of system (3)-(4) about (S, N ) leads to the following system, where we have used the notations S(t) and N (t) instead of S(t) − S and N (t) − N for the sake of simplicity.
The system (10)-(11) can be written The characteristic equation of system (10)-(11) associated with the steady state (S, N ) is defined by After calculations, this equation reduces to We recall that the steady state (S, N ) is locally asymptotically stable when all roots of (13) have negative real parts and the stability can only be lost if eigenvalues cross the imaginary axis, that is if pure imaginary roots appear. One can notice that λ = −δ < 0 is always an eigenvalue of (13). Therefore, we only focus on the equation We first analyze, in the next section, the stability of the trivial steady state E 0 . We establish necessary and sufficient conditions for the population's dying out. Then, in section 5, we concentrate on the behavior of the positive steady state E * .
then the unique real root of (15) is λ 0 = 0, and all other eigenvalues have negative real parts. One can easily check that λ 0 = 0 is a simple root of (15), since the first derivative of the Then the linear system is stable, but we cannot conclude to the asymptotic stability of the trivial steady state E 0 = (0, 0) of system (3)-(4) without further analysis. This is done in Proposition 4.2.
When condition (16) holds true, E 0 is in fact the only steady state of system (3)-(4) (see Proposition 2.1). In this case, we can show that E 0 is globally asymptotically stable.
We first show the following result.
Proof. Using (3), a classical variation of constant formula gives us, for t ≥ 0, Setting σ = θ − τ , this expression becomes Let ε > 0 be fixed. Since N is assumed to tend to zero when t tends to ∞, there exists T > 0 such that Then, for t ≥ T + τ , Using (17) and the fact that β(0) is a bound of β, we obtain, for t ≥ T + τ , Let t > 0 be such that Then, for t ≥ max{t, T + τ }, we obtain Thus, S(t) tends to zero as t tends to +∞, and the proof is complete.
We then prove the following result, dealing with the global asymptotic stability of E 0 .
Then all solutions (S(t), N (t)) of system (3)- (4) converge to the trivial solution (0, 0). Hence E 0 is globally asymptotically stable and the cell populations dye out.
Proof. Let (S(t), N (t)) be a solution of (3)-(4). We define, for t ≥ 0, Using (4), we can check that From condition (18) and the fact that β is decreasing, it follows that Thus, Y is decreasing. Since Y is a nonnegative function, we deduce that there exists y ≥ 0 such that lim t→+∞ Y (t) = y.
In particularly, Y is bounded, and consequently N is also bounded. We then deduce, with (4), that N ′ is bounded and, using a similar technic than the one used in the proof of Lemma 4.1, that S is bounded. Consequently, with (3), we obtain that S ′ is bounded.
Hence, V is a Lyapunov functional (see Hale and Verduyn Lunel [18]) on the set With assumption (18), G = C 2 . In the proof of Proposition 4.2, we did not directly use the properties of Lyapunov functionals, but the function Y is defined by Through Propositions 4.1 and 4.2, we obtained necessary and sufficient conditions for the global asymptotic stability of E 0 . Therefore, in the next section, we concentrate on the behavior of E * , the unique positive steady state of (3)-(4).

Local Asymptotic Stability of the Positive Steady State
We now turn our considerations on the stability of the unique nontrivial steady state of system (3)-(4), namely E * = (S * , N * ), where, from Proposition 2.1, S * is the unique solution of (7), and N * = (2e −δτ − 1)e δτ S * .
In particularly, 2e −δτ − 1 > 0. In this case, Proposition 4.1 indicates that the unique other steady state E 0 is unstable.

-(4) is locally asymptotically stable, and the system (3)-(4) undergoes a transcritical bifurcation.
When τ increases and remains in the interval [0, τ ), the stability of the steady state can only be lost if purely imaginary roots appear. Therefore, we investigate the existence of purely imaginary roots of (20).

Remark 4. Assumption (H 1 ) is not necessary for the existence of τ * , it just implies the uniqueness. This latter can be achieved with weaker conditions, as we will check on an example in section 6.
We assume, in the following, that (H 1 ) and (H 2 ) are fulfilled and τ ∈ [0, τ * ). System (22)-(23) is equivalent to Note that for τ ∈ [0, τ * ), B(τ ) < 0. Therefore, adding the squares of both sides of (25), purely imaginary eigenvalues iω of (20), with ω > 0, must satisfy So, in the following, we will think of ω as ω(τ ). Substituting this expression for ω in (25), we obtain From the above reasoning, values of τ ∈ [0, τ * ) solutions of system (27) generate positive ω(τ ), given by (26), and hence yield imaginary eigenvalues of (20). Consequently, we look for positive solutions τ of (27) in the interval (0, τ * ). Positive solutions τ ∈ (0, τ * ) of (27) satisfy where N 0 denotes the set of all nonnegative integers. We set Values of τ for which ω(τ ) = B 2 (τ ) − A 2 (τ ) is a solution of (25) are roots of the functions The roots of Z k can be found using popular software like Maple, but are hard to determine with analytical tools [9]. The following lemma states some properties of the Z k functions.
Therefore, provided that no root of Z k is a local extremum, the number of positive roots of Z k , k ∈ N 0 , on the interval [0, τ * ) is even. Moreover, if Z k has no root on the interval [0, τ * ), then Z j , with j > k, does not have positive roots.
Remark 5. The second statement in Lemma 5.3 implies, in particularly, that, if Z 0 has no positive root, then (25) has no positive solution, and equation (20) does not have pure imaginary roots.
In the following proposition we establish some properties of pure imaginary roots of equation (20), using a method described in [20].
In the following, we do not mention the dependence of the coefficients A and B (and their derivatives) with respect to τ . Now, from (30), we obtain Since ∆(λ, τ ) = 0, we deduce For τ = τ c , we obtain Then, Noticing that Since iω(τ c ) is a purely imaginary root of (20), then, from (26), Substituting this expression in (31), we obtain, after simplifications As we already noticed, if equation (20) has pure imaginary roots then necessarily B < 0. We then deduce (29) and the proof is complete.
Using this last proposition and the previous results about the existence of purely imaginary roots of (13), we can state and prove the following theorem, dealing with the asymptotic stability of E * . (i) If Z 0 (defined in (28)) has no root on the interval [0, τ * ), τ * defined in Lemma 5.2, then the positive steady state E * = (S * , N * ) of (3)-(4) is locally asymptotically stable for τ ∈ [0, τ ).
(ii) If Z 0 has at least one positive root τ c ∈ (0, τ * ) then E * is locally asymptotically stable for τ ∈ [0, τ c ) and a Hopf bifurcation occurs at E * for τ = τ c if Proof. First, from Lemma 5.1, we know that E * is locally asymptotically stable when τ = 0.
If Z 0 has no positive root on the interval (0, τ * ), then the characteristic equation (13) has no pure imaginary root (see Remark 5 and Lemma 5.3). Consequently, the stability of E * cannot be lost when τ increases. We obtain the statement in (i). Now, if Z 0 has at least one positive root, say τ c ∈ (0, τ * ), then equation (13) has a pair of simple conjugate pure imaginary roots ±iω(τ c ) for τ = τ c . From (32) together with Proposition 5.1, we have either By contradiction, we assume that there exists a branch of characteristic roots λ(τ ) such that λ(τ c ) = iω c and dRe(λ(τ )) dτ < 0 for τ < τ c , τ close to τ c . Then there exists a characteristic root λ(τ ) such that Re(λ(τ )) > 0 and τ < τ c . Since E * is locally asymptotically stable when τ = 0, applying Rouché's Theorem [12], we obtain that all characteristic roots of (13) have negative real parts when τ ∈ [0, τ c ), and we obtain a contradiction. Thus, In this case, a Hopf bifurcation occurs at E * when τ = τ c .
The result stated in (ii) leads, through the Hopf bifurcation, to the existence of periodic solutions for system (3)-(4).
In the next section, we apply the above-mentioned results of stability to a particular introduction rate β and we present some numerical illustrations.

Example and Numerical Simulations
We develop, in this section, numerical illustrations of the above mentioned results (mainly the ones stated in Theorem 5.1).
The parameter β 0 represents the maximal rate of introduction in the proliferating phase, θ is the value for which β attains half of its maximum value, and n is the sensitivity of the rate of reintroduction. The coefficient n describes the reaction of β due to external stimuli, the action of a growth factor for example (some growth factors are known to trigger the introduction of nonproliferating cells in the proliferating phase [8,28]). Then, from (9), the unique positive steady state of (3)-(4) exists if and only if δ < β 0 and 0 ≤ τ < τ = 1 δ ln 2β 0 δ + β 0 .

From (5)-(6), it is defined by
After computations, we can state that condition (24) is equivalent to Noticing that the function χ, defined in Lemma 5.2, is given by χ(y) = −nβ 0 θ n y n (θ n + y n ) 2 , one can check that (33) is equivalent to (H 2 ).
However, the function χ does not necessarily satisfy (H 1 ), which is too strong (as mentioned in Remark 4). For y ≥ 0, Consequently, χ is decreasing for y ≤ θ and increasing for y > θ. Taking y = β −1 (δ), we find that χ is decreasing on [0, β −1 (δ)] if and only if β 0 < 2δ. In this case (H 1 ) is fulfilled. If β 0 > 2δ, then χ is decreasing on the interval [0, θ] and increasing on [θ, β −1 (δ)], yet τ * is uniquely defined. Note that  Figure 1: The functions Z 0 (left) and Z 1 (right) are drawn on the interval [0, τ * ) for parameters given by (34) and n = 12. One can see that Z 0 has exactly two roots, τ 1 ≃ 4.52 and τ 2 ≃ 8.36, and Z 1 has no root. with Then, assuming that (33) holds true, we define the Z k functions, as in (28), for τ ∈ [0, τ * ) by We choose the parameters according to [21,29,30]: Notice that the value of θ is in fact normalized and does not influence the stability of system (3)-(4) since all coefficients actually do not depend on θ. The value of θ only influences the shape of the oscillations and the values of the steady states. Using Maple to determine the roots of Z n , we first check that Z 0 (and consequently all Z k functions) is strictly negative on [0, τ * ) for n ≤ 10. Hence, from Theorem 5.1, the positive steady state E * = (S * , N * ) of (3)-(4) is locally asymptotically stable for τ ∈ [0, τ ).
One can see on Figure 1 that Z 0 has two positive roots in this case, τ 1 ≃ 4.52 days and τ 2 ≃ 8.36 days, and that Z 1 is strictly negative, so all Z k functions, with k ≥ 1 have no roots. Consequently, there exist two critical values, τ 1 and τ 2 , for which a stability switch can occur at E * . For τ < τ 1 , one can check that the populations are asymptotically stable on Figure 2. In this case τ = 3.5 days and the solutions of (3)-(4) oscillate transiently to the steady state. Numerical simulations of the solutions of (3)-(4) are carried out with dde23 [31], a Matlab solver for delay differential equations. When τ = τ 1 , one can check that so condition (32) holds, and a Hopf bifurcation occurs at (S * , N * ), from Theorem 5.1. This is illustrated on Figure 3. Periodic solutions with periods about 15 days are observed at the bifurcation, and the steady state E * becomes unstable.
When τ increases after the bifurcation, one can observe oscillating solutions with longer periods (in the order of 20 to 30 days), as it can be seen in Figure 4. This phenomenon has already been observed by Pujo-Menjouet et al. [29,30]. It can be related to diseases affecting blood cells, the so-called periodic hematological diseases [19], which are characterized by oscillations of circulating blood cell counts with long periods compared to the cell cycle duration. Among the wide variety of periodic hematological diseases, we can cite chronic myelogenous leukemia [6,16,29], a cancer of white blood cells with periods usually falling in the range of 70 to 80 days, and cyclical neutropenia [10,19] which is known to exhibit oscillations around 3 weeks of circulating neutrophils (white cells), as observed on Figure 4. Eventually, one can note that when τ passes through the second critical value τ 2 , stability switches and the steady state (S * , N * ) becomes stable again (see Figure 5).

Discussion
We considered a nonlinear model of blood cell dynamics in which the nonlinearity depends upon the entire hematopoietic stem cell population, contrary to the common assumption used in previous works [5,6,10,21,22,29,30] dealing with blood cell production models. Then we were lead to the study of a new nonlinear system of two differential equations with delay (describing the cell cycle duration) modelling the hematopoietic stem cells dynamics.
We obtained the existence of two steady states for this model: a trivial one and a positive delay-dependent steady state. Through sections 4 and 5, we performed the stability analysis of our model. We determined necessary and sufficient conditions for the global asymptotic stability of the trivial steady state of system (3)-(4), which describes the population's dying out. Using an approach proposed by Beretta and Kuang [9], we analyzed a first degree exponential polynomial characteristic equation with delay-dependent coefficients in order to obtain the existence of a Hopf bifurcation for the positive steady state (see Theorem 5.1), leading to the existence of periodic solutions.
On the example presented in the previous section, we obtained long periods oscillations, which can be related to some periodic hematological diseases (in particularly, to cyclical neutropenia [10]). This result is in keeping with previous analysis of blood cell dynamics models (as it can be found in [21,29,30]). Periodic hematological diseases are particular diseases mostly originated from the hematopoietic stem cell compartment. The appearance of periodic solutions in our model with periods that can be related to the ones observed in some periodic hematological diseases stresses the interesting properties displayed by our model. Periods of oscillating solutions can for example be used to determine the length of cell cycles in hematopoietic stem cell populations that cannot be directly determined experimentally. Moreover, stability switches have been observed, due to the structure of the equations (nonlinear equations with delay-dependent coefficients). Such a behavior had been noted in previous works dealing with blood cell production models (see [29,30]), but it had never been mathematically explained.
We can note that our assumption that proliferating and nonproliferating cells die with the same rate may be too limitative, since Pujo-Menjouet et al. [29,30] already noticed that the apoptotic rate (the proliferating phase mortality rate γ) plays an important role in the appearance of oscillating solutions. However, by assuming that the two populations die with different rates, we are lead to a second order exponential polynomial characteristic equation, and the calculations are more difficult than the ones carried out in the present work. We let it for further analysis.