Oscillations in a white blood cell production model with multiple differentiation stages

TL;DR

This paper demonstrates the occurrence of a super-critical Hopf bifurcation in a multi-stage white blood cell production model, explaining persistent oscillations observed in certain hematopoietic diseases.

Contribution

It provides an explicit analytical expression for the bifurcation point and highlights the unique oscillatory behavior arising from the multi-compartment structure.

Findings

Identifies parameter ranges leading to oscillations

Connects model results with clinical observations

Shows oscillations are absent in simpler models

Abstract

In this work we prove occurrence of a super-critical Hopf bifurcation in a model of white blood cell formation structured by three maturation stages. We provide an explicit analytical expression for the bifurcation point depending on model parameters. The Hopf bifurcation is a unique feature of the multi-compartment structure as it does not exist in the corresponding two-compartment model. It appears for a parameter set different from the parameters identified for healthy hematopoiesis and requires changes in at least two cell properties. Model analysis allows identifying a range of biologically plausible parameter sets that can explain persistent oscillations of white blood cell counts observed in some hematopoietic diseases. Relating the identified parameter sets to recent experimental and clinical findings provides insights into the pathological mechanisms leading to oscillating blood cell counts.

Figures

Tables

Table 1. Table 1: Summary of existence and stability conditions for the steady states depending on the parameter values of a 1 subscript 𝑎 1 a_{1} and a 2 subscript 𝑎 2 a_{2} . ”stable” refers to local asymptotic stability.

	$a_1>a_2$
	$E_0$: stable $E_1:\nexists$ $E_2:\nexists$	$E_0$: stable $E_1:\nexists$ $E_2:\nexists$
		$E_0$: unstable $E_1$: stable $E_2:\nexists$
	$E_0$: unstable $E_1:\nexists$ $E_2$: exists
$a_1>\frac{1}{2},a_2>\frac{1}{2}$	$E_0$: unstable $E_1$: unstable $E_2$: exists	$E_0$: unstable $E_1$: stable $E_2:\nexists$

Table 2. Table 2: Parameter configurations leading to a Hopf bifurcation

Constellation	$a_1$	$p_1$	$d_1$	$a_2$	$p_2$	$d_2$	$d_3$
1		$\uparrow$		$\downarrow$			$\downarrow$
2	$\uparrow$	$\uparrow$					$\downarrow$
3		$\uparrow$		$\uparrow$		$\uparrow$
4		$\uparrow$			$\downarrow$	$\uparrow$
5		$\uparrow$				$\uparrow$	$\downarrow$
6				$\uparrow$	$\downarrow$	$\uparrow$
7					$\downarrow$	$\uparrow$	$\downarrow$
8	$\uparrow$		$\uparrow$			$\uparrow$
9			$\uparrow$	$\uparrow$		$\uparrow$

Equations

$$s(t)=\frac{1}{1+k{u}_3(t)},$$

$$\frac{d{u}_1}{dt}$$

$$\frac{d{u}_2}{dt}$$

$$\frac{d{u}_3}{dt}$$

$$\frac{d{\tilde{u}}_1}{dt}$$

$$\frac{d{\tilde{u}}_2}{dt}$$

$$\frac{d{\tilde{u}}_3}{dt}$$

$$\frac{d{v}_1}{dt}$$

$$\frac{d{v}_2}{dt}$$

$$\begin{array}{c}\frac{d{u}_1}{dt}\le \left(2\frac{a_1}{1+\frac{k}{M_2{M}_1}{u}_1}-1\right)u_1<0,\text{ for }{u}_1>\frac{M_2{M}_1}{k}\left(2a_1-1\right)=:K_1\end{array}$$

$$\begin{array}{c}\frac{d{u}_2}{dt}\le \left(2\frac{a_2}{1+\frac{k}{M_2}{u}_2}-1\right)p_2{u}_2+2K_1<0\\ \text{for }{u}_2>\max \left\{\frac{4a_2-1}{k}{M}_2,\frac{4}{p_2}{K}_1\right\}=:K_2\end{array}$$

$$\begin{array}{c}\frac{d{u}_3}{dt}\le 2p_2{K}_2-d_3{u}_3<0\\ \text{for }{u}_3>\frac{2p_2{K}_2}{d_3}=:K_3.\end{array}$$

$$E_1={\left(0,\frac{d_3}{p_2}{\bar{u}}_3^1,\frac{2a_2-1}{k}\right)}^{\mathrm{T}}$$

$$a_2>\frac{1}{2}.$$

$$E_2={\left(\left(1-\frac{a_2}{a_1}\right)\frac{p_2}{p_1}{\bar{u}}_2^2,\frac{d_3}{\left(2-\frac{a_2}{a_1}\right)p_2}{\bar{u}}_3^2,\frac{2a_1-1}{k}\right)}^{\mathrm{T}}$$

$$a_1>\frac{1}{2}\text{ and }{a}_2<a_1.$$

$$k=\frac{2a_1-1}{{\bar{u}}_3},d_3=\frac{2-\frac{a_2}{a_1}}{1-\frac{a_2}{a_1}}\frac{{\bar{u}}_1}{{\bar{u}}_3}{p}_1,p_2=\frac{d_3}{2-\frac{a_2}{a_1}}\frac{{\bar{u}}_3}{{\bar{u}}_2}.$$

