First, we consider the observed data associated with each animal (Koelle et al., 2017), as shown in Figure 1. Granulocytes in blood samples drawn from each animal at times points t _j , j = 1, 2, …, J are sequenced and clonal (barcode) abundances tabulated. The total abundance (number of mature cells of a given cell type), S ̂ ( t j ) , and the richness C ̂ s ( t j ) (the total number of different barcodes detected in each sample) are also recorded and plotted in Figures 1A, B. In this study, S ̂ ( t j ) denotes the total measured number of granulocytes (barcoded and unbarcoded) in a sample taken at time t _j . The fluctuations of S ̂ ( t j ) and C ̂ s ( t j ) across t _j may arise from varying sampled sizes across time points and/or fluctuations in the state of the animal.

An example of the abundances of each clone within the granulocyte population from Koelle et al. (Koelle et al., 2017) is shown in Figure 1C. In these experiments, tagged stem cells are transplanted back into a rhesus macaque at t = 0 so that initially each clone consists of a single cell. A series of 1 ≤ j ≤ J samples are taken at time t _j after implementation, yielding a set of mature cells. We denote the abundance of clone i (among granulocyte cells) in the sample taken at time t _j after transplantation as s ̂ i ( t j ) . The J measurements allow each clone of a particular mature cell type i to be characterized by a mean s ̂ i and variance σ ̂ i 2 defined by s ̂ i = 1 J ∑ j = 1 J s ̂ i t j σ ̂ i 2 = 1 J ∑ j = 1 J s ̂ i t j − s ̂ i 2 . (1)

Note that the total measured population of any cell type S ̂ ( t j ) = ∑ i = 1 C s s ̂ i ( t j ) . A scatter plot of σ ̂ i versus s ̂ i for all clones detected in a sample of granulocytes is shown in Figure 1D.

For each clone at ( s ̂ i , σ ̂ i ) we can evaluate the local density ρ ̂ , the number of clones within some size window. This density can be viewed as the concentration of data points shown in Figure 1D as a function of mean clone size s ̂ , and will be constructed using kernel density estimation (Rosenblatt, 1956; Parzen, 1962) of the data points in ( s ̂ , σ ̂ ) space. The unknown density function ρ ̂ is obtained by concatenating isotropic Gaussian kernel functions about each point and using an optimal, common bandwidth parameter, typically chosen as the value that minimizes the mean integrated squared error, or Kernel Density Estimation (KDE) (Silverman, 1986). The reconstructed density function can be thought of as a probability that a random clone arises in the volume (s, s + ds) × (σ, σ + dσ). For each clone at ( s ̂ i , σ ̂ i ) we can evaluate the local density ρ ̂ .

In the remainder of this paper we will develop a mathematical model that we can simulate to generate total populations, clonal populations, and their associated attributes (s _i , σ _i , ρ _i ). Note that while there are many existing mathematical models of hematopoiesis (Colijn and Mackey, 2005; Peixoto et al., 2011; De Souza and Humphries, 2019), they describe only time variations in total populations, rather than that of lower-population, individual clones. We will tune parameters of the model so that its predictions provide a reasonable match to the aforementioned measurements, paying particular attention to clone abundances and clone size variability.

Our mathematical model incorporates known and accepted features of hematopoiesis. Three main cell compartments are considered: hematopoietic stem cells (HSCs), transit amplifying progenitor cells, and peripheral mature cells. Although the stem cell population in bone marrow is large and can be described using a deterministic model, the initial populations within each clone are small and require a discrete stochastic description. We will then assume that a small sample of mature cells is drawn from the animals at times t _j , j = 1, 2, …, J and sequenced.

