## Abstract

Communication networks between cells and tissues are necessary for homeostasis in multicellular organisms. Intercellular (between cell) communication networks are particularly relevant in stem cell biology, as stem cell fate decisions (self‐renewal, proliferation, lineage specification) are tightly regulated based on physiological demand. We have developed a novel mathematical model of blood stem cell development incorporating cell‐level kinetic parameters as functions of secreted molecule‐mediated intercellular networks. By relation to quantitative cellular assays, our model is capable of predictively simulating many disparate features of both normal and malignant hematopoiesis, relating internal parameters and microenvironmental variables to measurable cell fate outcomes. Through integrated in silico and experimental analyses, we show that blood stem and progenitor cell fate is regulated by cell–cell feedback, and can be controlled non‐cell autonomously by dynamically perturbing intercellular signalling. We extend this concept by demonstrating that variability in the secretion rates of the intercellular regulators is sufficient to explain heterogeneity in culture outputs, and that loss of responsiveness to cell–cell feedback signalling is both necessary and sufficient to induce leukemic transformation in silico.

## Visual Overview

### Synopsis

Communication networks among cells, tissues, and organ systems are necessary for homeostasis in multicellular organisms. Intercellular communication networks are particularly relevant in stem cell biology, as stem cell fate decisions (self‐renewal, proliferation, lineage specification) are tightly regulated based on physiological demand and responsive to external perturbations. Hematopoiesis, the process by which blood cells develop, serves as a prototype for other stem cell systems. Motivated by the fact that hematopoietic stem cell transplantation is a curative therapy for a number of hematopoietic and immunological diseases, herein we explore the behaviour of intercellular regulatory networks as tools to regulate cell fate during *in vitro* human blood stem cell propagation.

We have integrated a diverse range of experimental and theoretical literature to develop for predictive purposes a novel mathematical model of hematopoiesis. Our resulting model links functional cellular assays to specific model outputs, defines cell‐level kinetic parameters such as cell cycle rates and self‐renewal probabilities as functions of culture variables, and simulates feedback regulation using cell–cell interaction networks (Schematically depicted in Figure 1A). Computational analyses of the system dynamics indicate that non‐cell autonomous parameters (cell–cell interactions) are dominant factors controlling stem cell growth. As a result of negative feedback interactions, our analysis suggest an antagonistic relationship between mature and primitive cell compartments, a finding empirically supported by numerous *in vitro* and *in vivo* observations.

To investigate the practical utility of the model, we cultured human umbilical cord blood stem cells under a range of conditions and dynamic perturbations (media exchange and mature cell depletion), and compared experimentally observed cell population outputs to model simulations. Our model simulations quantitatively recapitulate experimental results; culturing enriched progenitor populations at low initial cell densities, with frequent media exchange, and with progenitor re‐selection enhances stem and progenitor expansion as a consequence of reduced inhibitory feedback signalling. Using a protein array system, we identify a limited number of secreted molecules (TGF‐β1, CCL4, CCL5, and CXCL8) displaying dynamic profiles consistent with predicted intercellular regulators. Additionally, we show the variability in secretion rates observed for these putative regulators is sufficient to explain the experimentally observed distribution of cell population outputs.

Computational analyses of the model revealed that system dynamics are relatively robust to changes in kinetic parameters, but highly sensitive to topological alterations of the proposed cell–cell regulatory network. We extend this concept by demonstrating that deleting a single regulatory connection (negative feedback control of self‐renewal) is necessary and sufficient to reproduce the characteristic features of *in vitro* leukemic transformation (Figure 8).

In summary, our intercellular feedback model of blood stem cell growth predictively simulates the dynamic characteristics of both normal and malignant hematopoiesis *in vitro*. This model may therefore serve as a platform for further experimental interrogation of the regulatory mechanisms controlling stem cell fate *in vitro* and *in vivo*, and as a tool for the rational design of stem cell‐therapy bioprocesses.

We present a novel mathematical model of blood stem cell development incorporating cell level kinetic parameters as functions of inter‐cellular communication networks.

Negative feedback signalling between differentiated cells and stem and progenitor populations is a dominant factor regulating cell population outputs during in vitro culture of human blood stem cells, producing an antagonistic relationship between mature and primitive cell compartments.

Stem cell fate can be controlled non‐autonomously by the dynamic perturbation of cell‐cell signalling networks.

Dysregulated responsiveness to inter‐cellular signals is both necessary and sufficient to induce leukemic transformation in silico.

## Introduction

Systems‐biology research to date has primarily focused on elucidating the topological features and dynamics of intracellular networks, with the implicit assumption of cell populations as homogenous, autonomous units. Intercellular communication networks, represented with cell types as vertices, and functional (rather than molecular) interactions as edges, have been largely unexplored.

Communication networks among cells, tissues, and organ systems are necessary for homeostasis in multicellular organisms. For example, soluble factor‐mediated cell–cell networks have a dominant function in orchestrating immune reactions through ‘cytokine cascades’ in response to infection. Intercellular communication networks are particularly relevant in stem cell biology, as stem cell fate decisions (self‐renewal, proliferation, lineage specification) are tightly regulated based on physiological demand and responsive to external perturbations. Evidence suggests that stem cell dysregulation is fundamental to the progression of multiple cancers, degenerative diseases, and general aging phenomenon (Rossi *et al*, 2008). Stem cell‐based therapies are hence emerging as a foundational tool in regenerative medicine. One of the key challenges lies in controlling the emergent cellular and microenvironmental complexity that arises as stem cell populations develop *in vitro* and *in vivo* (Kirouac and Zandstra, 2006).

Hematopoiesis, the process by which blood cells develop, serves as a prototype for other stem cell systems. Hematopoietic stem cells (HSCs), at the apex of a developmental hierarchy, give rise to a series of increasingly differentiated and developmentally restricted progenitor cells, eventually producing all of the mature blood cell populations. *In vivo*, HSC fate decisions are regulated by cross‐talk with neighbouring cell populations either directly or through secreted factors (Wilson and Trumpp, 2006). Much experimental and theoretical work has been conducted to understand the structure and dynamics of the homeostatic control mechanisms *in vivo* linking blood cells in the circulation to stem cells in the bone marrow (Lajtha *et al*, 1962). Evidence suggests that mature blood cells suppress proliferation and differentiation of progenitors, and progenitors correspondingly suppress the expansion of stem cells through coupled negative feedback loops (Wichmann and Loeffler, 1985).

Hematopoietic cells are known to secrete and respond to a large number of regulatory proteins in lineage‐ and differentiation stage‐specific patterns (Billia *et al*, 2001; Majka *et al*, 2001). This results in complex and dynamic intercellular signalling networks, providing a mechanism by which cells interrogate and interpret their local microenvironment (their niche), propagate this information through signal transduction and gene regulatory networks, and respond by modulating cell fate decisions. A number of studies have attempted to reconstruct the intracellular molecular networks regulating stem cell fate (Muller *et al*, 2008), including mathematical modelling of genetic regulatory networks and intracellular feedback mechanisms (Glauche *et al*, 2007). However, intercellular regulatory mechanisms remain largely undefined.

Motivated by the fact that HSC transplantation is a curative therapy for a number of hematopoietic and immunological diseases, herein we explore the behaviour of intercellular regulatory networks as tools to regulate cell fate during *in vitro* human blood stem cell propagation. We have developed and used for predictive purposes a novel mathematical model of *in vitro* hematopoiesis by linking functional cellular assays to specific model outputs, by defining cell‐level kinetic parameters such as cell cycle rates and self‐renewal probabilities as functions of culture variables, and by simulating feedback regulation using cell–cell interaction networks. Our resultant model captures many facets of hematopoiesis, connecting internal model parameters and microenvironmental variables to measurable cell fate changes. We show that negative feedback signalling between differentiated cells and stem and progenitor populations is a dominant factor regulating culture output, and that stem cell fate can be controlled non‐autonomously by the dynamic perturbation of cell–cell signalling networks. We extend this concept by demonstrating that variability in the secretion rates of intercellular regulators are sufficient to explain variability in culture output, and that loss of responsiveness to cell–cell feedback signalling is both necessary and sufficient to induce leukemic transformation in silico. The development of quantitative models incorporating cell–cell regulatory networks should serve as an important tool to understand and control emergent cellular complexity *in vitro* and *in vivo* (Kirouac and Zandstra, 2008).

## Results

### A feedback‐based cell–cell interaction network model of hematopoiesis