$$\begin{array}{rl}& J\left({\bar{u}}_1^i,{\bar{u}}_2^i,{\bar{u}}_3^i\right)=\\ & \left(\begin{array}{ccc}2\frac{a_1}{1+k{\bar{u}}_3^i}-1& 0& -2a_{1}k\frac{1}{{\left(1+k{\bar{u}}_3^i\right)}^2}{\bar{u}}_1^i\\ 2\left(1-\frac{a_1}{1+k{\bar{u}}_3^i}\right)& \left(2\frac{a_2}{1+k{\bar{u}}_3^i}-1\right)p_2& -2p_2{a}_{2}k\frac{1}{{\left(1+k{\bar{u}}_3^i\right)}^2}{\bar{u}}_2^i+2a_{1}k\frac{1}{{\left(1+k{\bar{u}}_3^i\right)}^2}{\bar{u}}_1^i\\ 0& 2\left(1-\frac{a_2}{1+k{\bar{u}}_3^i}\right)p_2& 2p_2{a}_{2}k\frac{1}{{\left(1+k{\bar{u}}_3^i\right)}^2}{\bar{u}}_2^i-d_3\end{array}\right).\end{array}$$

$$J\left(0,0,0\right)=\left(\begin{array}{ccc}2a_1-1& 0& 0\\ 2\left(1-a_1\right)& \left(2a_2-1\right)p_2& 0\\ 0& 2\left(1-a_2\right)p_2& -d_3\end{array}\right).$$

$${\lambda }_1^0$$

$${\lambda }_2^0$$

$${\lambda }_3^0$$

$$J\left({\bar{u}}_1^1,{\bar{u}}_2^1,{\bar{u}}_3^1\right)=\left(\begin{array}{ccc}\frac{a_1}{a_2}-1& 0& 0\\ 2-\frac{a_1}{a_2}& 0& -\left(1-\frac{1}{2a_2}\right)d_3\\ 0& p_2& -\frac{1}{2a_2}{d}_3\end{array}\right).$$

$$\left(\lambda -\frac{a_1}{a_2}+1\right)\left[\lambda \left(\lambda +\frac{1}{2a_2}{d}_3\right)+\left(1-\frac{1}{2a_2}\right)d_3{p}_2\right]=0$$

$${\lambda }_1^1$$

$${\lambda }_2^1$$

$${\lambda }_3^1$$

$$p_2>\frac{1}{1-\frac{a_2}{a_1}}\left[\frac{1}{\gamma \left(a_1,a_2\right)}-\beta \left(a_1,a_2\right)\frac{d_3}{p_1}\right]p_1.$$

$$p_2<\frac{1}{1-\frac{a_2}{a_1}}\left[\frac{1}{\gamma \left(a_1,a_2\right)}-\beta \left(a_1,a_2\right)\frac{d_3}{p_1}\right]p_1,$$

$$\beta \left(a_1,a_2\right)$$

$$\gamma \left(a_1,a_2\right)$$

$$J\left({\bar{u}}_1^2,{\bar{u}}_2^2,{\bar{u}}_3^2\right)$$

$$=$$

$$0=$$

$$+\underset{=:b_2}{\underbrace{\left[\left(1-\frac{a_2}{a_1}\right)\left(1-\frac{a_2}{a_1}\left(1-\frac{1}{2a_1}\right)\frac{1}{2-\frac{a_2}{a_1}}\right)-\left(1-\frac{1}{2a_1}\right)\left(1-2\frac{a_2}{a_1}\right)\right]d_3{p}_2}}\lambda$$

$$+\underset{=:b_3}{\underbrace{\left(1-\frac{1}{2a_1}\right)\left(1-\frac{a_2}{a_1}\right)d_3{p}_2}}.$$

$$b_1{b}_2-b_3>0\iff$$

Full text

Keywords: Hopf bifurcation, hematopoiesis, oscillating blood cell counts, mathematical model, stem cells.

1 Introduction

This work is devoted to the study of Hopf bifurcations and emergence of oscillatory dynamics in a multi-compartmental model of healthy blood cell production (hematopoiesis). The model describes a multi-stage blood cell production process based on self-renewal and differentiation of stem and progenitor cells, which is needed for regeneration of mature white blood cells. Each maturation stage is treated as a homogeneous compartment and its time evolution is described by an ordinary differential equation with coefficients controlled by a nonlinear feedback signal that depends on the count of mature cells. The model was introduced in ref. [47] and then has been applied to study blood cell recovery after bone marrow transplantation [60, 59] and extended to model evolution and response to therapy of hematological diseases such as acute leukemias [61, 63, 57, 58] and myelodysplastic syndromes [68]. Although mathematical understanding of the underlying equations proved to be useful for model applications and interpretation in context of the patients’ data [6, 57], rigorous analysis of the underlying equations has been established only in the case of a two-compartment maturation structure in ref. [20]. In this work, we close the gap and provide analysis of a system involving an intermediate differentiation stage given by a three-compartment structure. While in case of the two-compartment model, the positive equilibrium is globally stable whenever it exists [20], our analysis shows that increasing the number of compartments may lead to the loss of stability of the positive equilibrium due to a super-critical Hopf bifurcation. This finding is of biological relevance, since it shows that the number of maturation stages may impact the system dynamics.

Periodic oscillations in a model with non-linear feedback mechanisms but without explicit delays have not been studied in the context of hematopoietic system so far. Our model shows that we can have cycling hematopoiesis as inherent property of the multistep maturation process, however, arising far away from the parameter regime corresponding to the healthy system. Periodic oscillations of blood cell counts are a rare but intriguing phenomenon that can be observed in humans and animals [23, 13, 11]. Cyclic neutropenia is the most frequent disease with oscillating blood cell counts. In cyclic neutropenia patients’ neutrophil counts show periodic oscillations with maxima that are significantly below the neutrophil counts of healthy individuals. Since neutrophils are responsible for immune defence, patients repeatedly suffer from infections [15]. The disease can be cured by transplantation of healthy bone marrow [53, 14]. Similarly, accidental transplantation of bone marrow from a patient with cyclic neutropenia transfers the disorder to a previously unaffected host as has been shown in ref. [36]. A mechanistic understanding of the disease, therefore, requires quantitative insights into blood cell formation and its regulations.

