Sequential Branch Points in Hierarchical Cell Fate Determination

Sequential Branch Points in Hierarchical Cell Fate Determination David Foster, Sui Huang, Stuart Kauffman arXiv:0903.4215v1 [q-bio.MN] 25 Mar 2009 Institute for Biocomplexity and Bioinformatics, University of Calgary, Biological Sciences Building, 2500 University Dr. N.W., Calgary, Alberta, Canada T2N 1N4 Abstract Multi-potent stem or progenitor cells undergo a sequential series of binary fate decisions which ultimately generates the diversity of differentiated cell. Efforts to understand cell fate control have focused on simple gene regulatory circuits that predict the presence of multiple stable states, bifurcations and switch-like transition. However, existing gene network models do not explain more complex properties of cell fate dynamics such as the hierarchical branching of developmental paths. Here, we construct a generic minimal model of the genetic regulatory network which naturally exhibits five elementary characteristics of cell differentiation stability, branching, exclusivity, directionality, and promiscuous expression. We argue that a modular architecture comprised of repeated network elements in principle reproduces these features of differentiation by repressing irrelevant modules and hence, channeling the dynamics in the high-dimensional phase space through a simpler, lower dimensional subspace. We implement our model both with ordinary differential equations (ODEs) to explore the role of bifurcations in producing the oneway character of differentiation and with stochastic differential equations (SDEs) to demonstrate the effect of noise on the system in simulations. We further argue that binary cell fate decisions are prevalent in cell differentiation due to general features of dynamical systems. This minimal model makes testable predictions about the structural basis for directional, discrete and diversifying cell phenotype development and thus, can guide the evaluation of real gene regulatory networks that govern differentiation. Key words: cell differentiation, genetic regulatory network, dynamics, modularity 1. Introduction The genetic regulatory network plays an essential role in controlling cell differentiation and maintaining the stability of cell types. Each cell type has its own stable gene expression profile, the set of genes in its genome which are actively expressed [1]. Thus, understanding cell differentiation requires a model of the behavior of many interacting genes which establish the expression profiles associated with each cell type. Gene regulatory interactions are often described as a network, with each node representing a gene and links between two genes indicating a regulatory interaction. The large scale characteristics of genetic regulatory networks have been studied in some detail [2, 3], but only in terms of network structure: the degree sequence [4], the prevalence of network motifs [5, 6], and so forth. In terms of GRN dynamics, most efforts have been “bottom-up” - devoted to simulating specific regulatory interactions or analyzing models of small genetic circuits [7, 8] and therefore focused on specific subsections of the network. Our goal is to work in a “top-down” fashion to construct a Preprint submitted to Journal of Theoretical Biology September 24, 2009 minimal dynamical model which captures many important features of cell differentiation. Most significantly, we capture the sequential branching differentiation of cells from the multipotent progenitor cell to more lineage-restricted cells and finally to terminal cells. To construct our model, we first note five elementary characteristics of cell differentiation which we wish to capture: stability, directionality, branching, exclusivity, and promiscuous expression. Of course, this is not meant to be an all encompassing or fully detailed description of cell differentiation. However, these simplified characteristics prove rich enough to create an interesting model. We review here our definition of each characteristic. Stability : Cell types should be stable states of the dynamics. A cell which is in specific cell type should remain in that cell type for a “long time” (relative to the time scale of gene activation, transcription and translation), unless altered by an external signal. Directionality: Progenitor cells should reliably differentiate into a variety of more specialized cell types [9] and should only very rarely to regain pluripotency by de-differentiating [10]. Branching: Differentiation should occur through discrete lineage specifications, in which a multipotent cell makes a binary fate decision and commits to one of two different lineages. Exclusivity: Once cells in our model have committed to a specific lineage, they should lose the ability to produce cell types from the other lineage. Each differentiation event both selects a specific lineage and forecloses the possibility of the other. Promiscuous expression: Multipotent progenitor cells should simultaneously express genes specific to all of their prospective lineages. For example, a hematopoitic stem cell expresses genes found in all of the lineages into which it can ultimately differentiate [11]. 2. Model specification To create a model which exhibits these characteristics, we take our inspiration from the purported network properties of real GRN. Biological networks exhibit modularity [12, 13, 14]. Modular networks can be divided into largely independent subnetworks, called modules. The nodes in a given module interact strongly with other nodes in the module, but only weakly with nodes in other modules. In our model, the network is composed of many modules, each responsible for a specific cell fate decision. These modules influence each other in an hierarchically organized fashion, as shown in Fig. 1, which ensures that differentiation events occur in a branching, sequential fashion. The dynamical state associated with each cell type is a product of the dynamics of the associated module. Arrows pointing from one module to another indicate influence of the former on the latter. Influence between modules acts as follows: a parameter in the regulated module depends on the expression level of a gene in the regulating module. In the top pairs of panels the system is in the “Multi” (multi-potent) state, and in the bottom pair, the system is in one of the two available progenitor cell types. On the dynamical side, the module A undergoes a change in its state (the specifics of which we will discuss shortly). In making the change, it exerts an influence on the two modules downstream, in this case turning off module C and allowing module B to enter a state which describes “Pro 1” (Progenitor 1). Notice that like module A, modules B and C are also linked to down-stream modules on which they exert influence, and that module A is itself influenced by a module “higher up” on the differentiation cascade. This allows the process to proceed sequentially through as many levels as are required. Our model suggests an important role for modularity: restricting the available flows in the state space of the dynamics. These restrictions ensure that only the appropriate pathways of differentiation are available to a given cell type. The modular structure of the network is what imposes the branching structure on cell differentiations. 2 Network Cell type x1 A x2 Multi B x3 C x4 x5 x1 A x2 x4 Pro 2 Multi B x3 Pro 1 x6 C x5 x6 Pro 1 Figure 1: Transition of a module in our network from the multipotent cell type to one of the progenitor cell types. Modules consist of two genes, each of which activates its own transcription and inhibits the transcription of the other gene. In the top panels, both genes in module A are active at a moderate level, represented by the gray coloration. This intermediate steady state is poised between the two subordinate lineages. Modules B and C are also both expressing their genes at a moderate level, as a result of the influence of module A. In the lower panels, the cell has differentiated from the multipotent cell type to progenitor type 1. Module A has switched from the intermediary stable state to the stable state in which gene x1 is highly active, whereas gene x2 has become inactive. As a result, module C has become completely inactive due to the signals from module A, whereas module B remains in the intermediate steady state, and is thus prepared to make further cell fate decisions. 3 3. Dynamics of a single module Previous work has studied the behavior of GRNs using a variety of different dynamical system models: Boolean networks [15, 16, 17], differential equations [18, 19], and chemical master equations [20, 21]. Although each model has its merits, for our current purposes differential equations provide the best combination of rich behavior and analytical tractability. To model a set of N interacting genes as differential equations, we associate the expression level (number of proteins produced by the gene which are present in the cell) of each gene with variable xi . The state of the system is given by the vector x = [x1 . . . xN ]. Since genes influence the transcription rates of other genes, the rate of transcription (or equally the instantaneous rate of change of the expression level) of each gene, ẋi , is given by: ẋi = Fi (x) (1) If we plot the possible expression levels of each gene as one axis in a graph, the state space of the dynamics is given by {R+ }N . For a given state, the growth or decay of the expression level of each gene defines a vector ẋ = [ ẋ1 . . . ẋN ], which points in the direction that the state of the system will move in state space under the dynamics of the system. This is the trajectory of the system. A stable fixed point is any point in state space in which the state of the system is “trapped” (under the deterministic dynamics, ignoring noise). Formally, this occurs for any point x⋆ such that ẋi⋆ = Fi (x⋆ ) = 0 for all i. In our model, each cell type is represented by a stable fixed point of the dynamics. Our goal is to construct the simplest plausible set of genetic regulatory interactions for which the state of the system moves, either by random fluctuations or external signals, through a sequence of branching, exclusive steps from the fixed point which represents a multipotent cell type, to states which represent a series of intermediary cell types, and finally to a terminal cell type. Huang et al [22] provide experimental evidence to show that the genetic circuit controlling hematopoetic (blood cell) stem cells consists of a pair of mutually repressing genes, each of which also acts as a self-activator. Such a feedback loop is alleged to be an over-represented motif in real genetic networks [4]. We have chosen the Huang switch for our cell-fate controlling module, as shown in Fig. 1. The differential equations describing this system are: b1 ψn1 a1 x1n + − k1 x1 x1n + θ1n x2n + ψn1 n n a2 x b2 ψ x˙2 = n 2 n + n 2 n − k2 x2 x2 + θ2 x1 + ψ2 x˙1 = (2) These equations model the dynamics of self activation and mutual inhibition with Hill functions. The parameters a1 , a2 give the maximum transcription rate due to self activation, and θ1 , θ2 describe the concentration of protein x1 or x2 at which the inflection point in the Hill function occurs. Parameters b1 , b2 give the basal rate of transcription, which is suppressed by a downward sloping Hill function with ψ1 , ψ2 controlling the concentration of protein at which its inflection point occurs. The cooperativity of the interactions is given by n, which controls the steepness of the Hill function. The rate of protein degradation is controlled by k1 , k2 . To simplify our analysis we impose certain symmetries on the equations, and select specific values for some parameters which produce the desired behavior. We make the maximum contributions of the basal expression rate and the self-activation the same (a1 = b1 ) and thus have 4 a single parameter which controls the overall maximum transcription rate of the gene. For simplicity, we also set the cross-over points of the two Hill functions (for the self-activation and inhibition) to the same value (θ1 = ψ1 = θ). Finally, we choose to make all parameter values a, θ, and k identical between the two equations (so k1 = k2 ), thus making the behavior symmetric between x1 and x2 . Thus we can write the new set of parameter values as: a1 , a2 , b1 , b2 = a θ1 , θ2 , ψ1 , ψ2 = θ = 1 n=4 k1 , k2 = k (3) We choose specific parameter values θ = 1 and n = 4. These parameter choices are not unique, and do not require careful tuning: as long as n ≥ 4 a wide range of values of θ, a, and k produce the desired tri-stable behavior. This leaves us with the new simplified equations:   4  x1 1   − kx1 + 4 x˙1 = a  4  x1 +4 1 x2 + 1   x 1   − kx2 x˙2 = a  4 2 + 4 x2 + 1 x1 + 1 (4) These equations have only two free parameters: a and k. These are the parameters which are influenced by other modules, following the schematic shown in Fig 1, and as detailed below. In order to better understand the influence of these parameters, we explore the bifurcation which occurs as they are change. A bifurcation is a dynamical event in which a small change of a parameter leads to a qualitative change in the dynamics of the system. Bifurcations can destabilize fixed points, or cause points which were previously unstable to become stable. While a fixed point may be stable under the deterministic dynamics of the system, in the presence of noise even a stable fixed point may not trap the system indefinitely. Thus bifurcations which change the stability of a fixed point are the dynamical equivalents of cell fate decisions in differentiating cells. For example, a parameter change may destabilize the fixed point representing a multipotent cell type, thus driving the state to a different fixed point, which represents the new progenitor cell type. For a = 2 and k = 1, the two variable system is “tri-stable”, meaning that there are three stable fixed points in the dynamcs. This state space is shown in the top left panel of Fig 2. The three stable fixed points are given by black dots: there is an intermediate fixed point, characterized by moderate amounts of x1 and x2 , and two exterior fixed points, characterized by large amounts of one variable, and small amounts of the other. Leaving the parameter a = 2 fixed and tuning k from 1 to 1.5 drives a bifurcation from tri-stablity to bi-stabliltiy, as shown in the subsequent three panels of Fig. 2. The intermediate fixed point is destabilized, leaving only the two exterior stable fixed points. The bifurcation also occurs if the overall transcription rate, a, is increased while k is kept constant –only the ratio of these parameters have a real effect on the qualitative dynamics. Their absolute values merely set the time scale. In their analysis of the tri-stable switch, Huang et al [22] argue that the “intermediate” fixed point represents the hematopoitic stem cell, and the two “exterior” fixed points represent the two committed lineages. They provide evidence that the differentiation from hematopoesis to the erythroid or myelomonocytic lineages is caused by the destabilization of the intermediate fixed point via a subccritical pitchfork bifurcation. This in turn forces the cell to shift to one of the 5 k=1 k = 1.166 k = 1.333 k = 1.5 [x2 ] [x2 ] [x1 ] [x1 ] Figure 2: Each panel shows the state space of the two gene toggle switch for the parameter values chosen (a = 2 and k = 1, 1.166, 1.333, 1.5). For k = 1.0 there are three stable fixed points in the dynamics (shown as black dots) and two unstable fixed points (shown as gray diamonds). The arrows give the trajectory of the system at each point in phase space, with longer arrows indicating faster movement and shorter arrows slower movement. As k changes from 1 to 1.5 the module undergoes a subcritical pitchfork bifurcation [23] (shown in more detail in Fig. 3). The two interior unstable fixed points converge on the central stable fixed point eventually destabilizing it when k ≈ 1.3, leaving only the two exterior stable fixed points. 6 [x1 ] a. b. c. k Figure 3: The x1 position of the fixed points as the parameter k changes from 1.0 to 1.5. The black circles represent stable fixed points while the gray diamonds represent unstable fixed points. As the module undergoes a subcritical pitchfork bifurcation, a cell at the intermediate stable fixed point moves along the trajectory shown by the arrow a). When k ≈ 1.3 the fixed point is destabilized and the cell is pushed by any random fluctuations from the now unstable fixed point to one of the exterior fixed points following either trajectory b) or the dotted trajectory. Once the cell has been drawn to one of the exterior fixed points, if the module reverses the bifurcation and once again becomes tri-stable, the cell will remain at its new stable fixed point, following trajectory c). So once the state of the system has undergone the bifurcation, it cannot easily be driven back to its starting position, thus making the differentiation event modeled by this bifurcation one-way. 7 two other fixed points. In Fig. 3 we show the bifurcation diagram corresponding to this event. The black dots indicate the stable fixed points, and the grey, the unstable fixed points. The x-axis shows the value of the bifurcation parameter, k. When k = 1.0 the circuit is tri-stabile: there are three stable fixed points. When k = 1.3, the intermediate stable fixed point is destabilized by the approach of the two unstable fixed points. Thus, the state of the system must migrate to one of the two remaining stable points. However, starting from one of the exterior fixed points, if we reverse the process and decrease the bifurcation parameter, the system will not return to the intermediate steady state–the exterior fixed points remain stable to small perturbations throughout. Thus, to reverse the cell differentiation, not only would the the intermediate fixed point need to be restabilized, the state of the system would also need to be driven from the exterior fixed point. This ensures that cell differentiation is directional: a simple, transient, one parameter perturbation can drive the system toward specificity, but there is no similarly simple way to de-differentiate back to multipotency. Noise is ubiquitous in biological systems and important for many functions [24]. The GRN as a whole, but especially the genes responsible for the process of development (and cell differentiation), must be organized for reliable performance. Thus we require that our model behaves robustly in the presence of noise. Feedback loops [25] and modularity [26] have been shown to contribute to robustness, and because our model makes use of these design characteristics, we find that robustness is easily attained. We simulate noise in the dynamics using a Langevin equation: ẋi = Fi (x) + gi (x)ηi (t) (5) Here, Fi (x) is the same function as in the deterministic dynamics. In the second term, an arbitrary function gm i (x) describes the influence of the random function ηi (t) on the variables. We choose to introduce only additive noise: gi (x) is simply a constant value for all i, with no dependence on x, as suggested elsewhere [27]. We can then simulate the evolution of the system forward by one time step of size ∆t with the equations: √ xi (t + ∆t) = xi (t) + ∆t ∗ ẋi (t) + ∆t ∗ g ∗ σi (6) The parameter g sets the amount of noise in the dynamics; we have chosen to make it the same for every gene. The parameter σi drawn at each time step from a gaussian distribution with mean 0 and standard deviation 1. The influence of the random variable is similar to a random walk. Since the variables in our system describe the concentration of protein products, we also require that no variable be allowed to drop below zero, since a negative concentration is non-physical. In the presence of noise, we find that the transition between stability and instability for the intermediate fixed point is no longer abrupt. As the module approaches the bifurcation point, the restoring force which pushes the state of the system back towards the intermediate fixed point weakens. We simulated 500 cells starting at each of the three stable fixed points for values of k = 1, 1.1, 1.2, 1.3 with g = .1 and ∆t = .001 for 40, 000 time steps (40 units of t, a fairly lengthy period in simulation time relative to the timescale of transcription and translation). The distribution of the final states for cells starting at the intermediate fixed point spreads as k increased (shown in Fig. 4). As the value of k increases, the variance of the gene expression increases. This behavior may soon be observable in vivo, using flow cytometry to measure the approximate expression level of a gene for each cell in a population of cells approaching differentiation. By observing the variance of genes as cells differentiate, we may gain a clearer picture of which genes are controlling cell fate decisions: the expression level of significant genes should become more varied across the population as the cells approach a cell fate decision. 8 k=1 k = 1.1 k = 1.2 k = 1.3 [x2 ] [x2 ] [x1 ] [x1 ] Figure 4: The behavior of the module in the presence of noise as it approaches the bifurcation. Each panel shows the final positions of 500 simulations (“cells”) starting from the central fixed point simulated for 40 units of t (in ∆t = .001 time steps). The noise magnitude is g = .1 and the overall activity a = 2.0 for various values of k. Black circles show the positions of the stable fixed points and gray squares show the end points of each simulation. As k increases, the variance of expression levels of x1 and x2 increases near the intermediate fixed point. For k = 1.0, 1.1 none of the simulated cells “escape” from the intermediate fixed point. For k = 1.2 approximately 50 cells (or 10%) escape, divided evenly between the two exterior fixed points. When k = 1.3 approximately 90% escape. In similar simulations at the same noise strenght, g = .1, if the “cells” start near one of the exterior fixed points, none of them transitioned to the intermediate fixed point for any parameter values tested. 9 Before the actual bifurcation point is reached, we see non-trivial numbers of cells have already moved to one of the exterior, ‘differentiated’ fixed points. This indicates that for values of k ≤ 1.2 the intermediate fixed point becomes metastable (weakly stable). Thus in the presence of noise, the state of the system will not remain near this fixed point indefinitely. As we discuss later, progenitor cell types may be metastable fixed points in which the state of the network rests for a short while before continuing its differentiation. This would explain the difficulty of culturing and maintaining such cell types in vitro. Finally, our model provides a potential dynamical explanation of the observed prevalence of binary (as opposed to trinary, etc) branches in cell differentiation. The sub-critical pitchfork bifurcation found in our modules arises from the topology of the regulatory interactions and is present for many choices of parameter values without careful tuning. This bifurcation leads naturally to binary cell differentiations. In Fig. 4 the distribution of final states does not spread evenly throughout the state space: it spreads preferentially along a pair of trajectories which correspond to the two sides of the unstable eigenvector of the intermediate fixed point. These trajectories, in turn, lead to one of the two exterior fixed points. In general, only one of the eigenvalues of the fixed point becomes positive at the bifurcation point. Thus, near the interior fixed point, flows along the other eigenvectors still push the state of the system toward the interior fixed point, concentrating flows along the unstable eigenvector. Trinary branchings, since they cannot depend on such simple mechanics of the state space, are comparatively much more fragile and require much more careful selection of parameters. We have found in our model that binary branching occurs ‘for free’ from a robustly occurring pitchfork bifurcation, with the simulated cells flowing away from their starting state along one of the two ends of the single, destabilized eigenvector. 4. Sequential Differentiation Now we are ready to build the full hierarchical model. In Fig. 5, we show a differentiation event in more detail: three modules are shown, A, B, and C. Module A influences B and C. In the top panel, module A is in the intermediate fixed point, but this fixed point is only meta-stable. This meta-stability is caused by the influence of the signal ∆A . If module A is part of a larger cascade, then the signal will originate from a gene in the module above A in the hierarchy; if A is at the top of the hierarchy then the signal may originate from outside the cell. Module A expresses moderate levels of x1 and x2 , which in turn influence the parameters kB , aB and kC , aC in modules B and C, through the functional relationships: kB = G(x1 ), kC = G(x2 ) aB = H(x1 ), aC = H(x2 ) (7) For simplicity, the same functions G and H are used throughout the network. In Fig. 6 we show how the piece-wise linear functions G and H respond to changes in the variable on which they depend. The specific functional form chosen is not important; sequential differentiation will occur as long as the functions reflect two basic requirements. First, if the value of x decreases dramatically, indicating that the gene x has been shut off, then transcription rate, a, of the downstream module, which is under the influence of x, should be dramatically decreased. To achieve this, for small values of x, H(x) should be small, effectively silencing the genes in the influenced module. Second, if the value of x increases substantially, then the fixed point in the influenced module should be destabilized, preparing it to undergo differentiation. Thus, the function G(x), 10 x1 kA = G(∆A ) = 1.25 x1 Module A kB = G(x1 ) = 1.0 x3 aB = H(x1 ) = 2.0 intermediate x2 kC = G(x2 ) = 1.0 x5 aC = H(x2 ) = 2.0 x4 x3 Module B x2 meta-stable intermediate x4 x1 kA = G(∆A ) = 1.25 Module B x4 x3 x6 x2 x2 x1 Module A meta-stable x1 differentiated kB = G(x1 ) = 1.25 x3 aB = H(x1 ) = 2.0 x5 Module C Noise drives differentiation: x6 x2 kC = G(x2 ) = 1.0 x5 aC = H(x2 ) = 0.2 inactive x4 Module C x6 x5 x6 Figure 5: The modular network and consequences of a single cell fate decision, in detail. In the top panel, the state of the system is poised in the meta-stable intermediate stable state of module A. Modules B and C are both in the intermediate stable state, expressing x3 , x4 , x5 , and x6 at a moderate level. In the lower panel, module A has undergone differentiation, driven by noise, switching to one of the two exterior stable states, x1 high, x2 low. This in turn has caused module B to become meta-stable by increasing the value of kB from 1.0 to 1.25. Module C is inactivated by the change in parameter aC from 2.0 to 0.2, essentially silencing the genes in the module. The k and a parameters in all three modules are functionally dependent on one of the genes in the module directly above them. In the case of module A, we denote the stimulus as a generic external signal ∆A to reflect the fact that it could come from another gene, or some external differentiation signal. 11 G(x) H(x) x Figure 6: The piecewise linear response functions G and H. In the top panel we show the function G, which controls the decay rate k, and in the bottom panel we show the function H, which controls the transcription rate a. Specifically, G = 1.0 for x ≤ 2.7 and G = 1.25 for x ≥ 3.0 and interpolates between these values for 2.7 ≤ x ≤ 3.0. Similarly H = 0.2 for x ≤ 1.0 and G = 2.0 for x ≥ 1.5 and interpolates between these values for 1.0 ≤ x ≤ 1.5. which controls k, should increase slightly in value in response to higher values of x. This higher value of G(x) will drive the bifurcation above. Returning to Fig. 5, moderate levels of x1 and x2 result in parameters in modules B and C which maintain stable, moderate expression levels of their genes x3 , x4 , x5 and x6 . Note that any modules downstream of B and C would also express their genes at an intermediate level. Since module A is in merely a meta-stable state, the state of the system is driven by noise after some time to one of the two exterior fixed points, in this case x1 high, x2 low. The result is shown in the lower panel: module C has been silenced by the decrease in expression of x2 (thus also shutting off all the genes which are in modules downstream from it, under the control of x5 or x6 ). Module B, under the influence of increased x1 , has become meta-stable, and is thus waiting to drive a differentiation event as the system “chooses” between expressing x3 and x4 . The suggested links that integrate the modules correspond to the simplest possible regulatory connections that generate the observed branching behavior. Whether real networks utilize similar connections remains to be seen, however, modules of mutually regulating transcription factors clearly are linked in real GRN. 12 Module A x1 x2 Meta-stable state Decided Module B x3 x4 Meta-stable state Decided Module C x5 x6 Inactivated t Figure 7: A simulation of three interacting modules with noise g = 0.1. At t = 5 the intermediate fixed point of module A is changed from stable to meta-stable by an external stimulus. For the next ≈ 10 units of time t the top module remains meta-stable, then finally falls into the x1 high exterior fixed point, as shown in the top panel. This silences module C (and thus genes x5 and x6 ) which are rendered inactive, as shown by the bottom panel. Finally, in the middle panel we see module B, which became meta-stable after module A made its cell fate decision and entered the exterior fixed point. Module B remains meta-stable for ≈ 25 units of t until it transitions to an exterior fixed point, completing the differentiation of the cell through this segment of the network. 13 5. Whole network results Fig. 7 shows the result of a simple three module simulation, with the noise parameter g = 0.1. The differentiation tree functions properly in the presence of noise: at the outset, all genes are expressed at a moderate level. At time t = 5, a signal is introduced which destabilizes the top module. After some time, the state of the the top module is decided in favor of x1 . As x2 decreases, we see that genes x5 and x6 are also silenced. The genes x3 and x4 continue to be expressed at moderate levels for some time, until through random fluctuations, the system exits the meta-stable state, resulting in increased x4 expression and decreased x3 expression. At the outset, when any outcome was possible and the cell multipotent, all of the genes were expressed. As the cell differentiates and loses its pluripotency, the genes which correspond to the lineages which are no longer accessible to the cell are silenced. As seen in the simulation described above, individual cells can remain in the meta-stable state for a long time before continuing their differentiation. We simulated 2000 iterations of the model, and measured the length of time the system spent in meta-stable states. We found an exponential distribution of waiting times, with some of our simulated cells spending as many as 290 units of t in the meta-stable state. For the purposes of comparison, we note that once the state escapes the region of state space close to the meta-stable fixed point and the transition to one exterior fixed points begins, the transition often occurs in less than 5 units of t. So while individual cells cannot be maintained in these meta-stable states indefinitely, they can be quite long lasting. We have not considered the proliferation of cells, as described by Hoffmann et al [27]. The proliferation rate is a function of the state of the GRN. Thus if the proliferation rate is high in the meta-stable and intermediate states, a population of cells could be maintained in the meta-stable state indefinitely. The number of progenitor cells found in vivo is the result of a balance between cells differentiating into the meta-stable state from more multipotent cell types, progenitor cells proliferating, and progenitor cells differentiating out of the progenitor cell type to a more lineage committed state through the stochastic dynamics described above. The metastability of the dynamical fixed point may provide an explanation for the limited capacity for self renewal of some progenitor cell types in vitro. 6. Conclusion Our model reproduces all of the desired behaviors of cell differentiation by using established network properties of real genetic nets. Cell differentiation occurs in a branching, exclusive fashion between a dynamically distinct cell types. A cell type higher in the differentiation tree exhibits promiscuous expression of all its accessible lineages. Cell fate decisions are directional: small parameter changes can easily drive differentiation down the tree, but not back up again. It is easy to add additional control mechanisms to this model. For example, we could add a necessity for external signals, in concert with the influence of “higher” modules, to drive each bifurcation. Our goal however is merely to create the minimal outlines of a dynamical system capable of capturing the general observed behaviors of cell differentiation. In addition, our model suggests a new role for modularity and an explanation for the prevalence of binary cell fate decisions. Modules allow for the easy implementation of repeated behaviors, in this case facilitating the repeated branchings observed in cell fate decisions. A hierarchical, modular network topology serves to channel trajectories in the state space through lower dimensional sub-spaces. This channels the state of the system through an orderly progression of fixed 14 points, representing the various cell types. Many of the important properties of cell differentiation can be captured quite naturally by a system of loosely coupled genetic modules. Further we found that binary cell fate decisions emerge effortlessly and robustly from the topology of our network modules due to the behavior of flows near the intermediate fixed point. The dynamical simplicity and robustness of binary cell fate decisions could explain their prevalence in cell differentiation. Building progressively more accurate whole-system models that recapitulate key functionalities provides an important complement to the traditional “bottom up” approach, which seeks maximal integration of all available details. A “top down” approach, focused on less detailed minimal models, allows us to focus on understanding the influence of the topology of the genetic regulatory network on the resulting dynamics. The details of our individual modules and their interconnections are of course over simplified and incomplete. However, emerging evidence from the topology of real genetic networks suggests that the overall modular architecture is likely to be correct. We hope that future experimental results will show to what degree the genetic regulatory network in general, and the cell differentiation network in particular can be treated as a collection of loosely coupled modules. As large-scale models improve through interaction with experimental evidence, we believe that the top down and bottom up approaches can eventually “meet in the middle”, fitting carefully specified subsections of the genetic network into a larger network framework. 7. Acknowledgements The authors would like to thank Jacob Foster for useful comments and discussions. References [1] Kauffman, S.A, 1969. Metabolic stability and epigenesis in randomly constructed genetic nets. J. Theor. Biol. 22, 437-467. [2] Davidson, E., et al., 2002. A Genetic Regulatory Network for Development. Science 295, 1669-1678. [3] Alon, U., 2009. Biological Networks: The Tinkerer as an Engineer. Science 301, 1866-1867. [4] Lee, T.I., et al., 2002. Transcriptional Regulatory Networks in Saccharomyces cerevisiae. Science 298, 799-804. [5] Shen-Orr, S., Milo, R., Mangan, S., Alon, U., 2002. Network motifs in the transcriptional regulation network of Escherichia coli. Nature Genetics 31, 64-68. [6] Milo, R., et al., 2002. Network Motifs: Simple Building Blocks of Complex Networks. Science 298, 824-827. [7] Chickarmane, V., et al., 2006. Transcriptional Dynamics of the Embryonic Stem Cell Switch. PLoS Comput. Biol. 2(9) e123. [8] Ackers, G.K., Johnson, A.D., Shea, M.A., 1982. Quantitative Model for Gene Regulation by lambda Phage Repressor. Proc. Natl. Acad. Sci. 79, 1129-1133. [9] Huang, S., Eichler, G., Bar-Yam, Y., Ingber, D.E., 2005. Cell Fates as High-Dimensional Attractor States of a Complex Gene Regulatory Network. Phys. Rev. Lett. 94, 128701. [10] Kal, T., Spradling, A., 2004. Differentiating germ cells can revert into functional stem cells in Drosophila melanogaster ovaries. Nature 428, 564-569. [11] Cross, M., Enver, T., 1997. The lineage commitment of haemopoietic progenitor cells. Current Opinion in Genetics & Development 7, 609-613. [12] Hartwell, L.H., Hopfield, J.J., Leibler, S., Murray, A.W., 1999. From molecular to modular cell biology. Nature 402, C47-52. [13] Ihmels, J., et al., 2002. Revealing modular organization in the yeast transcriptional network. Nature Genetics 31, 370-377. [14] Ravasz, E., Somera, A.L., Mongru, D.A., Oltvai, Z.N., Barabasi, A.L., 2002. Hierarchical Organization of Modularity in Metabolic Networks. Science 297, 1551-1555. [15] Kauffman, S.A., 2004. The ensemble approach to understand genetic regulatory networks. Physica A 340, 733-740. 15 [16] Kauffman, S.A., Peterson, C., Samuelsson, B., Troein, C., 2003. Random Boolean network models and the yeast transcriptional network. Proc. Natl. Acad. Sci. USA 100, 14796-14799. [17] Espinosa-Soto, C., Padilla-Longoria, P., Alvarez-Buylla, E.R., 2004. A Gene Regulatory Network Model for CellFate Determination during Arabidopsis thaliana Flower Development That Is Robust and Recovers Experimental Gene Expression Profiles. The Plant Cel 16, 2923-2939. [18] Gardner, T.S., Cantor, C.R., Collins, J.J., 2000. Construction of a genetic toggle switch in Escherichia coli. Nature 403, 339-342. [19] Chickarmane, V., Troein, C., Nuber, U.A., Sauro, H.M., Peterson, C., 2006. Transcriptional Dynamics of the Embryonic Stem Cell Switch. PLoS Comput. Biol. 2(9). [20] Ribeiro, A., Zhu, R., Kauffman, S.A., 2006. A General Modeling Strategy for Gene Regulatory Networks with Stochastic Dynamics. Journal of Computational Biology 13(9), 1630-1639. [21] Danos, V., Feret, J., Fontana, W., Harmer, R., Krivine, J., 2008. Rule-based modelling symmetries, refinements. Proceedings of FMSB 2008, Lecture Notes in Bio Informatics, to appear. [22] Huang, S., Guo, Y., May, G., Enver, T., 2007. Bifurcation dynamics in lineage-commitment in bipotent progenitor cells. Developmental Biology 305, 695-713. [23] Strogatz, S., 1994. Nonlinear Dynamics and Chaos. Perseus Publishing, Cambridge, 55-60. [24] Elowitz, M.B., 2002. Stochastic Gene Expression in a Single Cell. Science 297, 1183-1186. [25] Eldar, A., et al., 2002. Robustness of the BMP morphogen gradient in Drosophila embryonic patterning. Nature 419, 304-308. [26] Variano, E.A., McCoy, J.H., 2004. Networks, Dynamics, and Modularity. Phys. Rev. Lett. 92, 188701. [27] Hoffmann M., Chang H.H., Huang S., Ingber D.E., Loeffler M., et al., 2008. Noise-Driven Stem Cell and Progenitor Population Dynamics. PLoS One 3(8):e2922. doi:10.1371/journal.pone.0002922 16