Clonal dominance and transplantation dynamics in hematopoietic stem cell compartments

Hematopoietic stem cells in mammals are known to reside mostly in the bone marrow, but also transitively passage in small numbers in the blood. Experimental findings have suggested that they exist in a dynamic equilibrium, continuously migrating between these two compartments. Here we construct an individual-based mathematical model of this process, which is parametrised using existing empirical findings from mice. This approach allows us to quantify the amount of migration between the bone marrow niches and the peripheral blood. We use this model to investigate clonal hematopoiesis, which is a significant risk factor for hematologic cancers. We also analyse the engraftment of donor stem cells into non-conditioned and conditioned hosts, quantifying the impact of different treatment scenarios. The simplicity of the model permits a thorough mathematical analysis, providing deeper insights into the dynamics of both the model and of the real-world system. We predict the time taken for mutant clones to expand within a host, as well as chimerism levels that can be expected following transplantation therapy, and the probability that a preconditioned host is reconstituted by donor cells.


Author summary
Clonal hematopoiesis-where mature myeloid cells in the blood deriving from a single stem cell are over-represented-is a major risk factor for overt hematologic malignancies. To quantify how likely this phenomena is, we combine existing observations with a novel stochastic model and extensive mathematical analysis. This approach allows us to observe the hidden dynamics of the hematopoietic system. We conclude that for a clone to be detectable within the lifetime of a mouse, it requires a selective advantage. I.e. the clonal expansion cannot be explained by neutral drift alone. Furthermore, we use our model to describe the dynamics of hematopoiesis after stem cell transplantation. In agreement with earlier findings, we observe that niche-space saturation decreases engraftment efficiency. We further discuss the implications of our findings for human hematopoiesis where the quantity and role of stem cells is frequently debated.

Introduction
The hematopoietic system has evolved to satisfy the immune, respiratory, and coagulation demands of the host. A complex division tree provides both amplification of cell numbers and a variety of differentiated cells with distinct roles in the body [1][2][3]. In a typical adult human *10 11 terminally differentiated blood cells are produced each day [3][4][5]. It has been argued that the division tree prevents the accumulation of mutations, which are inevitable given the huge number of cell divisions [6][7][8]. At the base of the tree are hematopoietic stem cells (HSCs). These have the ability to differentiate into all hematopoietic cell lineages, as well as the capacity to self-renew [1,9], although the exact role of HSCs in blood production is still debated [10,11].
With an aging population, hematopoietic malignancies are increasingly prevalent [12]. Clonal hematopoiesis-where a lineage derived from a single HSC is overrepresented-has been identified as a significant risk factor for hematologic cancers [13][14][15]. To assess the risks posed to the host we need an understanding of how fast clones are growing, when they initiate, and if they would subvert physiologic homeostatic control.
The number of HSCs within a mouse is estimated at *0.01% of bone marrow cellularity [16,17], which amounts to *10,000 HSCs per host [3,16,18,19]. In humans this number is subject to debate; limited data has lead to the hypothesis that HSC numbers are conserved across all mammals [18], but the fraction of 'active' HSCs depends on the mass of the organism [20] (see also Refs [5,21] for a discussion).
Within an organism, the HSCs predominantly reside in so-called bone marrow niches: specialised micro-environments that provide optimal conditions for maintenance and regulation of the HSCs [22,23]. There are likely a finite number of niches within the bone marrow, and it is believed that they are not all occupied at the same time [16]. The number of niches is likely roughly equal to the number of HSCs, and through transplantation experiments in mice it has been shown that *1% of the niches are unoccupied at any time [16,24]. A similar number of HSCs are found in the peripheral blood of the host [16]. These free HSCs are phenotypically and functionally comparable to (although distinguishable from) bone marrow HSCs [19,25]. The HSCs have a residence time of minutes in the peripheral blood, and parabiosis experiments (anatomical joining of two individuals) have shown that circulating HSCs can engraft to the bone marrow [25]. It has also been shown that HSCs can detach from the niches without cell division taking place [19]. These findings paint a picture of HSCs migrating between the peripheral blood and the bone marrow niches, maintaining a dynamic equilibrium between the two compartments.
In this manuscript we construct a model from the above described processes, and we use this to answer questions about clonally dominant hematopoiesis. We first consider this in mice, where we use previously reported values to parametrise our model. The model is general enough that it also captures scenarios of transplantation into both preconditioned (host HSCs removed) and non-preconditioned hosts: the free niches and the migration between compartments also allows for intravenously injected donor HSCs to attach to the bone marrow niches and to contribute to hematopoiesis in the host. In the discussion we comment on the implications of these results for human hematopoiesis.