Cyclic neutropenia has been extensively studied using mathematical models. One hypothesis derived from mathematical models is that oscillations are caused by increased apoptosis /reduced proliferation of neutrophil precursors [3, 41] combined with a reduced entry of stem cells into the proliferative phase [10]. An alternative mechanism could be an increase of the death rates of stem cells [41, 44]. Other model-derived hypotheses for the origin of cyclic neutropenia include reduced maturation speed [70] or dysfunction of feedback mechanisms [67]. Experimental studies suggest abnormal responsiveness of cells to growth factors [24, 72] or increased apoptosis of progenitor cells [22] as possible reasons for the origin of periodic oscillations. One common feature of most models of cyclic neutropenia is that they include constant or distributed delays. Intuitively, a system with feedbacks the effects of which occur with a delay, can be supposed to oscillate if the feedback loop gain is large enough. However, whether oscillations indeed appear, depends on configuration of feedbacks. In ref. [16] it has been shown that a reduction of progenitor cell’s self-renewal is sufficient to explain oscillatory dynamics in a model with linear feedback regulation and without delays. However, a system of linear compartments effectively constitutes a delay distributed according to a convolution of negative exponential functions, i.e. non-central gamma type, [30](Section 17.8.7). Other mathematical modeling works studying cyclic neutropenia include [35, 31, 27, 21]. Different modeling approaches are reviewed in [26, 9].

The paper is organized as follows. In Section 2 we present a derivation of the considered model and its biological justification. In Section 3 we provide analytical results, including uniform boundedness of solutions and linear stability analysis. We provide criteria for the occurrence of a Hopf bifurcation and illustrate them by model simulations. In Section 4 we study systematically for which subsets of the biologically relevant parameter space Hopf bifurcations occur and we relate our findings to experimental results. Section 5 concludes with a short summary and a discussion of the obtained results.

2 Model motivation and formulation

Blood cells are continuously produced during the life of higher metazoans. This task is fulfilled by the hematopoietic (blood forming) system which is located in the bone marrow [55]. Hematopoietic stem cells (HSC) give rise to progenitor cells which subsequently produce mature cells [55]. Due to its vital importance hematopoiesis is a tightly regulated process. Complex nonlinear feedback mechanisms allow the organisms to adapt to environmental conditions and to efficiently respond to perturbations such as blood loss or infection [48]. Key processes during hematopoiesis are cell proliferation, self-renewal and differentiation. Proliferation denotes the division of one parent cell into two progeny. If progeny are of the same cell type as the parent, e.g., a progeny of a stem cell is again a stem cell, this process is referred to as self-renewal. The alternative scenario, where progeny are of a more mature cell type compared to their parent cell is referred to as differentiation [64, 65].

In this work, we study for which configurations of proliferation and self-renewal parameters periodic oscillations of blood cell counts can occur. We focus on a three-compartment version of the model describing dynamics of stem cells, progenitor cells and mature cells. Dynamics of each cell population is described by an ordinary differential equation. Denoting the number of cells per kg of body weight at time $t$ as $u_{i}(t)$, where $i=1$ corresponds to stem cells, $i=2$ to progenitor cells and $i=3$ to mature white cells, each cell type is characterized by the following parameters:

• •

Proliferation rate $p_{i}$, describing the frequency of cell divisions per unit of time. In accordance with biology we assume that mature cells do not divide [29].

• •

Fraction of self-renewal $a_{i}$, describing the fraction of progeny cells originating from division and returning to the compartment of their parent cell.

• •

Death rate $d_{i}$, describing the fraction of cells dying per unit of time. For simplicity, we assume that immature cells ($i=1,2$) do not die and that mature cells die at constant rates. This is a good approximation of reality [60, 29],

White blood cell production is regulated by negative feedback signals, such as G-CSF [48, 40]. Since signal dynamics take place on a faster time scale compared to cell divisions, a quasi-steady state approximation can be used to describe the signal concentration as a function of white blood cell counts [47, 62, 46].

[TABLE]

where $k>0$, see [47] for details. Rigorous proof of the corresponding quasi-steady state model reduction is presented in ref. [46].

Following the previous work [47, 60, 59], we assume feedback inhibition of the fraction of self-renewal by mature cells, i.e. $a_{i}(t)=a_{i}s(t)$.

The flux to division of healthy cells in compartment $i$ at time $t$ equals $p_{i}{u}_{i}(t)$. During division, a parent cell is replaced by two progeny cells. The outflux from mitosis at time $t$, therefore, equals $2p_{i}{u}_{i}(t)$, of which the fraction $2a_{i}(t)p_{i}{u}_{i}(t)$ stays in compartment $i$ (referred to as self-renewal). The fraction $2(1-a_{i}(t))p_{i}{u}_{i}(t)$ proceeds to maturation stage $i+1$ (process referred to as differentiation). Taking into account that mature cells do not divide and that the parent cell disappears as it gives rise to its progeny, we obtain the following system of differential equations

[TABLE]

where $p_{1},p_{2}>0,d_{3}>0$ and $k>0$. The initial conditions fulfill $u_{1}(0)>0$, $u_{2}(0)\geq 0$, $u_{3}(0)\geq 0$. A schematic of the model is depicted in Figure 1.