We first describe the initial conditions including the number of HSCs and HSPCs injected into each animal. As listed in Table 1, H is the total number of HSCs injected into each animal, among which H ₀ are untagged and H _i ≪ H ₀ contain barcode 1 ≤ i ≤ C _H (and are GFP+). The total tagged HSC population is H * ≡ ∑ i = 1 C H H i so that H = H ₀ + H*. The richness C _H ≲ C _L is the number of barcodes transferred into the animal, which is comparable to the richness of the barcode library C _L used in each experiment. Since H ₀ ∼ 10⁷ ≫ H _i , we will consider the probability distribution of only the tagged populations, which is described by the multinomial P H = H * ! ∏ i = 1 C H 1 C H H i 1 H i ! , (2) where H ≡ ( H 1 , H 2 , … , H C H ) and H* is the total number of GFP+ (barcoded) cells. Specifically, for animal ZH33 studied in (Koelle et al., 2017)) H ≈ 3 × 10⁷, H*/H ≈ 0.35, ∑ i = 1 C H H i ≈ 1.1 × 1 0 7 , C _H ≲ C _L ≈ 6 × 10⁵. Thus, the typical H _i ≈ H*/C _H ∼ 180.

A certain fraction η ₀ of the H HSCs home into the bone marrow, successfully engraft, and subsequently actively self-renewal and/or differentiate. Engrafted HSC populations are defined by h ( 0 ) = ( h 1 ( 0 ) , h 2 ( 0 ) , … , h C h ( 0 ) ( 0 ) ) , where the richness of engrafted HSCs in the bone marrow is C _h ≲ C _H . Transplantation efficiencies are typically single-digit percentages (Abbuehl et al., 2017; Radtke et al., 2020) and transplanted CD34⁺ cells contain significant numbers of progenitor cells. Thus, the fraction η ₀ ≪ 1; if η ₀ is sufficiently small (approximately ≲ 1/180), then we can safely assume that the initial clone populations in bone marrow are represented by very few cells. For simplicity, we approximate h _i (0) ≈ 1. Even if η ₀≰1/180, most barcodes will be represented by very few cells. We have verified that an initial condition in our model that allows for, say, some h _i (0) = 2, 3 does not qualitatively affect the mature cell populations.

The random selection of cells into the bone marrow can be thought of as a sampling (without replacement) process. Including the untagged population, the probability distribution of engrafted cells resulting from the injected tagged cell population H = ( H 1 , H 2 , … H C H ) is given by P h 0 | H = 1 H * h * 0 ∏ i = 1 C H H i h i 0 , (3) where H* and h * ( 0 ) = ∑ i = 1 C H h i ( 0 ) are the total initial numbers of barcoded injected cells and engrafted barcoded HSCs, summed over all clones. Note that the number of untagged transplanted cells h ₀(0) ≫ 1 is large so that we can approximate it by its deterministic value h ₀(0) ≈ η ₀ H ₀.

To extract the overall probability of initial condition h(0), we average Eq. 3 over the prior P ( H 1 , … , H C H ) and find P h 0 = ∑ H P h 0 | H P H = h * 0 ! ∏ i = 1 C H 1 C H h i 0 1 h i 0 ! , h * 0 ≡ ∑ i = 1 C H h i 0 . (4)

Besides the initial condition h _i (t = 0) = 1, the initial number of untagged HSCs h ₀(t = 0) is related to the transplantation efficiency and is generally unknown. Barcodes are associated with a GFP tag and the initial fraction of sampled cells that are GFP+ is ∼ 35 % . Since we assume a neutral model, it is reasonable to assume that the fraction of injected tagged cells H*/H is equivalent to the fraction of tagged cells in the engrafted population ∑ i = 1 C h ( 0 ) h i ( 0 ) / ∑ i = 0 C h ( 0 ) h i ( 0 ) ≈ 0.35 (although this ratio slowly decreases via extinction). The precise richness of HSC population in stem cell niche, C _h(t) < C _H is also unknown, but except for fluctuations, has a lower bound of C ̂ s , the total number of unique clones detected across all samples across all cell types. Thus, we take h ( 0 ) = ∑ i = 0 C h ( 0 ) h i ( 0 ) = h 0 ( 0 ) + C h ( 0 ) ≈ C h ( 0 ) / 0.35 .

