Towards an understanding of lineage specification in hematopoietic stem cells: A mathematical model for the interaction of transcription factors GATA-1 and PU.1

In addition to their self-renewal capabilities, hematopoietic stem cells guarantee the continuous supply of fully differentiated, functional cells of various types in the peripheral blood. The process which controls differentiation into the different lineages of the hematopoietic system (erythroid, myeloid, lymphoid) is referred to as lineage specification. It requires a potentially multi-step decision sequence which determines the fate of the cells and their successors. It is generally accepted that lineage specification is regulated by a complex system of interacting transcription factors. However, the underlying principles controlling this regulation are currently unknown. Here, we propose a simple quantitative model describing the interaction of two transcription factors. This model is motivated by experimental observations on the transcription factors GATA-1 and PU.1, both known to act as key regulators and potential antagonists in the erythroid vs. myeloid differentiation processes of hematopoietic progenitor cells. We demonstrate the ability of the model to account for the observed switching behavior of a transition from a state of low expression of both factors (undifferentiated state) to the dominance of one factor (differentiated state). Depending on the parameter choice, the model predicts two different possibilities to explain the experimentally suggested, stem cell characterizing priming state of low level co-expression. Whereas increasing transcription rates are sufficient to induce differentiation in one scenario, an additional system perturbation (by stochastic fluctuations or directed impulses) of transcription factor levels is required in the other case.