In order to simplify our calculations, we finally rewrite the model equations in dimensionless terms using the reparametrizations $\tilde{t}:=p_{1}\,t$, $\tilde{p}_{2}:=\frac{p_{2}}{p_{1}}$, $\tilde{d_{3}}:=\frac{d_{3}}{p_{1}}$, and $\tilde{u}_{i}\left(\tilde{t}\right):=u_{i}\left(\frac{\tilde{t}}{p_{1}}\right)$ for $i=1,2,3$:

[TABLE]

The parameter ranges and initial conditions remain unchanged. For convenience, we drop the ~-symbol in the remainder of this paper.

3 Model analysis

In this section we show uniform boundedness of solutions (Section 3.1), provide conditions for existence of non-negative steady states (Section 3.2) and linearized stability analysis (Section 3.3), and prove occurrence of Hopf bifurcation (Section 3.4).

3.1 Uniform boundedness of solutions

Theorem 1**.**

Solutions of (M1)-(M3) with positive initial values remain in the first octant and for sufficiently large times $t$ even in a compact subset $C$ which does not depend on the initial values.

Proof.

The proof follows the lines of the proof in ref. [6], adjusted to the three-compartment structure of the model. We consider the rescaled system (M1*)-(M3*) and define $C$ as the cuboid bounded by the planes $\{(u_{1},u_{2},u_{3})\in\mathbb{R}^{3}\mid u_{i}=0\}$ and $\{(u_{1},u_{2},u_{3})\in\mathbb{R}^{3}\mid u_{i}=C_{i}\}$ for $i=1,2,3$. Obviously, the orbits of solutions with nonnegative initial values remain in the first octant. To show uniform boundedness from above, we compute equations for the fractions $v_{1}:=\frac{u_{1}}{u_{2}}$ and $v_{2}:=\frac{u_{2}}{u_{3}}$ that lead to the following estimates,

[TABLE]

for all $u_{3}\geq 0$. Taking $M_{i}:=\max\{B_{i},v_{i}(0)\}$, we conclude $u_{2}\geq\frac{1}{M_{1}}\,u_{1}$ and $u_{3}\geq\frac{1}{M_{2}}u_{2}\geq u_{1}\,\frac{1}{M_{2}\,M_{1}}$. With these relations we obtain

[TABLE]

and

[TABLE]

as well as

[TABLE]

Taking $C_{i}=\max\{K_{i},u_{i}(0)\}$, we conclude about invariance of the set $C$.

∎

3.2 Existence of steady states

Existence and uniqueness of steady states has been systematically studied in ref. [62]. The following proposition summarizes the results.

Proposition 2**.**

The trivial steady state $E_{0}=\left(\bar{u}_{1}^{0},\bar{u}_{2}^{0},\bar{u}_{3}^{0}\right)^{\mathrm{T}}=\left(0,0,0\right)^{\mathrm{T}}$ of (M1) - (M3) exists for all parameter values. 2. 2.

There exists a semi-trivial steady state $E_{1}=\left(0,\bar{u}_{2}^{1},\bar{u}_{3}^{1}\right)^{\mathrm{T}}$ of (M1) - (M3) with positive components $\bar{u}_{2}^{1},\bar{u}_{3}^{1}$ given by

[TABLE]

if and only if

[TABLE] 3. 3.

There exists a strictly positive steady state $E_{2}=\left(\bar{u}_{1}^{2},\bar{u}_{2}^{2},\bar{u}_{3}^{2}\right)^{\mathrm{T}}$ of (M1) - (M3) given by

[TABLE]

if and only if

[TABLE]

Let us remark that for appropriate values of $d_{3}$, $p_{2},\text{ and }k$ the steady state $E_{2}$ can be any point in $\mathbb{R}^{3}_{+}$. We have $E_{2}=(\bar{u}_{1},\bar{u}_{2},\bar{u}_{3})$ if

[TABLE]

3.3 Linear asymptotic stability

In Proposition 3 we summarize the linear asymptotic stability of the steady states $E_{0}$ and $E_{1}$. In Theorem 4 we study the linear asymptotic stability of $E_{2}$ using the Routh-Hurwitz Criterion.

Proposition 3**.**

1. (i)

The steady state $E_{0}=\left(0,0,0\right)^{\mathrm{T}}$ is locally asymptotically stable if $\operatorname{max}\{a_{1},a_{2}\}<\frac{1}{2}$ and unstable if $\operatorname{max}\{a_{1},a_{2}\}>\frac{1}{2}$. 2. (ii)

The steady state $E_{1}=\left(0,\frac{d_{3}}{p_{2}}\,\frac{2\,a_{2}-1}{k},\frac{2\,a_{2}-1}{k}\right)^{\mathrm{T}}$ is locally asymptotically stable if $a_{1}<a_{2}$ and unstable if $a_{1}>a_{2}$.

Proof.

The Jacobian of system (M1*)-(M3*) at the steady state $\left(\bar{u}_{1}^{i},\bar{u}_{2}^{i},\bar{u}_{3}^{i}\right)$is given by

[TABLE]

(i) Consider the Jacobian matrix at $E_{0}$

[TABLE]

As $J\left(0,0,0\right)$ is a lower triangular matrix, we obtain the eigenvalues

[TABLE]

which implies (i).

(ii) Consider the Jacobian matrix at $E_{1}$

[TABLE]

We recall that for existence of $E_{1}$ it has to hold $a_{2}>1/2$. We obtain the characteristic equation

[TABLE]

and thus the eigenvalues