Self-renewal, death, and differentiation into progenitor cells all contribute to the stochastic dynamics of h _i . Although the total HSC population in the niche h ( t ) ≈ ∑ i = 0 C h ( t ) h i ( t ) is large and can be approximated deterministically, the HSC population of each clone h _i (t) may be small and must be treated stochastically. Under our neutral assumption, the intrinsic self-renewal rate r _h of HSCs does not depend on the barcode identity i. Since HSCs reside in niches in the bone marrow that place limits on growth, we assume the HSC proliferation rate follows a linearly decreasing form defined by a carrying capacity and the engrafted HSC population h(t) r h h t = r h 0 1 − h t K , h t ≡ ∑ i = 0 C h h i t , (5) where r _h(0) is the intrinsic proliferation rate of a single, isolated HSC. Note that the untagged HSCs are included through h ₀. Finally, we assume that HSCs die at rate μ _h and differentiate at rate α and that these rates, like the growth rate in Eq. 5, do not depend on barcode identity.

As shown in the Supplementary Material, the richness C _h(t) may progressively decrease from random HSC death and extinction and can be estimated by solving for the stochastic birth-death process (neglecting outflux from differentiation) and using generating functions to find E C h t ≈ C h 0 ψ t + ϕ t , (6) where ψ t ≡ e − ∫ 0 t r h t ′ − μ h d t ′ , ϕ t ≡ ∫ 0 t r h t ′ ψ t ′ d t ′ . (7)

In this expression, r _h(t) is approximated by r h ( h ̄ 0 ( t ) + h ̄ * ( t ) ) , where h ̄ 0 + h ̄ * ( t ) is given by the explicit solution to the deterministic birth-death process with carrying capacity K. Thus, given C _h(0), r _h(0), K, and μ _h determine how the expected richness decreases. Henceforth we treat C _h(t) in our model as the expected value E [ C h ( t ) ] derived from our stochastic birth-death model, i.e., we use C _h(0)/(ψ(t) + ϕ(t)) as the model for C _h(t).

We will also simulate the stochastic birth-death process for HSCs (see Supplementary Material for details), with the differentiation rate α that allows HSCs to differentiate into the progenitor/transit amplifying cell compartment (see Figure 2). Progenitor cells are further distinguished by their generation ℓ. Thus, ℓ not only measures generation number but also an effective differentiation state. Each HSC differentiation event leads to an ℓ = 0 progenitor cell. Since the number of HSCs of any one clone is small, the initial differentiation events from clone i follow a Poisson process with rate αh _i . The populations of the subsequent generations of progenitor cells in clone i are denoted by n i ( ℓ ) . Once the ℓ = 0 generation cells are generated, the number of progenitor cells quickly expand, so their dynamics will be described by a deterministic model as developed in (Xu et al., 2018) d n i ℓ t d t = Poisson α h i t − r n 0 + μ n 0 n i 0 t ℓ = 0 , 2 r n ℓ − 1 n i ℓ − 1 t − r n ℓ + μ n ℓ n i ℓ t 1 ≤ ℓ ≤ L − 1 , 2 r n L − 1 n i L − 1 t − ω + μ n L n i L t ℓ = L , (8) where Poisson(αh _i (t)) is the time-inhomogeneous point Poisson process describing HSC differentiation events. In other words, after a differentiation event at time t ₁, the probability density of the time Δt to the next differentiation event is given by α h i ( t 1 + Δ t ) exp − α ∫ 0 Δ t h i ( t 1 + s ) d s . We will use the values of h _i (t) from our stochastic simulations and sample from this inter-event time density to simulate realizations of differentiation events.

In Eq. 8, r n ( ℓ ) and μ n ( ℓ ) represent the proliferation and death rates of generation-ℓ progenitor cells, respectively, and ω is the terminal differentiation rate into mature blood. Each division can also be thought of symmetric differentiation producing successively more differentiated progenitor cells. To model the finite proliferative potential of progenitor cells, we set the maximum number of generations to ℓ = L, after which the L ^th generation cell can only terminally differentiate to mature blood.

Consider a single isolated differentiation event of an HSC at t = 0 belonging to a particular clone. The resulting progenitor cell population after this event is described by Eq. 8 without the Poisson(αh _i (t)) term but with an initial condition corresponding to a single ℓ = 0 cell: n i ℓ 0 = 1 ℓ = 0 , 0 otherwise. (9)