Introduction
The hematopoietic system consists of a variety of functionally different cell types, including mature cells such as erythrocytes, granulocytes, platelets, or lymphocytes, as well as several different precursor cells (i.e., premature cell stages) and hematopoietic stem cells (HSC) (Lord, 1997;Orkin, 2000). Most mature cell types have limited life spans ranging from a few hours to several months, which implies the existence of a source capable of replenishing these differentiated cells throughout the life span of an individual. This supply is realized by the population of HSC, which is maintained and even regenerated after injury or depletion throughout the whole life of the organism. This selfrenewal property is a major characteristic defining HSC (Loeffler and Roeder, 2002;Lord, 1997;Potten and Loeffler, 1990). A second major characteristic of HSC is their ability to contribute to the production of cells of all hematopoietic lineages, thus ensuring the supply of functionally differentiated cells meeting the needs of the organism. The process controlling the development of undifferentiated stem or progenitor cells into one specific functional direction (i.e., one specific hematopoietic lineage) is called lineage specification. It is generally accepted that the process of lineage specification is governed by the interplay of many different transcription factors (Akashi, 2005;Cantor and Orkin, 2002;Cross et al., 1994;Orkin, 1995Orkin, , 2000Tenen, 2003). Experimental results suggest that a number of relevant transcription factors are expressed simultaneously in HSC, although at a low level (Akashi et al., 2003;Hu et al., 1997). Some authors refer to this state of a low level co-expression as priming behavior (Akashi, 2005;Cross and Enver, 1997;Enver and Greaves, 1998). During differentiation the balanced co-expression of these potentially antagonistic transcription factors is assumed to be broken at some point (or even multiple points). Thereafter, the system is supposed to be characterized by an up-regulated level of some transcription factors, specific for a particular lineage, while other transcription factors are down-regulated. These observations suggest a transcription factor network, capable of switch-like behavior by changing from unspecific co-expression to different states of specific expression. However, the general underlying principles of the regulatory mechanisms are currently unknown. Particularly, it is unclear whether the assumption of a dynamically balanced low level co-expression state is justified or whether priming should rather be interpreted as the result of an inactive transcription factor network overlaid by stochastic fluctuations of transcription factor expression.
In this paper we propose a simple mathematical model describing different interaction scenarios of two transcription factors. Biologically, this simple two component network model is motivated by experimental observations on the transcription factors GATA-1 and PU.1, known to be involved in the process of lineage specification of HSC (Du et al., 2002;Oikawa et al., 1999;Rekhtman et al., 1999;Rosmarin et al., 2005;Tenen, 2003;Voso et al., 1994). The zinc finger factor GATA-1 is reported to be required for the differentiation and maturation of erythroid/megakaryocytic cells, while the Etsfamily transcription factor PU.1 supports the development of myeloid and lymphoid cells (reviewed by Cantor and Orkin, 2002;Tenen, 2003). For both, GATA-1 and PU.1, it has been demonstrated that they are able to stimulate their own transcription through an auto-catalytic process (Chen et al., 1995;Nishimura et al., 2000;Okuno et al., 2005;Tsai et al., 1991). Additionally, there are physical interactions between GATA-1 and PU.1 which induce a mutual inhibition and, therefore, favor one lineage choice at the expense of the other (erythroid/megakaryocyte vs. myeloid) (Du et al., 2002;Nerlov et al., 2000;Rekhtman et al., 1999Rekhtman et al., , 2003Voso et al., 1994;Yamada et al., 1998;Zhang et al., 1999Zhang et al., , 2000. In particular, two different mechanisms for the mutual inhibition of these two transcription factors have been suggested by experimental observations: On one hand, GATA-1 binds to the β3/β4 region of PU.1 (complex 1) and displaces the PU.1 co-activator c-Jun from its binding site, thereby, inhibiting the transcription initiation of PU.1 (Zhang et al., 1999). On the other hand, the inhibition of GATA-1 transcription is mediated by the binding of the N-terminal region of PU.1 to the C-finger region of GATA-1 (complex 2), thus blocking the binding of GATA-1 to its promoter (Zhang et al., 2000). That means, although both inhibition mechanisms are interfered through the formation of PU.1/GATA-1 heterodimers, the two complexes are structurally different. Whereas complex 1 (inhibition of PU.1 transcription by GATA-1) is known to bind to DNA, thus occupying a PU.1 promoter site, DNA-binding of complex 2 (inhibition of GATA-1 transcription by PU.1) has not been reported so far.
The mechanisms of antagonistic interdependence together with positive autocatalytic regulation provide a framework for the theoretical investigation of different scenarios of transcription factor interaction and their implications for the explanation of lineage specification control. Applying a mathematical model, which formalizes the described interactions, it is now possible to analyze different combinations of transcription factor activation and inhibition on a qualitative and quantitative level. The proposed model relies on principles suggested for the description of general genetic switches (e.g. Becskei et al., 2001;Demongeot, 2002, 2005;Gardner et al., 2000).
In this paper it is our objective to examine the following questions within the framework of this model structure: • Are the experimentally described interactions of the two transcription fac-tors sufficient to generate a switching behavior between a stable co-expression of two factors and the dominance of one of these factors?
• What are the conditions inducing such a qualitative change in the system behavior? • Is there evidence for a functional role of the (experimentally suggested) priming status?
To answer these question the following strategy is applied. Firstly, the model equations are derived on the basis of the described biological mechanisms of transcription factor interaction for GATA-1 and PU.1 (Section 2). Secondly, this model is analyzed with respect to the existence of steady state solutions and their dependence on the model parameters. According to our objective, to understand the mechanisms leading to switches between different stable system states, we focus our analysis particularly on the determination of bifurcation conditions, considering different scenarios of transcription factor interaction (Section 3). Finally, the obtained results are discussed in relation to the ongoing debate about lineage specification control in the hematopoietic system, specifically with respect to potential explanations of the experimentally suggested low level co-expression of transcription factors (priming) in undifferentiated progenitors and stem cells (Section 4).

Model description
Although our analysis is motivated by experimental observations of specific transcription factor interactions (GATA-1 and PU.1), the model may also be applied in the general context of two interacting transcription factors. In the following, the two transcription factors are denoted by X and Y .

General assumptions
The general design of the model structure is based on the following assumptions which are motivated by the experimental observations outlined in Section 1: • Both transcription factors, X and Y , are able to act as activator molecules: -If bound to their own promoter region, X and Y introduce a positive feedback on their own transcription. This process is referred to as specific transcription ( Fig. 1(a)). -X and Y are both able to induce an overall transcription which also effects potentially antagonistic transcription factors. Although such an interaction is most likely indirect, for the model we consider a mutual activation of X and Y by the opposing transcription factor, which we refer to as unspecific transcription ( Fig. 1(b)). We assume that transcription initiation is only achieved by the simultaneous binding of two X and Y molecules, respectively (i.e., binding cooperativity c = 2). This assumption is motivated by the result that a binding cooperativity c > 1 is a necessary condition for the existence of system bistability (see e.g. Becskei et al., 2001;Cinquin and Demongeot, 2005;Gardner et al., 2000). • There is a mutual inhibition of X and Y . Within this context, two possible mechanisms, based on the formation of two structurally different complexes of X and Y , are considered: -Joint binding of X and Y molecules to promoter sites ( Fig. 1(c)). Here, the DNA-bound XY -complex (Z 1 ) acts as a transcription repressor, which blocks the binding sites. This represents a mode of competition for free binding sites. -Formation of another XY -complex, called Z 2 , which neither binds to X nor Y DNA binding site ( Fig. 1(d)). In contrast to Z 1 , this represents a competition for free transcription factor molecules. Both inhibition mechanisms (including combinations of them) are considered for X as well as for Y .
To facilitate the analysis of the mathematical model we make the following simplifications: • Post-transcriptional regulation is neglected, i.e., the transcription of a gene is considered to ultimately result in the production of the corresponding protein (here, a transcription factor). • Time delays due to transcription and translation processes are neglected.
• Simultaneous binding of X/Y monomers together with a Z 1 -heterodimer, of two Z 1 -heterodimers, as well as of a X and a Y monomer at the same promoter are excluded from the analysis. • Interactions of X, Y as well as the promoter regions of the coding genes with further transcription factors are neglected.
Throughout the paper the following notations are used: x, y denote the molecule concentrations of X and Y , respectively. Z 1 denotes the DNA bound XYcomplex and Z 2 the structurally different XY -complex, which is not able to bind to promoter DNA. D x/y denotes free DNA binding sites within the promoter region of X and Y , respectively. In contrast, binding sites occupied by X or Y molecules or by the XY -complex Z 1 are denoted as D (c-e) Suggested mechanisms of transcription inhibition for X and Y by formation of XY -complexes: (c) A XY -complex, called Z 1 , bound to the promoter regions acts as a repressor; (d) The formation of a structurally different XY -complex (Z 2 ) competitively inhibits the DNA binding of X and Y molecules; (e) Combination of (c) and (d) as suggested for GATA-1 and PU.1 (Zhang et al., 1999(Zhang et al., , 2000 2

.2 Model equations
With these assumptions one can write down a set of chemical reaction equations which underly the system dynamics.
Herein we made the simplifying assumption that the DNA binding of X and Y always occurs as the binding of homodimers. That means, that the sequential binding of two monomers, as the second possibility of DNA binding, is not consider. The process of dimerization, as well as the DNA binding and dissociation, are regarded to be in quasi steady state.
Here and throughout the paper, the K i = k i /k i (i = 1, ..., 7) denote the equilibrium (dissociation) constants of the reactions, with k i andk i representing the forward and backward reaction rate constants, respectively. Finally, it is assumed that both transcription factor monomers, X and Y , decay with first order kinetics at rates k x 0 and k y 0 , whereas dimer-complexes are assumed to be stable.
The different mutual transcription inhibition mechanisms are illustrated in Figs. 1(c)-(e). First, we consider the formation the XY -complex Z 2 (see Fig.  1 Under the quasi steady state assumption Z 2 does not contribute to the mathematical description of the system dynamics. As shown in Fig. 1(c), there is also the possibility that X and Y form a structurally different heterodimer Z 1 , which is able to bind to the promoter regions, acting as a repressor for X and Y transcription, respectively: As with the promoter binding of X and Y , we collapse dimerization, which is assumed to be in quasi steady state, and DNA binding into one process, neglecting the sequential binding of monomers.
Under the posted quasi steady state assumptions, equations (1)-(7) lead to the following set of ordinary differential equations: Details of the derivation are given in Appendix A.

Symmetric system
To analytically derive steady state as well as potential bifurcation conditions, we restrict ourself in this section to the special case of a completely symmetric system, i.e.: (9) can be written in a dimensionless form as Equations (10) and (11) are a pair of coupled first order differential equations. The steady state (ẋ =ẏ = 0) is defined implicitly by The domain of these nullclines for x and y is restricted by the choice of parameters as outlined in Appendix B. The intersections of the nullclines correspond to the fixed points (x * , y * ) of the differential equations (10) and (11). Fixed points on the diagonal (x * , x * ) are traced under the simplifying condition x = y. In this case, equations (12) and (13) can be summarized by The first (trivial) fixed point of equation (14) is x * 1 = 0. Having eliminated this solution, the remaining quadratic equation yields two further non-trivial fixed points at Bifurcation points can be characterized by nullclines intersecting with equal slopes. The derivatives of equations (10) and (11) are evaluated to determine explicit conditions for bifurcation occurrence on the diagonal, considering s as the bifurcation parameter 1 . For simplicity the denominators in equations (12) and (13) are defined as P x = 1+x 2 +k u y 2 +k r xy and P y = 1+k u x 2 +y 2 +k r xy.
The partial derivative of equation (12) with respect to y leads to For bifurcation points on the diagonal (x = y), where the denominators P x and P y simplify to P * = 1 + x * 2 (1 + k u + k r ), equation (18) can be rewritten as Inserting P * in the form P * = x(s + uk u ) derived from equation (14) and neglecting the trivial solution the equality now reads To find the bifurcation points on the diagonal one needs to study the two distinct cases for x ′ = 1 and x ′ = −1 (see Appendix C).
Case I (x ′ = 1): Equation (20) satisfies the condition x ′ = 1 at This only coincides with the fixed points x * 2/3 derived in equation (15) if the expression under the radical in equation (15) vanishes, i.e., x * 2 = x * 3 . This is true for which corresponds to the condition defined in (16). This implies that the "birth" of the fixed points (x * 2 , x * 2 ) and (x * 3 , x * 3 ) coincides with the bifurcation condition x ′ = 1. The sequence of nullclines in Fig. 2(a)-(c) illustrates this behavior for an unspecific transcription rate u = 1. The nullclines in Fig. 2(a) do not intersect for s < s * 1 , i.e., there is no non-trivial fixed point along the diagonal. For s = s * 1 one common fixed point at exists, which marks the bifurcation point depicted in Fig. 2 The nullclines change qualitatively for a further increase in the bifurcation parameter s as shown in the sequence Fig.  2(d),(e).
In the case of a smaller unspecific transcription rate u = 0.4 the corresponding bifurcation is illustrated in Fig. 3(d). The two fixed points (x * 2 , x * 2 ) and (x * 3 , x * 3 ) generated at the diagonal are both unstable as depicted in Fig. 3(e),(f). The qualitative differences between the scenarios for small and large unspecific transcription u are more thoroughly investigated in the subsequent paragraphs.
For small u the condition for the occurrence of further bifurcations s * 2/3 ≥ s * 1 is violated. Numerical results indicate that two saddle-node bifurcations form fixed points off the diagonal at s < s * 1 as depicted in the sequence of nullclines in Figs. 3(b),(c). The saddle-node bifurcation on the diagonal is observed at s * 1 . For large s these scenarios show a comparable pattern of two up-regulated steady states with one high and one low expressed component and a further stable fixed point at (0, 0) (compare Figs. 2(h) and 3(f)). Fig. 4 comprise the above findings for k u = 0.8. The x-coordinate for the fixed points is shown depending on the bifurcation parameter s. For u = 1 (Fig. 4(a)) the birth of two fixed points through a saddle-node bifurcation can be seen at s * 1 , given by equation (22). Condition (24) defines the occurrence of the pitchfork bifurcation on the upper branch (x * 2 , x * 2 ) at s * 2 , whereas condition (25) is the equivalent for the lower branch (x * 3 , x * 3 ) at s * 3 . Note that the additional condition s * 2/3 ≥ s * 1 is fulfilled. The upper branch gives rise to three fixed points, one unstable (arising from the existing stable fixed point) at the diagonal at x * 2 and two new stable fixed points branching off this axis. For the lower case all three fixed points are unstable for s > s * 3 . The inset in Fig. 4(a) enlarges this bifurcation occurring at s * 3 . Fig. 4(b) illustrates the equivalent scenario for u = 0.4. The saddle-node bifurcation at s * 1 represents the formation of the unstable fixed points (x * 2 , x * 2 ) and (x * 3 , x * 3 ). In addition, two further saddle node bifurcations exist that also form stable fixed points. Since s * 2/3 < s * 1 these bifurcations do not occur on the diagonal. All branches in Fig. 4 that do not represent fixed points on the diagonal were determined numerically.  (22), (24) and, for the lower branch, numerical results. In the lower left region only one stable fixed point at (0, 0) exists. In the region marked with "two stable fixed points" one additional up-regulated stable fixed point exists besides the one at (0, 0). In "three stable fixed points" region two additional up-regulated stable fixed points exist. The dashed horizontal lines correspond to the sequences of nullclines in Figs. 2 and 3. numerical results. The sequence of nullclines given in Fig. 2 is illustrated by the dashed line at u = 1, with the dots referring to the subfigures for varying s. The dashed line at u = 0.4 gives a similar representation, with its correspondence in Fig. 3.

The bifurcation diagrams in
Figs. 2 and 3 both indicate that the basin of attraction for the fixed point at the origin (0, 0) is separated from the basins of attraction of the up-regulated stable states by a set of unstable fixed points. The sequences of graphs also illustrate that these unstable fixed points move towards the fixed point at (0, 0) for increasing s, thus continously reducing the size of its basin of attraction. However, this size characterizes the stability of the fixed point at the origin (0, 0) in response to external perturbations. Unlike the intermediate stable steady state, arising from the bifurcation at s * 1 depicted in Fig. 4(a), where a dynamically increasing s inevitably leads to one of the two up-regulated fixed points, the escape from the fixed point (0, 0) needs to be triggered by a perturbation that exceeds the size of its basin of attraction. Given the position of the unstable fixed point at the diagonal (x * 3 , x * 3 ) as a function of s in equation (15), an appropriate measure for the size of the basin of attraction is provided.

Asymmetric system
As indicated by Zhang et al. (1999Zhang et al. ( , 2000 the inhibition of PU.1 by GATA-1 and the converse are based on different mechanisms. The formation of the PU.1-GATA-1 complex, which we refer to as a Z 2 -complex, prevents free transcription factors from binding to their specific DNA binding sites. A competitive inhibition in this form affects both transcription factors, although Zhang et al. (2000) do not explicitly outline the consequences of binding of the PU.1-GATA-1 complex to the PU.1 binding site. On the other hand, GATA-1 prevents the binding of c-Jun to the DNA bound PU.1 protein and thus disables the transcription initiation of PU.1. This process explicitly targets the PU.1 binding sites and introduces a functional asymmetry of inhibition mechanisms.
The mathematical counterpart of this asymmetry is a specific binding rate K 6 > 0 while keeping K 7 = 0 (see equations (8), (9)). In terms of the dimensionless formulation in equations (10) and (11) this translates into two different rate constants k rx > 0 and k ry = 0. The additional binding mode (k rx > 0) can be interpreted as a reduction in the transcriptional activity of the X gene conferring a disadvantage relative to Y .
For any k rx > 0 there is a symmetry breaking which shifts the previously observed bifurcations off the diagonal and destroys the pitchfork bifurcation observed in the symmetric case for large u. The two up-regulated stable fixed points are not created instantaneously by the transformation of a previous stable state at the diagonal, but the initial stable point remains unchanged while a further (saddle node) bifurcation forms the second up-regulated stable point alongside with one unstable fixed point. This scenario is shown in the sequence of nullclines in Fig. 6(a)-(c). The parameter k rx regulates the distance between the up-regulated stable points and the extension of their basins of attraction. This is visualized in the bifurcation diagrams in Fig. 6 For small unspecific transcription rates u, where in the symmetric case the additional up-regulated stable states are created off the diagonal, no qualitative changes are introduced by the functional asymmetry.
The introduction of asymmetry is not necessarily based on different interaction mechanisms. It is plausible that auto-regulative transcription activation does not require identical transcription rates for the genes of interest. This can be described by relaxing the symmetry assumption of Section 2.1, which leads to gene specific transcription rates s x and s y . This asymmetry in transcriptional activity results in a qualitatively similar symmetry breaking as in the case of the mechanistic asymmetry: the pitchfork bifurcation occurring for large u is replaced by a remaining stable state alongside a saddle-node bifurcation forming the second up-regulated stable state (data not shown). The magnitude of the difference in the specific transcriptions rates s x and s y regulates the distance between the up-regulated stable states in the phase plane.
In a scenario where asymmetry of interaction mechanisms occurs alongside an asymmetry in the specific transcription rates, the effects on the system behavior combine, either amplifying or compensating each other.

Over-expression scenarios
Induced over-expression of a certain critical component is a common experimental method to study interaction dynamics between different transcription factors and has also been applied to the GATA-1 / PU.1 system (Nerlov et al., 2000;Rekhtman et al., 1999;Zhang et al., 2000). These experiments provide insight in the stability of the system, interaction time scales, and the role of co-factors and interaction mechanisms. We have applied an over-expression impulse of amplitude a oe and duration d oe to the model system given in equations (10), (11). Characteristics of the dynamic response are only valid under the outlined steady state assumptions. A qualitative overview of the simulation results is presented in Fig. 7. Starting from a fully symmetric system as studied in Section 3.1 where, for large s, the system is in one of the two up-regulated states (characterized by one high and one low expressed transcription factor) two modes of over-expression are applied: a short impulse over-expression of the lower expressed component and a long and steady overexpression of the same component. Not surprisingly, the model reacts to the over-expression with two distinct scenarios, depending on the intensity of the impulse. For a subcritical over-expression the system returns to the previous expression level (indicated in Figs. 7(a) and (d)), whereas for a supercritical situation the former expression state is reversed (indicated in Figs 7(b), (c), (e) and (f)). Translating this picture into the x vs. y phase plane, the supercritical over-expression corresponds to a change from one basin of attraction to another, induced by a crossing of the separatrix. Most available experimental techniques to artificially induce gene expression lead to a massive over-expression that significantly exceeds physiological levels, a scenario still underestimated by Figs. 7(c) and (f). A sensitively tuned expression experiment is more promising to elucidate critical intensities and time scales necessary to induce a permanent shift in the genetic expression patterns and thus to characterize the stability of the initial states.

Discussion
The presented model of transcription factor interaction is based on principles of coupled feedback regulations, which have previously been proposed for the description of general genetic switches (Becskei et al., 2001;Cinquin and Demongeot, 2005;Francois and Hakim, 2004;Gardner et al., 2000;Glass and Kauffman, 1973) and the modeling of prokaryotic gene regulation (McAdams and Arkin, 1998;Santillán and Mackey, 1998. Here, specific experimental knowledge of activation and inhibition mechanisms of two transcription factors (GATA-1 and PU.1), which play a key role in the myeloid/erythroid differentiation process of hematopoietic progenitor cells, is incorporated in this general framework.
Our model analysis particularly focuses on the investigation of the steady states of transcription factor expression and there dependence on parameter changes. In this context, we are able to analyze the experimentally suggested To facilitate the mathematical analysis, a number of simplifications have been made. We interpret the transcription factors described in the model (X and Y ) as representatives of a more complex factor formation rather than an explicit model of PU.1 and GATA-1 alone. Also, we are aware that most of the statements resulting from the model analysis are only semi-quantitative in the sense that for all model parameters, as there are DNA binding-, decay-, and transcription-rates, no experimentally determined estimates are available for the investigated system. In the same line of argumentation, details of the transcription/translation process, like the DNA binding sequence of transcription factor molecules and the delay induced by the processes of transcription and translation, have been excluded from the analysis. Although such phenomena can influence the dynamics of the system (Bundschuh et al., 2003;Vilar et al., 2002), these effects are speculative since detailed information about relevant rates and time scales are not available. The simplifications arising from the quasi steady state assumption outlined in Section 2 for dimerization and DNA binding impose further limitations on our model with respect to the exact description of the system dynamics (c.f. Pirone and Elston, 2004). However, these simplifications do not effect the steady state behavior, and, thus, do not alter the results derived in Section 3.
The functional role of the so called priming behavior is a question of particular biological relevance which is addressed by this model. It has been suggested that low level co-expression of multiple transcription factors, specific for different lineages, might be a characteristic of (hematopoietic) tissue stem cells (Akashi, 2005;Cross and Enver, 1997;Orkin, 2000). However, it is currently unclear whether priming corresponds to a stable state of low level co-expression or to a truly zero-expression overlaid by some random expression noise. Furthermore, there is a hypothesis that lineage specification induction might be a two stage process with a primary initialization of transcription factor network interaction (i.e., a transition from no expression to low level co-expression) and a secondary network-induced differentiation process (Enver and Greaves, 1998). This perspective immediately leads to the questions under which conditions such a two stage process can be established and whether such a sequence of different activation states of the transcription factor network requires (multiple) external induction signals or whether it represents a system inherent development.
The suggested model generates two characteristic modes of system stability depending on the magnitude of the specific transcription rate s: For small s only the trivial fixed point (0, 0) exists; for large s two additional up-regulated stable states are observed that are marked by the dominance of one factor over the other (dominated co-expression). These modes are maintained independently of a mechanistic or parametric asymmetry. Assuming a differentiation initiation by increasing the transcription rate s (e.g. by changes in chromatin structure (Berger and Felsenfeld, 2001;Rosmarin et al., 2005) or by alterations in activation/inhibition complexes (Hume, 2000)), the transition between the different stable states is the central mechanism characterizing lineage specification.
Within the proposed biological framework the trivial fixed point at (0, 0), which exists for all values of s, can be identified with the undifferentiated state of a cell where neither activation nor decision processes are observed. It should be mentioned that stability of this fixed point is specific for the outlined model and has not been observed for the general case of a toggle switch (c.f. Gardner et al., 2000). In logical extension, the two up-regulated stable fixed points, observed for large s, would be interpreted as expression states promoting one or the other lineage. These distinct states are characterized by a high auto-regulative expression of one dominating factor and a reduced expression of the antagonistic factor.
In Sections 3.1 and 3.2 it has been demonstrated that, for increased unspecific transcription rates u (but only in this case!), a further stable fixed point exists for intermediate s prior to the formation of the two distinct states of dominated co-expression. This particular fixed point is characterized by a balanced low level co-expression of the two antagonistic factors where no final commitment decision has been made. The resulting transition sequence between three distinct regions of multi-stability can be interpreted as a possible explanation for a two stage differentiation process mentioned above.
The induction of a system change from the stable trivial fixed point to the dominated or, if existent, balanced low level co-expression state, needs to be triggered either by a stochastic background expression or by an active impulse on the system. The unstable fixed point separating the trivial from the upregulated stable states is an indicator of the size of the basins of attraction. The observation that the unstable fixed point approaches the trivial one for increasing s indicates that the magnitude of the perturbation to introduce a transition from the zero-state to the co-expression states decreases in the same fashion: for a sufficiently large s even a small perturbation is able to initiate differentiation.
Concluding from these results, there are two different scenarios to explain the experimentally suggested priming behavior within the proposed model framework: (1) Priming might be considered as the existence of perturbations in the expression of transcription factors, imposed on a zero-expression state represented by the trivial fixed point at (0, 0), either in the form of stochastic background fluctuations (functional noise) or by active impulses. In this scenario, the perturbations are necessary components of the regulatory system to induce a differentiation process. It points to the potential role of stochastic effects in the context of decision making in stem cell differentiation as frequently suggested (see Kaern et al. (2005) for a review). (2) In contrast to this scenario, priming can also be explained by the balanced low level coexpression state, which becomes unstable for increasing specific transcription rates. Due to this parameter dependent loss of stability, this scenario would lead to differentiation without the need for external perturbations 2 . However, the balanced low level co-expression state is only existent if there is a certain degree of unspecific transcription.
Currently, our results do not allow to decide between the two scenarios. The introduction of artificial differentiation impulses of different intensities on uncommitted cells might be an appropriate way to tackle this question experimentally. Whereas, a low level co-expression priming (like in scenario (2)) would be unaffected by these perturbations, the system could be enforced to escape the priming status in scenario (1). Moreover the existence and the stability of the different stable system states depend sensitively on the model parameters. Due to the lack of available data on transcription and binding rates, we are currently not able to specify the biological relevant regimes more rigorously. Any experimental approximation of binding and transcription rates for the involved components supports the identification of the nature of priming.
The over-expression scenarios presented in Section 3.3 fail to explain experimental findings described by several authors (Nerlov et al., 2000;Rekhtman et al., 2003;Zhang et al., 2000). In spite of the induced up-regulation of one transcription factor it was observed that the transcription level of the antagonistic transcription factor remained more or less constant. These observations are in contrast to the model results presented here, in which the induced over-expression of the initially low expressed factor shifts the equilibrium to the opposing co-expression state. Retaining our model assumptions, a potential interpretation can be given as follows: One of the major functions of transcriptional regulators like GATA-1 and PU.1 is the activation of a set of lineage-specific genes which include further transcription and growth factors as well as functional components of the committed lineages (Tenen, 2003). In Sieweke and Graf (1998) and Tsai et al. (1991), the authors point to a continuously modulated set of cooperative lineage-inherent transcription factors changing with the state of differentiation. Such secondary complexes of transcription factors could in turn act as activators of the initial transcription factor, substituting for a simple auto-regulation and thus stabilizing the initial up-regulation pattern (Hume, 2000). In such a scenario our model would only account for the initial switching process. The experimentally observed stable transcription level of the antagonistic factor in over-expression experiments could be interpreted as a substitution of the auto-regulation by secondary transcription factor complexes.
Summarizing, the presented model is able to provide a quantitative explanation for possible mechanisms underlying lineage specification control in eukaryotic systems. It is able to generate parameter dependent changes in the system behavior, with alteration of the number of possible stable steady states. Specifically, the model explains states of stable co-expression as well as the situation characterized by an over-expression of one factor over the other. The conditions inducing shifts from one to another stable state (e.g. parameter choice, degree of system disturbances), however, depend in a sensitive manner on the assumed activation and inhibition mechanisms. Using the mathematical model, we were able to test several combinations of experimentally described feedback mechanisms with respect to their influence on the resulting stable states and provide possible explanations for the experimentally suggested differentiation priming of stem cells.

Acknowledgment
The authors thank Michael Mackey for his encouragement and critical discussions in the process of preparing this manuscript, and Michael Cross for his explanations of many biological details. Furthermore we acknowledge the Human Frontier Science Program which supported initiation of this work by short-term fellowship to I.R. at the Centre for Nonlinear Dynamics, McGill University, Montreal.

A Derivation of transcription factor dynamics
It is assumed that the transcription of transcription factors X and Y requires the existence of activator complexes, i.e., the binding of X or Y dimers to the promoter regions of X (D x ) and Y (D y ), respectively. As described in Section 2.2, we distinguish between a specific (see equations (1),(3)) and an unspecific (equations (2),(4)) transcription activation. Furthermore, there is the possibility that X and Y can act jointly as a repressor dimer Z 1 , inhibiting the DNA binding of the X and Y activator dimers (see equations (6), (7)).
The total amount of promoter sites for X and Y can be specified as the sum of unbound (free) and occupied (by repressor or activator molecules) promoter regions, i.e., obtained from assuming equations (1)-(4), (6), (7) to be in a quasi steady state, the fraction of promoter sites contributing to active X and Y transcription is given by K 1 x 2 + K 2 y 2 1 + K 1 x 2 + K 2 y 2 + K 6 xy (A.3) and D yy y + D xx y D tot y = K 3 y 2 D y + K 4 x 2 D y D y + K 3 y 2 D y + K 4 x 2 D y + K 7 xyD y = K 3 y 2 + K 4 x 2 1 + K 3 y 2 + K 4 x 2 + K 7 xy (A.4) respectively.
Taking the (first order) decay rates of X and Y into account, one immediately obtains equations (8), (9) by writing down the balance equations for X and Y .

B Domain of the nullclines
Under the equilibrium assumption equation (12) can be solved for: which describes a set of nullclines. There is an obvious singularity at x = u.
The solutions y 1/2 (x) are real for 0 < h(x) = k 2 r x 4 − 4k u x(u − x)(sx − 1 − x 2 ) with h(x) defined as the expression under the square root in the previous equation. For k r = 0 the roots of h(x) are located at where the singularity at x = u = x h 3 marks a boundary for distinct intervals in the domain: Fig. B.1(b), and

C Derivation of bifurcation condition
The nullclines of the symmetric system derived in (12) and (13) are interpreted as functions of x and y: x = f (x, y) = sx 2 + uk u y 2 1 + x 2 + k u y 2 + k r xy , (C.1) y = g(x, y) = sy 2 + uk u x 2 1 + k u x 2 + y 2 + k r xy . (C.2) To derive bifurcation conditions one has to determine the point of tangency of the nullclines f (x, y), g(x, y) at a steady state (x * , y * ), i.e.
From these statements, it follows that we have to consider the following equal-ities to find the bifurcation conditions for the symmetric system, restricting to symmetric steady states of the form (x * , x * ): df (x, y) dx (x * ,x * ) = dg(x, y) dx (x * ,x * ) = |1|. (C.4)