[TABLE]

As $\lambda_{2}^{1}$ and $\lambda_{3}^{1}$ have always negative real parts, the only condition that needs to hold for $E_{1}$ to be locally asymptotically stable is $a_{1}<a_{2}$. If $a_{1}>a_{2}$, $E_{1}$ is unstable.

∎

For $E_{2}$, we will show local asymptotic stability using the Routh-Hurwitz-Criterion.

Using the expression for $E_{2}$ from Proposition 2 we can prove the following

Theorem 4**.**

The positive steady state $E_{2}$ of system (M1) - (M3) is locally asymptotically stable if

[TABLE]

$E_{2}$* unstable if*

[TABLE]

where

[TABLE]

Proof.

We perform calculations for the transformed system (M1*)-(M3*).

Consider the Jacobian matrix at $E_{2}$, i.e.

[TABLE]

We obtain the characteristic equation

[TABLE]

We observe that under positivity conditions for $E_{2}$ ($a_{1}>a_{2}$ and $a_{1}>\frac{1}{2}$) the relations $b_{1}>0$ and $b_{3}>0$ hold true: $b_{3}$ is a product with positive factors only and therefore positive itself. The expression $b_{1}$ can be written as $b_{1}=\left(1-\frac{a_{2}}{a_{1}}\right)\,p_{2}+\left(1-P\right)\,d_{3}$ where $P$ is a product consisting of factors that are in $(0,1)$ under positivity conditions. Thus, $1-P$ is positive and, therefore, $b_{1}$ is, as a sum of two positive summands, positive as well.

We distinguish between the following parameter configurations. Details can be found in the book by Gantmacher [17].

[TABLE]

It remains to determine conditions so that the relations $b_{1}\,b_{2}-b_{3}>0$ and $b_{1}\,b_{2}-b_{3}<0$ respectively are satisfied, to complete the proof. Using $\beta\left(a_{1},a_{2}\right)$ and $\gamma\left(a_{1},a_{2}\right)$ as defined in the theorem, we can rearrange $b_{1}\,b_{2}-b_{3}$ as follows. Further details on how to proceed can be found in the appendix A.

[TABLE]

As the factors $\left(1-\frac{1}{2\,a_{1}}\right)$ and $\left(1-\frac{a_{2}}{a_{1}}\right)$ are positive under positivity conditions for $E_{2}$ and $d_{3}$ and $p_{2}$ are positive anyway, it suffices to find conditions which ensure that the remaining factor is of positive sign. We find

[TABLE]

Thus, the eigenvalues of the Jacobian matrix at $E_{2}$ have negative real parts if and only if relation (4) holds. This implies local asymptotic stability. In analogy there exist eigenvalues with positive real parts if and only if $b_{1}\,b_{2}-b_{3}<0$. Transforming back to the parameters of system (M1)-(M3) completes the proof of this theorem. ∎

Finally, we will see that the parameter range where $E_{2}$ exists and is unstable is bounded:

Proposition 5**.**

For $p_{1}=1$ the set

[TABLE]

i.e. the parameter range where $E_{2}$ exists, is positive and is unstable, is a bounded subset of the parameter space.

Proof.

First, we note that for all $\left(a_{1},a_{2},d_{3},p_{2}\right)\in A$ it holds that

[TABLE]

and thus

[TABLE]

where

[TABLE]

As the range of $a_{1}$ and $a_{2}$ is bounded anyway, it suffices to show that the sets $A_{a_{1},a_{2}}$ are uniformly bounded. To this end, we observe that for fixed $a_{1}$ and $a_{2}$ the boundary of each $A_{a_{1},a_{2}}$ consists of the part of the graph of the linear equation

[TABLE]

lying in the first quadrant and its axis intercepts (see Figure 2).

We will show that the intercepts of the $d_{3}$- and the $p_{2}$-axis $\frac{1}{\beta\left(a_{1},a_{2}\right)\,\gamma\left(a_{1},a_{2}\right)}$ and $\frac{1}{(1-\frac{a_{2}}{a_{1}})\cdot\gamma\left(a_{1},a_{2}\right)}$ respectively are uniformly bounded. For $p_{2}\in A_{a_{1},a_{2}}$ for all $a_{1},a_{2}\in(0,1)$ it holds:

[TABLE]

This is equivalent to $p_{2}<2$ for all $a_{1},a_{2}\in(0,1)$. Furthermore, we note that

[TABLE]

and thus, we also have for all $d_{3}\in A_{a_{1},a_{2}}$ $d_{3}\leq\frac{1}{\beta\left(a_{1},a_{2}\right)\,\gamma\left(a_{1},a_{2}\right)}<2\ \forall a_{1},a_{2}\in(0,1)$. This means that the set $A_{a_{1},a_{2}}$ are uniformly bounded and therefore, the set A consisting of all parameters for which $E_{2}$ is unstable is bounded. ∎

Table 1 summarizes existence and local asymptotic stability of the equilibria $E_{0}$ to $E_{2}$.

Corollary 6**.**

The steady state $E_{2}$ exists and is locally asymptotically stable if and only if

[TABLE]

In this case all other non-negative steady states are unstable.

3.4 Hopf bifurcation

