HSC Niche Dynamics in Regeneration, Pre-malignancy, and Cancer: Insights From Mathematical Modeling

Abstract The hematopoietic stem cell (HSC) niche is a crucial driver of regeneration and malignancy. Its interaction with hematopoietic and malignant stem cells is highly complex and direct experimental observations are challenging. We here develop a mathematical model which helps relate processes in the niche to measurable changes of stem and non-stem cell counts. HSC attached to the niche are assumed to be quiescent. After detachment HSC become activated and divide or differentiate. To maintain their stemness, the progeny originating from division must reattach to the niche. We use mouse data from literature to parametrize the model. By combining mathematical analysis and computer simulations, we systematically investigate the impact of stem cell proliferation, differentiation, niche attachment, and detachment on clinically relevant scenarios. These include bone marrow transplantation, clonal competition, and eradication of malignant cells. According to our model, sampling of blood or bulk marrow provides only limited information about cellular interactions in the niche and the clonal composition of the stem cell population. Furthermore, we investigate how interference with processes in the stem cell niche could help to increase the effect of low-dose chemotherapy or to improve the homing of genetically engineered cells.

In this section we develop a mathematical model describing the interaction of HSC with the bone marrow niche. The model is based on the following biological processes.
1. HSC require interaction with the bone marrow microenvironment to maintain their stemness [20,41,49]. HSC attached to the niche are in a quiescent state, i.e., they do not divide [21]. This assumption is supported by the finding that niche-specific molecules such as osteopontin or CXCR4 inhibit entry into the cell cycle [28,29,19,39,21,9].
2. HSC detach from the niche at a rate u H . Upon detachment HSC become activated. This asumption is in line with experiments showing that an interruption of CXCR4 or osteopontin related queues leads to cell cycle entry [28,29,39]. This transition can be written as where N H denotes a niche space that is occupied by an HSC, A H denotes an activated HSC and N E a vacant niche space. For simplicity, we assume that HSC activation occurs immediately after detachment.

Activated HSC can divide and differentiate into progenitors at a rate d A H . This transition is written as
corresponding to the loss of a stem cell, D H denotes a progenitor. 4. Activated HSC can divide at a rate r H . We assume that the offspring originating from division are inactive HSC, i.e., they have to reattach to the niche before they can divide again. Inactive HSC either differentiate at a rate d I H or reattach to the niche at a rate b H . We assume that inactive cells cannot be activated outside the niche. This implies that HSC differentiate in absence of the stem cell niche, as it is observed in cell culture experiments [20,41,49]. For the sake of simplicity we assume that activated HSC become inactive after performing one division. It is straightforward to extend the model to allow an arbitrary finite number of divisions before becoming inactive. On the population level it makes no difference whether we allow asymmetric stem cell divisions or not. In section S1. 3 we formally show that the model introduced here based on symmetric HSC divisions is equivalent to a model including both symmetric and asymmetric HSC divisions. The described processes can be written as where I H denotes an inactive HSC. 5. The progenitors, D H , are committed to differentiation and cannot rebind to the niche. Instead of explicitly modeling progenitor dynamics, we assume that the mature cell output is correlated with progenitor production. This is reasonable since progenitors can only perform a limited number of divisions and thus each progenitor gives rise to an approximately constant number of mature cells. 6. For the sake of simplicity we assume that the bone marrow niche has a constant finite capacity, denoted by K [48].
We neglect re-attachment of stem cells to the niche without division. Simulations and mathematical analysis suggest that this simplification does not impact model dynamics. Similarly, allowing more than one division of activated cells before inactivation has no impact dynamics in the stem cell niche [31]. The biological processes underlying the model are summarized in Figure 1 of the main text and in Figure S1. We assume that the cells are well-mixed. This view is supported by magnetic resonance imaging of acute myeloid leukemia patients who mostly show a non-focal diffuse marrow infiltration and by post-mortem biopsy studies of chronic myeloproliferative diseases [40,12]. This fits well to the concept that malignant cells spread through the blood stream to various bone marrow sites. Due to the high number of cells in the hematopoietic system [1,44] we describe population dynamics using ordinary differential equations (ODEs). We denote by A H the abundance (measured e.g., as cells per kg of body weight) of activated HSC, by I H that of inactive HSC, by N H that of occupied niche spaces and by N E that of empty niche spaces. The influx to the progenitor compartment per unit of time is denoted i D H . Based on the processes described above we obtain the following system of ODEs For brevity, we refer to the sum of I H and A H as 'free HSC'.