The subsequent populations at time t form a temporal “burst” of cells that are described by n i ( ℓ ) ( t ) which is the solution of Eq. 8 without the Poisson(αh _i (t)) term but using the initial condition in Eq. 9. If we assume that all progenitor generations carry the same division and death rates, r n ( ℓ ) = r n for 0 ≤ ℓ ≤ L − 1 and μ n ( ℓ ) = μ n for 0 ≤ ℓ ≤ L, we can find an analytical solution associated with a single isolated burst as n i L t = e − ω + μ n t L − 1 ! 2 r n r n − ω L ∫ 0 r n − ω t z L − 1 e − z d z . (10) We can evaluate all populations n i ( ℓ ) ( t ) for ℓ < L by solving Eq. 8 and using Eq. 10, as detailed in the Supplementary Material.

If we assume that mature cells do not appreciably proliferate ¹ , the mature cell population in clone i obeys d m i t d t = ω n i L t − μ m m i t , (11) where μ _m is the lineage-dependent turnover rate of mature cells. Using the solution to n i ( L ) ( t ) , we solve Eq. 11 to find (see Supplementary Material) m i t = ω ∫ 0 t n i L t ′ e − μ m t − t ′ d t ′ . (12)

The mature cell population burst (of a specific clone) arising from a single, isolated differentiation event is plotted in Figure 3A. Note that the expression for a mature cell burst given in Eq. 12 is derived from the specific initial condition Eq. 9; however, some low-ℓ progenitors are also initially transplanted (see below). Thus, m _i (t) will in general depend on the initial numbers of n ^(ℓ>0)(0).

Since Eq. 8 are linear, populations arising from a sequence of Poisson-distributed differentiation events can be constructed by adding those derived from single events occurring at times T _k . In this case, the resulting mature cell population at time t is given by m i t = ∑ k = 1 k max m i t − T k , T k max + 1 > t > T k max , (13) where m _i (t − T _k ) is the solution given in Eq. 12 (for the specific initial condition in Eq. 9). Two different sequences of bursts are shown in Figures 3B, C. In Figure 3B, we consider a clone with few HSCs such that αh _i (t) ≪ μ _m. This limit gives rise to differentiation events that occur rarely over the lifespan of mature cells, as depicted Figure 3B. In Figure 3C, we plot a sequence of more frequent mature cell population bursts that arise for more frequent differentiation events αh _i (t) ≫ μ _m from HSCs that are in higher population clones.

Recall that the calculations involving Eqs 8–12 are performed for each clone i, resulting in a series of time-dependent expressions for m _i (t) as per Eq. 12. These single event responses are then summed according to Eq. 13 to arrive at the total, time-dependent population of mature cells of a specific type and carrying the same barcode. These predicted whole-organism populations depend on the model parameters C h ( 0 ) , α , K , r n , μ h , μ n , μ m , L , ω . Since growth of transit amplifying progenitor cells is fast, we will henceforth assume r _n ≫ μ _n ≈ 0. All other rates are given in units of per day.

Since the clone abundances are derived from sequencing cells in a small fraction η (for rhesus macaque, η ∼ 10^–5 − 10^–4) of the animal’s blood taken at times t _j , we also need expressions for mature cell populations within a sample. Given the small sample sizes η, low population clones in the mature cell pool can easily be missed. For a given population m _i in the whole animal, the probability that s _i cells are captured in the sample is given by (Chao and Lin, 2012; Xu et al., 2020) P s i | m i , S , M = 1 M S ∏ i = 0 C h m i s i ≈ ∏ i = 0 m i s i η s i 1 − η m i − s i , (14) where, for a given cell type (for example granulocytes), S(t) = ∑ _i=0 s _i (t) is the total number of sampled cells (including untagged ones), M(t) = ∑ _i=0 m _i (t) is the total number of circulating mature cells (including untagged ones), and the sampling fraction is η = S/M (which we first assume is the same at each t _j ).

After computing m _i at time t _j using Eq. 12, we take the nearest integer value and use it for m _i in Eq. 14. We then draw a single value s _i (t _j ) from the binomial distribution, assuming η is given. Finally, we simulate our model to generate trajectories of M(t) and then determine the tagged sampled fraction S t ≈ H * H η M t . (15)