In this section, we will further investigate the change of the dynamical behavior when $p_{2}$ passes through $\frac{1}{1-\frac{a_{2}}{a_{1}}}\,\left[\frac{1}{\gamma\left(a_{1},a_{2}\right)}-\beta\left(a_{1},a_{2}\right)\,\frac{d_{3}}{p_{1}}\right]p_{1}$. We note that this is only possible if $\frac{1}{1-\frac{a_{2}}{a_{1}}}\,\left[\frac{1}{\gamma\left(a_{1},a_{2}\right)}-\beta\left(a_{1},a_{2}\right)\,\frac{d_{3}}{p_{1}}\right]p_{1}>0$ holds, i.e. for $d_{3}<d_{3}^{\text{max}}:=\frac{p_{1}}{\beta\left(a_{1},a_{2}\right)\,\gamma\left(a_{1},a_{2}\right)}$.

Theorem 7**.**

Let $d_{3}^{\text{max}}:=\frac{p_{1}}{\beta\left(a_{1},a_{2}\right)\,\gamma\left(a_{1},a_{2}\right)}$ and $d_{3}<d_{3}^{\text{max}}$. Then the steady state $E_{2}$ undergoes a Hopf bifurcation with bifurcation point $p_{2}=p_{2}^{*}:=\frac{1}{1-\frac{a_{2}}{a_{1}}}\,\left[\frac{1}{\gamma\left(a_{1},a_{2}\right)}-\beta\left(a_{1},a_{2}\right)\,\frac{d_{3}}{p_{1}}\right]p_{1}$, i.e. the Jacobian matrix $J$ at the positive steady state $E_{2}$ has two eigenvalues $\lambda_{1},\lambda_{2}$ for which the following relations hold

[TABLE]

Proof.

We consider the system (M1*)-(M3*). We recall that existence of $E_{2}$ requires $a_{1}>0.5$ and $a_{2}<a_{1}$. Let $P(x)=x^{3}+b_{1}\,x^{2}+b_{2}\,x+b_{3}$ be the characteristic polynomial of $J$. From the proof of Theorem 4 we know that $J$ has two purely imaginary eigenvalues unequal to zero if and only if $b_{1}\,b_{2}-b_{3}=0$, i.e., $p_{2}=p_{2}^{*}$. Thus, for $\lambda_{1,2}(p_{2}):=\mu(p_{2})\pm\omega(p_{2})$ it holds $\mu(p_{2}^{*})=0$, $\omega(p_{2}^{*})\neq 0$. It remains to show that $\mu^{\prime}(p_{2}^{*})\neq 0$.

Let us rewrite the characteristic polynomial $P(x)$ as

$P(x)=\left(x-\lambda_{1}\right)\,\left(x-\lambda_{2}\right)\,\left(x-\lambda_{3}\right)$, where $\lambda_{3}$ denotes the third eigenvalue. We can verify easily that the relation $b_{1}\,b_{2}-b_{3}=-P(\lambda_{1}+\lambda_{2}+\lambda_{3})$ holds. We will obtain an expression for $\mu^{\prime}(p_{2}^{*})$ by computing the derivatives of both sides of this equation with respect to $p_{2}$.

[TABLE]

Now we calculate the derivative with respect to $p_{2}$ of $b_{1}\,b_{2}-b_{3}$:

[TABLE]

By equating both expressions for the derivative we obtain:

[TABLE]

for $p_{2}^{*}>0$, since $a_{1}>a_{2}$ implies $\gamma(a_{1},a_{2})>0$. We note that $d_{3}<{d_{3}^{max}}$ implies $p_{2}^{*}>0$. since $d_{3}<d_{3}^{\text{max}}$. This completes the proof. ∎

Remark 8**.**

We note that the bifurcation point does not depend on $k$.

Remark 9**.**

It holds $b_{1}(p_{2})=-(\lambda_{1}(p_{2})+\lambda_{2}(p_{2})+\lambda_{3}(p_{2}))$ which equals $-\lambda_{3}(p_{2}^{*})$ if $p_{2}=p_{2}^{*}$. As we noted in the preceding proof, we have $b_{1}=-\lambda_{3}(p_{2}^{*})$ and therefore, we can explicitly compute the eigenvalues of the Jacobian matrix for the case $p_{=}p_{2}^{*}$.

[TABLE]

Similarly, $b_{3}(p_{2})=-\lambda_{1}(p_{2})\lambda_{2}(p_{2})\lambda_{3}(p_{2})$ and $b_{3}(p_{2}^{*})=-\lambda_{3}\left(p_{2}^{*}\right)\omega^{2}\left(p_{2}^{*}\right)$, which implies

[TABLE]

We note that the third eigenvalue $\lambda_{3}$ has a negative real part, since $\gamma(a_{1},a_{2})>0$. Consequently, for $p_{2}=p_{2}^{*}$ the orbits of the system exponentially approach a center manifold. On the center manifold, the orbits look essentially like the ones of the Hopf normal form as shown in ref. [38].

Figure 3 provides a numerical example for the super-critical Hopf bifurcation.

4 Biological implications

In this section we apply the mathematical results to gain insights into the origin of oscillating blood cell counts. We use our model to systematically study which parameters have to deviate from their physiological values to obtain persistent oscillations.

4.1 Numerical studies

Parameters of model (M1)-(M3) have been obtained using a combination of data from literature and patient data. It has been shown that the calibrated model can reproduce blood cell dynamics after bone marrow transplantation and during acute leukemia [58, 60, 59, 57, 61]. The parameters have been estimated as follows [58]:

$a_{1}=0.850$, $p_{1}=0.1/day$, $a_{2}=0.841$, $p_{2}=0.4/day$, $d_{1}=0.0$, $d_{2}=0.0$, $d_{3}=2.7/day$, $k=1.75\cdot 10^{-9}.$

In the following we refer to these parameters as reference values. To systematically check whether variation of one of these parameters can lead to a Hopf bifurcation, we search, where in the parameter space $p_{2}$ passes through the bifurcation point $p_{2}^{*}$. Biologically plausible intervals for the respective parameters are:

[TABLE]

These ranges are motivated as follows. The self-renewal fractions $a_{1}$ and $a_{2}$ correspond to the probabilities that a progeny cell belongs to the same maturation stage as its parent cell, therefore they assume values between zero and one. If $a_{1}\leq 0.5$, the stem cell population declines over time [62], which contradicts biological observations [62, 64]. For this reason we assume $a_{1}>0.5$. A proliferation rate of $1/day$ corresponds to more than one cell division per day, which is the biological upper limit required for genome duplication [50]. The value of $d_{3}$ corresponds to a neurophils half life between $5$ hours and 7 $days$ what is in agreement with measurements in humans [7, 54]. Of note biological studies suggest, that in cyclic neutropenia neutrophils half life does not significantly differ from reference values [23]. As neutrophils are considered as one of the cell types with shortest half life, we assume that $5$ hours is a reasonable lower bound also for the half life of immature cells [29]. For our numerical studies we subdivide the given intervals in 100 equidistant points.

There is biological evidence that oscillating blood cell counts may be related to death of immature cells [22]. For this reason we extend our numerical analysis to the case where stem and progenitor cell death rates, denoted as $d_{1}$ and $d_{2}$, can be positive. In this case we obtain the following expressions for $b_{1}$, $b_{2}$ and $b_{3}$.

[TABLE]

We numerically check at which locations of the parameter space $b_{1}b_{2}-b_{3}$ changes its sign from positive to negative. In the considered parameter space $b_{1}$ and $b_{3}$ are positive.

Remark 10**.**

Global boundedness of solutions also holds for the extended model with positive $d_{1}$ and $d_{2}$. The proof works analogously to the proof of Theorem 1 with $B_{1}:=\frac{1+p_{2}+d_{2}-d_{1}}{2(1-a_{1})}$ and $B_{2}:=\frac{2B_{1}+d_{3}+p_{2}}{2(1-a_{2})}$. The expression for the unique strictly positive steady state together with necessary and sufficient criteria for its existence are provided in ref. [62]. The same applies to the semi-trivial equilibria.

4.2 Results

Numerical simulations show that if only one parameter differs from its reference value, no Hopf-bifurcation is observed. The same holds if two parameters differ from their reference values. The minimal requirement for a Hopf-bifurcation to occur is that three parameters deviate from their respective reference values. All scenarios where Hopf bifurcation can occur are summarized in Table 2, more details are provided in Appendix B and Figure 5. For some parameters significant deviations from the reference values are required to observe Hopf bifurcation. One example for this is the neutrophil apoptosis rate, where a reduction to less than 25% of its original value is required. Taking into account the large variations of the respective cell parameters observed in healthy individuals under immune stimulation (5-10 fold), parameter changes of this order of magnitude can be considered biologically realistic [43].

Our findings are in line with the following results from literature:

The fact that Hopf bifurcations only occur if multiple parameters deviate from their reference values may explain why oscillating blood cell counts are rarely observed. 2. 2.

The death rates of immature cells are increased in most parameter configurations. This is in line with experimental findings [22]. However, increase of immature cell death is not necessary to obtain a Hopf bifurcation (see Constellations 1 and 2 in Table 2). 3. 3.

The model can reproduce neutrophil dynamics in cyclical neutropenia, see Fig. 4, if the feedback parameter $k$ deviates from its reference value. 4. 4.

For certain parameter values the model shows oscillations with a period of several months. Such oscillations have been observed in case of chronic myeloid leukemia [28, 56, 18]. 5. 5.

Many parameter configurations show increased death rates of stem and/or progenitor cells. This finding can explain why in some patients oscillations are induced by chemotherapeutic agents and other drugs inducing cell death [2, 33, 52].

The numerical studies provide the following insights into diseases with oscillating blood cell counts.

In some parameter configurations $a_{1}$ or $p_{1}$ or both are increased. Several studies show that increased stem cell self-renewal and proliferation are linked to malignancy [58, 57, 32, 34, 69]. This may explain why oscillating blood cell counts can be interpreted as a premalignant state [42, 15]. However, increased stem cell self-renewal and proliferation are not necessary for a Hopf-bifurcation to occur. This may explain why several studies do not see a relation between oscillating blood cell counts and malignancies [13]. 2. 2.

According to the considered model Hopf bifurcations can occur for decreased mature cell clearance $d_{3}$ but not for increased $d_{3}$. This is surprising, since for a long time increased mature cell death has been suspected as the origin of oscillating blood cell counts. Experiments, however, have shown that mature cell clearance is not increased in cyclic neutropenia [23]. In some patients with cyclic neutropenia a slight decrease of neutrophil clearance has been reported [23]. 3. 3.

Periodic auto-inflammatory syndromes, such as the PFAPA syndrome are associated with reduced neutrophil apoptosis [37, 66]. These syndromes are characterized by periodic inflammations with increased white blood cell counts [5]. Our results support the idea that decreased neutrophil apoptosis may contribute to the periodicity of the symptoms. 4. 4.

All detected parameter configurations involve changes in progenitor or mature cell parameters compared to healthy controls. It is controversial whether oscillating blood cell counts require functional deficits on the level of stem cells or whether alterations in the progenitor cell compartment are sufficient. Our model suggests that both scenarios can exist. In some cases stem cell properties ($a_{1}$, $p_{1}$, $d_{1}$) differ from their physiologic reference values (Constellations 1-5 and 8, 9), whereas in other cases they do not (Constellations 6 and 7). 5. 5.

