Visualizing hematopoiesis as a stochastic process

Hematopoiesis is an intricate system of checks and balances assuring that hematopoietic stem cells (HSCs) give rise to sufficient red blood cells, white blood cells, and platelets to maintain homeostasis and respond to physiological stress. HSCs are infrequent (estimated at 1 per 10⁴ to 10⁵ marrow cells in mice and 1 per 10⁸ marrow cells in humans),^1-4  reside in specified niches that protect HSCs by assuring their relative quiescence and genomic integrity, and give rise to multipotent progenitors (short-term repopulating cells) whose subsequent divisions and differentiation maintain blood cell production.

Conceptually, HSCs are independent decision-makers; a cell’s decisions to differentiate, replicate (self-renew), mobilize, or die are determined by its unique combination of specific genetic, epigenetic, and environmental inputs. Some of these cues are well defined and well studied, whereas others are unknown. An individual HSC’s behavior is not only unique but also dynamic, as its inputs can change over time. Thus, the fates of individual HSCs cannot be suitably predicted deterministically. Furthermore, direct observation of such fates is not possible, so indirect methods must be relied upon, rendering the study of HSCs difficult.

That HSC decision-making is complex and incompletely prescribed, yet outcomes (ie, differentiation, division, mobilization, or death) are clear and distinct, makes hematopoiesis an ideal system to be studied using stochastic modeling. Stochastic analyses have played a large role in modeling and understanding hematopoiesis, and have previously been applied to experimental observations in mice, cats, nonhuman primates, and humans. For example, we have applied this analytic approach to quantify mean rates of HSC differentiation, replication, and death; estimate the total number of HSCs and compute HSC frequencies among marrow cells^3-7 ; and show that the number of niches regulates the number of HSCs.⁸ 

However, understanding of mathematical subtleties and details of the assumptions and limitations of such models can be opaque, especially for hematologists who lack computational or statistical backgrounds. This motivated us to develop a tool that enables nonstatisticians to explore and understand how HSC decisions critically impact hematopoiesis. In this report, we introduce software that provides user-friendly capabilities for stochastic simulation and visualization of the hematopoietic process, allowing users to easily conduct virtual experiments and obtain estimates and uncertainty quantification from their outcomes.

This report presents the basic stochastic model underlying the hematopoiesis simulator proposed by Guttorp et al^9,10  and successfully applied to studies in mice, cats, baboons, and humans^3-5,7,11-13  and demonstrates its utility in several examples. We then describe the visualization tool, which makes available a flexible range of virtual experiments by allowing users to vary parameters and observe the effects of changes. Next, we analyze serial transplantation experiments in mice and studies where “expanded” HSCs are transplanted along with unmanipulated marrow cells in humans, providing new insights into these observations. The visualization tool, called “HematopoiesisSimulator” and an instruction manual are available for download free of charge for noncommercial use at https://github.com/guttorp/hemaviewer.

Mathematical models describe a natural phenomenon or practical problem quantitatively, and have been applied to hematopoiesis research since the 1960s.¹⁴  Stochastic models allow for random variations around some mean behavior as specified by the model parameters. By integrating randomness, relatively simply stochastic models can give rise to a broad range of observed outcomes. Stochastic modeling is therefore often more appropriate for describing very complex systems than deterministic models whose inputs must fully determine all observed outcomes. Deterministic models are more rigid, and when such a model is simplified, it often fails to adequately describe a complex system realistically. In systems featuring a rich ensemble of observed outcomes such as hematopoiesis, deterministic modeling is intractable due to the large number of inputs required to generate the full spectrum of observations. Additionally, stochastic models are more versatile than deterministic models in that they provide quantification of uncertainty or variability (reviewed in Ogawa,¹⁵  Guttorp,¹⁶  and Kimmel¹⁷ ).

As illustrated in Figure 1, HSC behavior in our model is described by 3 parameters: λ, the mean replication (self-renewal) rate; α, the mean apoptosis rate; and ν, the mean rate of differentiation to a multipotent progenitor (short-term repopulating) cell. An average rate for a specified outcome does not imply or require any particular biological mechanism (eg, symmetric division, asymmetric division, a feedback loop, or an age/replication history-dependent decision). We assume that fate decisions are Markovian, that is, that cells do not remember their past decisions, but rather act in a way dependent only on existing cues.¹⁶  The maximum number of HSCs (compartment 1 cells), that is, the maximum capacity of the HSC compartment, is denoted K.