The hematopoietic hierarchy can be divided into discrete cellular compartments, wherein compartment transitions are typically coincident with compartment size amplifying cell divisions. Taking advantage of differentiation‐state‐associated *in vitro* and *in vivo* assays, we have defined functional readouts as overlapping series of consecutive compartments. The functional readouts we consider are the immunodeficient (non‐obese diabetic (NOD)/Scid) mouse repopulating cell (SRC) assay for quantifying stem cells, the long‐term culture‐initiating cell (LTC‐IC) assay for quantifying primitive progenitors, and the colony forming cell (CFC) assay for quantifying committed progenitors (Coulombel, 2004). Hematopoietic cell populations are also broadly classified phenotypically based on their expression (Lin^{+}), or lack of expression (Lin^{−}), of blood lineage‐associated cell surface antigens. Frequency estimates in umbilical cord blood total nucleated cells (TNC) for cell compartment—assay relationships, as described in Materials and methods, are shown in Table I, and a schematic diagram of the model is depicted in Figure 1A. For clarity, all model parameters, variables, and simulated cell population outputs will henceforth be indicated with italics.

It is well established experimentally that self‐renewal probabilities decrease and cell cycling rates increase with differentiation (Young *et al*, 1996). We implemented Gaussian‐type functions to quantify proliferation rates (*u*_{i}) and the self‐renewal probabilities (*f*_{i}) as functions of compartment number (*i*) (see Materials and methods for details and Figure 1B and C for phase portraits).

We simulate a branching model of hematopoiesis by lumping differentiated (Lin^{+}) cells into three functional classes based on their functional feedback interactions with stem and progenitor cells; populations that secrete inhibitors, populations that secrete stimulators, and populations that secrete molecules with no net effect. The framework for the intercellular signalling network topology is based largely around *in vivo* experimental data on murine hematopoietic homeostasis (Wichmann and Loeffler, 1985) with the implicit assumption that the network will be at least partially reconstituted *in vitro*. We have designated compartment‐specific self‐renewal and proliferation rates as regulated by a balance between endogenously secreted inhibitors (negative feedback) and stimulators (positive feedback) (Figure 1D). We consider only soluble factor‐mediated signalling between blood cell populations, and thus limit our analysis to liquid suspension cultures; the additional complexities associated with stoma‐supported co‐cultures or *in vivo* hematopoiesis are not implicitly considered. The resultant model consists of 24 state variables [20 cell compartments (*X*_{i}) and 4 secreted regulatory molecules (SF1–4), defined in Figure 1A] and 16 internal parameters, their definitions and theoretically constrained ranges are given in Table II.

### Global parameter space analysis reveals a critical function for non‐cell autonomous parameters

To systematically explore the parameter space in an unbiased manner, functional outputs from theoretical 8‐day cultures (Madlambayan *et al*, 2005) were simulated for 1000 random parameter values log‐uniformly distributed within the constraints (‘LO’ and ‘HI’ values) defined in Table II. Figure 2A shows the statistical relationship among the *in vivo*, *in vitro*, and phenotypic assays, and Figure 2B shows the statistical relationship between individual parameter values and culture outputs. Pearson correlation coefficients (*r*) were calculated for each relationship and are shown on each plot.

Notably, there is a poor correlation between total nucleated cell (*TNC*) and stem cell (*SRC*) expansion (*r*=−0.029) and between progenitor (*CFC*) and stem cell (*SRC*) expansion (*r*=0.433). The model, thus, captures the recognition that these mature cell assays are poor surrogates for stem cell output (Zandstra *et al*, 1997). Our analysis suggests that under many conditions, the LTC‐IC assay is a better surrogate for stem cell content (*r*=0.703). Production of differentiated Lin^{+} cells correlates negatively with stem and progenitor cell output (*r*=−0.322 for *CFC* and *LTC‐IC*, *r*=−0.187 for *SRC*) likely due to the negative feedback‐biased architecture of the system.

As one might predict, the progenitor proliferation rate (*u*_{MAX}) is the best positive predictor of total cell (*TNC*), progenitor (*CFC*), and primitive progenitor (*LTC‐IC*) expansion (*r*=0.695, 0.452, and 0.274, respectively). However, surprisingly, this parameter (*u*_{MAX}) is inversely correlated with stem cell (*SRC*) expansion (*r*=−0.161). Non‐intuitively, the strongest positive correlate to stem cell (*SRC*) expansion is secretion rate of the proliferation inhibitor *SF1* (*sr1*) (*r*=0.227), a finding that suggests a negative relationship between total cell and stem cell growth. The secretion rate of the self‐renewal inhibitor *SF2* (*sr2*) is the best negative correlate with *CFC*, *LTC‐IC*, and *SRC* expansion (*r*=−0.14, −0.372, and −0.479, respectively). This analysis suggests that non‐cell autonomous parameters (soluble factor‐mediated cell–cell interactions) may be dominant factors controlling stem cell growth. It is notable that while the parameter correlations are low, particularly with regard to *SRC* output (the highest *r*^{2}=0.0515), most are statistically significant (∣*r*∣>0.062 corresponds to *P*<0.05 for *n*=1000 simulations). The wide distribution of outputs observed is indicative of highly non‐linear parameter interactions, necessitating the use of stochastic global optimisation methods to estimate parameter values. Thus, we next devised an appropriate reverse engineering strategy using earlier published data.

### Parameter training identifies self‐renewal inhibitor (SF2) exposure time as a key regulator of stem and primitive progenitor cell output

Our laboratory has previously developed a bioprocess for the clinical‐grade expansion of umbilical cord blood (UCB)‐derived HSCs (Madlambayan *et al*, 2006). This system allows us to dynamically perturb cell–cell signalling networks, and thereby interrogate our theoretical model. Specifically, Lin^{−} cell negative selection (S) reduces the secretion rates sr1 to sr4, and media exchange (E) reduces the concentration of factors SF1 to SF4 (Figure 3A). Eight‐day culture outputs for combinations of Lin^{−} cell selection (S) versus no selection (NS), and media exchange (E) versus no exchange (NE) at culture day‐4, as reported earlier (Madlambayan *et al*, 2005), are shown in Table III. Performing Lin^{−} cell selection and media exchange in combination (S/E) synergistically enhanced primitive progenitor (LTC‐IC) and stem cell (SRC) output, while the effect of the described culture manipulations on other system outputs (%Lin^{+}, TNC, or CFC expansion) was not statistically significant (*P*>0.1).

To estimate parameter values, we used the data displayed in Table III to define an objective function based on the sum of squared residuals between experimental and simulated culture outputs (see equation (25)). As the objective landscape was expected to be multi‐modal, a hybrid stochastic method was developed to estimate the unknown parameters by minimising the objective function. A genetic algorithm was first used to globally search the parameter space, followed by the Nelder–Mead Simplex algorithm for local optimisation. The estimated parameter vector is displayed in Table II (‘Est’ column), and the results are summarised in Figure 3B. Coefficients of variation (CVs) for the parameter estimates (Table II ‘CV’ column) were calculated using a computational bootstrapping approach, and indicate that a number of parameters are non‐identifiable (i.e. system dynamics can be simulated using multiple parameter combinations, indicated in Table II ‘ID’ column; ‘+’=identifiable, ‘−’=non‐identifiable). Such ‘sloppy’ parameter estimates are in fact a consistent feature of systems biology models, and may reflect underlying robustness of biological networks (Gutenkunst *et al*, 2007).

Importantly, the model captures the experimental results, and explains the synergistic and targeted effect of the culture manipulations (S/E) on stem and primitive progenitor populations through reduced exposure (concentration × time) to the proliferation and self‐renewal inhibitors *SF1* and *SF2* (Figure 3C). To further investigate the relationship between individual parameters and system dynamics, parameter sensitivity analysis was performed.

### Phase portraits and parameter sensitivity analysis demonstrate antagonistic relationship between mature and primitive cell compartments

Parameter sensitivities represent emergent systems‐level properties, which cannot be deduced by analysing individual components or reactions (Savageau, 1971). Parameter sensitivity coefficients (*S*) were calculated for each system output (*O*) ∼ parameter (*P*) combination (equations (26) and (27)), producing a 5 × 16 sensitivity coefficient matrix, depicted as a heat map in Figure 4. Parameters and system outputs were organised through unsupervised hierarchical clustering, and distances (1−*r*) represented by associated dendograms.