Materials and methods
Our model, shown schematically in Fig 1, contains two compartments for the HSCs. The bone marrow (BM) compartment consists in our model of a fixed number, N, of niches. This means that a maximum of N HSCs can be found there at any time, but generally the number of occupied niches is less than N. The peripheral blood (PB) compartment, however, has no size restriction. The number of cells in the PB and BM at a given time are given by s and n, respectively.
The dynamics are indicated by arrows in Fig 1. Our model is stochastic and individual based, such that events are chosen randomly and waiting times between events are exponentially distributed. Simulations are performed using the Gillespie stochastic simulation algorithm (SSA) [26]. HSCs in our model are only capable of dividing when attached to a niche; outside the niche, pre-malignant cells are incapable of proliferating due to the unfavourable conditions. Upon division, the new daughter HSC enters another niche with probability ρ, or is ejected into the PB. Here ρ depends on the number of free niches, i.e. ρ = ρ(n), and should satisfy ρ(N) = 0, such that a daughter cell cannot attach if all niches are occupied. Likewise, the migration of a HSC from the PB to the BM should depend on the number of empty niches. We choose the attachment rate as a(N − n)/N per cell. In general, cells can die in both compartments. However, we expect the death rate in the PB, δ, to be higher than the death rate in the BM, δ 0 , as the PB is a less favourable environment.
For our initial analysis, we assume there is no death in the BM compartment (δ 0 = 0), and new cells are always ejected into the PB (ρ = 0). These assumptions are relaxed in our detailed analysis, which can be found in the Supporting Information (S1 Supporting Information). The bone marrow (BM) compartment has a fixed total of N niches. At a given time, n of the niches are occupied, and N − n remain unoccupied. The peripheral blood (PB) compartment has no size restriction, and at a given time contains s HSCs. A HSC in the BM can detach at rate d and enter the PB, while a cell in the PB can attach to an unoccupied niche with rate a(N − n)/N. Here (N − n)/N is the fraction of unoccupied niches. HSCs may die in the PB or BM with rates δ and δ 0 . Reproduction (symmetric division) of HSCs occurs at rate β. The new daughter cell attaches to an empty niche with probability ρ, otherwise it is ejected into the PB. Dynamics are concretely described by the reactions in Eq (1).
A two-compartment model has been considered previously by Roeder and colleagues [27,28]. The rate of migration between the compartments is controlled by the number of cells in each compartment, as well as a cell-intrinsic continuous parameter which increases or decreases depending on which compartment the cell is in. This parameter also controls the differentiation of the HSCs. Further models of HSC dynamics, for example [29][30][31][32][33], have not considered the migration of cells between compartments. For example, Dingli et al. consider a constant-size population of HSCs in a homogeneous microenvironment [31,32]. Competition between wildtype and malignant cells then follows a Moran process. In our model the BM compartment size is fixed, but cell numbers can fluctuate.
To initially parametrise our model we consider only one species of HSCs: those which belong to the host. In a steady-state organism, the number of HSCs in the PB and BM are close to their equilibrium values, which are labelled as s Ã and n Ã , respectively. These values have been reported previously in the literature for mice, and are provided in Table 1. Other previously reported values include the total number of HSC niches N, HSC division rate β, and the time that cells spend in the PB, which we denote as ℓ. Using these values we can quantify the remaining model parameters δ, d, and a. These results are discussed in the next section.
When considering a second population of cells, such as a mutant clone or donor cells following transplantation, we may want to impose a selective effect relative to the host HSCs. We therefore allow the mutant/donor cells to proliferate with rate β 2 = (1 + ε)β, where ε represents the strength of selection. For ε = 0, the mutant/donor cells proliferate at the same rate as the host HSCs. In the S1 Supporting Information we consider the general scenario of selection acting on all parameters. Our analysis delivers an interesting result: the impact of selection on clonal expansion is independent of which parameter it acts on (provided δ 0 = ρ = 0).
For clonality and chimerism we use the same definition: the fraction of cells within the BM compartment that are derived from the initial mutant or the donor population of cells. Typically, experimental measurements of clonality and chimerism use mature cells rather than HSCs. However, this is beyond the scope of our model so we use HSC fraction as a proxy for this measurement. We are therefore implicitly assuming that HSC chimerism correlates with mature cell chimerism. The literature on the role of HSCs in native hematopoiesis is split [10,35] (also reviewed in [36]). For the division rate of HSCs in mice we use the value β = 1/39 per day. This is the average division rate of all HSCs within a host deduced from CFSE-staining experiments [34], but again there is some disagreement in reported values for this quantity [34,37,38]. These differences arise from the interpretation of HSC cell-cycle dynamics.
More concretely, our model consists of four sub-populations: n 1 is the number of host or wildtype cells located in the BM, and s 1 is the number of cells of this type in the PB. Likewise, n 2 and s 2 are the number of mutant/donor cells in the BM and PB, respectively. The cell numbers are affected by the processes indicated in Fig 1 (with δ 0 = ρ = 0). The effect of these events Table 1. Parameter values from empirical murine observations. These are equilibrium values in healthy mice.