S1.2 Model for competition of multiple stem cell clones
In this section we extend the model to account for clonal competition.
Clinically relevant conditions such as clonal hematopoiesis, AML or MPN are related to the presence of multiple stem cell clones [44,18,6,32,26,36,38]. In this section we extend the model to take into account such scenarios. For the sake of simplicity we consider two stem cell populations with different properties. In case of clonal hematopoiesis these two stem cell populations correspond to two HSC clones with slightly different properties. In case of malignancy the two clones correspond to healthy (HSC) and malignant or leukemic stem cells (LSC). In a setting after gene therapy the cell populations correspond to host-derived wild-type and transplant-derived transduced cells. Extension of the model to more than two clones is straightforward but will not be considered in this work. In the model both stem cell types reside in the same niche and compete for niche spaces, similar to the work of Wang et al. [44], Stiehl et al. [38] and Ashcroft et al. [3].
The model is visualized in Figure S1. For the sake of convenience we denote one cell type by H and the other by L.
Homogeneous stem cell population. HSC detachment from the niche occurs at rate uH , division of activated stem cells followed by differentiation at rate dA H , differentiation of inactive stem cells at rate dI H , division of active stem cells at rate rH and attachment of inactive stem cells to the niche at rate bH . Counts of niche-bound cells are denoted NH , counts of inactive cells as IH and counts of activated stem cells as AH .
Competition of two stem cell clones. Subscript H denotes the first (e.g., healthy) clone, while L denotes the second (e.g., leukemic) clone. Figure S1: Schematic representation of the model. (a) Homogeneous stem cell population, corresponding to a healthy scenario. (b) Competition of two stem cell populations, corresponding to AML, MPN, clonal hematopoiesis or following transplantation of donor cells or genetically engineered cells. In panel (b), we use the abbreviations introduced in panel (a). Empty niches are denoted by N E . The index H refers to healthy cells, the index L to donor-related, mutated or malignant cells.
In analogy to the healthy system described in section S1.1, both clones are subjected to the following biological processes.
1. Stem cells attached to the niche are quiescent and upon detachment they become activated. This assumption is motivated by the observation that malignant cells also enter cell cycle if they are mobilized from the niche or if niche-related signaling is disrupted [32,34,33,46]. HSC and LSC detach from the niche at rates u H and u L respectively: where N L denotes a niche space that is occupied by a LSC and A L denotes an activated LSC.
2. Activated stem cells divide and differentiate at rates d A H and d A L . This can be written as corresponding to the loss of stem cells. Here, D L denotes a LSC-derived progenitor.
3. Activated cells divide into two inactive cells at rates r H and r L respectively. Inactive cells can either differentiate (at rates d I H and d I L ) or reattach to the niche (at rates b H and b L ): where I L denotes an inactive LSC.
Denoting by A L the abundance (measured e.g., as cells per kg of body weight) of activated LSC, by I L that of inactive LSC, by D L that of LSC-derived progenitors and by N L that of niches occupied by LSC we obtain the following system of ODEs: In the model both stem cell clones compete for the same niche spaces. Depending on the cell properties (proliferation rate, attachment rate, detachment rate etc.) one clone can out-compete the other clone from the niche.

S1.3 Equivalence of models with and without asymmetric HSC divisions
In this section we show that it has no impact on system dynamics whether HSC can divide asymmetrically or not.
In the model derivation we have assumed that HSC always divide symmetrically (point 4 in section S1.1). In the following we show that up to a reparametrization, this is equivalent to a model allowing symmetric and asymmetric HSC divisions.
In a scenario where asymmetric HSC divisions are possible, the processes where S r is the rate of symmetric self-renewing divisions, S d the rate of symmetric divisions resulting in differentiated cells, and a is the rate of asymmetric divisions. The full mathematical model with explicit asymmetric divisions is then Defining d A H = S d + r H − S r = S d + a 2 yields exactly the system of equations as (S6). This demonstrates the equivalence of the model with and the model without asymmetric divisions.
Later, we will see that r H > d A H is a necessary condition for the existence of a positive homeostatic equilibrium (see equations (S30) and (S26)). We note that the inequality r H > d A H is equivalent to S r > S d , since d A H = S d + r H − S r . Therefore, this condition can be interpreted such that the rate of symmetric selfrenewing divisions, S r , must exceed the rate of symmetric differentiating divisions, S d , regardless of the absence or presence of asymmetric divisions.

S1.4 Remarks on model design
In this section we motivate why our model distinguishes between activated and inactive stem cells.
We aim to derive a model that is in line with the following experimental and clinical observations: (i) Stem cells in the niche are quiescent i.e., they do not divide (ii) Stem cells outside the niche can divide and/or differentiate (iii) The system has a homeostatic equilibrium, i.e., an equilibrium where cell counts are positive (iv) The system can recapitulate bone marrow transplantation, i.e., if a small number of HSC is added to empty niches the HSCs expand.
(v) In absence of a stem cell niche the number of stem cells declines over time.
In the following we consider different candidate models which are summarized in Figure S2.

Candidate Model 1
Let us consider a niche of fixed capacity K. We distinguish between non-dividing stem cells in the niche (N ) and active stem cells outside the niche (A). Stem cells in the niche die/differentiate at rate µ and detach from the niche at a rate u. Stem cells outside the niche are activated to divide at rate r and then divide into two offspring. Cells outside the niche are removed at rate d due to death or differentiation. The attachment of cells is proportional to the product of free niche spaces K − N and cells outside the niche A. The proportionality factor is referred to as b. The model is illustrated in Figure S2 (A). These assumptions lead to the following equations.
We linearize around the trivial equilibrium A = 0, N = 0 and obtain the following characteristic polynomial If d − (2r − 1) > 0 the characteristic polynomial has two roots with negative real parts, i.e., the trivial steady state is not unstable. This contradicts (iv). For simplicity, we neglect the case d − (2r − 1) = 0, which is improbable to be biologically realized since it is prone to be disturbed by noise. Therefore, we consider d − (2r − 1) < 0. In absence of the niche, the dynamics of stem cells is governed by which implies exponential growth of HSC. This contradicts to (v). Therefore, this model cannot satisfy (i)-(v).