Equations 5–15 represent a hybrid stochastic-deterministic model since the self-renewal process of a small number of HSCs in each clone and the final sampling step (Eq. 14) are modeled as discrete stochastic processes, while proliferation and differentiation of higher population progenitor and mature cells are treated deterministically via Eqs 8–12. The values s _i (t _j ) obtained from the model in Eqs 8–14 are used to generate the predicted clone population mean and standard deviations, s _i and σ _i , according to Eq. 1. Both s _i , σ _i are then used to determine the density of points ρ _i . These three values for each clone are then compared to their corresponding values constructed from data s ̂ i ( t j ) , as we detail in the next section. These predictions, along with the predicted richness C _s(t) and total sampled granulocyte population S(t) provide the basis for comparing with measured data and parameter inference.

The published experimental data in (Koelle et al., 2017) provide data for four rhesus macaques ZH33, ZG66, ZH19, and ZJ31 over a period of 49, 42, 36 and 20 months, respectively. Samples were taken after autologous transplantation with myeloablative conditioning at times t _j , 1 ≤ j ≤ J. Possible sampled cells are of five types: T cells, B cells, monocytes, granulocytes and NK cells. Barcoded myeloid (granulocytes and monocytes) and B-cells reach a noisy equilibrium after approximately 1 month, whereas for T-cells the time frame is longer, between 5 to 17 months. Furthermore, since granulocytes comprise a majority of white blood cells, we apply our mathematical and statistical model to the granulocyte lineage. At each sampling time t _j , the experimental data from (Koelle et al., 2017) reveals how many cells are sampled from each clone and what type each cell is. Across sampling time points, these sampled populations contain information on the overall abundance, how these abundances fluctuate in time, and the density of the number of clones detected as a function of abundance and abundance variability. Since ZH33 has the longest follow-up period, we use data from this macaque to compare experimental data with our mathematical predictions.

Validation of our model will rely on matching predictions with available data in the form of C ̂ s ( t j ) , S ̂ ( t j ) , and ( s ̂ i , σ ̂ i , ρ ̂ ( s , σ ) ) . Since the model contains many parameters and the data is noisy and “sparse,” the model will likely overfit. Therefore, we carry out the parameter estimation by hand in stages, imposing limit on parameter values that are physiologically feasible.

First, we compare C _h(t) (Eq. 6) with the richness C ̂ s ( t j ) shown in Figure 1A to provide a constraint among μ _h, K, C _h(0), and r _h(0). We assume that the C _h(0) associated with granulocytes is slightly greater than the total richness across all samples after the first two, C h ( 0 ) ≳ C ̂ s > 2 . This is equivalent to assuming that granulocyte richness after about 2 months arises solely from barcoded HSCs. The cumulative post-two-month richnesses C ̂ s > 2 for animals ZH33, ZG66, ZH19, and ZJ31 are 2,335, 2007, 4,007, and 30732, respectively. Although the sample specific C ̂ s ( t j ) quickly decreases for t > t ₁, our model prediction for C _h(t) follows Eq. 6 and decays more slowly. By estimating C _h(0) and using Eq. 6, we generate the predicted C _s(t) and S(t) by simulating our full model and comparing them to C ̂ s ( t j ) and S ̂ ( t j ) . This allows us to further constrain the parameters C _h(0), r _h(0), K, and μ _h. Note that no analytic formula exists for C _s(t).

We first set K ≈ 100C _h(0) as the niche carrying capacity since smaller or larger values of K cannot provide the correct average clone sizes or approximately matching values of C _s. This comparison allowed us to obtain rough constraints and approximations to some parameters, particularly μ _h, r _h(0), and C _h(0). Discrete sets of values that are consistent with C ̂ s ( t j ) were selected and further pruned by using the remaining data.

The sampled richness C ̂ s ( t j ) in ZH33 exhibits a sharp decrease after the first sample time, without a corresponding collapse in the sampled abundances of granulocytes S ̂ ( t j ) . Our model explains this phenomenon by the initial condition; namely, the initially transplanted population of CD34⁺ cells contains some partially differentiated HSPCs. Progenitor cells of barcode i (described in our model by the populations n i ( ℓ ) ) are initially transplanted so that some n i ( ℓ ) ( t = 0 ) > 0 particularly for small ℓ (cells with a low degree of differentiation).