Once an HSC commits to differentiation (leaves compartment 1 and enters compartment 2), the resulting clone contributes to hematopoiesis for a finite period of time. The mean time interval during which contributing clones support blood cell production is denoted 1/μ, and we assume that all contributing clones on average contribute equally.

Models such as this are termed stochastic hidden 2-compartmental models because we only can observe probabilistic behaviors of the second compartment (contributing clones). For example, the observations may consist of sampling committed progenitor cells in marrow or sampling granulocytes in blood, and not the direct observation of HSCs. In previous work, we demonstrated the capability of the stochastic model to analyze hematopoiesis in cats,^6,11,12  mice,^3,5,18  primates,⁷  and humans.⁴  The estimated parameters are consistent with other in vivo or model-based parameter determinations.^4,5  Our studies suggest that K is very similar across animals of different size and lifespans,⁵  as experimental data from mice³  and cats⁵  are consistent with clinical observations of human marrow transplantation.^4,5  Results also suggest that μ is invariant between species, although its value is difficult to infer with precision.^13,19 

Simulating hematopoiesis requires only the specification of the easy-to-conceptualize parameters K, λ, α, ν, and μ, the number of HSCs of the start of the simulation, and the number of contributing clones at the start of the simulation (Table 1).

The Java-based simulator enables visualization of the stochastic model over time (see Figure 2). The program allows researchers to vary the 4 rate parameters λ, α, ν, and μ independently for 2 types of HSCs (termed a and b). Thus, competitive transplantations where the 2 cell types have different behaviors can be simulated. In addition, researchers designate the initial number of HSCs (R0, ie, the number [N] of reserve [compartment 1] cells at time 0), the initial number of contributing clones (C0, ie, the number of contributing [compartment 2] cells at time 0), and the initial percentage of type a cells in each compartment. Total time period of simulation and sampling schedules can be altered as well. In addition, the simulator offers the option to have the HSC reserve increase without bounds or to set the capacity of the HSC reserve and adjust replication fates so that N remains ≤K, where K can be assigned. To cap the population of the HSC reserve, the user can designate how the model behaves when the total number of HSCs becomes greater than or equal to the capacity K. This cap can be set as a hard limit, neglecting any replication events that are set to occur when K is reached. Alternatively, the event may happen with reduced probability so that K acts as a soft limit to the population. Data from Chen et al⁸  suggest that the number of niches regulates the number of HSCs to maintain homeostasis. Conceptually, when an HSC replicates, 1 daughter cell remains in the niche and the second mobilizes and transits the marrow space or circulation. If it encounters an open niche, it engrafts; otherwise, it either dies or differentiates. Taking λ = 0 when N > K (ie, the HSC replication being neglected) can be understood as the mathematical equivalent to a replication event occurring and 1 daughter cell dying. The soft limit option further allows the user to designate the probability by which the homeless HSC dies vs differentiates.

After parameters are chosen, the simulator produces graphical outputs and summary statistics for both the HSC reserve (compartment 1) and the contributing clones (compartment 2). The time series of the percentage of type a cells in each compartment is also plotted. The cellular composition graph and statistics of the 2 compartments in the process of simulation can be viewed and a certain time revisited by using a slider (labeled “Plotting Delay”). The simulator also offers control over the time interval between each sample. The simulation results, along with the parameter values can be exported into (and later loaded from) a file, or simply saved as a screenshot. Moreover, the simulator provides the functionality to perform many replicate simulations. The multiple runs functionality makes it easy to assess variability of outcomes and the sensitivity of the model to parameters. The virtual serial transplantation experiment in the next section illustrates one such example. Complete details about the software are provided in the supplemental Instruction manual.