Candidate Model 2:
We consider a more general form of candidate model 1. We assume that proliferation and differentiation of cells outside the niche and clearance of stem cells in the niche are regulated by a feedback signal that depends on a weighted sum of A and N . The weights are denoted by α, β ≥ 0. The model is illustrated in Figure S2 (B). This results in the following equations.
The term f (αA + βN )A denotes the difference of the influx of unbound stem cells due to proliferation and the outflux due to differentiation and death. The term µ(αA + βN ) ≥ 0 describes the loss of nichebound stem cells due to differentiation and death. Both are assumed to be continuously differentiable. For f (αA + βN ) = (2r − 1) − d, µ(αA + βN ) = µ we obtain candidate model 1. This model covers regulations of proliferation and differentiation by Hill-type functions such as in [24,11]. For α = β = 1 the feedback depends on the total stem cell count, for α = 0, β = 1 the feedback depends on the number of stem cells in the niche and for α = 1, β = 0 it depends on stem cells outside the niche.
We linearize around the trivial equilibrium and obtain If f (0) < 0, the characteristic polynomial has two roots with negative real parts and the trivial equilibrium is not unstable. This contradicts to (iv). For simplicity, we neglect the case f (0) = 0, which is improbable to be biologically realized since it is prone to be disturbed by noise. Therefore, assume f (0) > 0. Then we obtain the following dynamics in absence of the niche (N = K = 0): The linearization of the right hand-side around A = 0 is f (0), which is positive. This implies that the stem cell count cannot converge to zero in absence of the niche, which violates (v). Therefore, this approach cannot satisfy (i)-(v).

Candidate Model 3:
We consider a scenario where the number of niches is not constant but depends on the stem cell counts. Denote by N E the number of empty niches, by N the number of stem cells residing in niches, i.e., the total number of niches is N E + N . Denote by A the number of stem cells detached from the niche.
Denote the production of empty niches by f (N E , N, A) ≥ 0. It might be positive even if N E = 0, since niche cells may be produced by non-niche cells. Denote the decay of empty niches by µ(N E , N, A) ≥ 0. The attachment of stem cells to the niche is proportional to the product of free niche spaces N E and cells outside the niche A. The proportionality factor is referred to as b. The detachment rate is denoted by u(N E , N, A) ≥ 0. Death or differentiation of HSC in the niche occurs at rate g(N E , N, A) ≥ 0. Death of niche cells which are attached to HSCs occurs at rate h(N E , N, A) ≥ 0. Death or differentiation of stem cells outside the niche occurs at rate d(N E , N, A) > 0. We assume that the stem cell proliferation rate r may be regulated by a feedback depending on the weighted sum of stem cells and set r ≡ r(αA + βN ) for α, β ≥ 0. The model is illustrated in Figure S2 (C). This results in the following system of equations. All coefficients are assumed to be continuously differentiable functions.
For the total number of stem cells it holds If it holds r(0) < 0, a small amount of transplanted cells will not expand. We therefore assume r(0) > 0. For simplicity, we neglect the case r(0) = 0, which is improbable to be biologically realized since it is prone to be disturbed by noise. In absence of the niche (N = N E = 0), stem cell dynamics are governed by which is unstable in A = 0, which means that stem cells can survive without the niche. This is a contradiction.

Conclusion
Since these straightforward models violate biological assumptions, we introduce the inactive stem cell state.
The reason why we introduce the inactive state is the observation that HSCs cannot expand or be maintained in absence of the stem cell niche [20,41,49]. In a model without the inactive state, where stem cells are either attached to the niche (and quiescent) or detached from the niche (and active) we would observe a self-maintaining population of active stem cells, even if the niche is removed. The inactive state implies that in absence of the stem cell niche all stem cells eventually differentiate. The inactive state in the model corresponds to early differentiation events (which are still reversible if the cell rebinds to the niche). This state may be linked to changes of gene expression occurring if stem cells expand outside the niche, as they have been described in [14] or as they occur in the case of niche impairment [43].