As shown in Figure 4, a fraction of the initial clones are HSPCs. Once these HSPC clones generate a burst of mature cells, they disappear from the animal without being renewed since there are no corresponding HSCs carrying the same barcode. Thus, the HSPC contribution to the overall sampled richness C ̂ s ( t j ) largely disappears after about 2 months. However, the total mature granulocyte population S ̂ ( t j ) does not suffer a decline since HSPCs lost due to terminal differentiation are replaced by HSC differentiation. The subsequently sampled mature cell richness then reflects the richness C _h(0) of the initially transplanted HSCs. We propose this partial HSPC transplantation as a mechanism for the observed rapid decrease in C ̂ s observed in some animals. The shape of ρ ̂ ( s ) (see Figure 7D) can inform our estimate of the initial progenitor population n ^(ℓ)(t = 0). Maxima in ρ ̂ ( s ) can be accounted for by offspring of initially transplanted progenitor cells of different stages ℓ, with n ^(ℓ)(0) generating smaller clones for larger ℓ (fewer remaining generations to expand).

Next, we consider the small clones, predominantly arising at short times, and find their average value at the first sample taken at t ₁. These small clones also yield the highest density values ρ ̂ ( s ) , but mostly disappear at longer times. Therefore, we assume they predominantly arise from initial progenitor cells. We can then generate the prediction m(t = t ₁) from our model assuming an initial condition n ^ℓ (0) (and assuming no HSC contribution by setting α = 0). This approximation provides a constraint on the deterministic progenitor cell parameters r _n, L, ω, η. In these experiments, the typical η ∼ 10^–5, so we find that L = 22, and collect a set of feasible values for r _n, ω, and η that provide a good starting point for estimating the other parameters in the model. Note that r _n, L, ω, η can largely compensate each other at this level of comparison. In other words, sets of different ranges of values of one parameter will fit equally well provided some other parameters are also correspondingly adjusted.

Next, we least-squares minimize S ̂ ( t j ) − S ( t j ) , where the model prediction S(t) is given by Eq. 15. This further helps fix α. Once a set of parameters that allow for a reasonable match of C _s(t _j ) and S(t _j ) have been defined, as shown in Figures 5A, B, we then refine the fitting to the mean−standard deviation ( s ̂ i , σ ̂ i ) scatter plot shown in Figure 1D.

To compare our modeling results to the remaining experimental data ( s ̂ i , σ ̂ i , and ρ ̂ i ) and to fine-tune the estimates of the other parameters, the most natural intuition would be to compute the Euclidean distance (or any other relevant distance metric) between model predictions and experimental datasets and tune parameters so as to minimize this distance. However, since a number of clones go extinct in our model, the data size (for s ̂ i , σ ̂ i , and ρ ̂ i ) varies in time.

Thus, before comparing predicted clone size distributions to measured results, we first cluster the data according to the values of s ̂ i and σ ̂ i . Recall that ρ ̂ ( s , σ ) ≈ ρ ̂ ( s ) (since σ is highly correlated with s) is still determined from KDE using the raw, unclustered data ( s ̂ i , σ ̂ i ) . Clustering is performed using k-means to partition the data into multiple regions (in s ̂ - and σ ̂ -space) such that the Euclidean distance between a point and the center of its cluster is smaller than its distance to all other cluster centers. The goal is not to cluster the ( s ̂ i , σ ̂ i ) points according to any real feature, but to simply reduce the dimensionality of the problem and to control the number of effective data points before applying least squares comparisons. Although there are no obvious features in the ( s ̂ i , σ ̂ i ) data, k-means clustering of the distribution of points does yield an optimal number of clusters k* via the “elbow” method where the curvature of the sum of square errors (distortion score) is maximal (Yuan and Yang, 2019). After implementing k-means clustering using the Python yellowbrick package, we find that the optimal number of clusters is typically k* ≈ 50 ± 3 depending on the initial randomization and partitioning process. Subsequent results, however, are insensitive to the precise numbers of clusters used as long as k* ≈ 50.