Whenever a Hopf bifurcation occurs, the counts of stem, progenitor and mature cells oscillate. This is in line with bone marrow examinations showing oscillations of immature cells [23]. Furthermore it fits to the clinical observation that in many patients not only white but also other blood cell counts oscillate [39, 23]. Of note the observation of oscillating stem cell counts does not imply that stem cell parameters have to deviate from their reference values. Alterations of progenitor cell parameters can be sufficient (Constellations 6 and 7). 6. 6.

For many parameter configurations blood cell counts oscillate within physiological bounds. This may suggest that oscillating blood cell counts not necessarily lead to clinical symptoms. The occurrence of oscillating blood cell counts in healthy patients is so far controversial [51, 12, 13]. 7. 7.

If the value of $k$ remains unchanged, all considered parameter configurations lead to neutrophil counts that are too high to be compatible with those observed in cyclic neutropenia [13]. However, changes in the parameter $k$ lead to the clinical observed low neutrophil counts. Examples are depicted in Figure 4. This observation is in line with the experimental finding that response of immature cells to feedback signals is altered in patients with cyclic neutropenia [24]. To observe low neutrophil counts $k$ has to be higher than its reference value. This means that the effect of cytokines is reduced in patients with cyclic neutropenia, as it has been observed experimentally [24]. 8. 8.

The bifurcation point is independent of the value of the feedback parameter $k$. This means that changes in the feedback signal that affect all cells simultaneously cannot produce oscillations. Consequently, substitution of feedback signals (such as G-CSF) cannot prevent oscillations. This is in line with clinical observations [25]. This finding supports the idea that alterations in the feedback mechanism may be responsible for the low neutrophil counts in cyclic neutropenia but not for the oscillations.

5 Discussion

In this work we prove the occurrence of Hopf bifurcation in a mathematical model of white blood cell formation. The model describes time evolution of a cell population structured by three maturation stages (stem, progenitor and mature cells) that are regulated by a nonlinear feedback mechanism. We show that for appropriate parameter choices a super-critical Hopf bifurcation occurs and a stable limit cycle emerges. This constitutes a major difference to the two-compartment version of the model that distinguishes only between mature and immature cells and has a globally stable steady state [20]. This finding demonstrates that the number of maturation stages can impact dynamical features of the system.

The considered model has been applied to clinical data and shows a good agreement with reality [47, 60, 58, 57, 59, 61]. It predicts that a Hopf bifurcation can occur for biologically plausible parameters and thus provides possible explanations for disease mechanisms that lead to oscillating blood cell counts. Our systematic numerical study suggests that sustained oscillations only occur if multiple parameters deviate from their physiological reference values. This finding is in line with the observation that oscillating blood cell counts are rarely observed.

Biological data and theoretical results have linked the occurrence of oscillating blood cell counts to increased death rates of immature cells. This finding is supported by our model, however, our results suggest that increased immature cell death is not necessary to obtain a stable limit cycle. Alterations of mature cell death rates together with stem or progenitor cell self-renewal and proliferation are also sufficient. The impact of perturbed self-renewal on oscillating blood cell counts has been discussed in the context of a linear model [16].

One of the most common diseases exhibiting periodic oscillations of white blood cells it cyclic neutropenia. Our model suggests that the low cell numbers detected in these patients can only be reproduced if the response of immature cells to feedback signals is reduced. This result is in line with in vitro experimental findings. Interestingly, the occurrence of Hopf-bifurcation is independent of the parameters that describe the feedback signal. This suggests that alterations in the feedback that simultaneously affect all cell types cannot lead to oscillating cell counts. However, in presence of parameter configurations that lead to oscillations, an additional alteration of the feedback signal can impact their amplitude. This result leads to the hypothesis that the experimentally detected reduced response of immature cells to cytokines in patients with cyclic neutropenia is the pathogenic mechanism leading to low cell counts; however it does not contribute to the occurrence of oscillations. This hypothesis is further supported by clinical data showing that administration of growth factors such as G-CSF increases white blood cell counts but does not lead to cessation of the periodic oscillations [49, 25].

The finding that multiple different parameter configurations can result in oscillating blood cell counts may explain parts of the heterogeneity among cyclic neutropenia patients. It furthermore suggests that different detected mutations [45, 19, 1, 4, 71, 8] may lead to different pathogenic mechanisms that result in similar symptoms.

In summary, we have proven the occurrence of a Hopf bifurcation in a non-linear three-compartment model of white blood cell formation. The Hopf bifurcation is a unique feature of the three-compartment setting and does not occur in the 2-compartment version of the model. We identify biologically plausible parameter sets that lead to a stable limit cycle and relate them to clinical and experimental findings. This quantitative approach can help to understand the pathogenic mechanisms and the clinical heterogeneity of different diseases that lead to oscillating blood cell counts.

Appendix A Supplementary calculations to the proof of Theorem 4

In the following we provide the calculations leading to equation (3).

[TABLE]

Appendix B Parameter configurations leading to Hopf bifurcation

[TABLE]

[TABLE]

Biologically plausible parameter regions where the Hopf-bifurcation exists are visualized in Figure 5. The reported values for $a_{1}$ and $a_{2}$ correspond to the self-renewal fraction in presence of maximal stimulation. The self-renewal fraction at time $t$ is given by $a_{1}s(t)$ and $a_{2}s(t)$ with $s(t)<1$.

Acknowledgments

This work was supported by research funding from the German Research Foundation DFG (Collaborative Research Center SFB 873, Maintenance and Differentiation of Stem Cells in Development and Disease, subproject B08).