The system outputs organise into two primary clusters, one containing differentiated (*Lin*^{+}) cells, total cells (*TNC*), and mature progenitors (*CFC*), and the other containing primitive progenitors (*LTC‐IC*) and stem cells (*SRC*). This clustering pattern emphasises the modular nature of cell compartments, and the relative independence of primitive and mature cell outputs. The parameters correspondingly fall into two general categories; those that positively influence stem cell (*SRC*) expansion while repressing total cell (*TNC*) growth, and vice versa. This demonstrates an antagonistic relationship between mature and primitive cell compartments, a result of the (negative) feedback interactions.

To highlight the effects of individual parameters on system dynamics, select phase portraits are shown in Figure 4B, wherein system outputs [*TNC* and *SRC*] are plotted against time, for individual parameter values [progenitor proliferation rate (*u*_{MAX}), self‐renewal inhibitor and stimulator secretion rates (*sr2* and *sr4*, respectively)] varied over the ranges given in Table II. In all cases, there is a clear divergence in the population dynamics for total cells (*TNC*) versus stem cells (*SRC*); while the total cell (*TNC*) numbers display near‐exponential growth throughout the 8‐day cultures, primitive progenitors (*LTC‐IC*) and stem cell (*SRC*) numbers plateau mid‐culture and then begin to decline. This bi‐phasic growth curve has been observed in many blood progenitor cultures, regardless of the choice of cytokine cocktail (Zandstra *et al*, 1997; Mobest *et al*, 1999).

As expected, the most sensitive parameter controlling stem cell growth is the probability of self‐renewal (*f*_{MAX}). This is consistent with most current strategies for inducing *in vitro* stem cell expansion through the identification and use of specific modulators of self‐renewal (Sauvageau *et al*, 2004). However, our analysis suggests for the first time that the system dynamics are predominantly regulated by the non‐stem cell autonomous negative feedback signals (*SF1* and *SF2*) rather than stem cell‐autonomous factors. Following this, sustained stem cell amplification should depend on removal or blocking of these signals. In fact, factors inducing total cell expansion [i.e. increasing progenitor proliferation rate (*u*_{MAX}) through stimulatory growth factors such as interleukin (IL)‐3] may actually repress stem cell growth through enhanced accumulation of cells secreting inhibitory factors. We next sought to use this insight in a new series of experiments to predictively manipulate culture variables to control stem and progenitor cell output.

### Dynamic culture perturbations reveal intercellular feedback control of stem and progenitor cell output

Our model simulations predict that cell–cell feedback would produce non‐linear cell population dynamics in response to culture perturbations. In a first series of experiments, we measured the *in vitro* production of undifferentiated (Lin^{−}) cells, differentiated (Lin^{+}) cells, and progenitors (CFC) in response to consecutive culture perturbations. UCB‐derived Lin^{−} cells were cultured for 8 days unmanipulated (NS/NE), or subject to the Lin^{+} cell selection and media exchange procedure at culture day‐4 (S/E‐d4), or at day‐4 and day‐6 consecutively (S/E‐d4,d6) (Figure 5A and B). The relative frequency of differentiated (%Lin^{+}) cells at culture day‐8 was unaffected by the number of times the Lin^{+} cells were depleted. Although the production of undifferentiated (Lin^{−}) cells across conditions was unaffected, total output of differentiated (Lin^{+}) cells was increased in response to the S/E procedures (*P*=0.079, *t*‐test). Additionally, progenitor (CFC) output was enhanced in response to consecutive perturbations (*P*=0.075, *t*‐test). These non‐intuitive results can be attributed to the dynamic regulation of progenitor cell proliferation and differentiation by the intercellular feedback control loops.

The wide sample‐to‐sample variability in the expansion potential of blood stem cells necessitates a large number of replicate experiments to achieve statistical significance (Koller *et al*, 1996). We examined what sources of exogenous (culture variables) or endogenous (biological) variability could account for the variability typically observed in culture. To test this, 8‐day unmanipulated (NS/NE) culture simulations were run with individual parameter values and select combinations normally distributed (CV=1) around the values given in Table III. Experimental (*n*=9) and simulated (*n*=100) output distributions were compared using a ranking metric based on the sum of squared residuals between total cell (TNC versus *TNC*), progenitor (CFC versus *CFC*), and primitive progenitor (LTC‐IC versus *LTC‐IC*) population expansion means and CVs (Supplementary Figure S1). Variabilities in the secretion rates *sr1*–*sr4* most closely simulated the distribution of experimental cell population outputs. As shown in Figure 5C, the distribution of culture outputs [CV for TNC=0.64, CFC=0.67, and SRC=0.73] is closely recapitulated in silico [CV for *TNC*=0.60, *CFC*=0.83, and *SRC*=0.76], providing evidence that variability in culture outputs may be attributable to variability in intercellular signalling dynamics. That is, microenvironmental noise may have a function in guiding cell fate.

### Intercellular feedback signalling regulates LTC dynamics

We next asked whether our model simulations are capable of accurately predicting cell population dynamics for extended (16‐day) liquid cultures. Although in unmanipulated cultures, stem and progenitor cell numbers rapidly plateau and begin to decline within a week, our model simulations predicted that by performing consecutive culture manipulations, stem and progenitor cell growth would be enhanced. Liquid cultures were maintained as described with Lin^{−} cell election and media exchange (S/E) performed every 4 days. As shown in Figure 6A, total cells (TNC) continue to proliferate throughout the 16‐day culture, whereas progenitor (CFC) expansion plateaus after day‐12, and primitive progenitor (LTC‐IC) numbers begin to decline after a maximal expansion of approximately 12‐fold at day‐8. Cultures subjected to media exchange alone every 4‐days (NS/E) were predicted to decline faster (Figure 6A, dashed lines). The model simulations fit the experimental data for media exchange efficiencies (defined as the percentage of media replaced) of approximately 90% or greater (Supplementary Figure S2), and suggest that the culture decline was due to increased secreted factor (*SF1* and *SF2*) feedback from the accumulating differentiating cells (Figure 6B). To interrogate the predicted dynamic correlation between inhibitory protein secretion (*SF1* and *SF2*) and culture decline, we profiled conditioned media obtained at culture days 4, 8, 12, and 16 for a set of 31 ligands using a liquid chip cytokine array. Although the sample‐to‐sample variability in protein concentrations detected was large, the chemokines MIP‐1β (CCL4), IL‐8 (CXCL8), and RANTES (CCL5), as well as TGF‐β1 display dynamics qualitatively similar to that of the theoretical proliferation and self‐renewal inhibitors SF1 and SF2 (Figure 6C). The remaining 27 ligands were not consistently detected above background (Supplementary Table SI). Many additional factors are likely to be involved; however, technical limitations of proteomic technology currently preclude systematic measurement for the small sample volumes obtained from these primary human cell cultures. Although the specific functional effects of combinations of these factors on stem cell fate have yet to be elucidated, these data show that a number of ligands are endogenously secreted with dynamics qualitatively similar to our model predictions.

Model simulations predicted that increased frequency of culture manipulations at later time points (day‐8+) would improve primitive progenitor (LTC‐IC) expansion, while having no effect on the output of total cells (TNC) or committed progenitors (CFC). To test this prediction, we conducted a series of 12‐day liquid cultures as described (S/E‐d4,8), or with an additional S/E procedure at day 10 (S/E‐d4,8,10). As predicted, the additional culture manipulation had no effect on total (TNC) or progenitor cell (CFC) output (not shown); however, primitive progenitor (LTC‐IC) expansion was enhanced from 11±4 to 19±3‐fold (*P*=0.10, *t*‐test) (Figure 6D). Primitive progenitor cell output can, thus, be specifically manipulated by dynamically perturbing intercellular feedback. We next asked whether our intercellular feedback model could provide mechanistic insight into the effects of culture variables on cell population outputs.

### Culture strategies to manipulate intercellular feedback signalling: plating density and input cell population enrichment

Thus far, all cell cultures described have been initiated at Lin^{−} cell densities of 10^{5} cell/ml. However, it has been reported in a number of studies that low cell densities and progenitor enrichment are associated with greater stem cell expansions (Sandstrom *et al*, 1995; Kohler *et al*, 1999; Xu *et al*, 2000), although the underlying mechanisms are largely undefined. To study the effect of plating density on culture outputs, both *in vitro* and in silico cultures were initiated with UCB Lin^{−} cells at densities of 10^{4}, 5 × 10^{4}, 10^{5}, and 5 × 10^{5}/ml. As shown in Figure 7A (upper panel), lower cell seeding densities induce greater expansions of total cells (TNC), progenitors (CFC), and primitive progenitors (LTC‐IC).

