Because PCML is believed to arise in the stem cell compartment in the bone marrow, here we consider a stem cell model with a proliferating and a resting G₀ phase [5]. A cell entering the proliferating phase can either die by apoptosis at a rate γ or divide a fixed time τ after entry (point of cytokinesis). Immediately after division, the two newly born cells, the daughter cells, enter the resting phase. They can either remain in this resting phase, or exit from this phase by differentiating into one of the committed cell lines (red cells, white cells, platelets) at a rate δ, or exit and re-enter the proliferating phase at a rate β.

If we denote the density of proliferating cells at time t by P(t), and the density of resting cells at time t by N(t), the conservation equations are given by (see [6–8] for further details)

Note that the solution to Eq. (2) is independent of the behavior of the solution to Eq. (1) but the converse is not true. Knowing the behavior of the solutions of Eq. (2) it is easy to obtain the behavior of the solutions of Eq. (1). Thus we concentrate our study on Eq. (2).

A solution to Eq. (2) is a continuous function $N:[-\tau ,+\infty )\to {\mathrm{R}}_{+}$, for all t>0. We denote the continuous function $\varphi :[-\tau ,0)\to {\mathrm{R}}_{+}$, φ(t)=N(t) for all t∈[−τ,0], as the initial condition for N. Using the method of steps, it is easy to prove (see [9]) that for every ϕ∈C([−τ,0]), there is a unique solution to Eq. (2), where C([−τ,0]) is the space of continuous functions on [−τ,0].

The purpose of this paper is to apply the technique used in [11] to this HSC model in order to study the effect of various parameter changes on the dynamics and to thereby hopefully obtain some insight into the HSC origin of PCML.

Our objective is to derive some understanding of the possible mechanisms giving rise to the long-period oscillations seen in PCML by looking at the influence of each parameter in our HSC model on the stability of the steady states and the period of the solutions that result once stability is lost.

In this section, we briefly summarize the stability results analytically proved by Mackey and Pujo-Menjouet in [4]. Before analyzing the period of the solutions, it is necessary to give the regions where stability occurs and the ones where the solutions become unstable through the Hopf bifurcation. To that end, we determine the steady states of Eq. (2), and then simplify the problem using dimensionless variables and determine the regions of stability.

3.1 Steady states

3.2 Local stability

3.2.1 Influence of the sensitivity n

When n varies from 1 to 10 (Fig. 1), we observe that oscillations appear (for n about 10), then increase in period and amplitude up to a ‘saturation’ limit (Fig. 2) due to the fact that the Hill function β given by (3) approaches the Heaviside step function when n→∞, as explained in Section 4. Consequently, to increase the oscillation period and amplitude in a larger way, the other parameters involved in the model play a crucial role.

Note that for the next simulations in this section, we choose n=12 not only because we are focused on long period oscillations but also because in Section 4 we are interested in the linear case where n tends to infinity. Due to the saturation effect, it is then sufficient to take n=12 as a good parameter value for comparison.

3.2.2 Influence of the loss rates δ and γ

Consider first the loss rate δ from the resting phase. One can observe five different stages as δ increases.

• 1. When δ=0, the solution converges rapidly to the nontrivial steady state (Fig. 3). Then as δ increases, small oscillations start to appear.
• 2. If we increase δ slightly, one notices an abrupt rise in the amplitude and period of the oscillations occurring in a narrow interval of change in δ. We denote by δ_g the value at which this transition occurs where the amplitude and period are at their maximum. From our calculations (shown in Fig. 3), this value is about δ_g≈0.46567 day⁻¹.
• 3. Let δ become larger than δ_g. Then, the amplitude and period decrease slowly as δ continues to increase.
• 4. This phenomenon persists until the solution starts to exhibit a succession of small-amplitude oscillations just before a ‘big amplitude jump’. These big jumps are, however, smaller than the amplitudes seen in the previous case, and continue to decrease as δ increases. This behavior is shown in Fig. 4 when δ is close to 0.5. Due to this small oscillation effect the period again increases. We denote these solutions as having a ‘tail-jump’ shape.
• 5. When a too large value of δ is reached, the solutions converge to the trivial steady state x≡0 due to a high loss rate in the resting phase leading to the extinction of the cell population.

This behavior is subtle. It is rather easy to observe the five different stages but it is more difficult to explain them by considering the nonlinear equation (7). However, if we transform this nonlinear equation (7) into a simpler linear system and if we consider the explicit expression of the Hopf period given by T_H in (12), it is then possible to get a better understanding of what happens. This will be done in Section 4.