S2 Parameterization
In this section we fit the model to experimental data and set the model parameters for further simulations.
The quantification of stem cell numbers is challenging and different approaches have been developed. The estimates of stem cell number in an individual mouse are approximately 5000 − 15000, depending on age, sex, strain and methodology [25,10]. The percentage of empty niches can be quantified using unconditioned transplants, as described in [5], which suggests that in equilibrium 0.1 − 1% of niches are empty. The same study suggests that approximately 1−5% of HSC enter blood stream each day. Parabiosis experiments in mice suggest that 5 − 10% of HSC travel in bloodstream within 7 weeks [47]. A quantitative study of murine steady state hematopoiesis is given in a work by Busch et al. [8]. This work estimates a lower bound of 17,000 HSC per mouse (0.006% of bone marrow nucleated cells), 30% of which are active.
Quantifications in the human system are even more difficult. One possibility is to use CD34+CD38-ALDH+ cells [44]. This leads to approximately 10 7 HSC per kg of body weight. The number of human HSCs actively contributing to blood cell formation at a given time is estimated by 50,000-200,000 based on somatic mutations [23]. HSC mobilization studies in humans show that the number of circulating HSCs can increase more than 50 fold [17], implying that less than 2% of HSC are circulating under equilibrium conditions. After bone-marrow transplantation the average time HSC spend in the blood stream is approximately one day [15].
In the following we calibrate our proposed model to the murine system. The number of niches is set to K = 15000, based on estimations in the work of Catlin et al. [10]. Model dynamics remain qualitatively unchanged if different values for K are chosen.
Stem cell niche dynamics can be experimentally studied based on transplantation of stem cells in unconditioned hosts [4]. The work from [4] provides evidence that consecutive transplantation of small HSC doses (12700/7 cells on 7 consecutive days) leads to an increased chimerism compared to the transplantation of the total cell dose (12700 cells) at a single time point. This observation suggests that transplanted stem cells home to niche spaces that are unoccupied under steady state conditions. To re-equilibrate the system releases random (i.e., donor and host-derived) stem cells from the niche. In case of sequential transplantations this process is reiterated and donor-derived stem cells home to the niche at the expense of the randomly released host-derived cells. As a proof of concept, we demonstrate that our model can reproduce these experimental observations.
For this purpose we fit the model to the data of [4], using simulated annealing. The experimental data was extracted from Figure 6 of Bhattacharya et al. [4] based on visual inspection. Simulated annealing was performed using the MATLAB implementation by Vandekerckhove [42] with the following settings: Initial temperature: 10 3 , Stop temperature: 10 −3 , Max tries: 30. We assume that host and transplant cells have identical parameters, i.e. we set  [4]. As in the work of Stiehl et al. [37], we do not consider all transplanted cells to be stem cells, but reduce the number of cells added by one order of magnitude. Simulations indicate that if the transplanted cell count was much further reduced, the final chimerism in the model was smaller than in the experiments from Bhattacharya et al. [4], regardless of model parametrization. The transplanted stem cells added are assumed to be an equal mixture of activated and inactive cells, i.e., A L (0) = I L (0)¸where t = 0 is the time of transplantation. This assumption is made for the sake of simplicity since the impacts of stem cell harvesting and stem cell processing on HSC properties during transplantation are not well understood [45]. Based on our numerical analysis this assumption has no relevant impact on the obtained results. For comparison with experimental data we calculate the chimerism at time t as The parameters obtained from the simulated annealing are given in Table S1. The model can fit the experimentally observed short-term (one week) and long term (16 weeks) dynamics, however with different parameter sets. The results are shown in Figures S3 and S4. Qualitative model dynamics are similar for both parameter sets. All depicted simulations are based on the parameter set that matches long-term dynamics. Using the other parameter set leads to the same qualitative results. Our investigations did not find any parameter-set fitting both data-sets simultaneously. However, both parameter sets from the fitting match the data qualitatively, i.e., for both parameter sets multiple transplant doses imply an increased chimerism compared to the single dose scenario and the steady state value for the single dose scenario is identical for both parameter sets. A reason for the differences in the fitted parameters might be that in the model HSCs immediately home to the bone marrow after transplantation. Mechanisms related to homing that occur on the time-scale of a single day are not included in the model and could potentially explain the discrepancy. In addition, biological explanations for increasing uptake of subsequent transplants over time could also play a role, but are not considered in the model, since no short-term data of the multi-dose transplantation is available.
The number of HSC contained in CD34 + cell grafts and the fraction that homes to the bone marrow cannot be exactly quantified. For simulations of clinical scenarios we assume, in agreement with clinical data, that a standard human graft contains 3 · 10 3 HSC per kg of body weight and that the remaining transplanted CD34 + cells are progenitors [37]. We assume that HSC are equally distributed between the activated and the inactive state. Numerical simulations support that the qualitative results are robust with respect to changes in the number of transplanted activated and inactive HSC.
In formula (S43) we compute the steady state production of new progenitors, which is equal to the steady state division rate of stem cells under equilibrium conditions (since in equilibrium in average one dividing stem cell gives rise to one stem cell and one progenitor). It turns out that the progenitor production and thus the stem cell division rate is equal to u H , which in our estimates is, depending on the data used for the fit, between 1:25 and 1:250 per day. The stem cell birth/differentiation rate estimated by Busch et al [8] is 1:110 per day and is between our fits. Figure S3: Simulated HSC chimerism (best fit) after single dose transplant in an unconditioned host. Black points depict the short-term data from Bhattacharya et al. [4]. The red line shows the fit of the model to the short-term data (R 2 = 0.86). Parameter values are provided in Table S1. The dotted line (long-term parameters) is shown for the purpose of comparison, it results from the fit of the model to the chimerism after 112 days, which is depicted in Fig. S4. Figure S4: Simulated HSC chimerism (best fit) after 7 small transplants on consecutive days compared to a single dose transplant. Hosts are unconditioned. Black dotted lines correspond to the chimerism after 16 weeks (112 days) as measured in Bhattacharya et al. [4]. Parameter values are provided in Table S1. "LT" represents the fit of the model to the chimerism data at day 112, which results in the "long-term parameter set". "ST" is depicted for comparison and refers to model dynamics using the "short-term parameter set" obtained from fitting the model to the data depicted in Fig. S3. Table S1: Best fit of parameter values found using simulated annealing. We assume that host cells and transplanted cells have the same parameters. Therefore u denotes both u H and u L in the simulation, and similarly for the other parameters. In this section we summarize the equilibrium states of the model.
The model possesses the following four equilibrium states. The formal derivation is presented in section S4.1.
• Stem cell free equilibrium: This equilibrium corresponds to the absence of stem cells and will lead to the death of the organism. This equilibrium also corresponds to the state after conditioning for bone marrow transplantation before the stem and progenitor cells are infused or to a situation with graft failure. This equilibrium exists for all parameter values.
• Homeostasis: This equilibrium corresponds to the homeostatic state of a healthy organism. All cells are wild-type cells. As shown in section S4.1, this equilibrium only exists if activated stem cells divide more often than they differentiate and if generation of inactive stem cells outweighs their differentiation.
• Out-competition of wild-type cells: This equilibrium corresponds to a state where all wild-type cells have been out-competed. In case of stem cell transplantation this means that all HSC are derived from the graft. In case of malignancy, this state corresponds to the out-competition of HSC by leukemic cells which leads to the death of the organism. In case of clonal hematopoiesis of indeterminate potential (CHIP), this corresponds to monoclonal hematopoiesis. This equilibrium exists only if activated mutated (or respectively transplanted) stem cells divide more often than they differentiate and if generation of inactive stem cells outweighs their differentiation.
• Coexistence in equilibrium: In this equilibrium all population counts are positive. This corresponds to coexistence of wild-type and other cells. This type of equilibrium exists only for very specific parameter choices which is biologically unlikely to be realized.
Detailed calculations are provided in supplementary section S4.1.