Accessibility
A Wolfram Mathematica notebook containing the analytical details can be found at https:// github.com/ashcroftp/clonal-hematopoiesis-2017. This location also contains the Gillespie stochastic simulation code used to generate all data in this manuscript, along with the data files.

Steady-state HSC dynamics in mice
By considering just the cells of the host organism, we can compute the steady state of our system from Eq (2), and hence express the model parameters δ, d, and a in terms of the known quantities displayed in Table 1. These expressions are shown in Table 2, where we also  (Table 1), we find disparate dynamics in our model. At one extreme, the average time a cell spends in the BM compartment (1/d) can be less than two hours (for s Ã = 100 cells and ℓ = 1 minute). Thus under these parameters the HSCs migrate back-and-forth very frequently between the niches and blood, and the flux of cells between these compartments over a day (s Ã /ℓ) is significantly larger than the population size. In fact, under these conditions 144,000 HSCs per day leave the marrow and enter the blood. With slower turnover in the PB compartment (ℓ = 5 minutes, but still s Ã = 100), the average BM residency time of a single HSC is eight hours, and 28,800 HSCs leave the bone marrow per day. At the other extreme, if the PB compartment is as small as reported in Ref. [19] (s Ã = 1 cell), then the residency time of each HSC in the bone marrow niche is between 8 and 290 days (for ℓ = 1 and 5 minutes, respectively). Under these conditions the number of cells entering the PB compartment per day is 1,440 and 288, respectively. For an intermediate PB size of s Ã = 10, the BM residency time is between 17 and 90 hours (for ℓ = 1 and 5 minutes, respectively), and the flux of cells leaving the BM is a factor ten greater than for s Ã = 1.