Consider next the influence of the apoptosis rate γ from the proliferating phase. Numerically, we found that γ influences the behavior of the solutions in the much same way as δ. Namely, as γ increases, the amplitude and period of the oscillations increase and again there is a transition point at which there is a sudden jump in both, and then a slight decrease in the amplitude and period. We have estimated the value at which this happens to be γ≈0.2 day⁻¹ (Fig. 5). A succession of ‘tail-jump’-like shapes then occurs until finally the solution approaches the trivial steady state x≡0.

3.2.3 Influence of the cell cycle time τ

The influence of τ on the behavior of the cell population corresponds to the predictions based on the local stability analysis in Section 3.2. Thus, when τ is small and close to 0 (so the cell cycle time is short and the cells rapidly proliferate), the solution is stable and rapidly converges to the nontrivial steady state. As τ increases, it reaches the point τ=τ_crit, the nontrivial steady state becomes unstable and there is a Hopf bifurcation to limit cycle behavior. Further increases lead to complex solution behavior that are still under investigation and will be the object of future work. Finally, if τ is increased past the value τ_max, the cell cycle time becomes so large that the population is extinguished and the nontrivial steady state ceases to exist. In this case the single remaining steady state is the trivial one and the cell population dies out. This is clearly due to the parameter $\kappa =2{\mathrm{e}}^{-\gamma \tau }$, because as τ increases, κ tends toward 0 and so the right-hand side of Eq. (7) is ruled only by the negative term. The solution reaches then the trivial steady state. All these behaviors are shown in Fig. 6.

3.2.4 Influence of the maximal reentry rate β₀

The maximal reentry rate β₀ plays an important role in determining the amplitude and period of the periodic solutions. As shown in Fig. 7, when β₀ is sufficiently large, the more β₀ increases, the larger the amplitude and the period of the oscillations are. However, when β₀ is close to 0, the behavior is more complicated. Indeed, in this case, one can observe from Eq. (7) that the behavior of the solutions is ruled by the two loss death rates γ and δ, and it is not so straightforward to understand what is happening by considering the full equation (7). To give a more accurate interpretation of this complex behavior we transform the nonlinear equation (7) into a piecewise linear one, by considering the extreme case of n→∞. This is the subject of the next Section 4.

In [4], the explicit solution to Eq. (2) was computed for n→∞ to obtain the oscillation period and amplitude analytically. In this section we are not interested in the analysis of [4]. Our purpose is to go further in the interpretations of the numerical results presented in Section 3.2 and explain the complicated influence of the loss rates observed in the previous section.

Before showing the numerical results, note that when n→∞, the feedback function β can be approximated by the Heaviside step function β(x_τ)=β₀[1−H(x_τ−θ)], where:

Since the numerical simulations in Section 3.2 were carried out for n=12, when n tends to infinity one observes only slight differences in the numerical simulation for the same set of parameters. As explained in Section 3.2.1, this is due to a ‘saturation’ limit influenced by the Hill function β. For this reason, we do not show the effect of τ in this section, as well as the effect of large values of β₀. However, the transformed problem is useful in the sense that, thanks to its simplicity, we are able to interpret the different behaviors induced by the loss rates.

Indeed, concerning the loss rate δ from the resting phase, we can explain the five different stages shown in the previous section.

• (i) If δ is very low, i.e., δ≅0, the right-hand side of the first and second equations of the system (14) are close to 0 and the solution x does not oscillate, or does so only slightly, and rapidly reaches the non-trivial steady state.
• (ii) For larger values of δ, the amplitude and period in oscillations reach their maximum very rapidly; these maximum values are due to the fact that the difference between −αx and Γx_τ (or −δx and Γx_τ) in the right-hand side of the third (or fourth) equation of (14) is at its optimum. In other words, all the conditions are met to obtain the maximal amplitude and period oscillations. This value of δ is difficult to find analytically, but is approximated by δ≈0.46567 day⁻¹ in Section 3.2.2. The rapid evolution in such a short interval of δ is not well understood analytically but is confirmed in [4]. Indeed, the authors were able to give an explicit but quite complicated form of the amplitude and to show numerically that in this interval the increase of the amplitude to its maximum is remarkably fast. However, because the explicit form found in [4] is complicated we were not able to provide any simple interpretation.
• 3. After this abrupt increase, the amplitude and period decrease slowly due to a smaller and positive difference between −αx or −δx and Γx_τ. We can see this phenomenon if we study the Hopf period given by T_H in (12). It is clear (even analytically) that T_H decreases as δ increases (up to a critical value of δ as explained in the next case). This is shown in Fig. 8. Because the influence of τ and γ on the Hopf period T_H is similar, we do not show the effect of the increase of γ in this figure. Thus, with our chosen parameters, to obtain maximal oscillations, δ has to be very close to an ‘optimal’ value that is neither too small nor too large and τ and γ have to also stay at a low level.
• 4. When δ is too large, the difference between −αx or −δx and Γx_τ decreases, but is still positive. Then the solution x as well as x_τ have some difficulty in going to 1. This means that some very small oscillations occur due to the fact that the difference −αx or −αx−Γx_τ (or −δx−Γx_τ) is positive, then negative following the different values taken by x and x_τ. This phenomenon leads to very long period oscillations with very small amplitudes and the ‘tail-jump’-like shape. This period tends to infinity as δ tends to a value δ_crit. In Fig. 8, if we increase τ with τ=1 day and τ=3 days), one can observe a shorter interval where oscillations occur, i.e. δ_crit becomes smaller as τ gets bigger. The cell population tends to be extinguished for lower values of δ when τ is larger.
• 5. Then, when δ>δ_crit, δ becomes too large, so the difference −αx−Γx_τ (or −δx−Γx_τ) remains negative and the oscillations vanish with the solution x converging to the trivial solution x≡0.