In a second series of simulations, cultures were initiated in silico using cell populations with different amounts of stem cell enrichment [from *i*=20 (∼TNC, Lin^{−} plus Lin^{+}) to *i*=8 (∼LTC‐IC)], and 8‐day culture outputs were simulated (Figure 7B, upper panel). Although the corresponding experiment was not performed, experimental data for similar studies initiated with mononuclear cells (TNC), Lin^{−}, Lin^{−}CD34^{+}, and Lin^{−}CD34^{+}CD38^{−} enriched populations have been published (Conneally *et al*, 1997; Kohler *et al*, 1999). Consistent with earlier reports, expansion of stem cells (*SRC*), progenitors (*LTC‐IC*, *CFC*), and total cells (*TNC*) directly correlate with the degree of progenitor cell enrichment.

Examination of the simulated secreted factor concentrations (Figure 7A and B, lower panels) under these conditions yields mechanistic insight. Cultures seeded at lower cell densities and/or with more enriched populations are cumulatively (concentration × time) exposed to lower levels of endogenously produced proliferation and self‐renewal inhibitors (*SF1* and *SF2*), resulting in greater stem and progenitor cell expansions.

Our intercellular feedback model is thus capable of simulating many features of *in vitro* hematopoiesis. As the model contains a large number of (non‐identifiable) parameters, we performed analysis to determine whether our results could be attributable to parameter ‘over‐fitting’. Intracellular network connections were systematically re‐wired producing a set of nine structurally altered models (S1–9) with conserved or decreased number of free parameters. We then calculated the Akaike information criterion (AIC) for each, a metric used to rank alternative models based on their ability to explain data with a minimum number of free parameters (Landaw and DiStefano, 1984). Our original model performs significantly better than 8/9 altered models; however, deletion of the positive proliferation feedback connection (SF3) produces an incrementally higher ranking (model ‘S9’), as the data are equivalently described, but with two fewer parameters (Supplementary Figure S3). As our culture media is supplemented with high concentrations of proliferation‐inducing cytokines, endogenous proliferative signals may be obscured. However, in other *in vitro* or *in vivo* environments, positive feedback loops may have a significant function (Eaves and Eaves, 1988). This connection was hence maintained to enable broader model applicability. We next queried whether the model could be extended to hematopoietic pathologies.

### Loss of responsiveness to endogenous self‐renewal inhibitors simulates stem cell leukemic transformation

Progenitor enriched (Lin^{−} or CD34^{+}) UCB cells have been used to generate *in vitro* models of leukemic stem cells (LSCs) by transformation with leukaemia‐associated gene fusions such as *MLL‐ENL* (Barabe *et al*, 2007), *MLL‐AF9* (Wei *et al*, 2008), and *TLS‐ERG* (Warner *et al*, 2005). The resulting cell lines are immortalised, yet growth factor‐dependent and growth factor‐responsive, display partially blocked myeloid differentiation patterns, and induce acute myeloid leukaemias *in vivo* on transplantation into immunodeficient mice.

We used data from Warner *et al* (2005) as a training set, and the reverse engineering strategy described earlier, to determine model parameter sets that most closely simulate long‐term LSC cultures. The simplest notion of tumourgenesis, as applied in chemo‐ and radiation therapies, is that cancerous cells have a higher intrinsic proliferation rate. To simulate this, the progenitor proliferation rate (*u*_{MAX}) was increased by 50%. We could account for a partial differentiation block through reducing the Lin^{+} maturation rate (*u*_{+}) by 50%. Enhanced self‐renewal associated with LSCs was simulated by increasing the probability of self‐renewal (*f*_{MAX}) by 15%. As malignant transformation involves multiple mutations, we also considered the effect of modulating all three parameters (*u*_{MAX}, *u*_{+}, and *f*_{MAX}) simultaneously. Finally, we considered a situation wherein stem cell responsiveness to endogenous self‐renewal inhibitors (*SF2*) are diminished. Simulated culture outputs are shown in Figure 8. Simulations were run using a range of parameter values; however, the results were qualitatively unaffected.

Remarkably, modulating any of the three parameters (*u*_{MAX}, *u*_{+}, or *f*_{MAX}) independently had little effect on the long‐term *in vitro* growth kinetics. Modulating the three parameters simultaneously had additive effects; however, the cell line was not immortalised as the stem cells (*SRC*) were gradually depleted as a result of accumulating differentiated cells and associated feedback inhibition. When the stem cells were simulated to be unresponsive to the endogenous self‐renewal inhibitor (*SF2*), they display all of the characteristic features of *in vitro*‐derived LSCs; sustained proliferation (immortalisation), reduced differentiation, and constitutive self‐renewal. The dynamic oscillations result from the weekly re‐plating, as the transformed cells are still responsive to the endogenous proliferation inhibitor (*SF1*). Modulating this single topological feature of the model (which may represent a combination of several distinct molecular events) is, therefore, necessary and sufficient to induce leukemic transformation in silico. This is consistent with the parameteric and structural analysis of our model; system dynamics are robust against alterations in kinetic parameters, but sensitive to alterations of the regulatory network structure.

## Discussion

The cellular network model described predictively recapitulates many disparate features of *in vitro* hematopoiesis, explaining the enhanced stem and progenitor cell expansions observed at low cell densities, progenitor enrichment, frequent media exchange, and progenitor re‐selection as a consequence of reduced inhibitory feedback signalling. The model suggests scenarios in which stem cell fate is non‐cell autonomously controlled by intercellular signalling dynamics. This basic concept could be extended to malignant hematopoiesis, as dysregulated responsiveness to intercellular signals was sufficient to reproduce characteristic features of *in vitro* leukemic transformation. Additionally, biological variability in the secretion rates of endogenous regulators was sufficient to explain the large sample‐to‐sample variability observed in culture outputs. Although the model was developed specifically for blood stem cell growth *in vitro*, the underlying principles may be more broadly applicable to other stem cell systems both *in vitro* and *in vivo*.

Although many individual aspects of our model have been considered earlier, to our knowledge this is the first mathematical model of stem cell growth that systematically considers (1) quantitative functional readouts as lumped and overlapping cell compartments, (2) both proliferation rates and self‐renewal probabilities as functions of differentiation state, and (3) proliferation rates and self‐renewal probabilities as functions of endogenously secreted positive and negative regulatory molecules. Although some of the structural assumptions may be disputed (i.e. the hierarchy may exist as a functional continuum rather than discrete stages (Quesenberry, 2006) and cells may alter their functional status (undergo compartment transitions) without proliferating (Kent *et al*, 2008)), these alterations have minimal effects on model outputs. Additionally, the use of differential equations limits the model's applicability to situations in which the stochastic effects of low stem cell numbers are negligible. Modelling the growth of limiting numbers of stem cells requires the use of single cell‐based (stochastic) models (Roeder *et al*, 2008). Further, our analysis has been limited to liquid suspension cultures; extension to stoma‐supported culture or *in vivo* hematopoiesis would require additional complexities such as direct cell–cell contact, blood–stroma interactions, and/or hormonal regulatory mechanisms to be incorporated.

Earlier models of *in vivo* hematopoiesis incorporating cell–cell feedback have set kinetic parameters (or probability functions) as directly responsive to neighbouring cell densities (Wichmann and Loeffler, 1985; Roeder and Loeffler, 2002). The limited number of models of *in vitro* hematopoiesis have either ignored intercellular feedback (Varma *et al*, 1992) or considered it indirectly by setting proliferation kinetics as cell density dependant, assumed to be a result of nutrient depletion as in microbial fermentations (Peng *et al*, 1996). Hence, our model is the first to explicitly consider a mechanism (secreted regulatory molecules) whereby cells are capable of sensing local cell densities and differentiation status, and applying this regulatory feature to *in vitro* culture.