The portability of Java code removes the platform restriction on the simulator. The simulator is composed of 2 core parts: the simulation engine and the presentation layer. The simulation engine is responsible for carrying out the stochastic simulation. It is built into its own jar file and can be invoked by other programming languages such as R (www.r-project.org), a language and environment for statistical computing and graphics, via the rJava interface.²⁰  The presentation layer accepts user input, sends a request to the simulation engine, and renders the output graphically. In addition to the ready-to-use HematopoiesisSimulator tool, we release all source code and example R programs for plotting and data postprocessing so that advanced users may freely interface and modify all software functionality to best fit their own purposes. An open-source implementation along with thorough documentation are also available in our GitHub repository (https://github.com/guttorp/hemaviewer).

As an initial application of the visualization tool, we simulated serial transplantation experiments in mice. In typical serial transplantation experiments, investigators transplant murine marrow cells into an irradiated mouse, allow hematopoiesis to reconstitute, then 3 to 4 months later, transplant the same number of marrow cells from the first recipient mouse into a second irradiated recipient mouse; and 3 to 4 months later, into a third irradiated recipient, etc. HSC proliferative capacity or robustness is then defined as the number of transplantations that are feasible before host animals fail to reconstitute their hematopoiesis with the serially transplanted cells. It is assumed that HSCs exhaust themselves (are no longer capable of additional divisions) and that this exhaustion explains why cells do not persist after the third or fourth serial transplantation. This assay is a standard method to quantitate HSC robustness.^21-23 

A virtual experiment using the simulator tool and parameter values optimized for murine hematopoiesis (λ = once per 2.5 weeks, α = once per 20 weeks, ν = once per 3.4 weeks, and μ = 6.9 weeks)³  reveals that this conclusion is premature. Intrinsic in the reasoning in the paragraph above is the requirement that the HSC compartment has fully reconstituted by 3 to 4 months after transplantation. However, stochastic simulation (see the paragraph that follows and Figure 3) shows that this assumption is incorrect and, as a consequence, that the “exhaustion” observed in HSCs during serial transplantation can reflect HSC dilution, as was argued by some investigators in the early 1980s.^24,25 

Figure 3 shows cumulative outcomes of 1000 independent trials of a serial transplant experiment, where 10 × 10⁶ nucleated marrow cells are transplanted into irradiated mice. As a mouse is estimated to have 2.8 × 10⁸ marrow cells,²⁶  10 × 10⁶ cells are ∼4% of the total marrow cell number. After 12 weeks (3 months), 10 × 10⁶ marrow cells from the first recipients are “transplanted” into a second irradiated recipient, and so forth. Thus, hematopoietic reconstitution in each successive recipient begins with ∼4% of the number of HSCs present in the previous donor at the end of the 12-week period. We observe that most virtual mice fail to reconstitute after 3 serial transplantations; and as seen in studies of real mice, none survive after the fourth transplantation. The simulation study thus demonstrates that this outcome may be explained solely by the fact that fewer and fewer HSCs are present among the 10 × 10⁶ marrow mononuclear cells transplanted at later time points; the findings do not require that HSCs age or lose proliferative potential after repeated transplantation. If a longer time interval, that is, 24 weeks (6 months) between transplantations is allowed, then 4 serial transplantations are feasible, as shown in the rightmost plot in the bottom panel of Figure 3. Similar outcomes are seen if 20 × 10⁶ (8%), not 10 × 10⁶ (4%), of marrow cells are transplanted to each recipient. This too is consistent with experimental observations. As also shown in Figure 3, if 20 × 10⁶ marrow cells are transplanted each 24 weeks, then two-thirds of the recipients survive after the fourth transplantation.

There are several implications in these findings. First, they show that serial transplantation studies assay HSC number. Components of HSC robustness such as a high HSC replication rate and the slowness of HSCs to differentiate are also measured, but not the differentiation potential of HSCs (data not shown). Second, after the transplantation of 10 × 10⁶ marrow cells, there is a 27-month wait (ie, longer than the average lifetime of a mouse) before its HSC compartment is fully reconstituted, even though the circulating blood cell count normalizes after 3 to 4 weeks. That a transplanted mouse throughout its subsequent life would be especially sensitive to chemotherapy or radiation that injures HSCs is thus anticipated. The same issue pertains to larger animals and humans as can be demonstrated using the visualization tool. A 70-kg man transplanted with 2 × 10⁸ nucleated marrow cells per kilogram (ie, ∼150 HSCs in a 70-kg man) would not fully reconstitute the stem cell reserve (11 000-22 000 HSCs) for roughly 22 years, and would not have large numbers of HSCs (eg, over 1000 HSCs) for ∼12 years (simulations not shown).

A second example where simulating human hematopoiesis can help interpret clinical observations concerns hematopoietic reconstitution after cord blood transplantation when 2 units are transplanted. In adults, 1 strategy is to “expand” 1 cord blood unit with cytokines, such as notch ligand, to assure a more rapid initial reconstitution of hematopoiesis. A second unmanipulated cord blood unit is simultaneously transplanted to assure that blood cell production is maintained. The clinical observation is that after 1 year, some cells derived from the first unit can persist.²⁷  To determine whether this implies that HSCs were present in the expanded cell population that was infused, that is, that the culture conditions not only expanded progenitors but also maintained some HSCs, we used the visualization tool to simulate from the stochastic model beginning with no HSCs, but only with cells in compartment 2. As shown in Figure 4, persistence of short-term repopulating cells is indeed a feasible explanation for the clinical observations. This observation reveals that finding residual contribution from the expanded cord blood unit at 1 year does not imply that the unit contained HSCs. By effectively expanding progenitors, the concomitant graft of HSC (the second unit) could be allowed sufficient time to replicate, reconstitute the stem cell reserve, and support hematopoiesis in the long run. This too can be shown via simulation by adding HSCs from the second unmanipulated cord blood unit into the simulations as population b. Multipotent progenitors from the expanded unit assure hematopoietic reconstitution during the first months after transplantation, whereas HSCs from the unexpanded unit support long-term hematopoiesis in these simulations, as in practice.

Visualizing the outcomes of stochastic simulations illustrates what can happen under a broad range of hypothetical experiments. It allows evaluation of the plausibility of possible explanations intuitively through data visualization and rigorously through quantification of their likelihood. Visualizing hematopoiesis as a stochastic process thus provides a powerful, informative method for analyzing and interpreting observational data.

Because experiments can be simulated under different experimental conditions, the visualization tool (Figure 2) can also help researchers design experiments to assure an interpretable outcome before initiating the study, saving both cost and time. For example, various methods for competitive transplantation studies in mice might be tested by examining the simulated outcomes when 2 types of HSCs with different behaviors, type a and type b, are transplanted. Results of the simulation could then guide selection of an optimal experimental design. By simulating serial competitive transplantation using the approach in Figure 3, what conditions might distinguish a defect in stem cell number from a defect in stem cell quality could be determined. Studies might also be simulated in cats or nonhuman primates, the levels of variability between results compared, and sample sizes estimated, that is, determining that many studies would be needed for an answer up to a desired precision.

As another example, human hematopoiesis can also be simulated (ie, Figure 4) and a variety of hypotheses about clonal progression in the myeloproliferative disorders or the emergence of a paroxysmal nocturnal hemoglobinuria clone after aplastic anemia can be tested. What if a single HSC were able to replicate on average 2 times faster than other HSCs: would its progeny dominate blood cell production? And if so, how fast? What if the replication rate of aberrant HSCs was on average only 1.2 times faster than other HSCs? If the multipotent progenitors derived from this HSC survived 2 times longer than the multipotent progenitors derived from normal HSCs, would the clone then dominate? What if HSC differentiation were inhibited (eg, the chance that an HSC differentiated was halved) but a single HSC that was able to differentiate normally remained: would its progeny dominate hematopoiesis and, if so, with what time course and with what likelihood? What if both HSC numbers were decreased and differentiation impaired but a single HSC that was able to differentiate normally remained? Such questions can be tested as virtual patient outcomes could be generated; then, whether they resemble observed patient data could be assessed. Figure 5 displays how the simulator tool can be used to easily gain intuition about the first 2 questions, and in effect, serve to conceptualize possible initiations and evolutions of clonal hematopoiesis of indeterminate potential.^28,29 

The visualization tool makes simulation accessible and lets the investigator or clinician envision hematopoiesis as a dynamic and competitive process. Conceptualizing and visualizing hematopoiesis using this versatile tool allows for testing of hypotheses, assessment of their plausibility, and quantification of outcomes.

The authors are grateful to Kenny Harrington, Angie Zhu, Youyi Fang, Andrew Hau, Arman Bilge, and Zijin Wang for programming help.

This work was in part supported by National Institutes of Health, National Heart, Lung, and Blood Institute grant R01 HL46598.

Contribution: J.X. designed and performed the simulation studies, analyzed data, and wrote the manuscript; Y.W. helped with programming and the user interface; and P.G. and J.L.A. conceived and designed the visualization tool, planned experiments, reviewed and analyzed data, and wrote the manuscript.

Conflict-of-interest disclosure: The authors declare no competing financial interests.

Correspondence: Janis L. Abkowitz, University of Washington, Box 357710, Seattle, WA 98195; e-mail: janabk@u.washington.edu.