A slightly similar interpretation can be given for the death rate γ related to the proliferating phase, and so we do not go into the details. However, the length of the Hopf period T_H is more difficult to study analytically when γ varies. If we analyze it numerically one can note that the period length decreases as γ increases up to a certain value. Then due to the effect of the ‘tail-jump’ oscillations as explained above we can see that the period length will increase up to infinity when γ tends to a value we denote by γ_crit. For γ>γ_crit, the periodic solutions disappear, and the solution x converges to the trivial solution x≡0, which means that the cell population is extinguished.

We give two numerical illustrations of the Hopf period T_H when γ varies. One figure is with two different values of δ and the other is with two different values of τ. In the first figure, we increase δ (Fig. 9 with δ=0.05 day⁻¹ and δ=0.1 day⁻¹), which implies an increase in the period length. In the second figure, which is similar to Fig. 8, if we increase τ (Fig. 10 with τ=1 day and τ=3 day), we can see a shorter interval where oscillations occur, i.e. γ_crit becomes smaller as τ gets bigger. The cell population tends to be extinguished for lower values of γ when τ is larger.

The influence of the proliferating phase duration τ has been fully explained analytically in Section 3.2 and illustrated in Section 3.2.3. The transformed problem does not bring further insight.

However, the system (14) is very useful to interpret the influence of the maximal reintroduction rate β₀. As noted in Section 3.2.4, when β₀ is sufficiently large, then α∼β₀ and Γ∼2β₀ and the more β₀ increases, the more the amplitude and period of the oscillations increase. However, when β₀ decreases and is close to 0, α∼δ and $\Gamma \sim 2{\mathrm{e}}^{-\gamma \tau }$. In that case, the population dynamics become more complicated and depend on δ and γ. The behavior is then under the influence of these loss rates and so it is similar as the one presented in Section 3.2.2 and above (see Figs. 11 and 12).

The aim of this note was to show how short cell cycle durations in the HSC (order of one to two days) can give birth to long period oscillations (order of 40 to 80 days) as observed in the disease PCML. We also wanted to see how these oscillations could occur and in which way each of the five HSC parameters (n, β₀, τ, γ, and δ) influence the oscillation amplitude and period.

Our analysis of Eq. (7) showed that instability occurs in a very small interval [0,τ_n] of τ. It is in this interval that our attention has been focused, and we showed analytically that periodic solutions can arise in this interval. We performed numerical investigations to see if large amplitude and long period oscillations could exist and under which circumstances. We observed three different groups in the parameters.

• • The first one consists of n and β₀ and corresponds to the influence of the Hill function β(N). Changes in n and β₀ can produce oscillations, but due to a saturation effect increases in n do not imply long period and large amplitude. However, if β₀ is large, long period and large amplitude oscillations are possible.
• • The second group consists of the loss rates γ and δ. It appears that under changes in these parameters, it is possible to obtain large amplitude and long period oscillations. However, the evolution of these oscillations is very sensitive to slight changes in both γ and δ. We showed that the behavior can be totally different in small intervals of change of these parameters, and thus the HSC population seems to be highly sensitive to changes in its loss rates (mortality or differentiation).
• • The third group, consisting of τ, is interesting in the sense that this is the main focus of our study. Indeed, because our goal was to show that long period oscillations could occur when the cell cycle duration τ is small, our attention was concentrated only in small variations of τ (no more than 3 days). Like the loss rate, the behavior of the solution showed a high sensitivity to changes in τ, and we observed that the interval of τ over which we can obtain long-period oscillations with large amplitude is quite narrow, i.e. τ can be neither too small nor too large.

All these five parameters can be combined together to give rise to the maximal values of amplitude and period oscillations, but this optimum point is difficult to find due to the strong influence of each parameter.