The process of cellular differentiation is central to our understanding of the nature of multicellular living systems, their stability in a changing environment, and how such systems fail in diseases, such as cancer [
1,
2]. This developmental process of individual cells in a multicellular organism committing to their specialized phenotypic characteristics is temporally coordinated by a complex dynamical system comprised of large numbers of interacting genes and their products [
3-
5]. Not surprisingly, dynamical systems theory has been used to study cell differentiation [
6-
8].
Despite its tremendous importance, there is very little accumulated knowledge of the process of differentiation from a systems perspective and of the role of molecular programs involved in this process. Even for an agent that causes differentiation to a common recognizable state, we do not know whether the cells, as manifestations of the underlying dynamic bio-molecular systems, always follow common or different molecular paths (or system state trajectories). In the latter case, we also do not know which of those paths is the most stable and least reversible.
Since a cell's phenotype and behavior are largely determined by the activities of the genes and proteins constituting a genetic network, it follows that the rules of interactions between these elements translate directly into rules of cellular behavior. That is, the enormous state space of a genetic network (i.e., the space of all possible configurations activities of the constituent elements) becomes reduced into a relatively small number of trajectories and steady states (attractors) of the dynamical system. Kauffman postulated that these attractor states in model networks correspond to the cell types in multicellular organisms [
9,
10], and the process of differentiation corresponds to a trajectory (in the state space) leading into one of the attractors. The cellular fate is thus determined by the attractor in which the genetic network eventually ends up; this can, to a large extent, be controlled by appropriate external stimuli that place the system into different initial states. It is important to note that many trajectories ensuing from such different initial states can flow to a common attractor and thus constitute its basin of attraction.
Consider that small molecule chemicals, such as dimethyl sulfoxide (DMSO) and a host of others can induce cell differentiation in a variety of cell systems along with concomitant cellular properties [
11-
15]. This rather amazing fact implies the pre-existence of cellular fates that need only be selected by means of external stimuli rather than created by specific molecular events. This 'selection' of cell fates occurs by means of the inherent nature of the dynamical system to flow to an attractor when placed in some initial transient state and thus, differentiation is a process of selecting a particular attractor in a genetic regulatory network. This view has been supported experimentally by Huang
et al. using genome-wide mRNA expression profiling [
16] as well as by means of analyzing cell fates in response to generalized physical stimuli, such as cell distortion [
17]. For a more extended discussion on this topic, see [
10].
The homeostatic stability of a differentiated cell is a consequence of the underlying stability of the attractor – 'nearby' states, which may occur as a consequence of natural environmental variation, simply flow back to the attractor. It is known that normal cells have a balanced state of proliferation and differentiation, resulting in homeostatic stability [
18,
19]. A block of normal differentiation and abnormal reversal of differentiation (sometimes called de-differentiation) [
20] are believed to be some of the hallmarks of cancer [
21]. Accordingly, therapeutic strategies have been designed to facilitate cancer cells to reenter the differentiation program, often termed differentiation therapy [
22,
23].
The success of such therapeutic strategies depends on our ability to systematically determine appropriate molecular 'lever points', the perturbations of which place the biomolecular system into states that are poised to differentiate. Indeed, such a strategy corresponds to placing the system in a state by means of a stimulus, such as a therapeutic agent, and allowing the system to naturally flow toward an attractor that corresponds to the desired cellular endpoint [
24-
26]. To identify such targets for intervention, it is necessary to characterize the underlying molecular mechanisms, such as transcriptional regulatory networks, governing the process of differentiation. Systems biology approaches, which are predicated on global measurements and data integration, are now beginning to reveal transcriptional machinery underlying complex biological processes [
27-
30]. The rationale behind our study was to explore whether the aforementioned systems-level view of cell fates as attractors and differentiation as a route toward an attractor, when coupled with computational systems biology approaches, is informative for elucidating the transcriptional regulatory mechanisms governing differentiation.
To this end, we have selected a well-established differentiation model, human promyelocytic leukemia cells (HL60) originally isolated by Dr. Steven Collins from a 37-year-old female acute promyelocytic leukemia (APL) patient [
31]. The HL60 is a multi-potent cell line that can be stimulated to differentiate using a variety of chemical agents, including DMSO [
32], all-
trans-retinoic acid (ATRA) [
33], 1,25
α-dihydroxyvitamin
D3 [
34], 12-O-tetradecanoylphorbol 13-acetate (TPA) [
35], and granulocyte macrophage colony-stimulating factor (GM-CSF) [
36]. With the addition of ATRA, the HL60 cells differentiate into neutrophils, while displaying the early differentiation marker, CD11b, which begins to be expressed within one day of treatment [
37]. Although there are others, CD11b is an early differentiation marker, which allows one to capture the initial stage of the process. The CD11b+ differentiated HL60 cells can be stained with fluorescent-labeled anti-CD11b antibody and easily recognized by commonly used flow-cytometry methods and isolated by flow-sorting for further culturing and experimentation, as we have done here. The HL60 system was also used by Huang
et al. [
16] to demonstrate the correspondence between cell fates and high-dimensional attractor states of the underlying genetic network.
One could reason as follows. If we could place the HL60 into a state from which the system would dynamically flow towards the "neutrophil" attractor, as demonstrated by Huang
et al., then the genes that show altered behavior along the time-course trajectory relative to unstimulated controls could be hypothesized to be implicated in the neutrophil differentiation process. This, of course, may be the case, though the interpretation is confounded by the possibility that the genes exhibiting altered behavior are responsive to the particular mechanisms activated by the stimulus used, such as ATRA. Indeed, Huang
et al. also confronted this conceptual difficulty when they compared trajectories from ATRA-treated and DMSO-treated HL60 cells, finding that certain genes may behave differently simply as a result of different stimuli activating different biological pathways, while many other genes dynamically converge towards a common attractor, despite the system flowing from distinct starting states corresponding to ATRA and DMSO treatments [
16]. To identify genes that are not stimulus dependent, but are involved in the process of neutrophil differentiation, one could then use only one treatment, but in a way that allows one to alter cellular fate, namely, terminal differentiation into neutrophils or reversion back to the undifferentiated state.
The HL60 was shown to exhibit such behavior in two separate studies both demonstrating that this differentiation process contains at least two steps in which a precommitment stage precedes the decision to differentiate. Yen
et al. observed that with continuous exposure of ATRA at a high concentration, the HL60 proceeds through differentiation, but upon removal of the stimulus, the HL60 falls back to the undifferentiated state [
38]. By analogy, such a precommitment stage corresponds to a gradually sloping plateau between a valley and a mountain such that a ball sitting on this plateau would roll down into the "undifferentiated" valley in the absence of additional energy necessary to make it over the "terminally differentiated" mountain. More recently, Chang
et al. reported a population of "primed," undifferentiated CD11b- cells upon exposure to a low dose DMSO [
39]. Though these cells are negative with respect to the CD11b marker, thus considered to be "undifferentiated," upon encountering a second dose of DMSO stimulation, they exhibited an increased rate of differentiation, suggesting that the first low dose DMSO had placed them in a "primed" intermediate differentiated state.
We thus decided to determine two different treatments, both with ATRA, but with different concentrations and incubation times such that the two cell populations corresponding to these treatments would be poised at the same stage of differentiation (precommitment), but so that one population follows through to the terminally differentiated neutrophil attractor, while the other reverts back by dynamically flowing towards the undifferentiated state. The genes that would exhibit different behavior between these two trajectories would then be potentially implicated in the differentiation process.
To identify two such precommitment states, we used the percentage of CD11b+ cells at the end of a particular treatment as a measure of the stage of differentiation. We performed 80 ATRA treatments consisting of 8 doses (0.0005 μM to 1 μM) and 10 durations (4 to 13 days) in triplicate and measured percentages of CD11b+ cells, relative to an isotype antibody control, using FACS analysis. Consider loci in the two-dimensional dose × duration stimulus space, where all points on a particular locus correspond to a constant fraction of CD11b+ cells. Thus, two cell populations on the same locus can be said to be at the same "stage" of differentiation at least as it pertains to CD11b. We chose two such populations, one with a higher dose and a shorter duration and the other with a lower dose and a longer duration, such that the cells treated with the higher dose proceeded with differentiation into neutrophils while the cells treated with the lower dose reverted back to the undifferentiated state, despite both populations exhibiting the same percentage of CD11b+ cells at the end of their respective treatments. The cells were live-sorted, cultured in fresh media, and profiled every 24 hours with microarrays for five days in triplicate. This additional Fluorescence Activated Cell Sorting (FACS) step mitigates the confounding effects of cellular heterogeneity due to subpopulations that do not initiate the differentiation program (i.e. CD11b- cells). In this manner, the gene expression programs of the two cell populations, one differentiating and one reverting, could be analyzed using computational approaches.
We defined a criterion to identify genes whose behavior over time exhibits a divergence between the two treatments. It is these genes that are hypothesized to be involved in the differentiation process. We analyzed the promoters of these genes and found that they are overrepresented with known transcription factors functionally linked to myeloid differentiation, cell cycle, and development. This study points to the utility of incorporating a systems-level view of global dynamics, as distinguished from the dynamics or kinetic behavior of the individual elements of a system, into a computational analysis framework that can be used for studying transcriptional regulatory mechanisms governing a complex biological process such as differentiation.