Stem Cell Proliferation and Quiescence—Two Sides of the Same Coin

The kinetics of label uptake and dilution in dividing stem cells, e.g., using Bromodeoxyuridine (BrdU) as a labeling substance, are a common way to assess the cellular turnover of all hematopoietic stem cells (HSCs) in vivo. The assumption that HSCs form a homogeneous population of cells which regularly undergo cell division has recently been challenged by new experimental results. For a consistent functional explanation of heterogeneity among HSCs, we propose a concept in which stem cells flexibly and reversibly adapt their cycling state according to systemic needs. Applying a mathematical model analysis, we demonstrate that different experimentally observed label dilution kinetics are consistently explained by the proposed model. The dynamically stabilized equilibrium between quiescent and activated cells leads to a biphasic label dilution kinetic in which an initial and pronounced decline of label retaining cells is attributed to faster turnover of activated cells, whereas a secondary, decelerated decline results from the slow turnover of quiescent cells. These results, which support our previous model prediction of a reversible activation/deactivation of HSCs, are also consistent with recent findings that use GFP-conjugated histones as a label instead of BrdU. Based on our findings we interpret HSC organization as an adaptive and regulated process in which the slow activation of quiescent cells and their possible return into quiescence after division are sufficient to explain the simultaneous occurrence of self-renewal and differentiation. Furthermore, we suggest an experimental strategy which is suited to demonstrate that the repopulation ability among the population of label retaining cells changes during the course of dilution.


Simulation algorithm and model equations
The described model of HSC organization is mathematically represented as a single-cell based, stochastic process. I.e., the development of each individual cell in the system is simulated according to a set of defined rules including stochastic decisions. These rules are applied at discrete time steps (∆t =1 hour) to simultaneously update the status of all model cells.
The status of each model cell is given by its cell specific affinity a, its membership to a signaling context m ∈ {A, Ω}, and its position in the cell cycle c. The cell cycle for an activated cell in signaling context Ω is modeled as a sequence of G1 phase, S phase, and G2/M phase, in which the duration of latter phases (τ S and τ G2/M ) is fixed. In contrast to former versions of the model with fixed G1 phase [1,2], this duration τ G1 is now randomly chosen from an exponential distribution with meanτ G1 upon division of the parental cell. Cells changing into signaling context A exit from G1 and enter a cell cycle inactive phase (G0) until they are activated again and transit into signaling context Ω. To realize an update step, the actual total number of cells with a > a min in signaling context A (N A (t)) and Ω (N Ω (t)) is determined. Based on these numbers, the status of each model stem cell is updated as follows: (1) If the cell resides in signaling context A, it changes to signaling context Ω or stays in A with probabilities ω and 1 − ω, respectively. If it stays in A, its cell specific affinity a is increased by the factor r (regeneration coefficient), until a has reached its maximum value a max = 1. If the cell changes to signaling context Ω, its position in the cell cycle will be set to the begin of S-phase.
(2) If the cell resides in signaling context Ω, it changes to context A or stays in Ω with probabilities α and 1 − α, respectively. Herein, a change to signaling context A is only possible in G1-phase of the cell cycle for cells with a > a min . If the cell stays in signaling context Ω, a is decreased by the factor 1/d (with regeneration coefficient d) and the cell cycle position is incremented. In case of cell cycle completion (i.e., c(t) = τ G1 + τ S + τ G2/M ), c(t) is set to 0 and a new identical cell is generated (cell division). Individual durations of G1 phase for both daughter cells are randomly chosen from an exponential distribution. If a has reached the minimal value a min , the cells cannot reenter into signaling context A and proceed development in Ω.
The transition probabilities α and ω depend on the actual affinity of the cell a(t), on the fixed parameters a min and a max , and on the transition characteristics f α and f ω , respectively: The transition characteristics f α and f ω itself, depend on the total number of cells (N A ,N Ω ) in the respective target context. They are modeled by a general class of sigmoid functions: The parameters ν 1 , ν 2 , ν 3 , and ν 4 determine the shape of f α/ω .Ñ A/Ω is a scaling factor for N A/Ω . It is possible to uniquely determine ν 1 , ν 2 , ν 3 , and ν 4 by the more intuitive values f α/ω (0), Applied parameter choices of these functions are given in Supplementary Table 1. Within the model a representative population of HSCs includes all cells in signaling context A as well as a proportion of cells in signaling context Ω. The coverage of this HSC pool is indicated by the blue box in Figure 3 of the main document. The position of the right sided boundary of this box is defined in terms of transit time t trans that a cell remains within this box after having passed the threshold value a min . Motivated by the observation that even the most sophisticated protocols for the purification of HSCs yield at maximum 50 % of repopulating stem cells [3], the simulated HSC population intentionally includes a significant number of cells with restricted repopulation potential, i.e. cells a < a min . These cells can by definition not enter signaling context A, and are thus not able to contribute to system repopulation.
The estimated turnover times are averages over the time between two division events of an individual cell. For activated cells (i.e. cells are only in signaling context Ω for their whole time of existance) the individual time beween two division corresponds to the sum τ G1 + τ S + τ G2/M . Cells with intermediate residence in signaling context A are denoted as quiescent cells. This temporary state of inactivation potentially prolongs the time between the parental division and the next division event, leading to the increased average turnover of these cells.
The mathematical model is implemented in C++ and has been tested for UNIX-derived operating systems. The source code (including parameter files for the described scenarios) can be obtained from the authors.