S3.2 HSC dynamics in absence of niche cells
In the absence of niche cells it holds K = 0 implying N H (0) = N E (0) = 0. Under this assumptions, as shown in Figure S5, A H and I H converge to zero. This is in line with experimental observations from cell culture experiments [20,41,49].

S3.3 Response to perturbations of homeostatic cell counts
After perturbations, e.g., by blood loss, immune reaction, radiation, chemotherapy or bone marrow transplantation the hematopoietic system comes back to its equilibrium state. This property is captured by our model. When slightly perturbing homeostatic cell counts, the system returns to the homeostatic state. In supplementary section S4.4 we prove this property rigorously. An exemplary simulation is provided in figure S6.

S3.4 Impact of model parameters on equilibrium HSC counts and progenitor production
In this section we study how the cell counts of the healthy equilibrium and the production rate of progenitors depend on model parameters.
Changes of parameter values affect homeostatic HSC counts. The parameter dependence of the steady state cell counts is mathematically studied in section S4.2 of this supplement and summarized in Table 1 of the main text. Figure S7 illustrates how sensitively homeostatic cell counts depend on perturbations of parameter values. The figure shows the impact of a 5% perturbation of the parameter values from section S2. It implies that inactive HSC counts have a high sensitivity to changes of the rates of self-renewing divisions and activated cell differentiation (a 5% perturbation of the parameter implies a more than 40% perturbation of cell counts). Homeostatic cell counts are relatively insensitive to perturbations of all other parameters, i.e., a 5% perturbation of the respective parameter leads to less than 10% change of homeostatic cell counts. Especially, the model is  relatively insensitive to perturbations of the rates of attachment, detachment and inactive cell differentiation. Therefore, we consider perturbations of 100% in Figures 2-4 of the main text, where we study the transition dynamics in response to changes of these parameters.

S4.1 Steady states
In this section we calculate the steady states summarized in section S3.1.
The existence and uniqueness of solutions follows from the Theorem of Picard-Lindelöf. The system (S17a) -(S17i) has three equilibria that always exist, and a fourth which only exists given certain conditions. We present all four equilibria here.

Stem cell free equilibrium
The system always has a trivial steady state

Monoclonal equilibria
There exist two equilibria at which one clone has become extinct and the other clone is present at positive quantities. We denote these equilibria by E We note that for I * H > 0 and A * H > 0 to hold, it is necessary that r H > d A H . This corresponds to the biological requirement that proliferation rate of activated stem cells has to be larger than their differentiation rate. It is intuitive that this is required for maintenance of the cell population, since in a scenario where differentiation outweighs proliferation the cell number will decline to zero.

Substituting this relation in 0 =İ
Trivially, this holds for I * H = 0, which implies A * H = 0 and in turn N * H = 0. This is the stem cell free steady state.
, which thus yields an expression for N H in the non-trivial steady state: . This allows us to find the equilibrium of activated stem cells, A * H since: Figure S7: Impact of niche processes on homeostatic stem cell counts. The panels show how equilibrium stem cell counts (A-C, E-G) and progenitor production (D, H) change if a given parameter is increased or decreased by 1 % (A-D) or 5% (E-H) compared to the fitted parameters provided in Table S1. Black boxes denote an increased parameter value while gray boxes denote a decreased parameter value. Note that the impact of some parameter changes is too small to be visible on the depicted scale. The number of niche-bound stem cells is independent from the detachment rate.
The monoclonal equilibrium E * H with positive N H , A H and I H , is then given as: Note that the number of empty niches in the steady state is simply given as N *

This equilibrium has three positive components if and only if r H > d A H and Kb
Analogously the monoclonal equilibrium of the second clone (denoted by subscript L) has three positive components if and only if r L > d A L and Kb L Coexistence equilibria Imposing additional conditions on the parameters, both cell clones can coexist in an equilibrium state of the system. The equilibrium will be denoted as E * B (with the B denoting Both as both clones are present). Observe that the relation between the steady states of A H and I H shown in equation (S26) holds by the same argument as in the single clone case. As such we have the relations and