Clonal dominance in mice
Clonal dominance occurs when a single HSC generates a mature lineage which outweighs the lineages of other HSCs, or where one clone of HSCs outnumbers the others. The definition of when a clone is dominant is not entirely conclusive. Previous studies of human malignancies have used a variant allele frequency of 2%, corresponding to a clone that represents 4% of the population [39,40]. For completeness we investigate clonality ranges from 0.1% to 100%.
In the context of disease, this clone usually carries specific mutations which may confer a selective advantage over the wildtype cells in a defined cellular compartment. The de novo emergence of such a mutant occurs following a reproduction event. Therefore, in our model with ρ = 0, after the mutant cell is generated it is located in the PB compartment, and for the clone to expand it must first migrate back to the BM. This initial phase of the dynamics is considered in general in the next section of transplant dynamics, where a positive number S of mutant/donor cells are placed in the PB. We find (as shown in the S1 Supporting Information) that the expected number of these cells that attach to the BM after this initial dynamical phase is We then apply a fast-variable elimination technique to calculate how long it takes for this clone to expand within the host [41,42]. This procedure reduces the dimensionality of our system, and makes it analytically tractable. A full description of the analysis can be found in the S1 Supporting Information, but we outline the main steps and results of this procedure below. We first move from the master equation-the exact probabilistic description of the stochastic dynamics-to a set of four stochastic differential equations (SDEs) for each of the variables via an expansion in powers of the large parameter N [43]. We then use the projection method of Constable et al. [41,42] to reduce this system to a single SDE describing the relative size of the clone. This projection relies on the weak-selection assumption, i.e. 0 ε ( 1. The standard results of Brownian motion are then applied to obtain the statistics of the clone's expansion. In particular, the probability that the mutant/donor HSCs reach a fraction 0 < σ 1 of the occupied BM niches is given by where z 0 is the initial clone size can be found explicitly from Eq (3), such that z 0 = n 2 /N. We also have ξ = n Ã /N, and Λ is a constant describing the strength of deterministic drift relative to stochastic diffusion. Concretely, we have The mean time for the clone to expand to size σ (i.e. the mean conditional time) is written as Here B is another constant describing the magnitude of the diffusion, and is given by Although a general closed-form solution to Eq (6) is possible, it is too long to display here. Instead we use an algebraic software package to solve the second-order differential equation. A similar expression to Eq (6) can be obtained for the second moment of the fixation time, as shown in [44] and repeated in the S1 Supporting Information. The first scenario we consider is the expansion of a neutral clone (ε = 0); i.e. how likely is it that a single cell expands into a detectable clone in the absence of selection? It is known that the time to fixation of a neutral clone in a fixed-size population grows linearly in the system size [45]. Interestingly and importantly, in intestinal crypts this fixation is seen frequently because N ¼ Oð10Þ [46]. In the hematopoietic system, however, it likely takes considerably longer than this due to the relatively large number of stem cells. Solving Eq (6) with ε = 0 gives the mean conditional expansion time as From this solution we find that it takes, on average, 5-45 years for a neutral clone to reach 1% clonality (*100 HSCs). Expanding to larger sizes takes considerably longer, as highlighted in Fig 2. Therefore, clonal hematopoiesis in mice is unlikely to result from neutral clonal expansion; for a clone to expand within the lifetime of a mouse it must have a selective advantage. Neutral results for human systems are considered in the discussion. When the mutant clone has an advantage, there is always some selective force promoting this cell type. Therefore the probability of such a clone expanding is higher than the neutral case, as seen from Eq (4). In Fig 2 we illustrate the time taken for a single mutant HSC to reach specified levels of clonal dominance for different selective advantages. Advantageous clones (β 2 /β > 1) initially grow exponentially in time [Fig 2(a)], and are much faster than neutral expansion (β 2 /β = 1). These clones can reach levels of up to 90% relatively quickly, however replacing the final few host cells takes much longer. The advantage that a mutant clone must have if it is to represent a certain fraction of the population in a given period of time can be found from Fig 2(b). For a single mutant to completely take over in two years, it requires a fold reproductive advantage of β 2 /β % 2 [dashed lines in Fig 2(b)]. This means that the cells in this clone are dividing at least twice as fast as the wildtype host cells. To achieve 1% clonality in this timeframe, the advantage only has to be β 2 /β % 1.2. For the clone to expand in shorter time intervals, a substantially larger selective advantage is required. For example, 100% clonality in six months from emergence of the mutant requires β 2 /β % 5.5, i.e. the dominant clone needs to divide more than five times faster than the wildtype counterparts.
As shown in the S1 Supporting Information, Eqs (4) and (6) are equivalent to the results obtained from a two-species Moran process. This suggests the two-compartment structure is not necessary to capture the behaviour of clonal dominance. However, the consideration of multiple compartments is required to understand transplantation dynamics, as covered in the next section.

Transplant success in mice
We now turn our attention to the scenario of HSC transplantation. As previously mentioned this situation is analogous to the disease spread case, with the exception that the initial 'dose' of HSCs can be larger than one. We first consider the case of a non-preconditioned host. We then move onto transplantation in preconditioned hosts, where all host cells have been removed.
Engraftment in a non-preconditioned host. Multiple experiments have tested the hypothesis that donor HSCs can engraft into a host which has not been pretreated to remove some or all of the host organism's HSCs [16,19,34,[47][48][49][50][51][52][53]. These studies have found that engraftment can be successful; following repeated transplantations mice display a chimerism with up to 40% of the HSCs deriving from the donor [48][49][50].
In this scenario we start with a healthy host organism and inject a dose of S donor HSCs into the PB compartment, in line with the experimental protocols mentioned above. These donor cells can be neutral, or may have a selective (dis)advantage. Injecting neutral cells reflects the in vivo experiments described above, while advantageous cells can be used to improve the chances of eliminating the host cells. Transplanting disadvantageous cells would reflect the introduction of 'normal' HSCs into an already diseased host carrying advantageous cells. We do not consider this scenario further here, as the diseased cells are highly unlikely to be replaced without host preconditioning.  (6). Shaded regions are the predicted standard deviations, using the formula presented in the S1 Supporting Information. Here ℓ = 3 minutes, s* = 100, and the remaining parameters are as in Table 1. We can separate the engraftment dynamics of these donor cells into two different regimes: i) the initial relaxation to a steady state where the total number of HSCs is stable, and ii) longtime dynamics eventually leading to the extinction of either the host or donor HSCs. We focus on these regimes separately. Upon the initial injection of the donor HSCs, the PB compartment contains more cells than the equilibrium value s Ã . This leads to a net flux of cells attaching to the unoccupied niches in the BM until the population relaxes to its equilibrium size. Once the equilibrium is reached, the initial dynamics end, and the long-term noise-driven dynamics take over (discussed below). The challenge for this first part is to determine how many of the donor HSCs have attached to the BM at the end of the initial dynamics.
We identify two distinct behaviours which occur under low and high doses of donor HSCs. If the dose is small (S ⪡ N À n Ã ), then the number of donor HSCs that attach to the BM is given by Eq (3), and is proportional to the dose size S. To obtain this result we have assumed that the number of occupied niches remains constant, such that each donor cell has the same chance of finding an empty niche. However, if the dose of donor HSCs is large enough then all niches become occupied and the BM compartment is saturated; attachment to the BM can only occur following a detachment. Using this assumption, the initial dynamics can then be described by the linear ODEs where s(t), the total number of cells in the PB compartment, is found from _ s ¼ bN À ds. A derivation of Eq (9) can be found in the S1 Supporting Information.
The predicted chimerism, and the accuracy of these predictions, at the end of the initial phase are shown in Fig 3. The efficiency of donor cell engraftment decreases in the large-dose regime (S > N À n Ã ). This is simply because the niche-space is saturated, so HSCs spend longer in the blood and are more likely to perish. If the lifetime in the PB (ℓ) is short, then we have more frequent migration between compartments, as highlighted in Table 2. Hence smaller ℓ leads to higher chimerism. The approximation from Eq (9) becomes increasingly accurate for larger doses. The two approximations break down at the cross-over region between small and large doses. In this regime the number of occupied niches does not reach a stable value.
If the donor cells have a selective (dis)advantage, then the deterministic dynamics predict the eventual extinction of either the host or donor cells. However, the selective effect is usually small and only acts on a longer timescale. Therefore the initial dynamics are largely unaffected by selection, and we assume neutral donor cell properties when we model the initial dynamics.
The inefficiency of large doses can be overcome by administering multiple small doses over a long period. In this way we prevent the niches from becoming saturated and fewer donor cells die in the PB. Hence we should be able to obtain a higher level of engraftment when compared to a single-bolus injection of the same total number of donor HSCs. These effects have been tested experimentally [19,[48][49][50]. Parabiosis experiments are also an extreme example of this; they represent a continuous supply of donor cells [25]. As shown in Fig 4, our model captures the same qualitative behaviour as reported in the experiments: Multiple doses lead to higher levels of chimerism at the end of the initial phase of dynamics. This effect is highlighted more when the total dose size is large. Using our analysis we know how efficient each dose is, and what levels of chimerism can be achieved. Hence our model can be used to optimise dosing schedules such that they are maximally efficient.
We note here that engraftment efficiency is only important when donor cells are rare and there is no danger to life. This is the case, for example, in experimental protocols when tracking small numbers of cells. It makes sense to here use the multiple dosing strategy. Transplantation following preconditioning, however, provides a more viable approach to disease  Table 1.
https://doi.org/10.1371/journal.pcbi.1005803.g004 treatment where patient survival needs to be maximised. In this case the dose size should be increased, but it should be considered that there are diminishing returns in engraftment when the dose size is large enough to saturate all open niches. This dose size can be read from Fig 3. The long-term dynamics are handled in the same way as the clonal dominance results above; we use Eqs (4) and (6) to show how the number of donor cells injected into the PB affects the probability that they expand (as opposed to die out), and the time it takes for the host cells to be completely displaced. One key result is that a dose of just eight HSCs with an advantage of β 2 /β = 1.1 has over 50% chance to fixate in the host. However, the time for this to happen is *16 years. With a reproductive advantage of β 2 /β = 1.5 the success rate is *95% for the same dose, and the time taken now falls to *4 years. Further results are found in the S1 Supporting Information.
Engraftment in a preconditioned host. HSC transplantation procedures are often preceded by treatment or irradiation of the host-referred to as host preconditioning. This greatly reduces the number of host HSCs in the BM compartment. For this section we assume complete conditioning such that no host HSCs remain, i.e. myeloablative conditioning. Following the pre-treatment, a dose of donor HSCs of size S is injected into the PB compartment. We then want to know the probability that these cells reconstitute the organism's hematopoietic system. For this section we assume that donor HSCs have identical properties to the wildtype cells (i.e. no selection). We further assume that all donor HSCs have the potential to reconstitute the hematopoietic system in the long-term-in experiments this is not always the case as not all cells which are sorted as phenotypic HSCs (as defined by surface markers) are functional, reconstituting HSCs (see e.g. [54]). Because of this assumption, we only show results for the injection of a single donor HSC into the conditioned host (S ¼ 1). Higher doses lead to a greater probability of reconstitution. A further assumption is, that the host maintains (or is provided with) enough mature blood cells during the reconstitution period to sustain life.
We consider two approaches for estimating the probability of hematopoietic reconstitution. As a first-order approximation, the probability that a single HSC in the PB compartment dies is ψ = δ/(δ + a). Here we have assumed that all niches are unoccupied, such that the attachment rate per cell is a(N − 0)/N = a. For a dose of size S, the reconstitution probability is φ ¼ 1 À c S . Hence we have This prediction, Eq (10), is shown as dotted lines in Fig 5, which, however, does not agree with results from the model. The second approach considers all possible combinations of detachments and reattachments, as well as reproduction events. This leads to a reconstitution probability, given a dose of S donor cells, of which is derived in the S1 Supporting Information. This result, Eq (11), is shown as solid lines in Fig 5, and is in excellent agreement with the reconstitution probability observed in simulations. From these results we can conclude that, in our model, HSCs migrate multiple times between the PB and BM before they establish a sustainable population. It is also the case that in this model a single donor HSC is sufficient to repopulate a conditioned host in *90-99% of cases across all the parameter ranges reported in Table 1.

Discussion
We have introduced a mathematical model that describes the back-and-forth migration of hematopoietic stem cells between the blood and bone marrow within a host. This is motivated by the literature of HSC dynamics in mice. The complexity of the model has been kept to a minimum to allow us to parametrise it based on empirical results. The model is also analytically tractable, permitting a more thorough understanding of the dynamics and outcomes. For example, on long timescales we find the that the two-compartment model is equivalent to the well-studied Moran model. Meanwhile, analysis of the reconstitution of a preconditioned mouse shows that in our model HSCs migrate multiple times between the BM and PB compartments before establishing a sustainable population. Given these dynamics we first investigate clonal dominance, where a clone originating from a single mutant cell expands in the HSC population. In mice we find that a selective advantage is required if the clone is to be detected within a lifetime: A clone starting from a single cell with a reproduction rate 50% higher than the wildtype can expand to 1% clonality in one year. A cell dividing twice as fast as the wildtype reaches > 10% clonality in the same timeframe. Such division rates can be reached by MPN-initiating HSCs [55]. The requirement of a (reconstitution), and we assume no further extinction events occur once this upper limit has been reached. Dotted lines are the 'first-order' prediction of Eq (10). Solid lines are the predictions of Eq (11) which account for detachments, reattachments, and reproduction events. Remaining parameters are as in Table 1.
The model also captures the scenario of stem cell transplantation. Engraftment into a nonpreconditioned host is analogous to clonal dominance, except that the clone is initiated by multiple donor cells. For small doses of donor HSCs, the number of cells that attach to the BM is directly proportional to the size of the dose. For larger doses the BM niches are saturated, leading to lower engraftment efficiency. Donor chimerism can be improved by injecting the host with multiple small doses, as opposed to a large single-bolus dose of the same size. This agrees with results that have been reported in the empirical literature [19,[48][49][50]. Following preconditioning of a mouse to remove all host cells, we find that a single donor HSC is sufficient to repopulate a host in *90-99% of cases. This result rests on the assumption that the donor stem cell was, in fact, a long-term reconstituting HSC, which may not be the case in experimental setups.
In the S1 Supporting Information we consider the effects of death occurring in the BM niche (δ 0 6 ¼ 0), and the direct attachment of a new daughter cell to the bone marrow niche (ρ 6 ¼ 0). We find that death in the niche increases the migration rate of cells between the PB and BM compartments, which can greatly reduce the attachment success of the low-frequency mutant/donor cells. However, the direct attachment of daughter cells to the niche has no effect on the initial attachment of donor/mutant cells, and on the level of chimerism achieved in the initial phase of the dynamics.
Broadening the scope of our investigation, clonality of the hematopoietic system is a major concern for human health [13][14][15]39]. Clinical studies have shown that 10% of people over 65 years of age display clonality, yet 42% of those developing hematologic cancer displayed clonality prior to diagnosis [13]. Our model, and the subsequent analysis, can be applied to this scenario. However, the number of HSCs in man is debated, with estimates of *400 [20,58], Oð10 4 Þ [18], or Oð10 7 Þ [5]. Estimates as high as Oð10 9 Þ can also be obtained by combining the total number of nucleated bone marrow cells [59] with stem cell fraction measurements [1,60,61]. In Fig 6 we summarise how neutral and advantageous clones starting from a single HSC expand in human hematopoietic systems. We find that 4% clonality [39,40] can be achieved in a short period of time for even neutral clones [Fig 6(a)]. If the human HSC pool is Oð10 3 Þ or smaller, we would expect clonal hematopoiesis and the associated malignancies to be highly abundant in the population, perhaps more-so than they currently are [13][14][15]39]. On the other hand, for a system size of N = 10 6 it takes thousands of years for a single neutral HSC to expand to detectable levels, making neutral expansion extremely unlikely to result in clonal hematopoiesis. Therefore, for clonal hematopoiesis to occur in a pool of this size or larger [5] the mutants would require a significant fitness advantage over the wildtype HSCs. We also consider a range of parameters, and even relax the α = % = 0 condition, in S1 and S2 Figs. We find no significant differences in the predictions of our model.

