Cellular differentiation and attractor dynamics

In this work, differentiation is viewed as a set of transitions between attractor basins produced by a dynamical genetic regulatory network. This model of differentiation was pioneered by Kaufman and extended by many others [6, 3, 7, 5]. Borrowing from early work by Waddington[8], the landscape created by these attractor basins has been termed an epigenetic landscape [1]. A conceptual model of such an epigenetic landscape is shown in Figure 2. In this view, each cell type occupies an attractor basin at a particular level of a potential energy landscape. A cell can be moved out its attractor basin in response to an external signal or to gene expression noise. Once it crosses the barrier that delimits the basin, it moves down to another attractor basin lower in the epigenetic landscape. There are at least two possible paths leaving each attractor basin, with each downhill path leading to a different basin that represents a distinct, more specialized cell type. Once a cell descends into a new basin, the large potential energy barrier between the new lower basin and upper starting basin makes it unlikely for a more specialized cell to make the transition back to a progenitor cell. This process of cells moving out of an attractor basin in response to external signals or to gene expression noise and descending into attractor basins of lower potential energy that correspond to more differentiated cells is repeated at each level of the lineage tree.

This potential energy barrier that must be crossed to move between attractor basins is called the epigenetic barrier. Schmulevich et al. [9] proposed a method of quantifying this barrier termed the mean first passage time (MFPT), defined as the average number of state transitions needed to move from one attractor basin to another during the noisy operation of a Boolean regulatory network. The MFPT provides a measure of the probability of a particular transition between two attractor basins, with low MFPTs indicating a high likelihood of the transition, and high MFPTs indicating a low likelihood for this transition. Details on the calculation of MFPT values and all other aspects of the procedures are given in Methods; this section will only provide an overview.

The forward and reverse MFPT values between two attractor basins (simply called attractors from this point forward), att₁ and att₂, provide information on the directionality of the transition. Directionality is a key element of differentiation, as under normal circumstances, cells transition from less to more mature states, but not in the reserve direction. For the pair of attractors att₁ and att₂, we define a directional transition to occur if att₁ → att₂ (reaching att₂ from att₁) has a significantly larger MFPT than the MFPT of att₂ → att₁.

Another important aspect of cellular differentiation captured by MFPT is the probability of making a transition between any pair of different cell types. This is important in shaping the structure of a lineage tree and in stabilizing cell types. For example, progenitor cell types should not differentiate into cell types off the normal lineage path, and terminally differentiated cells must be prevented from dedifferentiation or transdifferentiating into other cell types. Therefore, the MFPT should be high in both directions for unfavored transitions between attractors. We term this separation, with high separation occurring when the MFPTs of att₁ → att₂ and att₂ → att₁ are both large.

Given the directionality of differentiation and the large separation of the majority of cell types within a linage tree, a plot of the distribution of MFTPs of the forward (for example, att₁ → att₂) and reverse (att₂ → att₁) transitions between all possible pairs of cell types within a lineage tree is expected to show clustering in the regions of directionality and separation. This is shown in Figure 3. In this plot of forward and reverse MFPTs between all possible attractor pairs produced by a set of gene regulatory networks, the quadrant with a low forward MFPT and high reverse MFPT represents attractor pairs (cell types) that are linked with a strong directional transition. In contrast, the quadrant with high MFTPs in both the forward and reverse directions represents well separated attractor pairs. This region of high separation represents low probability transitions between cells types, such as transdifferentiation or differentiation off the normal lineage pathway. Using this reasoning, if adding a small multistable switch to a larger GRN enhances the directionality of transitions between attractors, then in a plot like the one shown in Figure 3, there should be an increase in frequency of attractor pairs in the regions labeled directional. Similarly, if a multistable switch added to a gene regulatory network increases the separation between pairs of attractors, then there should be an increase in the region of Figure 3 labeled separate. This is the basis of the approach followed in this work.

An important point to note is that in a MFPT representation of biologically realistic lineage trees, the proportion of attractor pair transitions in the separate region will far exceed the proportion in the directional quadrant. This is because the topology of actual linage trees leads to there being significantly fewer directional transitions than well separated transitions. Intuitively, this stems from the ideas that the number of favored transitions between different cell types is much smaller than the number of theoretically possible transitions, and that most of the theoretically possible transitions are unfavored events such as dedifferentiation and transdifferentiation. Mathematically, the possible number of well separated transitions is on the order O(b^2h) while the number of directional transitions is of the order O(b ^h ), where b is the branching factor of differentiation tree (number of children for each node) and h is the height of the tree measured as the number of cell type transitions between a stem cell and a terminally differentiated cell. This expected difference in the proportions of separate and directional attractor pair transitions is important when interpreting the effects of adding multistable switches to random Boolean genetic regulatory networks (see below).