Substitution of the first relation inİ
Using the analogous equation for the second clone we find that either I * L = 0 or In the case I * L = 0 we also have A * L = 0 due to the relation in equation (S32). This in turn implies that N * L = 0 sinceṄ L = −u L N L when I L = 0 = A L . In this case we either obtain the stem cell free steady state E * 0 or the monoclonal steady state E * L . The analogous argumentation applies to the case I * H = 0, where we obtain either the stem cell free steady state E * 0 or the monoclonal steady state E * H . Equations (S33) and (S34) imply the following necessary condition for coexistence of two clones in an equilibrium state: Notably this condition does not depend on the rates of detachment from the niche (u H and u L ).
The condition (S35) is sufficient for existence of a manifold of coexistence equilibria. Indeed, set and consider the following manifold E * B for α ∈ R. We now show the following: For α ∈ (0, K − N * E ) the positive coexistence equilibria are given by E * B .
This is the case, because for α ∈ R holds that

S4.2 Parameter dependence of steady state cell counts
In this section we calculate how equilibrium cell counts change under small perturbations of model parameters. This was studied numerically in section S3.4.
We can use the partial derivatives of the steady state cell counts to find the parameter-dependence of the steady state values. This method leads to the results discussed in Section 3.1 and in Table 1 of the main text.
For brevity, we omit all but one of these calculations. Apart from the following example, all the cases are straightforward.
As such, > 1. Therefore, we obtain two biologically relevant parameter regimes, which are illustrated in Fig. S8: Figure S8: The dependence of the homeostatic count of activated HSC A * H on the parameter r H is relatively complex. Note the singularity at

S4.3 Parameter dependence of steady state progenitor production
In this section we calculate how equilibrium progenitor production changes under small perturbations of model parameters. This was studied numerically in section S3.4.
The production of progenitors per unit of time is described in the model by In the steady state situation, where dI H dt = dA H dt = 0, the production of progenitors should be in an equilibrium with the outflux of progenitors due to mature cell production. Hence determining how the steady state production of progenitors is influenced by changes of models parameters gives insights about the effect of the respective parameter changes on steady state mature cell production. In the following we study how perturbations of the parameters impact the steady state production of progenitors. The steady state production of progenitors is i * Thus, in steady state, the production of progenitors scales with u H N * H . Calculating the partial derivatives of u H N * H with respect to the parameters and assuming N * H > 0 reveals the dependence on parameter perturbations discussed in the main text: i * D H increases with K, u H , r H and b H , but decreases when d A H or d I H increases. In other words, increasing differentiation of HSC decreases the production of progenitors and in turn mature blood cells.

S4.4 Stability analysis
Linear stability of the homeostatic state In this section we mathematically show the linear asymptotic stability of the homeostatic equilibrium in absence of mutated cells. It was studied numerically in section S3.3.
Proposition: In absence of a second clone the homeostatic steady state is linearly stable whenever it exists

Perturbations of the homeostatic state and clonal fitness
In this section we investigate under which conditions the healthy equilibrium can be destabilized by adding a second stem cell clone. Based on this we derive the expressions for the clonal fitness.
We assume that < K which is necessary and sufficient for existence of the two monoclonal steady states E * H and E * L .
In the following we show that if we add a small number of malignant cells to the homeostatic state these Figure S9: Phase plane illustrations of equilibria. The circles and the thick black line denote equilibria. The white circle corresponds to the stem cell free equilibrium E * 0 . E * H and E * L correspond to monoclonal equilibria. The thick line E * B corresponds to the manifold of equilibria where two clones coexist. The arrows depict system dynamics in the presence of a single clone. If both clones have equal fitness, the system will move asymptotically toward a point on the thick line, which depends on the initial conditions. In panel a the two clones have equal fitness, whereas panel b depicts a scenario where the mutated clone (denoted by L) has a competitive advantage.
Proposition: Consider the model (S17a)-(S17i) of two competing clones. Let The homeostatic steady state is linearly asymptotically Proof: Consider the linearization of the model around the healthy equilibrium.
The linearization of the model around the homeostatic equilibrium E * H is given by where A 1 is the linarization of the model for one HSC clone around the homeostatic steady state, A 2 is a real 3 × 3 matrix and The characteristic polynomial of A is The latter holds since A is a block triangular matrix. The eigenvalues of A 1 correspond to the eigenvalues of the linearization of the model for one HSC clone around the homeostatic steady state. We know from the previous section that they are negative. We now have to investigate the sign of the real parts of the roots of χ A4 (X). It holds . Together with the expression for N * E,E H , we obtain The characteristic polynomial expands to We now check the Routh-Hurwitz Criterion. We write χ A4 (X) := X 3 + a 2 X 2 + a 1 X + a 0 . The polynomial has only roots with negative real parts if (i) a 2 > 0 and (ii) a 2 a 1 − a 0 > 0 and (iii) a 0 > 0.
(i) is satisfied, since it is a sum of positive terms. (ii) expands to which is a sum of positive terms and hence positive. As shown above (iii) holds if and only if Taken together this implies that the homeostatic state cannot be destabilized by a second clone if However in the opposite case the homeostatic state is destabilized. In the first case, leukemic cells introduced to the system die out whereas in the second case they grow. This implies that the relation between decides whether an additional clone introduced in the system can persist or not.
These findings provide the theoretical foundation for the following notion of fitness. If we define the fitness of the HSC by and the fitness of LSC by , then the clone with the larger fitness can destabilize the homeostatic state of the clone with lower fitness but not vice versa.
S5 Peripheral blood malignant cell burden can change faster or slower compared to the malignant cell burden in the stem cell niche Simulations indicate that the malignant cell burden in the stem cell compartment and in the peripheral blood can differ significantly. Figure S10a provides an example for a scenario where the malignant cell burden in the stem cell niche increases earlier and is higher compared to the progenitor compartment, while figure S10b provides an example for a scenario where the malignant cell burden increases earlier and is higher in the progenitor compartment compared to the stem cell compartment.