Figure 6A compares ( s ̂ i , σ ̂ i ) from experiments and from our hybrid multicompartment model. Figure 6B shows the distortion score (blue) and the convergence time (green) as a function of the number of clusters. The optimal number of clusters k* = 48 arises at the elbow of the distortion score curve. Figure 6C shows the clustered data (for animal ZH33) against the clustered model predictions. The radius of each circle, w _k , k = 1, …, 53 denotes the fraction of data points (fraction of the total number of observed clones) assigned to cluster k and is thus a coarse-grained representation of the local data density.

Thus, we have clustered measured data according to P = { p 1 , … , p k * } ≡ { ( s ̂ 1 , σ ̂ 1 , w ̂ 1 ) , … , ( s ̂ k * , σ ̂ k * , w ̂ k * ) } , where ( s ̂ k , σ ̂ k ) denotes the center values of cluster k and w ̂ k is the area of the k ^th cluster. For a fixed set of all model parameters, we generate predictions Q = { q 1 , … , q g * } ≡ { ( s 1 , σ 1 , w 1 ) , … , ( s ℓ * , σ ℓ * , w ℓ * ) } (in this expression and in the following, ℓ, ℓ* denote matrix indices and not progenitor cell generations). The optimal number of clusters derived from our hybrid stochastic-deterministic model, ℓ* will in general be different from the number of clusters k* derived from data, but is typically also about ℓ* ≈ 50.

In order to compare the clustered data P to the model prediction Q, which are matrices with different numbers of rows, we use the Wasserstein metric, or Earth mover’s distance (EMD) (Villani, 2009; Kolouri et al., 2017) to define a distance between them. Let d _k,ℓ denote the distance between cluster p _k and cluster q _ℓ so that the matrix a d = {d _k,ℓ } catalogues all possible cluster-cluster distances. We aim to compute a flow map f = {f _k,ℓ } that yields the minimal distance between clusters p = {p _k } and clusters q = {q _ℓ } by finding min f ∑ k = 1 k * ∑ ℓ = 1 ℓ * f k , ℓ d k , ℓ (16)

Subject to the constraints f k , ℓ ≥ 0 , ∑ ℓ = 1 k * f k , ℓ ≤ w k , ∑ k = 1 ℓ * f k , ℓ ≤ w ℓ , (17)

For all 1 ≤ k ≤ k* and 1 ≤ ℓ ≤ ℓ* and ∑ k = 1 k * ∑ ℓ = 1 ℓ * f k , ℓ = ∑ k = 1 k * w k = ∑ ℓ = 1 ℓ * w ℓ . (18)

After finding the optimal flow f, we evaluate the EMD as E M D P , Q = ∑ k = 1 k * ∑ ℓ = 1 ℓ f k , ℓ d k , ℓ ∑ k ′ = 1 k * ∑ ℓ ′ = 1 ℓ * f k ′ , ℓ ′ . (19)

Model parameters are varied until our model-derived predictions best match the clustered data by minimizing the EMD. In this paper, we consider only the granulocyte lineage since it is the most abundant and reliably measured with minimal complex dynamics and regulation.

Finally, note the fluctuations in S(t _j ) which are too large to be captured by the intrinsic stochasticity in our model, as indicated in Figure 5B. These “unknown” fluctuations can arise from a number of processes, including variable sampling fractions η(t _j ) at each time point t _j and fluctuating animal state due to infections, stress, inflammation, etc. These may influence total mature cell populations month to month. Some of these effects can be effectively accounted for by adjusting the mean sample fraction η = 10^–5 by an amount Δη(t _j ) at each time point. Using these values of Δη(t _j ) to match the data S(t _j ), and then readjusting the parameters, we find good comparison between the experimental measurements and our model, as shown in Figures 7A–C.

To further show consistency, we then plotted the predicted density and compare it with the data-derived (using KDE) density ρ ̂ in Figure 7D. A few large, highly variable clones remain not well reproduced by our model and are discussed in the next section. The parameter estimation procedure was applied to the three other animals in (Koelle et al., 2017), showing reasonable, consistent matching (see Supplementary Material).