We investigated how the addition of the four multistable switches shown in Figure 4 influenced the attractor landscape produced by randomly generated Boolean regulatory networks. A conventional node-and-edge diagram of each multistable switch used in biological literature is depicted in the figure, followed by a more informative logic circuit representation. The first logic circuit (Figure 4.a) is usually referred to as a bistable switch (BS) or toggle switch[10]. We call the second logic switch (4.b) a mutual inhibition switch (MI00). Note how the less informative node-and-edge diagrams for these two distinct logic circuits are identical. The next two multistable switches extend mutual inhibition with the addition of one (MI+0) or two (MI++) positive feedback loops. MI++ is sometimes referred to as tristable switch.

Multistable switches in myeloid differentiation

An important example of cellular diversification is the well studied system of hematopoiesis. During hematopoiesis, multipotent stem cells (hemocytoblasts) differentiate into either myeloid or lymphoid progenitors [11]. A sub-tree of the myeloid lineage tree is illustrated in Figure 1. This figure shows that common myeloid progenitor (CMP) cells produce two pluripotent cell types (megakaryocyte-erythrocyte progenitor (MEP) cells and granulocyte-monocyte progenitor (GMP) cells) that in turn produce terminally differentiated erythrocyte, megakaryocyte, monocyte and granulocyte cells.