S6 Multiple stem cell doses increase chimerism in the unpreconditioned host but not in the preconditioned host
To reduce complications resulting from low circulating blood cell counts after bone marrow transplantation different approaches have been suggested. One approach involves splitting the transplant in multiple doses that are infused on subsequent days [4,22]. In unpreconditioned mice the latter results in a higher chimerism [4]. This observation supports the concept that there exist empty niches under steady state conditions that are filled by subsequent transplantation of small cell doses [4]. Our calibrated model is able to reproduce this observation as discussed in supplementary section S2.
It is a clinically relevant question, whether the engraftment of transplanted cells in humans can be improved if the transplant is split into multiple doses. This question not only arises in the context of bone marrow transplantation in the case of hematologic malignancy, it is also relevant when genetically engineered HSC are used to cure metabolic diseases [27,2,16]. Furthermore, a transplant splitting approach could be useful if it is combined with a reduced intensity conditioning (RIC) regime. This could allow to further reduce the complications of stem cell transplants in the elderly.
To obtain insights into the question how preconditioning affects the advantageous impact of splitting a transplant into multiple doses, we simulate the following scenario: Starting at homeostasis, one third of the niche-bound HSC were replaced by a leukemic clone with competitive advantage. This scenario corresponds to the setting of a newly diagnosed acute myeloid leukemia (AML) patient. Preconditioning was modeled as pretransplant ablation of a prescribed fraction of healthy and leukemic cells. We then simulated transplantation of healthy HSC either as single dose or as multiple subsequent doses adding up to the same cell-count. To evaluate the outcome, we use the residual disease 100 days after transplantation.
The number of transplanted healthy HSC were chosen as 3.75 · 10 4 for the first dose and 2 · 10 4 for the subsequent doses to mirror the methodology of Landau et al. [22]. Note that this is two orders of magnitude lower than the cell-counts described by Landau et al. [22], since we do not assume all transplanted cells be to HSC.
The results are shown in Figure S11. Our simulations indicate that the benefit of splitting a transplant into multiple doses is greater for the unpreconditioned host. In the unpreconditioned host a reduction in the residual disease at day 100 of up to 8% was achieved when the transplant was split in multiple smaller doses. However, with preconditioning this effect becomes negligible. This finding is in line with results from a clinical trial [22].
The model suggests that subsequent transplantation of multiple cell doses could be beneficial compared to a single dose transplant only in absence of relevant preconditioning.  Figure S10: Time evolution of the malignant cell burden in the stem and progenitor compartment. In both panels, the niche is invaded by a clone with a fitness advantage due to a higher proliferation rate (referred to as clone 2). In panel (b) the advantage is smaller compared to panel (a). In addition in the scenario depicted in panel (a) the invading clone has a decreased detachment rate compared to the wild-type clone, while in panel (b) the detachment rate of the invading clone is increased compared to the wild-type cells. The simulations show that depending on stem cell parameters the malignant cell burden in the progenitor compartment can increase before or after the malignant cell burden in the stem cell compartment. Parameters: (a), r 2 = 1.2r 1 and u 2 = 5u 1 ; (b) r 2 = 4r 1 and u 2 = 0.4u 1 . The remaining parameters are as in Table S1. Figure S11: Residual disease after bone marrow transplantation with and without preconditioning. Before preconditioning, 33% of the stem cell niche are occupied by a leukemic clone with a competitive advantage. Without treatment, the leukemic cell mass increases to approx. 65% within 100 days, shown as a dashed blue line. Before stem cells are transplanted a variable fraction of leukemic and healthy cells is ablated due to preconditioning. Preconditioning is assumed to be equally effective for healthy and malignant stem cells. To model the transplant either one dose of free healthy stem cells or multiple doses at subsequent days were added to the system. The total amount of transplanted cells is equal in both cases. The cell counts were based on the experimental work of Landau et al. [22] as described in the text. The effect of splitting the transplant in multiple doses decreases in case of preconditioning.
S7 Mobilizing agents could decrease the residual disease in the unpreconditioned host but lead to an increase in the preconditioned host.
As shown in the previous section, our simulations indicate that the chimerism obtained after stem cell transplantation to an unpreconditioned host crucially depends on the amount of free niche spaces. This finding led us to question how the chimerism could be affected by stem cell mobilizing agents.
Mobilization agents force stem cells to detach from the niche. This could improve the efficiency of a bone marrow transplant following reduced or no preconditioning. One question arising in this context is whether no or reduced intensity conditioning (RIC) accompanied by mobilization agents and stem cell transplantation could be a treatment option for hematological cancers in the elderly not eligible for conventional transplantation regimens that require high dose chemotherapy. On the one hand, the use of mobilization agents could increase the number of donor cells homing to the niche and thus reduce the malignant cell burden. On the other hand detached stem cells are more prone to proliferate [34,32,33,46]. For this reason stem cell mobilization could trigger the expansion of the malignant cell population. To better understand the relative contribution of these two processes to disease evolution we simulated a scenario where a malignant stem cell clone with a competitive advantage initially occupies one third of the bone marrow niche. First we consider a scenario without preconditioning. A transplant of healthy stem cells is introduced at day 0 and the detachment rate of healthy and malignant stem cells is increased at the day of transplantation due to mobilization agents. Later the detachment rate returns to its original value. The model was simulated using the parameters and transplant doses given in section S2 of this supplement. We assume that mobilizing agents act immediately when administered. This is a good approximation of reality, since the peak effect of mobilizing agents is reached after 8-12 hours [30]. This is short compared to the pre-clinical phase of malignant diseases which takes multiple years [13]. Similarly, we neglect the time transplanted cells require to home to the bone marrow, which is approximately one day [4,15].
The simulations suggest that stem cell mobilization for short periods of time leads to an increased engraftment of healthy cells compared to the setting without mobilization. However, the number of malignant cells increases if the duration of mobilization treatment is longer than a critical time t c . The stronger the mobilization effect, the earlier the malignant cell expansion dominates and the smaller the value of t c . Figure S12 shows after which duration of mobilization treatment the malignant cell burden starts to increase and how this depends on the strength of the mobilization effect.
To consider the effect of mobilizing agents combined with preconditioning, we simulated the same scenario as above, with a malignant stem cell subpopulation occupying one third of the stem cell niche. On day 0, preconditioning was simulated by removing a certain percentage of all stem cells. We assume that preconditioning acts equally on healthy and malignant stem cells. After preconditioning, a transplant of free healthy cells was added, accompanied by mobilizing agents that were effective for 0.5 days.
In figure S13 we show the changes to the residual disease after 20 days in comparison to the same scenario without mobilization. The simulations suggest that following efficient preconditioning (i.e. preconditioning that ablates large numbers of stem cells), increased mobilization leads to an increase in residual disease. However, if preconditioning is performed with reduced intensity, the increased mobilization can lead to a decrease in residual disease.
The effect of mobilizing agents may also be of relevance when genetically engineered cells have to be introduced to the stem cell niche of an unpreconditioned host, as in the case of gene therapy. Next we addressed the question whether the number of genetically engineered cells homing to the bone marrow could be increased by using mobilizing agents prior to transplantation. Figure S14 shows that in case of no preconditioning the number of genetically engineered stem cells homing to the niche increases if mobilizing agents are administered before transplantation, however this effect becomes negligible in the context of preconditioning.
Taken together our simulations suggest that mobilizing agents could improve the homing of transplanted cells to the niche only in hosts with reduced or no preconditioning. Figure S12: Impact of mobilizing agents on the nonpreconditioned host. In an unpreconditioned host stem cell mobilizing agents can increase the number of graft-related cells homing to the niche. This allows for a reduction of the malignant cell burden in absence of preconditioning. However, if the effect of the mobilizing agents lasts longer than a critical time t c , malignant cells detached from the niche expand and increase the disease burden. To assess this effect, we compare the following two hypothetical therapy approaches: (i) bone marrow transplantation without preconditioning, (ii) bone marrow transplantation without preconditioning accompanied by mobilizing agents for up to 1.3 days. We compare the malignant cell burden after 100 days. When the disease burden after strategy (ii) is not higher compared to strategy (i) we define the effect of mobilization as beneficial. This is indicated by the green area. At the beginning of the simulations a leukemic clone occupies one third of the niche. The horizontal axis quantifies by how much the mobilizing agent increases the rate of stem cell detachment from the niche relative to wild type HSC, the vertical axis quantifies how long the mobilizing effect remains in the system. Figure S13: Impact of mobilizing agents after preconditioning. Before transplantation a fraction of HSC was ablated by preconditioning ("0%" indicates no preconditioning, "20%" indicates ablation of 20% of healthy and malignant stem cells from the niche.) After ablation an HSC transplant was added. Post transplantation HSC detachment was increased by mobilizing agents for 0.5 days (A), 1 day (B) or 2 days (C). The x-axis quantifies by which factor mobilization agents increase the HSC detachment rate. The y-axis shows the increase of malignant cells due to the mobilizing agent (residual disease after therapy with preconditioning, transplantation plus mobilization minus residual disease after therapy with preconditioning and transplantation). Negative values indicate that the mobilizing agent has a beneficial effect. For two days of mobilization there is no beneficial effect, regardless of the strength of the mobilization. Preconditioning of 0% is also shown in Fig. S12. The curve for 0% preconditioning in (A) corresponds to a horizontal line in Fig. S12 intersecting the y-axis at 0.5 days. Red color in Fig. S12 corresponds to positive values in (A), green color to negative values. Analogous relations hold for (B) and (C). In (D) the HSC detachment rate is increased by a factor of 20 for the duration indicated on the x-axis. The y-axis shows therapy outcome after 100 days. The preconditioning of 0% in (D) corresponds to a vertical line in Fig.  S12 intersecting the x-axis at 20. Again, negative values correspond to green, positive to red color in Fig. S12. Figure S14: Impact of mobilizing agents on the homing of genetically engineered HSC in unpreconditioned hosts. In the nonpreconditioned host, the number of genetically engineered stem cells (referred to as clone 2) homing to the bone marrow can be increased by the use of mobilizing agents prior to transplantation. In the simulation depicted in black HSC mobilization has been increased by a factor of 2 during the day before transplantation, which is a standard duration in clinical practice [35]. At the time of transplantation, the drug-induced mobilization was assumed to have stopped. For comparison, the blue dotted lines depicts a scenario with no increase in HSC mobilization.