Modeling large fluctuations of thousands of clones during hematopoiesis: The role of stem cell self-renewal and bursty progenitor dynamics in rhesus macaque

In a recent clone-tracking experiment, millions of uniquely tagged hematopoietic stem cells (HSCs) and progenitor cells were autologously transplanted into rhesus macaques and peripheral blood containing thousands of tags were sampled and sequenced over 14 years to quantify the abundance of hundreds to thousands of tags or “clones.” Two major puzzles of the data have been observed: consistent differences and massive temporal fluctuations of clone populations. The large sample-to-sample variability can lead clones to occasionally go “extinct” but “resurrect” themselves in subsequent samples. Although heterogeneity in HSC differentiation rates, potentially due to tagging, and random sampling of the animals’ blood and cellular demographic stochasticity might be invoked to explain these features, we show that random sampling cannot explain the magnitude of the temporal fluctuations. Moreover, we show through simpler neutral mechanistic and statistical models of hematopoiesis of tagged cells that a broad distribution in clone sizes can arise from stochastic HSC self-renewal instead of tag-induced heterogeneity. The very large clone population fluctuations that often lead to extinctions and resurrections can be naturally explained by a generation-limited proliferation constraint on the progenitor cells. This constraint leads to bursty cell population dynamics underlying the large temporal fluctuations. We analyzed experimental clone abundance data using a new statistic that counts clonal disappearances and provided least-squares estimates of two key model parameters in our model, the total HSC differentiation rate and the maximum number of progenitor-cell divisions.

Using forms given in Eq. (A3), since both φ and ψ are independent of h, we can define Thus, the probability distribution P (h, t) can be written as (A5)

B Alternative model of progenitor aging
An alternative model to the one we have analyzed allows younger-generation progenitor cells (ℓ < L) to differentiate into peripheral blood. Since each generation can differentiate with rate ω, the progenitor cell dynamics are slightly modified from those in our main model: Poisson(αh(t)) − (r n + µ n + ω)n (0) (t), ℓ = 0, 2r n n (ℓ−1) (t) − (r n + µ n + ω)n (ℓ) (t), Moreover, the dynamics of the mature peripheral blood population obey The solution to Eqs. (B1) and (B2) following a single differentiation event is These results can be applied to the model and analyzed and simulated using the same procedures as described in the main text. However, certain parameters have to be re-interpreted. For example, using the same value of ω = 0.16 will significantly increase the effective death rate for progenitor cells of each generation. Fortunately, as we will show later, this alternative mechanism should not affect our main conclusion as the parameter-fitting results are not sensitive to the exact shape of cell bursts.