To construct a GRN that simulates the dynamics of the myeloid differentiation, we extracted a set of regulatory gene expression levels of all cell types in Figure 1 from three datasets of distinct experiments available at ArrayExpress database (http://www.ebi.ac.uk/microarray-as/ae/): E-GEOD-5606, E-GEOD-8407, and E-GEOD-18483. Motivated by Krumsiek et al. [11], we picked 11 transcription factors that play important roles in myeloid differentiation: GATA-1, GATA-2, FOG-1, EKLF, Fli-1, SCL, C/EBPα, PU.1, cJun, EgrNab, and Gfi-1; note that the EgrNab, represents an integration of Egr-1, Egr-2 and Nab-2. Using these genes and their expression profiles, we utilized a search tool to infer a GRN for myeloid differentiation as a Boolean network (manuscript in preparation). This network includes 4 well-known gene interactions that represent multistable switches [11], [1]: a) An MI++ switch between GATA-1 and PU.1; b) An MI++ switch between GATA-2 and PU.1; c) A bistable switch between Fli-1 and EKLF; and d) A bistable switch between Gfi-1 and EgrNab. We computed the MFPT between attractors of this network that represent the cell types of the myeloid lineage tree. The pairwise forward and reverse MFPT values between all pairs of attractors of this network are depicted in the Figure 5 (red circles); we also included the MFPT values for the attractors of the original network proposed by Krumsiek and colleagues (green diamonds) that contains only four attractors as the terminally differentiated cell types. This figure shows that the majority of transitions in myeloid differentiation fall in either the separation or directionality regions shown in Figure 3.

Multistable switches in random networks

We showed that the myeloid differentiation network, with its multistable switches, generates directional transitions and well separated attractors. How general is this result? We extended our study to examine the role of these switches in a large space of cellular differentiation networks.

Density plots were generated for 9 different types of networks: RBN networks without any added multistable switch and RBNs with one or two identical copies of each of the four types of multistable switches. Figure 6 shows these density plots. Each plot shows the forward and reverse MFPT between all attractor pairs generated by 5000 networks of a single type.

The MFPT density distribution produced by RBNs without any added multistable switch (Figure 6a) shows no clustering in the directional or separate regions of the plot. Instead, the forward and reverse MFPTs of most of the transitions are equal and of intermediate values and therefore fall in the mid-range of the diagonal. Adding a single multistable switch of any type to the RBN had a modest effect of increasing the density of attractor pairs in the separate region. Adding two multistable switches of the same type to the RBN had a much stronger effect on increasing the frequency of well separated attractor pairs. This is reflected in an increased density in the separate region of the MFPT plots. The particular kind of multistable switch had little impact on this effect; instead, the critical element was adding two rather than one multistable switch to the RBN.

There was a modest increase in the density of attractor pairs in the directional regions of the MFPT plot when two identical multistable switches were added. However, as discussed above, a major clustering of MFPT values in the directional region is not expected in networks that produce lineage trees. The modest increase in directionality gained by adding multistable switches is likely to be significant. In contrast to the effect on separation, there was a difference between the multistable switch types in increasing directionality: The MI++ switch type did not increase directionality, but all three of the other types did. To better illustrate these enrichments in directional and separate regions, Figure 7 shows the difference between the MFPT distribution of networks with two embedded multistable switches and the base-line random network distribution.

Cell differentiation and attractor dynamics

Boolean networks [6] have proved effective in representing GRN structure and dynamics in many systems, including Drosophila development [12, 13], angiogenesis [14], eukaryotic cell dynamics [15], and yeast transcription networks [16]. Each gene in a network is represented as a node whose regulation by other genes is modeled using updating rules as logic functions. An expressed gene is assigned the value true and a non-expressed gene the value false.

A Boolean network with n genes has 2 ⁿ possible states, denoted as Ŝ. At each step in the simulation, the next state ŝ_t+1∈ Ŝ is determined by applying each gene's logic function (representing the regulatory interactions) to the current value of the genes in ŝ _t . Let this computation be defined as ŝ_t+1← D(ŝ _t ) where D(ŝ _t ) is the deterministic mapping function that finds the next state of the network given the current state. As the network is executed by repeated applications of D(ŝ _t ), the state will reach a previously visited state, and thus, since the dynamics are deterministic, enter into an attractor which represents a fixed point of the system. Attractors can be single states, called point attractors, or consist of more than one state that the network continuously transitions between, called cyclic attractors. Let â = D^∗(ŝ) be the resulting network attractor state reached when starting at ŝ and applying the logic functions until the attractor state â is reached.

In this work, cell types are considered attractors in the state space of possible gene expression profiles [17] and cell differentiation is modeled as the process of transitioning from one attractor to another [18].

Network construction

A random Boolean regulatory network is generated by randomly connecting a varying number of nodes, then instantiating each node with a randomly generated logic function. To replicate networks found in natural systems, we created only networks that operate in the critical domain, rather than ordered or chaotic. Critical networks implement maximal information flow [3] and have the lowest attractor basin entropy [19]. Evidence that GRN's tend to be critical is given in [20]. To generate critical networks, the parameters are set according to s = 2qp _N (1 − p _N ) where s is the sensitivity of the network to perturbations in gene values, p _N is the probability of the output of each Boolean function being 1, and q is the count of inputs to each Boolean function [21]. When s = 1 a single bit change is on average propagated to one other node and the network is in the critical domain. In an ordered network, s < 1 and perturbations tend to die out, while in a chaotic network, s > 1 and perturbations tend to grow. In this work s was fixed at 1 and p _N was adjusted depending upon the value of q.

The attractors of each random Boolean regulatory network are determined, then Markov chain analysis is performed to determine the transition probabilities between all possible states. This allows determination of the MFPTs between each pair of attractors [9]. The MFPTs allow the construction of a graph whose nodes are attractors and weighted edges are the MFPT value between different nodes. Figure 8.a shows a sample graph. MFPT graphs for cellular differentiation are expected to have a small MFPT value for forward edges (moving from less to more specialized cell types), large values for reverse edges, and large values in both directions for transitions between attractors at the same level of tree (level is the number of transitions from the root). In [5] a method was introduced that applied successively higher MFPT thresholds to prune edges from this complete MFPT tree as a means to identify separation among subsets of close attractor states as illustrated in Figure 8(b). The effects of changing the threshold from low to high was proposed as a possible mechanism for cellular differentiation with the low threshold representing pluripotency and the process of raising the threshold as type specialization as attractors become more and more isolated. This model proposes that cells differentiate by actively controlling their sensitivity of expression noise and can account for the observation that terminally differentiated cell states tend to be more stable than pluripotent states.

Network search

We perform a uniform Monte Carlo search over the space of critical random Boolean networks. For each network we find the attractors and compute the MFPT between all possible attractor pairs (extended from code posted at http://code.google.com/p/pbn-matlab-toolbox[9]). Using the MFPT values, for each type of multistable switch added to the network, we draw a density plot where the x-axis is the forward MFPT and the y-axis is the reverse MFPT (we consider the edge with lower MFPT as forward). The acquired density plots are used to determine the distribution of directional, non-directional, separated and non-separated probability transitions between attractors in each network.

Network types

We investigated nine types of networks. Approximately 5 ∗ 10⁴ networks of each type were explored to find 5, 000 networks of each type with five or more attractors. The different network types come from the use of the 4 multistable switches that are shown in Figure 4. The first switch (Figure 4 a) is a bistable switch, a small local circuit with feedback loops. This is a common switch in biological networks and it controls binary branch points between two mutually exclusive cell lineages[1, 10]. The truth table of the functions in this switch is [1,1,0,1] for binary numbers [00,01,10,11] respectively. The other three switches all encode mutual inhibition between two genes. The first is MI00 and is based on the network synthesized in [10]. MI00 includes two incoherent feedback loops. The final two switches extend mutual inhibition with the addition of one MI+0 or two MI++ positive (coherent) feedback loops. These two switches were explored in [22], where it was shown that the positive feedback loops can introduce additional shallow attractor basins in continuous ODE network models.

The four switches were used as described above to construct nine different types of networks: no motif, one BS, one MI00, one MI0+, one MI++ and then four more network classes each with two of the same switch. Note that when a motif defined in Figure 4 is embedded, two nodes of the original RBN are selected randomly, their logic functions replaced and inputs and outputs rewired. For illustration, consider how a MI0+ motif is embedded into a RBN. Starting with a RBN (see Figure 4 (a)), two nodes are selected randomly and their truth tables are changed to [0,1,0,0]. Then, the 2⁰ input of the second node is wired to the output of first node and, conversely, the 2⁰ input of first node is wired to the output of the second node. The small or-gate and not-gate are not considered in wiring, because they were previously considered in the truth tables of their respective nodes.

Mean first passage time

where d _ij is equal to 1 if there is a deterministic transition from x to y in the network dynamics, otherwise it is 0; p _e is the probability of a single bit flip resulting from noise and h _xy is the Hamming distance between two states; n is the total number of nodes in the network.

A low MFPT between two states indicates that starting from the first state, the second state is easily reached by gene expression noise. Figure 9 shows F _k , kF _k , and MFPT for the transition between two arbitrary attractors. As this figure shows, the a to b transition has a lower MFPT compared to the other.

At each network state update D(ŝ) there is a probability that the state will change as a function of the Hamming distance (h) between the current state and the subsequent state ŝ_t+1← D(η(ŝ _t , r)). MFPT models uniform gene expression noise by considering probabilistic bit flips at every possible state of the network and deriving the distribution of passage times from analysis of the corresponding Markov process. Statistically, the probability distribution of bit flips can be seen as a binomial distribution, thus the probability of r bit flips, η(ŝ _a , r) is $\left(\begin{array}{c}\hfill h\hfill \\ \hfill r\hfill \end{array}\right)p^r{\left(1-p\right)}^{h-r}$, where p is the probability of a single bit flip and h is the total number of bits.

Evaluating epigenetic barriers: MFPT vs. transitory bit flips

Villani et al. [5] studied noise-driven network transitions in RBNs. They introduced a measure of the probability of network transitions as the likelihood of attractor transition under expression noise. In this measure, for each pair of attractors {a _i , a _j }, P(i, j) is the portion of single one-step bit flips (transitory perturbations) in the nodes of all states of attractor a _i which will result in a transition from a _i to a _j under noise-free dynamics. The measure of likelihood of network transition under noise is similar to MFPT, but it does not consider gene expression variability throughout the network. MFPT better models global expression noise by considering probabilistic bit flips at every possible state of the network and deriving the distribution of passage times from analysis of the corresponding Markov process.

Since one-off bit flips consider noise only as a single bit changes and only when the network has reached its attractor states, it could serve as an efficient yet heuristic measure of the MFPT. To test this idea, a study was performed on a set of small critical networks where for each network and each pair of attractors, P(i, j) was compared with MFPT(i, j). Figure 10 depicts the relationship between MFPT and P for 100 arbitrary Boolean networks that have 5 or more attractors. Each point represents the epigenetic barrier between two attractors measured in MFPT and P. Since the networks studied in these experiments are small and do not have many attractors, many points are located in the line P = 0. The regression line in this figure shows that as MFPT increases P tends to decrease. P and MFPT are modestly correlated for these small networks and it is unclear how well one-off bit flips can accurately estimate MFPT when network size grows. Since the networks in our experiments are small, we only consider MFPT because of its realism in modeling expression noise.

Evaluating epigenetic barriers: MFPT vs. Hamming distance

An intuitive idea is that MFPT between attractors has a direct relationship to the Hamming distance that separates these attractors. However, we found that this is not the case. Instead, network dynamics, not the Hamming distance, is the main contributor to the MFPT between attractors. As an example of the limitations of Hamming distance, consider that the MFPT(a _i , a _j ) and MFPT(a _j , a _i ) can be different, but that the Hamming distance between these attractors is the same. However, even though there is not a strong relationship between MFPT and Hamming distance, a weak correlation between the average of the forward and reverse MFPT between attractors and their Hamming distance can be detected. This is depicted in Figure 11 which shows MFPT versus the Hamming distance obtained from 100 RBNs containing 8 nodes. As the Hamming distance increases, the upper-bound of MFPT values also increases (r = 0.1027 for Hamming distance and average MFPT). In Figure 11, the box represents the central 50% of the points and the red bar shows the median of data.