Limitations
Our model has been kept to a minimal level of biological detail to allow for parametrisation from experimental results. This has the added benefit of analytic tractability. The model is constructed under steady-state conditions, which is the case for neutral clonal expansion. However, in the case of donor-cell transplantation following myeloablative preconditioning, we are no longer in a steady state. Here we expect some regulatory mechanisms to affect the HSC dynamics, including a faster reproductive rate and a reduced probability of cells detaching from the niche. There are also possibilities for mutants to exploit or evade the homeostatic mechanisms [63]. Different mechanisms of stem cell control have recenty been considered for hematopoietic cells [64], as well as in colonic crypts [65].
The steady state assumption is also unable to capture the different dynamics associated with ageing. For example, in young individuals the hematopoietic system is undergoing expansion. In our model there is no distinction between young and old systems. In S3 Fig we demonstrate the impact of a (logistically) growing number of niches. Such growth means clonal hematopoiesis is likely to be detected earlier, and therefore would increase our lower bound estimate on the number of HSCs in man. Telomere-length distributions have been used to infer the HSC dynamics from adolescence to adulthood, and have suggested a slowing down of HSC divisions as life progresses [66]. Faster dynamics in early life would lead to a higher incidence among young people, which again increases our lower bound estimate.
It is also not entirely clear how to extrapolate the parameters from the reported mouse data to a human system. Here we have taken the simplest approach and appropriately scaled the unknown parameters. However, hematopietic behaviour may differ between species. For example, results of HSC transplantation following myeloablative therapy in non-human primates have shown that clones of hematopoietic cells persist for many years [67,68]. This could  (6), and shaded regions represent the calculated standard deviation (details in the S1 Supporting Information). Remaining parameters are β = 1/ 40 week −1 [62], ℓ = 60 minutes, s* = 0.01N, and n* = 0.99N. Here ℓ, s*, and n* are extrapolated from the murine data, where ℓ human % 10ℓ mouse , which follows the same scaling as the HSC division rate, β. Further parameter combinations are shown in S1 and S2 Figs. References refer only to the source of parameters; no part of this figure has been reproduced from previous works.
https://doi.org/10.1371/journal.pcbi.1005803.g006 be due to single HSCs remaining attached to the niche and over-contributing to the hematopoietic system, or due to clonal expansion of the HSCs to large enough numbers such that a contributing fraction will always be found in the BM. Both of these mechanisms are features of our model: the time a cell spends in the BM is much longer than the time in the PB and can be increased further by tuning the model parameters, namely by decreasing s Ã or increasing ℓ. Changes to these parameters seems to have little effect on our predictions of clonal expansion, as shown in S1 and S2 Figs. Clonal extinctions are also a feature of our work, and have been identified in non-human primates [68].
A more general point to discuss is the role of hematopoietic stem cells in blood production. In our model we are only considering HSC dynamics, however it has been proposed that downstream progenitor cells are responsible for maintaining hematopoiesis [11] in mice. Hence, myeloid clonality would also be determined by the behaviour of these progenitor cells. On the other hand, an independent study found that HSCs are driving multi-lineage hematopoiesis [10], suggesting we are correct in our approach. Again we also expect there to be variation between species in this balance of HSC/progenitor activity. With little quantitative information available, we have assumed that HSCs are the driving force of steady-state hematopoiesis across mice and humans.