Numerous molecular genetic studies have identified a negative relationship between stem cell proliferation and self‐renewal—genetic mutations that constitutively enhance HSC cycling generally result in a loss of long‐term repopulating cells, and vice versa (Orford and Scadden, 2008). Our findings suggest that this antagonistic relationship may be mediated (at least partially) by non‐cell autonomous effects—rapid cycling results in the accumulation of progenitors that inhibit self‐renewal through feedback signalling. This is supported by recent demonstrations that normally quiescent stem cells are induced to rapidly proliferate in response to progenitor depletion *in vivo* (Wilson *et al*, 2008). Data from various *in vitro* culture systems are also consistent with our model, indicating that soluble endogenous factors limit stem cell amplification, and progenitor re‐enrichment and media exchange/dilution enhance stem and progenitor cell growth (Eaves and Eaves, 1988; Gammaitoni *et al*, 2004; Flores‐Guzman *et al*, 2006). The antagonistic relationship between mature and primitive cell output hence has direct implications for bioprocess design studies.

Our analysis suggests two complementary approaches to control the intercellular signalling networks established in culture. If the molecular regulators are known in sufficient detail, they may be targeted through the use of blocking antibodies or receptor antagonists against inhibitory factors, the exogenous addition of competing stimulatory factors, or by direct modulation of intracellular signalling pathways using small molecules. The specific inhibitory/stimulatory factors may change over time as the culture evolves, likely necessitating dynamic media supplementation. Conversely, a global culture manipulation strategy may be used, wherein media exchange or dilution and cell subpopulation‐specific selection procedures are performed to indirectly control the intercellular signalling network dynamics.

Although there are likely to be many exogenous and endogenous sources of culture variability, we demonstrate that variabilities in the secretion rates of the endogenous regulators are sufficient to reproduce the variability observed in culture outputs. It has been reported that donor‐to‐donor variability in culture output is reduced though the use of stromal feeder cells (Koller *et al*, 1996), and wide distributions in stem and progenitors numbers is not observed between individuals *in vivo* under normal physiological conditions. The lack of additional regulatory controls in liquid cultures may result in the amplification of initial biological differences, which would be compensated for by stromal cells *in vitro* and *in vivo*. Understanding the emergence of such biological noise in culture will be important for the design of robust cell therapy bioprocesses.

A main limitation of the model in its present form is that the secreted regulators *SF1‐4* remain largely theoretical. We have, nonetheless, identified a limited number of candidate regulatory factors, that is, TGF‐β1, MIP‐1β (CCL4), IL‐8 (CXCL8), and RANTES (CCL5). It is notable that TGF‐β1 is a well‐established inhibitor of stem and progenitor expansion *in vitro* and *in vivo* (Yamazaki *et al*, 2009), and IL‐8 suppresses myeloid colony formation *in vitro* (Broxmeyer and Kim, 1999). Many molecules with stimulatory effects on proliferation and self‐renewal have been reported as being expressed by hematopoietic cells. Importantly, it is unlikely that any one or even a few of these molecules could be added or blocked to simulate model predictions. Instead, we propose that it is the net balance of opposing signals, integrated by the cell through intracellular signalling pathways, which regulates cell fate decisions.

Our simulations show that while perturbing kinetic constants has little effect on long‐term system dynamics, modifying the network structure (deleting the self‐renewal negative feedback loop) is both necessary and sufficient to induce leukemic transformation. This is consistent with experimental data, as leukemic cells are known to display abnormal growth‐factor response characteristics, and specifically to be resistant to TGF‐β1 inhibition (Lin *et al*, 2005). A number of other groups have modelled leukaemia using both differential equation‐ and stochastic single cell‐based models. Although all models rely on a basic assumption that LSCs have a competitive advantage against normal HSCs *in vivo*, the cell‐level parameters responsible for this competitive advantage have yet to be conclusively defined. Either LSCs are assumed to be unresponsive to an intrinsic limit in the stem cell reserve size (Michor *et al*, 2005) or to have an increased amplification rate (Roeder *et al*, 2006), both possibilities being conceptually consistent with our model. Unresponsiveness to inhibitory signals would allow the LSCs to expand above their normal physiological limit, resulting in an increased amplification rate and selective advantage against normal HSCs.

The large uncertainties associated with some of the parameter estimates indicate that the model is parametrically non‐identifiable, an outcome consistent with the structure of the cellular network (i.e. ratios of competitive interactions may be conserved rather than their characteristic values). Non‐identifiability appears to be a universal feature of many systems‐level models, and does not necessarily limit their predictive value (Gutenkunst *et al*, 2007). Consistent with recent models of intracellular signalling in mammalian (Chen *et al*, 2009) and bacterial systems (Barkai and Leibler, 1997), our analysis indicates that key system behaviours are determined primarily by network structure rather than the precise tuning of kinetic parameters.

In summary, our intercellular feedback model of blood stem cell growth predictively simulates the dynamic characteristics of both normal and malignant hematopoiesis *in vitro*. This model may, therefore, serve as a platform for further experimental interrogation of the regulatory mechanisms controlling stem cell fate *in vitro* and *in vivo*, and as a tool for the rational design of stem cell‐therapy bioprocesses.

## Computational methods

The hematopoietic hierarchy can be divided into a number of discrete compartments, from long‐term repopulating HSCs (LT‐HSCs) to fully differentiated mature cells. Each compartment can be viewed as representing a cell population at a distinct state of maturation, with unidirectional transition between compartments (differentiation) associated with cell cycling. Although recent evidence shows that LT‐HSCs may undergo functional transitions (differentiation) before mitosis (Kent *et al*, 2008), differentiation is generally co‐incident with population amplifying cell divisions. For simplicity, lineage specification is considered only after differentiation (Lin^{+}). A cell population balance can be constructed around each compartment (*i*) in which the number of cells in the compartment (*X*_{i}) is dependant on the number of cells entering from the earlier compartment (*X*_{i−1}), the cell proliferation rate (*u*_{i}), and the probability of self‐renewal (*f*_{i}). This is a deterministic model, which nonetheless can be viewed as incorporating stochastic elements whose impacts are negligible at the population level.

The cellular growth rate for compartment *i* is given by the equation:

A system of ordinary differential equations (ODEs) is therefore constructed, which describes the growth of each cellular compartment for a total of *n* compartments, with compartment 1 (*X*_{1}) representing LT‐HSCs (*X*_{0} is non‐existent, hence set as=0), and terminally differentiated mature cells represented by compartment *n* (*X*_{n}). Specific compartments can be ascribed to experimentally measurable cellular assays. The functional measures considered are long‐term NOD‐SRC, LTC‐IC, and CFCs, which readout stem cells, primitive progenitors, and mature progenitors, respectively. We additionally characterise the cells phenotypically as undifferentiated Lin^{−} or differentiated Lin^{+}. Starting with a stem cell (*i*=1), the number of cells in compartment *i* that can be generated through symmetric differentiation divisions is given by:

And subsequently the total number of cells (*X*_{T}) in all compartments, up to and including differentiation stage *n* generated from a stem cell is given by:

Hence the total number of compartments (*n*) readout in a given functional assay can be estimated through rearrangement of equation (3):

### Model assumptions

We assume that the SRC assay quantifies LT‐HSCs with absolute accuracy, although this assumption may be disputed. Estimates derived directly from equation (4) based on experimentally measured frequencies in UCB TNCs are presented in Table I. Only the last compartment within a given population is considered, and the first compartment is set to 1 (LT‐HSC). Although LT‐HSCs will not readout in the CFC assay, and it is unclear whether LT‐HSCs readout in the LTC‐IC assay (Coulombel, 2004), placement of the first compartment has negligible impact on resulting calculations. Similarly, accounting for the self‐renewal potential associated with HSCs and immediate descendants will increase the calculated total number of compartments (*n*); however, a few additional divisions in the early stages of the hierarchy will have negligible effect on the population balance calculations. The simplifying assumption that *X*_{i+1}=2*X*_{i} will likely not hold true across all compartments for a system in steady state, and this amplification coefficient may in fact be a dynamic parameter. However, in lieu of experimentally defined *in vivo* kinetic parameters, it is a reasonable estimate supported by our analyses (see Supplementary text and Supplementary Figure S4) and other experimental data. It is notable that our calculation of *n*=20 is consistent with earlier theoretical and experimental‐based estimates ranging from 17 to 30 (MacKey, 2001; Shochat *et al*, 2002).

The self‐renewal probabilities (*f*_{i}), proliferation rates (*u*_{i}), and cell death rates must be specified for each compartment to solve the cell population balance ODEs defined in equation (1). As these internal variables are inaccessible to experimental measurement, we chose to estimate parameters based on generalised functions and a reverse engineering strategy to estimate specific values within experimentally or biological constraints.