C Mean extinction time for a clone
As a function of the initial number h of HSCs in a clone, the mean extinction time (MET) T (h) under the steady-state approximation r h = µ h obeys [3,4] [ with an absorbing boundary condition T (0) = 0. By iterating Eq. (C1), we find which can be again iterated to obtain To solve for T (1), we invoke a reflecting boundary condition T (H ss ) − T (H ss − 1) = 1/(µ h H ss ) [5], where to find Upon using Eq. (C5) in Eq. (C3), we find which is the MET for a discrete system.
We can also approximate T (h) by considering h as a continuous variable, and replace the summations in Eq. (C6) by integrations to find a simpler, more insightful approximation to T (h): where we have used The continuous approximation to the MET matches the exact result quite well (relative error 5%) for all values of h.

D Effective parameters and symmetric HSC differentiation
There are differing reports on the measured death rates for circulating granulocytes. We have used the most recently reported value µ m = 1 per day for humans. The effect of changing the value of µ m → µ ′ m on our analysis is a reinterpretation of L e . By rewriting Eq. (13) as A + where one additional round of progenitor doubling compensates for the doubled loss rate of mature granulocytes. One may argue that the change in µ m can also be compensated for by doubling A + ss , which would have a different effect on the burstiness of the model compared to doubling L e . However, when re-fitting the data with µ ′ m = 2 or 0.2, we observed that (A + ss ) * did not change much, with most of the effect of modifying µ m absorbed by changes in L * e . Similarly, uncertainties in other parameters can also be subsumed into L e . For example, setting µ (L) n = ω > 0 implies that only half of the generation-L progenitors contribute to the peripheral blood. For a model with µ (L) n = 0 to generate an equivalent effect, we can halve the number of mature cells by using an effective maximum generation parameter L ′ e = L e − 1. This indicates that the intrinsic clone size fluctuations demonstrated in the experimental data strongly constrain A + ss . Another possible modification of our mechanistic model is to allow for the possibility of symmetric HSC differentiation. The effect of symmetric differentiation can again be subsumed into the parameter L e without qualitatively affecting our analysis. Assume a proportion 0 ≤ q ≤ 1 of HSC differentiations are symmetric, producing on average 1+q generation-0 progenitor cells. After L e rounds of proliferation, the 1+q generation-0 progenitors produce on average (1 + q) × 2 Le mature cells. This is equivalent to an exclusively asymmetric differentiation model (q = 0) with L ′ e = L e + log 2 (q + 1). We also expect symmetric differentiation to slightly increase the speed of coarsening since each HSC differentiation is also accompanied by the HSC's death and clones represented by a single HSC would disappear under symmetric differentiation. However, given the small rate α of HSC differentiation, the large number C h of clones, and the insensitivity of our results to the distribution h i , the data cannot quantitatively resolve the symmetric-asymmetric modes of HSC differentiation.

E Alternative objective functions and statistical insights
We developed our data analysis based on the statistics of the quantity y i , the time averaged relative clone sizes for those clones exhibiting z absences across their longitudinal samples. While reasonable parameter estimates were obtained from fitting to data, we also considered alternative objective functions. Specifically, we looked at the standard deviation σ i = 1 J J j=1 (f i (t j ) − y i ) 2 quantifying the temporal fluctuations of the relative sizes of each clone i. The way we construct an alternative objective function is similar to the way we constructed Y z . Recall for Y z , we calculated the average abundance across only those clones with the same z i = z absences across time. However, unlike z i which takes a finite set of discrete values {1, 2, ..., J − 1}, σ i is a continuous variable so we have to artificially bin their values. Instead, we bin clones with similar y i and study the average of their associated σ i 's. Since the distribution y i is non-linear with a long tail, we evaluated ln y i to obtain the near-linear distribution shown in Fig. E1(a), sorted ln y i into equal-width bins, and calculated the average of the associated σ i s. Dividing the values of ln y i into bins labeled by k, we compute in analogy with the definition of Y z . The objective function can be straightforwardly defined as It is also unclear how to set upper and lower bounds on the range of y i for comparison (in contrast to the natural bound on 1 ≤ z ≤ J − 1) because an unconstrained set of clones will be sensitive to the underlying h i distribution (an undesirable property). In Fig. E1(b), we fit the data from animal RQ5427 using MSE σ and find L * e ≈ 24.4, consistent with our previous estimate using Y z . While it is also possible to choose σ i as a measure of clone population fluctuations, we list several advantages ofẑ i over σ i for the current dataset. Note that the number of disappearances z i of each individual clone is defined on a finite set of integers (unlike the continuously measured σ i ), making it easier to bin clones with the same z values. Different clones i will exhibit different time-averaged abundances y i but may have the same value of z i . As shown in Fig. 4 in the main text, the largerẑ i is, the smaller the corresponding lnŷ i tends to be. The robust correlation between z i and y i encodes the level of fluctuations for a clone of certain size. For a given y i , the larger z i , the "burstier" the dynamics, implying a smaller number of tagged HSC differentiations per unit time (a smaller A + ss ). Another advantage of using z i statistics emerges when fitting model results to the pattern of the measured data in Fig. 4 in the main text. Average sizes y i (and the underlying h i ) associated with clones having 1 ≤ z ≤ 7 all contain at least one absence. This constraint naturally controls the upper and lower bounds of h i in a particular z bin (1 ≤ z ≤ 7), based on the burstiness of the model. Exact knowledge of the configuration {h i } is not required for fitting these y i data.
Thus, dividing clones into z bins provides us with a natural way to exclude unconstrained clone sizes. In other words, the theoretical values of y i (and the underlying h i ) associated with bin z i = 0 can be arbitrarily and unreasonably large, and such a possibility should be excluded. Similarly, all y i below a threshold size generate z i = J (clones that never appeared in the sampled blood) and do not provide any statistical power. This advantage of using z i can also be confirmed by visual inspection of Fig. 9(b) in the main text. Several very large clones do not follow the general statistical pattern and show extremely large variances. Without manually filtering out these clones, our fitting in Fig. 1(b) results in a larger L * e = 24.4 than the L * e = 23.4 obtained in the main text using Y z statistics. Finally, another option for comparing model with data is to use correlation functions. In this approach, the sampling gap ∆t j varies between 5 and 11 months, so the usual autocorrelation function with equal time gaps cannot be rigorously defined. We use the one-sample-gap autocorrelation function and bin values of ln y i in analogy to Eq. (E1) to define and construct an autocorrelation-based objective function Since the inter-sample intervals ∆t j are larger than a typical burst size ∆τ b ≈ 32 days, cells in different samples likely originate from different HSC differentiation events. Thus, the fluctuations of clone sizes are uncorrelated from sample to sample, as shown in Fig. E1(c). Randomly distributed between -1 and 1, the values of R i are centered about the line R = 1 2−J , corresponding to the majority of clones that have z i = J −1 (only 1 non-zero sample). Data fitting using R i and MSE R is ill-conditioned and cannot resolve L * e , as shown in Fig. E1(d).

F Simulation of the forward model
To generate predictions, we first choose values of θ model = {λ, C h , r n , L e } and simulate our model, including sampling, to find s i (t j ). To simulate each realization of our model we

Sample a fraction ε(t j ) =Ŝ
M + (tj ) of the total peripheral cell count M + (t j ) = i m i (t j ). Here, S + (t j ),M + (t j ), and the times t j are defined by the experiment. We used the Python package numpy.random.binomial. The cell counts of each clone are s i (t j ). Use the simulated total tagged cell counts in the samples S + (t j ) = i s i (t j ) to normalize si(tj ) S + (tj ) = f i (t j ). Up to this point, we have generated a data matrix f i (t j ) of size C h × J. The simulated, model-derived configurations f i (t j ) are then compared with experimentally measured valuesf i (t j ). The parameter L e that minimizes the mean-squared error will be chosen as the least-squares estimate L * e .

G Robustness to samping frequency and threshold
The robustness of our inference of L * e to sampling frequency is demonstrated for animal RQ5427 by excluding some time samples. In Figs. G1(a-h), we plot the MSE function by including only the first j = (8, 7, . . . , 1) time samples of the data. In this data set (animal RQ5427), the MSE remains meaningful, and the reconstruction of L * e is unchanged as long as at least four or five time samples are used. This conclusion is independent of which sampling time points are excluded. Since the system is well-approximated by a statistical steady state, the key determinant for robust inference is the number of samples included in the analysis. Robustness to a larger threshold of clone sizes is also demonstrated by eliminating clones whose average abundances are under a certain threshold in both the experimental and simulated data. In Figs. G2(a-h), we plot the MSE corresponding to the clone frequency thresholds 1.16 × 10 −5 , 2.03 −5 , 3.41 × 10 −5 , 8.84 × 10 −5 , 1.66 × 10 −4 , 3.30 × 10 −4 , 6.78 × 10 −4 , 1.46 × 10 −3 , respectively. Using these thresholds, the numbers of clones retained in the analysis are 482, 428, 375, 322, 268, 215, 159, and 107, corresponding to 90%, 80%, 70%, 60%, 50%, 40%, 30%, and 20% of the 536 total number of clones detected in animal RQ5427. Figs. G2 show that as long as 200 clones are included (a-f), the MSE yields a clear LSE L * e = 23.4 − 23.6. Only at very high thresholds, where only 20-30% of the clones are retained, does the minimum of the MSE shift to slightly higher values L * e ≈ 23.8, 24.3 as shown in Figs. G2(g-h), respectively. Thus, we conclude that the inference of L * e from the data is fairly insensitive to sampling threshold provided a reasonable number of clones (typically 200) are included in the analysis.   Figure G2: MSEs for animal RQ5427 using successively higher clone detection thresholds. Unique reconstruction of L * e is robust (a-f) even if only 40-50% of the clones are counted. (g-h) At even higher thresholds, the LSE for L * e increases only very slightly.