Simulation of dilution kinetics
For the simulation of label dilution kinetics each cell is additionally characterized by an individual label content b(t). As the dynamics of label uptake during BrdU administration are potentially biased by the cytotoxic effects of BrdU itself [4] and because we currently have no quantitative information about this process, the model does not consider this process. Instead, upon start of label dilution a certain fraction F 0 of all HSCs receives an initial label b(t = 0) = b 0 (individual cells are chosen randomly within this population). Upon division the daughter cells receive half the parental label content. Cells become undetectable if the individual value b(t) falls below a detection threshold b t .
Without loss of generality we have chosen b 0 = 0.5 as the initial fraction of label cells and b 0 = 0 for the non-labeled cells. The number of divisions N that are necessary to label dilution below the level of detectibility is regulated by adjusting the detection threshold b t . Parameters are provided in Supplementary Table 2.
As outlined in the main document, we additionally studied the more general case in which the intial label is not fixed to b 0 = 0.5 but distributed according to a given probability density. Details are provided in Section 4 below.

Simulation of competitive retransplantations.
For the simulation of the competitive retransplantation experiments, without loss of generality, a fraction of F 0 = 71% of all HSCs is initially labeled. At different time points during label dilution, the HSC population is seperated in L+ (b(t) > b t ) and Lsubpopulations (b(t) < b t ). 20 cells from each pool are retransplanted with 20 competitor cells chosen randomly from the HSC population of a corresponding model system for which the labeling routine had not been applied. Over time the 40 transplanted cells repopulate an initially empty model system and reestablish a homeostatic situation in which a fraction of cells is donor derived (originating either from the L+ or L-population) and the remaining fraction of cells originates from the pool of competitor cells. The fraction of donor derived cells represents the experimentally observable engraftment level. Engraftment levels at each time point in Figure 5 of the main document are estimated from 100 indivdual transplantation studies. The variability of the engraftment levels (indicated by the standard deviation) results from two aspects: First, the cells are chosen randomly from the selected L+ and L-populations. Second, the engraftment is modeled as regulated, but intrisically stochastic process. Very likely, both aspects do also account for the observed variabilty in the corresponding experimental results. To demonstrate the influence of the initial labeling routine an additional scenario is studied in which only HSCs with a > 0.9 are label with b 0 = 0.5. All other cells remain b 0 = 0. Due to the construction of the model, in which the cells in signaling context A accumulate at a = 1, the particular labeled population is almost completely composed of quiescent cells. Instead of assigning the same initial label content to all cells we have additionally studied the case that b 0 is distributed around a mean value µ =b 0 = 0.5. As we have defined b in the range [0, 1], we use Beta distributed random variables for this approach. The Beta distribution with parameters p > 0 and q > 0 has density where Γ is the gamma function. The distribution is symmetric aroundb 0 = 0.5 for p = q.
Densities for different values of p = q are shown in Figure 1A. Label dilution kinetics are obtained for the corresponding assignment of the initial label b 0 according to these probability densities. The detection threshold is set to b t = 0.02. Results are shown in Figure 1B. Even for a highly heterogenous spread of the initial label content b 0 the dilution kinetics are closely similar to the case where all cells are initially labeled with b 0 = 0.5 (corresponding to p = q = 10 6 ).   (*) For Figure 5B, F 0 = 100% of the cells are labeled in a restricted region.