Stimulation with high concentrations of hematopoietic growth factors such as stem cell factor (SCF), Flt3 ligand (FL), and thrombopoietin (TPO) are known to inhibit apoptosis, and cultures are maintained such that concentrations of nutrients (i.e. amino acids) and waste products (i.e. lactate) are not limiting. It is, therefore, reasonable to assume the cellular death rate to be negligible.

It has been documented that the proliferation rate (*u*_{i}) of hematopoietic cells varies with the stage of maturation such that the proliferation rate (*u*_{i}) of mature progenitors (CFC)>primitive progenitors (LTC‐IC)>stem cells (SRC) (Mobest *et al*, 1999). The *in vitro* conditions under study do not allow for the functional maturation of differentiating cells, and this can be accounted for by a setting a lower proliferation/differentiation rate for the Lin^{+} cells. We define the proliferation rate of Lin^{−} cells as a function of compartment number by a Gaussian‐type function:

where *n*_{MAX} is the compartment with the maximum proliferation rate (analogous to the mean), *u*_{MAX} is the proliferation rate of this compartment, and *D*_{GR} is the growth rate decay term (analogous to the variance), which defines the steepness of the change in cycle rates between successive compartments. Physiological limits bound all the three parameters between 1 and 10.

By definition, only the LT‐HSC population has the capacity for self‐renewal probability (*f*_{i}) to exceed 50% *in vivo*; however, ST‐HSCs must by definition also have some capacity to self‐renew (<50% *in vivo*), and downstream progenitor populations have also been documented to undergo limited self‐renewal divisions (Marley *et al*, 2003). Self‐renewal probabilities (*f*_{i}) should, therefore, diminish with the stage of differentiation, as described theoretically in Roeder and Loeffler (2002). We define the probability of self‐renewal (*f*_{i}) as a function of compartment number by a Gaussian‐type function with maximum set at *i*=1:

where *f*_{MAX} is the maximal self‐renewal probability of the LT‐HSC compartment (*X*_{1}), (limited to [0 1]) and *D*_{SR} is the self‐renewal decay term (analogous to the variance), which defines the steepness of the change in self‐renewal probabilities between successive compartments, bounded between 1 and 10 based on physiological limits.

### Explicit time‐dependent terms

When purified hematopoietic progenitor cells are placed in culture, there is an initial lag time before the cells enter cycle (G_{0} to S/G_{2}/M/G_{1} transition), after which cells cycle with approximately constant doubling time. A lag phase is therefore introduced into the system through the explicit time‐dependant function:

where *t* is the culture duration (days), τ_{D} is the time for 50% of the cells to enter cycle, and *kt* is the Hill coefficient defining the rate at which cells are induced to cycle (approaching a step function at τ_{D} as *kt* → ∞). On the basis of experimental observations (Ko *et al*, 2007), setting τ_{D}=2 days and *kt*=4 produces a reasonable kinetic response. For simplicity, this lag‐phase term is applied equally to all the cell compartments, whereas in reality there would likely be differences among HSCs, progenitors, and mature cells in their rate of entry to cycle on stimulation (Punzel *et al*, 2002).

### Non‐linear terms

It was our goal to incorporate as few structural assumptions as possible into the model so as to minimise any systemic bias. For models containing few compartments (i.e. 3 or 4), it is possible to use random methods to structure the intercellular regulatory relationships (positive feedback, negative feedback, feed‐forward), followed by optimisation algorithms to select the top performing network topologies (Socolovsky *et al*, 2007). However, for a 20‐compartment model this approach would be computationally unfeasible due to the combinatorial explosion of possible topologies. Hence, we used the wealth of experimental and theoretical literature on hematopoietic regulation to define a candidate intercellular network topology.

Murine HSC transplantation studies demonstrate that numbers of stem cells, progenitors, and total cells in the bone marrow and circulation are regulated through feedback control mechanisms. Feedback control is inferred from observations that there exist non‐autonomous established ‘set points’ for cell numbers in each compartment, such that cellular expansion *in vivo* correlates negatively with cell dose transplanted (Iscove and Nawa, 1997). More recently, *in vivo* imaging studies have demonstrated directly that normally quiescent LT‐HSCs are induced to rapidly enter cycle and self‐renew in response to chemical or radiation‐induced progenitor depletion, returning to quiescence after regeneration of the bone marrow (Wilson *et al*, 2008; Xie *et al*, 2009). Appropriate hematopoietic regeneration requires independent regulation of stem cell proliferation and self‐renewal, and evidence suggests that these processes are in fact regulated independently by distinct cell compartments *in vivo*; bone marrow progenitors (CFU) inhibit self‐renewal independent of proliferation effects (Blackett and Botnick, 1981), whereas mature cells in the circulation normally suppress stem and progenitor cell proliferation in the bone marrow (Cheshier *et al*, 2007). However, under appropriate circumstances (i.e. infection), circulating white blood cells can directly induce the proliferation of HSCs through the secretion of inflammatory cytokines (Essers *et al*, 2009).

Mature cell populations in fact display lineage‐specific functional effects on stem and progenitor cells; platelets (Foss *et al*, 2008) and NK cells (Fardoun‐Joalland *et al*, 1994) both secrete factors that enhance progenitor expansion, whereas macrophages (Xu *et al*, 2000) and red blood cells (Cheshier *et al*, 2007) secrete factors that inhibit progenitor expansion. Stem and progenitor cell fate decisions are regulated by a balance of antagonistic stimulatory versus inhibitory soluble factors (Cashman *et al*, 1990; Jacobsen *et al*, 1994) produced by distinct cell populations in the microenvironment (Wright *et al*, 1979). Importantly, these feedback mechanisms seem to be (at least partially) reconstituted *in vitro*, as progenitor output in both stroma supported (Oh *et al*, 1994) and liquid suspension cultures (Kohler *et al*, 1999) is enhanced by frequent media exchange, cell‐density reduction, and progenitor re‐selection (Gilmore *et al*, 2000; Flores‐Guzman *et al*, 2006).

The intercellular regulatory network topology described in Wichmann and Loeffler (1985) consisting of three inter‐related negative feedback control loops is capable of simulating *in vivo* hematopoietic response to multiple perturbations. The topology consists of (1) auto‐regulatory HSC feedback, (2) intramedulary short‐range feedback, and (3) blood–bone marrow feedback. Although the mechanisms of feedback control are suggested to be soluble factors, the model is formulated by setting proliferation and self‐renewal rates as directly responsive to cell densities, implying a linear relationship between cells and secreted factor concentrations. We expanded on this concept by considering proliferation rates (*u*_{i}) and self‐renewal probabilities (*f*_{i}) as regulated by a balance of secreted inhibitory and stimulatory factors (SF1‐4). Mature (Lin^{+}) cell subpopulations secrete factors (*SF1*) that inhibit the proliferation of progenitors (Lin^{−}), and progenitors (Lin^{−}) in turn secrete factors (*SF2*) that inhibit the probability of self‐renewal of stem and primitive progenitor cells. These negative feedback loops are coupled as *SF2* secretion is induced by SF1. We additionally consider two factors secreted by mature (Lin^{+}) cell subpopulations, SF3 and SF4, which stimulate progenitor proliferation and self‐renewal, respectively (positive feedback). This regulatory architecture is described by the following set of ODEs; secretion of *SF1* through *SF4* is given by:

As factors *SF1* through *SF4* are theoretical, the corresponding terms *L*_{1} through *L*_{4}, representing ligand concentrations producing half‐maximal dose response, are arbitrary and hence set to 1 for simplicity. Equations (13) and (14) then reduce to:

We can then set practical bounds on secretion rates sr1 though sr4 based on dynamic dose–response relationships. For cell densities in the range of 5 × 10^{5} cell/ml, secretion rates of 10^{−7} to 10^{−5} pg/cell/day result in 50% maximal response reached between 0.5 and 50 days, and we similarly set Ls=10^{−2} to 10^{1}. Hill coefficients (*k1–4* and *ks*) for the terms are set between 0.1 and 10, covering the range of relevant sigmoidal dose–response curves. Combining equations (5), (7), and (14), and (6) and (15) to define *u*_{i} and *f*_{i}, respectively, as functions of compartment number (*i*), time (*t*) and secreted factor concentrations (SF1–4) generates the functions:

Such that (0⩽*u*_{i}⩽*u*_{MAX}) and (0⩽*f*_{i}⩽*f*_{MAX}). It should be noted that *SF1* through *SF4* most likely represent groups of ligands rather than single factors. Additionally, secretion of *SF1*, *SF3*, and *SF4* by Lin^{+} cells can be equivalently interpreted as a single homogenous population secreting all three factors, three distinct populations each secreting a single factor, or more likely, a combination thereof. Within the heterogeneous Lin^{+} population, there are likely subpopulations biased towards *SF1*, *SF3*, and SF4‐type secretion profiles; however, this additional complexity is not considered. Proliferation rates and self‐renewal probabilities of the progenitor populations will hence be dependent on the exogenous growth factor(s) provided (setting *u*_{MAX}, *u*_{+}, *n*_{MAX}, *f*_{MAX}) and the balance between inhibitory and stimulatory factor accumulation in the culture microenvironment (*SF1*:*SF3*, *SF2*:*SF4*).

### Parameter estimation

The resulting model contains 16 undefined parameters (*P*_{1} through *P*_{16}) (see Table II). Data from Madlambayan *et al* (2005) (see Table III) was used as a training set to estimate parameters through a reverse engineering strategy. The proportion of differentiated (Lin^{+}) cells produced (%Lin^{+}) and cell population fold‐expansion values (*TNCX*=total nucleated cell, *CFCX*=colony forming cell, *LTC‐ICX*=long term culture‐initiating cell, *SRCX*=scid repopulating cell) represent observable functions (*O*_{i}) of system variables (*X*_{i});

where *t*_{0} and *t*_{f} represent culture initiation and analysis time points, respectively. The observable variables may be represented as terms in the observation vector (*O*):

Which defines the system output (*O*). A cost function (*J*) to be minimised was defined as a weighted (*w*) sum of squared residuals between experimental (Oe) and simulated (Os) system outputs:

A total of 20 system outputs (%Lin^{+}, TNC, CFC, LTC‐IC, and SRC expansion × 4 culture conditions) are therefore used to fit the 16 internal parameters. Each term in the equation *J* is normalised (division by Oe^{2}), and a weighting factor (*w*_{i}) is included (set to 1 for all terms initially) introducing a user‐defined bias based on which terms are deemed more or less important.

Because of the highly non‐linear and multi‐modal nature of the cost function (*J*), a hybrid stochastic method was used to solve this non‐linear programming problem. A genetic algorithm was first used as a global method to search the entire parameter space, followed by the Nelder–Mead Simplex algorithm to search for a local minimum in the vicinity of output from the genetic algorithm. These were implemented using the MATLAB *gatool* and *fminsearch* functions, respectively. Multiple settings for the genetic algorithm were tested and the following were selected, producing a solution that converges within 150 generations; population size=160 [10 × number of parameters (Np)], initialised as a random log‐uniform distribution; elite count=8; crossover fraction=0.5; mutation rate=0.5. After model validation studies, objective functions were created based on the new data, and a generalised pattern search algorithm (*psearchtool*) was used to adjust parameter values within 95% confidence intervals (±2 s.d.) of the estimates defined in Table III.

### Model analysis—local parameter sensitivities

To understand the relationship between system responses and variations in individual parameter values, local parameter sensitivity analysis was performed. The sensitivity coefficient (*S*) is defined as:

which is defined for each system output (*O*_{i}) and system parameter (*P*_{j}). Individual parameters were altered individually by 1% (Δ*P*=0.01) from their estimated values, and resulting changes in system outputs (Δ*O*) were determined. The resulting expression essentially denotes the percentage change in output *O*_{i} resulting from a 1% change in parameter *P*_{j}. This analysis produces a 5 × 16 sensitivity coefficient (*S*) matrix.

To quantify relationships between parameters and system outputs, this matrix was converted to a heat‐map image for visualisation, and 2D, unsupervised hierarchical clustering used to organise rows (*P*_{j}) and columns (*O*_{i}). Clustering and visualisation was performed using dChip v2006 software (www.dChip.org). Sensitivity values were normalised (mean=0, s.d.=±1, across rows) and rows/columns clustered through the centroid‐linkage method with a distance metric of 1−*r* (Pearson correlation coefficient) as described in Eisen *et al* (1998).

### Model analysis—parameter identifiability

Parameter estimation accuracies are central to measuring identifiability of mechanistic models. The Fisher information matrix and Cramer–Rao theorem are commonly used to estimate the lower bound of parameter estimation errors. However, this approach assumes that the model is linear with respect to parameters, whereas our model is highly non‐linear. We, therefore, implemented a bootstrapping approach as described in Kremling *et al* (2004), which is more computationally intensive but requires no underlying assumptions. On the basis of experiments and results described in Figure 2, we generated a synthetic data set of 50 experiments from the distributions reported in Table III. We then ran our parameter estimation algorithm on each individual data set, resulting in 50 estimated parameter vectors (*P*). We filtered the results for parameter vectors (*P*) producing cost functions (*J*) <5, and characterise the distribution of individual parameter estimates (*Pi*) by their CV. Non‐identifiable parameters are defined as those for which it is not possible to determine with 95% confidence (estimated as ±2 s.d.) that their values are non‐zero (Zak *et al*, 2003).

### Model analysis—structural discrimination

The large number of non‐identifiable free parameters persuaded us to examine the sensitivity of model simulations to structural alterations of the regulatory network. The regulatory architecture is based on theoretical and experimental evidence rather than an exhaustive search of all possible topologies due to computational limitations. We chose to systematically ‘re‐wire’ intracellular network connections by perturbing the sign of interactions (stimulatory↔inhibitory) while conserving the number of free parameters. The following perturbations are defined as A, B, and C:

A: SF2 effect on *f*_{i}—inhibitory → stimulatoryB: SF4 effect on *f*_{i}—stimulatory → inhibitoryC: SF1 effect on sr2—stimulatory → inhibitory

These perturbations were implemented alone and combinatorially, and the resulting structurally altered models defined as S1 through S7:

S1: AS2: BS3: CS4: ABS5: ACS6: BCS7: ABC

Perturbing the sign of interaction for the effect of SF1 and SF3 on progenitor proliferation rate (*u*_{i}) is equivalent to simply deleting one of the feedback connections; hence, we considered two additional structural alterations that have the effect of reducing the number of free parameters:

S8: delete SF1 functional effects (k1=−∞)S9: delete SF3 functional effects (sr3=0, k3=0)

We ran the parameter estimation algorithm on each of the structurally altered models (S1–S9) in addition to our original model [control (*C*)] using an extended objective function (*J*) incorporating experimental data reported in Figures 3B, 6A, D, and 7A. The Akaike Information Criterion (AIC) was then computed for each:

where Np=number of parameters and Ne=number of independent experimental measurements. Alternate models are then ranked, with the lowest AIC corresponding to the model best able to describe the data with minimum free parameters (Landaw and DiStefano, 1984).

### Software

All model simulations and computational analysis was performed using MATLAB R2008a software (The Mathworks, Natick, MA, USA), and differential equations were solved using the non‐stiff numerical solver *ode23* with default error tolerances. To ensure that results were not affected by the propagation of rounding errors, select simulations were performed using the non‐stiff solver *ode45* and the stiff solver *ode15s* with different error tolerances, and the results were unaffected.

## Materials and methods

### Cell sample collection and processing

UCB samples were obtained from consenting donors according to procedures accepted by the ethics boards of Mt Sinai hospital (Toronto, ON, Canada), Joseph Brandt hospital (Burlington, ON, Canada), and Credit Valley Hospital (Mississauga, ON, Canada). Mononuclear cells were obtained by first mixing the UCB sample with 10% pentastarch (Bristol‐Myers Squibb Canada, Montreal, QC, Canada) at a 1:5 volumetric ratio. The sample was then centrifuged for 10 min at 50 *g*, and the upper (leukocyte rich) plasma layer was removed and centrifuged for 10 min at 400 *g* to obtain a cell pellet. Red blood cells were depleted by suspending the cells for 10 min in red blood cell lysis buffer (0.15 M NH_{4}Cl, 0.01 M KHCO_{3}, 0.1 mM EDTA). Lineage depleted (Lin^{−}) cells were isolated from the mononuclear cell fraction using the StemSep system (Stem Cell Technologies, Vancouver, BC, Canada). This process depletes cells expressing cell surface antigens CD2, CD3, CD14, CD16, CD19, CD24, CD56, CD66b, and CD235a.

### Cell culture