Conclusion
In conclusion, this simple mathematical model encompasses multiple HSC-engraftment scenarios and qualitatively captures empirically observed effects. The mathematical calculations provide insight into how the dynamics of the model unfold. The analytical results, which we have verified against stochastic simulations, allow us to easily investigate how parameter variation affects the outcome. We now hope to extend this analysis, incorporating further effects of disease and combining this model with the differentiation tree of hematopoietic cells.
Supporting information S1 Supporting Information. Supporting mathematical details. Contains detailed derivations of all equations presented in the manuscript, including the details of the projection method. The analysis is carried out for unrestricted parameters, including selective effects on all parameters and permitting death in the BM space, as well as direct attachment of new daughter cells to the niche. (PDF) S1 Fig. Clonality in man: More parameter combinations and death within niches. Time taken until a clone initiated from a single cell represents 4% [39,40] of the human HSC pool, as a function of the total number of niches in the system. Colours represent the selective advantage of the invading clone. Solid lines correspond to death only within the niches (α = 0), while dashed lines represent equal death rates in both compartments (α = 1; see S1 Supporting Information for details). Lines are generated using mathematical formulae in the S1 Supporting Information. Remaining parameters are β = 1/40 week −1 [62], n Ã = 0.99N, and % = 0. Some predictions are missing when d 0 and/or a 0; these parameter regimes are incompatible with our model.  [39,40] of the human HSC pool, as a function of the total number of niches in the system. Colours represent the selective advantage of the invading clone. Solid lines correspond to the daughter cell entering the PB compartment after reproduction (% = 0), while dashed lines represent daughter cells remaining in the BM (% = 1; see S1 Supporting Information for details). Lines are generated using mathematical formulae in the S1 Supporting Information. Remaining parameters are β = 1/40 week −1 [62], n Ã = 0.99N, and α = 0. Some predictions are missing when d 0 and/or a 0; these parameter regimes are incompatible with our model. (TIFF) S3 Fig. Predicted and simulated incidence curves of clonal hematopoiesis in man: Constant size and under growth. The cumulative probability density function (CDF) of times to reach 4% clonality [39,40] when starting from a single neutral mutant in a normal host. Shaded regions are incidence curves from simulations using either a constant niche count of N = 10 4 , or a logistically growing number of niches with _ N % rNð1 À N=KÞ, where K = 10 4 , r = 0.3 per year, and N(0) = K/20. These parameters represent a maturation period of *20 years to reach N % K. Lines are predicted incidence curves which assume normally-distributed times to clonality, using the mean and variance formulae as described in the S1 Supporting Information, and constant population size as indicated in the legend (minimum and maximum number of niches). Remaining parameters are β = 1/40 week −1 [62], n Ã = 0.99N, s Ã = 0.01N, and ℓ = 60 minutes. Finally, we only consider here the conditional incidence time, which have been normalised by the fixation probability. This probability is 20 times larger for the neutral mutant in the growing model when compared to the fixed number of niches. (TIFF)