UCB Lin^{−} cells were seeded at 10^{5} cells/ml (unless otherwise noted) in serum‐free Stem Span media (Stem Cell Technologies, Vancouver, BC, Canada) supplemented with 100 ng/ml SCF (Amgen, Thousand Oaks, CA, USA), 100 ng/ml FL (Amgen), and 50 ng/ml TPO (R&D Systems, Minneapolis, MN, USA), 1 μg/ml low‐density lipoproteins (Calbiochem, La Jolla, CA, USA), and penicillin‐streptomycin at 100 U/ml and 100 μg/ml, respectively (Invitrogen, Carlsbad, CA, USA). The cell suspension either injected into a cell culture bag of appropriate volume (2, 7, or 12 ml) through the self‐sealing rubber septum using a sterile syringe attached to a threaded cannula, or placed into wells of a 24‐well tissue culture plate (Corning Inc, Corning, NY, USA). Cultures were carried out using a culture bag‐based bioprocess (described below) for experiments requiring Lin^{+} cell depletion, otherwise 24‐well tissue culture plates were used. Cultures were maintained on an orbital shaker at 37°C in a humidified atmosphere of 5% CO_{2} in air.

### Subpopulation selection and media dilution

The bioprocess used in these studies for depletion of *in vitro*‐generated Lin^{+} cells was described in Madlambayan *et al* (2006). It consists of two gas‐permeable cell culture begs connected through a magnetic selection element, used to remove the Lin^{+} cells (or any other antibody‐labelled cell subpopulation). The system is completely closed, sterile, autoclavable, and disposable (single use), making it attractive for clinical applications. Cell selection in the bioprocess was performed in a similar manner to the StemSep system using the reagents provided with the kit, as described in Madlambayan *et al* (2006). The cell culture bag was flushed with Lin^{+} antibody cocktail and magnetic colloid (dextran‐coated iron particles) as per manufacturer's instructions. This effectively attached the dextran‐coated iron particles to the Lin^{+} cells. The cell culture was subsequently allowed to flow the selection element that was placed in a magnetic field, allowing the iron‐labelled Lin^{+} to be retained in the element, whereas non‐labelled Lin^{−} cells flow through to the secondary culture bag. For flow rate control, a peristaltic pump was attached upstream of the cell culture bag using a septum/cannula connection, and used to drive the cell solution through the selection element at a flow rate of 1.3 ml/min.

Media dilution was performed on the enriched Lin^{−} cells by removing the secondary culture bag from the selection element and placing it into a 50 mL conical centrifuge tube. Paper batting was used to stabilise the culture bag during centrifugation. The bag was centrifuged for 7 min at 200 *g*, after which a cell pellet was visible at the bottom. Conditioned media was then removed through the self‐sealing septum using a sterile syringe and fresh media was added through the same septum. The cell culture bag was then placed back into the incubator.

### Phenotypic analysis

Staining for Lin^{+} marker expression on culture‐generated cells was accomplished by suspending 5 × 10^{4} cells in 100 μl ice cold Hank's balanced saline solution containing 2% (v/v) human UCB serum (HBSS‐HS). The cells were then incubated with Lin^{+} antibody cocktail followed by magnetic colloid (dextran‐coated iron beads) as described, washed twice in HBSS‐HS, and finally stained with saturating amounts of FITC‐labelled anti‐dextran antibody (Stem Cell Technologies) for 30 min on ice. For isotype controls, the Lin^{+} antibody incubation step was not performed. All samples were washed in HBSS‐HS and stored on ice before analysis either on a FACSCanto (BD Biosciences, San Jose, CA, USA) or Coulter EPICS XL (Beckman Coulter, Fullerton, CA, USA) flow cytometer.

### Progenitor cell assays

Cells were assayed for CFC frequency by plating 500 cells into 1.5 ml methylcellulose‐based medium (MethoCult H4434, Stem Cell Technologies) containing 1% methylcellulose in Iscove's Modified Dulbecco's Medium, 30% fetal bovine serum, 1% bovine serum albumin 10^{−4} M 2‐mercaptoethanol, 2 mM l‐glutamine, 50 ng/ml SCF, 10 ng/ml granuloctye‐macrophage colony stimulating factor, 10 ng/ml IL‐3, and 3 U/ml erythropoietin. After 14 days of incubation at 37°C in a humidified atmosphere of 5% CO_{2} in air, duplicate cultures were visually scored for CFC content (colony number and lineage composition). Cells were assayed for LTC‐IC frequency by plating 10^{3} freshly isolated Lin^{−} cells, or 5 × 10^{3} to 2 × 10^{4} culture‐generated cells in triplicate onto irradiated (6000 rad) murine stromal cells (M2‐10B4) on collagen‐coated six‐well plates in MyeloCult H5100 medium containing 10^{−6} M hydrocortisone (Stem Cell Technologies). After 5 weeks of culture at 37°C with weekly half‐media exchanges, the contents of each well were harvested using 0.25% Trypsin with 0.38 g/l EDTA in HBSS (Invitrogen, Carlsbad, CA, USA), and plated into methylcellulose‐based media. LTC‐IC content was determined by enumerating CFCs present after 14 days of incubation.

### Conditioned media proteome analysis—luminex liquid chips and ELISA

Conditioned media samples were assayed in triplicate using the Biosource Human Cytokine 30‐Plex detection kit (Invitrogen, Carlsbad, CA, USA). These kits use Luminex (Luminex Co, Austin, TX, USA) microspheres as a fluid platform for multiplex sandwich ELISA. The ‘microspheres’ consist of 5 μm polysytene beads bar‐coded through unique ratios of APC: APC‐Cy7 dye. Each colour‐coded microsphere contains primary capture antibody against an individual ligand, which in combination with secondary PE‐conjugated detection antibody, can be used to quantify the concentration proteins in a test sample through flow cytometry. Briefly, 50 μl antibody‐coated microspheres per test were suspended and aliquoted into wells of a 96‐well filter plate (Millipore, Billerica, MA, USA) and washed in wash buffer. A measure of 50 μl of incubation buffer plus 100 μl of culture media sample (or appropriately diluted standard) was then added to each well, and the plate was covered and incubated on an orbital shaker at room temperature for 2 h. After sample incubation, the wells were washed and 100 μl of detection antibody mixture was added to each well and incubated on an orbital shaker at room temperature for 1 h. After incubation, wells were washed and incubated with 100 μl streptavidin‐PE for 30 min. After this final incubation, wells were washed and resuspended in 5 ml polystyrene falcon tubes (BD Biosciences, San Jose, CA, USA) for analysis using a BD FACSCanto flow cytometer (BD Biosciences, San Jose, CA, USA). APC versus APC‐Cy7 gates for each of the 30 cytokine‐specific bead regions were defined based on manufacturer's datasheets, and cytokine signals quantified by comparing PE fluorescence levels to standard dilution curves prepared using the Human Cytokine 30‐Plex Standard mixture provided. Cytokine‐specific standard curve equations were developed by fitting data to a four‐parameter logistic curve:

where PE is the raw PE‐fluorescence intensity, *x*_{i} is the cytokine concentration, and *A*, *B*, *C*, and *D* are the unknown parameters, fit using the Nelder–Mead Simplex non‐linear optimisation algorithm.

TGF‐β1 is not one of the molecules included in the kit, hence conditioned media samples were separately analysed using a human TGF‐β1 Quantikine ELISA Kit (R&D Systems, Minneapolis, MN, USA) as per manufacturer's instructions.

## Acknowledgements

We thank F Notta and J Dick for providing the raw data from Warner *et al* (2005) and for helpful discussions. Support for this work was provided by Insception Biosciences, the Guggenheim Foundation, the Canadian Stem Cell Network, and the Natural Sciences and Research Council of Canada. DCK is an NSERC Post Graduate scholar and PWZ is the Canadian Research Chair in Stem Cell Bioengineering.

## Conflict of Interest

The authors declare that they have no conflict of interest.

## Supplementary Information

Supplementary text [msb200949-sup-0001.doc]

Supplementary Figures S1–4 [msb200949-sup-0002.pdf]

Supplementary Figure Legends S1‐4 [msb200949-sup-0003.doc]

Supplementary Table S1 [msb200949-sup-0004.doc]

## References

This is an open access article under the terms of the Creative Commons Attribution‐NonCommercial‐NoDerivs License, which permits use and distribution in any medium, provided the original work is properly cited, the use is non‐commercial and no modifications or adaptations are made.

- Copyright © 2009 EMBO and Nature Publishing Group