INTRODUCTION

Many developmental patterning processes are controlled by morphogen gradients (Freeman and Gurdon,2002; Tabata and Takei,2004). In most cases, these gradients are established within fields of cells by diffusion of a secreted signaling molecule from a localized source. Distinct cell fates are induced at defined ranges of the morphogen concentration along the gradient. Conversion of a graded morphogen signal into discrete cellular territories implies the existence of sharp thresholds of morphogen concentration that cells can distinguish and interpret to choose between alternative cellular fates. The question of how such thresholds arise along morphogen gradients has long been investigated (Lewis et al.,1977; Meinhardt,1982; Goldbeter and Wolpert,1990; Houchmandzadeh et al.,2002; Bollenbach et al.,2005; Howard and Ten Wolde,2005; Melen et al.,2005) and is of key importance in patterning processes controlled by morphogen gradients, such as segmentation of the body axis (Pourquié,2003).

The segmented or metameric aspect of the body axis is a basic characteristic of many animal species ranging from invertebrates to human. The vertebrate body is built on a metameric organization, which consists of a repetition along the antero-posterior (AP) axis of functionally equivalent units, each comprising a vertebra, its associated muscles, peripheral nerves, and blood vessels. The segmented distribution of vertebrae derives from the earlier metameric pattern of the embryonic somites, which are epithelial spheres generated in a rhythmic fashion from the mesenchymal presomitic mesoderm (PSM). The segmental pattern was proposed to be established in the PSM by a mechanism involving an oscillator (the segmentation clock), which is thought to set the periodicity of the process, and a traveling wavefront defined by antagonistic gradients of the signaling molecules fibroblast growth factor (FGF) and retinoic acid (RA), which control the spacing mechanism of somite boundaries (Pourquié,2003).

Segmentation is first manifest as a pre-pattern of striped gene expression, initially seen at a defined level of the PSM called the determination front (Fig. 1) (Dubrulle et al.,2001). Position of the determination front is set by the antagonistic gradients of FGF and RA signaling that arise from the posterior and anterior ends of the PSM, respectively (Fig. 1) (Dubrulle et al.,2001; Sawada et al.,2001; Diez del Corral et al.,2003; Moreno and Kintner,2004). The posterior-to-anterior gradient of FGF signaling is initially set up as an fgf8 mRNA gradient, which is then translated into a ligand and then a signaling gradient across the PSM (Dubrulle and Pourquié,2004). This FGF gradient is thought to be responsible for maintaining the cells of the posterior PSM in an immature state and to prevent them from activating their segmentation program (Dubrulle et al.,2001). An RA gradient is established in opposite orientation in the PSM to control the activation of segmentation genes (Diez del Corral et al.,2003; Moreno and Kintner,2004; Vermot and Pourquié,2005). Retinoic acid is produced by RALDH2, the RA biosynthetic enzyme, which is expressed in an AP gradient in the somites and PSM (Blentic et al.,2003; Diez del Corral et al.,2003; Novitch et al.,2003; Vermot et al.,2005). Furthermore, RA is degraded by the enzyme CYP26 of the cytochrome P450 family, which is expressed in the posterior PSM (Fig. 1) (Blentic et al.,2003). Thus, the RA gradient is viewed as a classical morphogen diffusion gradient, with a source in the anterior PSM and a sink in the posterior end of the embryo.

Figure 1.

Opposite gradients of FGF and RA signaling in the developing PSM. Schematic representation of the expression domains of fgf8 and raldh2 and of their respective targets Cyp26 and mesp2. Gene expression domains are shown in gray.

The FGF and RA signaling pathways are coupled through mutual inhibition. Thus, FGF signaling activates cyp26 and represses raldh2 expression in the PSM, whereas RA signaling restricts FGF signaling to the posterior PSM either by restricting fgf8 mRNA expression or by activating the dual specificity phosphatase MKP3, which in turn antagonizes FGF signaling (Diez del Corral et al.,2003; Moreno and Kintner,2004). At the determination front level, cells respond to a periodic signal from the segmentation clock, by activating the mesp genes in a segment-wide domain (Fig. 1) (Morimoto et al.,2005). mesp gene expression is repressed by high FGF signaling (Delfini et al.,2005) and requires RA signaling for their transcription, thus explaining their restricted expression at the determination front level (Moreno and Kintner,2004).

In this report, we propose a model in which the mutual antagonism of FGF and RA gradients generates a sharp threshold associated with a phenomenon of bistability. This phenomenon, which involves the coexistence between two stable steady states, occurs in a spatial window within the PSM. We suggest that the abrupt bistable steady-state switch that can occur in this window can explain the coordinated segmental gene activation that takes place in the presumptive segment at the determination front level in response to the periodic signal of the segmentation clock.

RESULTS AND DISCUSSION

Bistability Arising From Mutual Inhibition of RA and FGF Signaling

We consider a relatively simple model for the mutual inhibition of RA and FGF (see Model and Simulation Procedures section). This model is one in a family of possible models, which differ by the details of the inhibitory processes. The results should not depend qualitatively on the particular implementation of the inhibitions as long as the mutual inhibition of RA and FGF is sufficiently strong. Thus, in the model considered (Fig. 2), we assume that FGF promotes the expression of the enzyme CYP26 that hydrolyzes RA (this effect will be measured by an activation constant KA), while RA impairs the synthesis of FGF from fgf8 mRNA (this effect will be measured by an inhibition constant KI). Similar results should be obtained when assuming that RA promotes the degradation of fgf8 mRNA or that FGF inhibits expression of raldh2. Experimental evidence exists in favor of these two processes, which could further contribute to mutual inhibition of RA and FGF (Diez del Corral and Storey,2004). Moreover, we consider that RA synthesis and FGF synthesis are governed by two opposite gradients along the PSM. Thus, the concentration of the RA-synthesizing enzyme decreases from a maximum value RALDH20 in the anterior extremity down to zero at the posterior extremity of the PSM, while the gradient in fgf8 mRNA decreases from a value of M0 down to zero in the opposite direction. The kinetic eqs. [1]–[4] govern the time evolution of the four variables of the model, namely the concentrations of RA, FGF, CYP26, and cyp26 mRNA, while eqs. [5] define the gradients in RA and FGF synthesis (see Model and Simulation Procedures section).

Figure 2.

Scheme of the regulatory interactions between RA and FGF signaling. Synthesis of RA is catalyzed by the enzyme RALDH2, while RA is hydrolyzed by the enzyme CYP26. The inhibitory effect exerted on RA by FGF arises from the induction of cyp26 expression by FGF. The inhibition exerted by RA on FGF occurs through impeding the rate of fgf8 mRNA translation. As shown by the model built according to this regulatory scheme, bistability readily arises from the mutual inhibition between RA and FGF.

We wish to explore the dynamics resulting from mutual inhibition of RA and FGF in a given point along the PSM axis. To that end, we first determined the maximum rate of RA synthesis and the amount of fgf8 mRNA in this point according to the values predicted by the gradients, and then used these values to integrate numerically the kinetic equations that govern the time evolution of the four variables of the model. We then determined the steady states to which the system evolves in the course of time, as a function of parameters KI and KA, which measure, respectively, the strength of inhibition of FGF signaling by RA and of RA signaling by FGF. Inhibition of FGF signaling by RA is enhanced when the inhibition constant KI diminishes, while the inhibition of RA signaling by FGF increases with the activation constant KA.

To characterize the steady state by a single quantity reflecting the relative importance of FGF and RA signaling, we use the ratio ρ of relative saturation of FGF and RA receptors (see Model and Numerical Simulations section). Thus, values of ρ larger or smaller than unity will reflect predominance of FGF over RA signaling, or of RA over FGF. In Figure 3, we show the steady states of the coupled RA-FGF signaling system as a function of parameters KI and KA. In Figure 3A, ρ is determined as a function of KI while KA is kept at an intermediate value KA = 0.2 nM. Similarly, in Figure 3B, ρ is obtained as a function of KA for KI = 0.2 nM.

Figure 3.

Mutual inhibition of FGF and RA generates bistability. The curves show the ratio ρ of FGF to RA receptor occupancy at steady state, given by eqn. [13], as a function of KI (A) and KA (B). A value of ρ larger (smaller) than unity indicates the predominance of FGF (RA) over RA (FGF) signaling. The results have been obtained in a particular point of space corresponding to the value x = 15, which defines a pair of values of νs1 and MF according to eqs. [5]. Parameter values are: M0 = 5 nM, RALDH20 = 7.1 nM, ks1 = ks2 = ks3 = kd3 = kd4 = 1 min−1, kd1 = 1 nM−1min−1, kd2 = 0.28 min−1, kd5 = 0, n = m = 2, V0 = 0.365 nM/min, VsC = 7.1 nM/min, L = 50, Kr1 = Kr2 =1 nM. Moreover, KA = 0.2 nM (A) and KI = 0.2 nM (B). FGF and RA receptor occupancy were obtained after determining the steady-state values of variables RA, MC, C, and F through numerical analysis integration of eqs. [1]–[4] by means of the AUTO program. The stable steady states were also determined by integrating these equations numerically using the Berkeley Madonna program, starting from various initial conditions for F, in the range 0–100 nM, and setting the other initial conditions as RA = MC = C = 0.1 nM. When bistability occurs in the region bounded by the two vertical dashed lines, the evolution toward one or the other stable steady state depends on the initial conditions. Then, above a threshold initial value of F, the system evolves to branch 1 (high FGF/low RA), while below this threshold the system evolves to branch 2 (low FGF/high RA). An unstable steady state (dashed line) separates the two stable steady states. In the absence of bistability, regardless of the initial conditions the system evolves to the same, unique, stable steady state. The above numerical values considered are semi-arbitrary and were chosen to exemplify bistability. The phenomenon occurs over a wide range of parameter values, generally in a window as a function of a given control parameter, as shown here and also in Figures 4–9. Bistability is intimately associated with the pattern of mutual inhibition that links RA and FGF signaling. The hatched vertical lines delimit the bistability window.

Figure 4.

Bistability as a function of KI and KA. A: Domain of bistability in the parameter plane formed by the activation constant KA and inhibition constant KI, which measure, respectively, the regulatory effects of FGF and RA. The remaining panels show the steady-state values of (C) RA, (D) FGF, and of (B) the ratio ρ of FGF to RA receptor saturation as a function of KA for the four increasing values of KI indicated by horizontal lines in A. In the bistability range, the middle steady state is unstable. The diagrams were established by means of the program AUTO (see Model and Simulation Procedures section) for the parameter values of Figure 3.

Figure 5.

Bistability as a function of rate of RA synthesis and of fgf mRNA. A: Domain of bistability in the parameter plane formed by the rate of RA synthesis by RALDH2, νs1 (in nM/min) and by the concentration of fgf mRNA, MF (in nM). The remaining panels show the steady-state values of (C) RA, (D) FGF, and (B) the ratio ρ of FGF to RA receptor saturation as a function of MF for the three increasing values of νs1 indicated by horizontal lines in A. The diagrams were established by means of the program AUTO (see Model and Simulation Procedures section) for the following parameter values: ks2 = ks3 = 0.1 min−1, kd1 = 0.2 min−1, kd2 = kd3 = kd4 = 0.01 min−1, n = m = 2, V0 = 0.001 nM/min, VsC = 0.1 nM/min, Kr1 = Kr2 = 1 nM, KI = KA = 1 nM. Here, νs1 and MF, which appear in eqs. [1] and [4], are directly treated as bifurcation parameters without using their definition by eqs. [5].

Figure 6.

Bistability as a function of kd5. A: Domain of bistability in the parameter plane formed by the activation constant KA and inhibition constant KI, for different values of parameter kd5, which measures the rate of basal RA degradation, not regulated by FGF. B: Ratio ρ of FGF to RA receptor saturation for increasing values of kd5. The curves show that the bistability domain progressively shrinks and eventually disappears as kd5 increases. The diagrams were established as a function of KA by means of the program AUTO for the parameter values of Figure 3.

Figure 7.

Gradients of mutually inhibiting RA and FGF signaling define a bistability window along the PSM. A: Spatial distribution of RA. In a domain bounded by two critical values of position x, two branches of stable steady states coexist (filled squares or circles). These branches are separated by a branch of unstable steady states (dashed line). B: Spatial distribution of FGF. The stable steady states (filled squares or circles) correspond to the steady-state levels of RA indicated by the same symbol in A. C,D: Based on the levels of RA and FGF in A and B, the curves show, respectively, the occupancy of the RA and FGF receptors, defined by eqn. [12]. E: Ratio of FGF to RA receptor saturation, determined according to eqn. [13]. A value of ρ larger (smaller) than unity indicates the predominance of FGF (RA) over RA (FGF) signaling. The blue and red colors refer, respectively, to FGF or FGF receptor saturation, and to RA or RA receptor saturation. In A–E, data (points) have been obtained for eqs. [1]–[4] by use of the program AUTO that generates the stable or unstable steady states as a function of position x, which is considered as a bifurcation parameter. Data have also been obtained by numerical integration of eqs. [1]–[4] as a function of position x, using different initial conditions for F in the range 0–100 nM, with the initial values RA=100 nM and Mc = C = 0.1 nM. A–E were established in the absence of diffusion of RA or FGF for the parameter values of Figure 5, with M0 = 5 nM, RALDH20 = 4 nM, ks1 = 4 min−1, L = 100 (this length corresponds to a PSM length of 1 mm). For these parameter values, the range of bistability along the PSM extends from x = 6.2 to 61.8. The linear gradients in RALDH2 and MF obey eqs. [5], with RALDH2 = RALDH20 (1-(x/L)). The vertical dashed lines in A and B corresponding to points x1 and x2 delimit the bistability window.

Figure 8.

Bistability in the presence of exponential gradients in RALDH2 and fgf mRNA. The ratio ρ of FGF to RA receptor saturation is shown as a function of position x in the presence of nonlinear gradients. A:RALDH2 decreases exponentially from the value RALDH20 = 4 nM to 0 as x increases, while fgf mRNA (MF) decreases exponentially from the value M0 = 5 nM to 0 in the opposite direction. The data are obtained using eqs. [1]–[4], for the parameter values of Figure 7. The exponential gradients are defined as MF = M0 exp[a1(x−L)/L] and νs1 = ks1RALDH20 exp (-a2x/L), with a1 = a2 = 5. The value of MF or RALDH2 at the boundary is divided by 2 in B,C, respectively, while in D the two values are halved as compared with case A. E: The domain of bistability for cases A–D is compared to make clearer the slight movement of the bistability range in the different situations considered. The vertical dashed lines corresponding to points x1 and x2 delimit the bistability window.

Figure 9.

Effect of diffusion of RA and FGF on the occurrence of bistability. To investigate the effect of diffusion of RA and FGF, we extended the model by distinguishing between the intracellular and extracellular concentrations of the two chemical species and by including diffusion terms for extracellular RA and FGF. The set of four kinetic eqs. [1]–[4] must then be replaced by the set of six kinetic eqs. [6]–[11] (see Model and Simulation Procedures section). The curves show the ratio ρ of FGF to RA receptor saturation functions defined by eqs. [13] and [14] for the 6-variable model. A: Bistability in the extended, 6-variable model is illustrated in the absence of diffusion by setting the diffusion coefficients of RA and FGF, DRA and DF, equal to zero. B: Bistability can similarly be obtained in the presence of diffusion of RA and FGF. The diffusion coefficient of RA, DRA, was given the experimentally determined value of 6 × 10−6 cm2/min (Eichele and Thaller,1987). Given that diffusion of FGF is much slower (Scholpp and Brand,2004), the diffusion coefficient DF was given the value 6 × 10−9 cm2/min. In a narrow range of position x, the system evolves to either one of two different steady states depending on initial conditions. Outside this range, the system evolves to the same steady state regardless of initial conditions. Parameters values are as in Figure 7, with KA = 0.005 nM, kd5 = 0.01 min−1, kt1 = 0.1 min−1, kr1 = 0.1 min−1, kt2 = 0.001 min−1, kr2 = 0.001 min−1, ke1 = 0.1 min−1, ke2 = 0.2 min−1, θ = 10, Kr1= 1nM, Kr2 = 0.01 nM (Olwin and Hauschka,1989). Gradients in fgf mRNA and RALDH2 are as in Figure 7.

When the inhibition exerted by RA on FGF signaling predominates, a single steady state exists in which the RA response is high compared to FGF. The converse situation occurs when FGF inhibition of RA signaling predominates. When the two inhibitions are very weak, FGF and RA both reach levels close to those obtained in the absence of inhibitory coupling. A markedly different picture is obtained when the two inhibitory effects are sufficiently strong and of comparable magnitude. Then, depending on the initial conditions, the system evolves either to a steady state with high FGF and low RA, or to a steady state with low FGF and high RA. This phenomenon of coexistence between two stable steady states is known as bistability.

Bistability naturally gives rise to a threshold phenomenon. Indeed, the two stable steady states are separated by an unstable steady state (the set of unstable steady states obtained as a function of KI or KA is represented by a dashed curve in Fig. 3). Stability of a steady state means that when the system moves away from it as a result of some perturbation, it will eventually return to this steady state. However, in the case of bistability, there exists a threshold associated with a critical perturbation beyond which the system will evolve to the alternative stable steady state. Each of the two stable steady states possesses its own basin of attraction defined by the set of initial conditions from which the system will evolve to this particular steady state. The threshold, which reflects the existence of the unstable steady state, corresponds to the minimal perturbation that drives the system out of the basin of attraction of one stable steady state into the basin of attraction of the other stable steady state. As a result of this threshold-generating mechanism, sharp transitions in the levels of FGF and RA at steady state may occur in the region of bistability. At a given value of KI and KA, these transitions are associated with all-or-none switches between the two stable steady states.

Besides depending on KA and KI, the bistability phenomenon depends also on the other parameters that govern the levels of RA and FGF. These parameters include the rate constants for synthesis and degradation of RA by RALDH2 and CYP26, the rate of expression of cyp26 and half-life of its mRNA, and the corresponding parameters for FGF. To determine the robustness of the bistability domain to variations of the system parameters, we constructed bifurcation diagrams in the KA-KI and RALDH2 (νs1)-MF parameter planes. Thus, Figure 4A shows the domain of bistability as a function of KA and KI. Four horizontal cuts corresponding to increasing values of KI are made through and above the bistability domain. Figure 4B–D shows for each of these four values the steady-state value of RA, FGF, and of the ratio ρ of FGF to RA receptor saturation as a function of KA. The comparison of the curves indicates that bistability disappears predominantly in a horizontal manner through progressive shrinking of the bistability domain, although a vertical component can also be seen, since the upper and lower branches of steady states progressively come closer to each other. A similar conclusion can be drawn from the data shown in Figure 5 where bistability is investigated as a function of RALDH2 (νs1) and MF. Figure 5 further indicates that the domain of bistability widens as these two parameters progressively increase. The disappearance of bistability through shrinking of the bistability domain can be seen directly from the bifurcation diagrams both when KA or KI increases (Fig. 4A–D) and when RALDH2 or MF decreases (Fig. 5A–D). Together, these data indicate that for defined sets of strength of mutual inhibition of the RA and FGF gradients, the concentration of RA and FGF in the PSM can reach two stable steady states depending on the initial conditions. The set of conditions where the two steady states co-exist defines the bistability window.

The inclusion of two parallel routes for RA degradation in the model allows us to further illustrate the role of mutual inhibition in the occurrence of bistability. In Figure 6, the domain of bistability is determined as a function of parameter kd5, which measures RA degradation not catalyzed by CYP26 (i.e., independent from FGF). An increase in kd5 has the effect of masking the contribution of CYP26 to the degradation of RA. At large values of kd5, the effect of FGF-controlled RA degradation via CYP26 becomes negligible, so that mutual inhibition of RA and FGF effectively ceases to operate. The results in Figure 6 show, accordingly, that the domain of bistability progressively shrinks and eventually vanishes as kd5 increases.

Bistability Becomes Localized in the Presence of Gradients

In the absence of gradients (i.e., in spatially homogenous conditions) and for appropriate values of KA and KI, bistability would occur everywhere in space along the PSM. However, we know that RA and FGF establish mutually antagonist gradients along the PSM. These gradients allow the spatial positioning of the bistability window in the embryo. To illustrate the role of these antagonistic gradients, we first consider the simplest case where they are linear. Thus, the level of fgf8 mRNA goes from the highest value, M0 at the posterior end of the PSM down to zero at the anterior end, while the rate of RA synthesis by RALDH2 goes from a maximum value of ks1RALDH20 at the anterior end to a zero value at the posterior end. In Figure 7, we examine how the system of antagonistic FGF and RA signaling behaves in the PSM in the presence of two opposite gradients of RALDH2 and fgf8 mRNA controlling, respectively, RA and FGF synthesis. We will first describe the behavior of the system in the absence of diffusion and will then determine how the results are affected when taking into account the diffusion of RA and FGF.

In theory, sharp transitions in the levels of RA and FGF could occur for two critical values of position x along the PSM AP axis, which define the bistability domain bounded by the two vertical lines at positions x1 and x2 in Figure 7A,B. Within this domain, at a given position x, two stable steady-state levels of RA (Fig. 7A) and of FGF (Fig. 7B) coexist, separated by an unstable steady state. Upon moving continuously from a posterior-to-anterior position along the PSM (0 marks the level of the newly formed somite while 100 corresponds to the posterior tip of the PSM), starting from branch 1 (high FGF/low RA), the system will shift abruptly to branch 2 (low FGF/high RA) when the domain of bistability ends at position x1 (Fig. 7A,B). As the system shifts from branch 1 to branch 2, the level of FGF receptor saturation falls abruptly while the level of RA receptor occupancy rises sharply (Fig. 7C,D). The ratio of FGF to RA receptor occupancy allows us to capture this switch in a single curve (Fig. 7E).

To address the effect of the shape of the gradients, we investigated the occurrence of bistability as a function of position x when the gradients in RALDH2 and fgf mRNA become exponential. Again, as expected, a domain of bistability can be demonstrated in these conditions (Fig. 8A). The width of this domain depends on the characteristics of the gradients, such as the maximum value at the boundary and the steepness of the exponential decline in RALDH2 and fgf mRNA. The domain of bistability is not much affected when the values of RALDH2 or MF at the boundaries are decreased by 50% (Fig. 8B–D). As shown in Figure 8E, which compares the data shown in Figure 8A–D, only a slight displacement of the bistability range toward the anterior or posterior end of the PSM is observed in these conditions, depending on which one of the gradients is decreased. This is consistent with the lack of alteration of the segmentation process observed in heterozygote mutants in which the amount of fgf8 or RA is lowered by half (Meyers et al.,1998; Niederreither et al.,1999). The occurrence of bistability is therefore a robust phenomenon, which persists over a wide range of conditions.

Extending the Model to Incorporate Diffusion of RA and FGF

To incorporate into the model the effect of diffusion of RA and FGF along the PSM, it becomes necessary to distinguish the concentrations of RA and FGF in the intracellular and extracellular medium. The dynamics of the system is then governed by six instead of four kinetic equations, and eqs. [6]–[11] replace eqs. [1]–[4]. The evolution of extracellular RA and FGF is now governed by the partial differential equations [10] and [11]. Besides diffusion, the new system of equations incorporates transport of RA and FGF between the intracellular and extracellular medium, as well as degradation of extracellular RA and FGF (see Model and Simulations Procedures section).

As shown in Figure 9A, when diffusion terms in eqs. [10] and [11] are neglected (i.e., when setting DRA = DF = 0), bistability can occur in the extended model much as in the previous, 4-variable version. The results obtained for bistability in the 4-variable model in Figures 3–8 can thus be recovered in the extended model. The effect of diffusion can be tested by assigning finite values to the diffusion coefficients of RA and FGF, with the constraint that RA diffuses much faster than FGF. Whereas the diffusion coefficient of RA in embryos was estimated to be of 6 × 10−6cm2/min (Eichele and Thaller,1987), FGF8 was indeed shown to diffuse considerably slower in the zebrafish embryo (Scholpp and Brand,2004). The results of numerical integration indicate (Fig. 9B) that the bistability phenomenon is preserved in the presence of RA and FGF diffusion, but occurs only over a portion of the domain of bistability observed in the absence of diffusion (compare Fig. 9A and B). The transitions between the two branches remain very abrupt but are no more of an all-or-none type. Because of the smoothing effect of diffusion, intermediate values between the two branches of stable steady states can indeed be obtained over a narrow range of position x near the boundaries of the bistability domain (Fig. 9B).

Relation to the Determination Front and the Segmentation Clock

As new somites form at the anterior extremity of the PSM, new cells are constantly added at its posterior end as a result of the axis elongation process. Consequently, the relative position of cells in the PSM becomes continuously more anterior (whereas their absolute axial position does not change). Newly formed PSM cells are first exposed to the highest concentrations of FGF8 and, therefore, start on branch 1 (the only branch of steady state accessible at that end) (Figs. 7, 10). These cells then continue to move along branch 1 until they reach the bistability window at x2 (Figs. 7A,B, 10). Once in the bistability domain, cells that reside on branch 1 (high FGF/low RA) are poised to respond to an external stimulus by an all-or-none switch to branch 2 (low FGF/high RA). This switch can occur before the system reaches x1 (the anterior point of the PSM, where the bistability domain ends), in response to a suprathreshold increase in RA signaling or decrease in FGF signaling, or from a combination of both.

Figure 10.

Bistability arising from mutually antagonistic gradients of FGF and RA signaling, and its role in segment determination. Sequential stages leading to segment determination during a time window corresponding to one oscillation of the segmentation clock in the PSM (one somite formation). a,c,e,g: Position of the FGF (green) and RA gradients (purple) in the PSM. The purple and green curves indicate the level of ρ corresponding to the ratio of FGF signaling over RA signaling at steady state in the PSM schematized below (Fig. 7E). In the green part of the curve, FGF dominates over RA (branch 1), while in the purple part, RA dominates over FGF (branch 2). b,d,f,h: Progression of the segment-specification process. Posterior elongation of the embryo leads to the posterior displacement of the bistability window (between the vertical dashed lines in x1 and x2; Fig. 11) but its relative position in the PSM remains constant. As a result of this posterior displacement, new cells from the posterior PSM constantly enter the bistability window, at the level of x2, in green (d, f). This model implies that the size of the bistability window has to be equal to or larger than the segmental domain. It is more likely that the domain is larger rather than fitting exactly the segmental domain and consequently, the segment that was determined during the previous oscillation cycle lies in the bistability window (hatched black, limited posteriorly by x3). The clock signal (vertical arrow) is assumed to be provided synchronously to cells in the bistability window but only cells found between x3 and x2 can respond because cells anterior to x3 have already responded to the clock signal during previous oscillation cycles. This periodic signal triggers the bistable transition from branch 1 (high FGF signaling) to branch 2 (high RA signaling) in cells located between x2 and x3. In response to this transition, cells synchronously activate genes such as Mesp2 in a domain defining the future segment (hatched purple). This also results in defining the future posterior boundary of the newly specified segment, which corresponds to position x3.

Figure 11.

Displacement of the domain of bistability along the PSM as a function of the concentration of fgf8 mRNA. Positions x1 and x2 denote, respectively, the lower and upper limit of the region of bistability (see Fig. 7A,B) along the PSM. When x < x1, a single steady state exists and is characterized by a high level of RA and low level of FGF, while for x > x2, the steady state is characterized by a high level of FGF and a low level of RA. When the concentration M0 of fgf8 mRNA at the posterior boundary of the PSM decreases, the drift in x1 and x2 indicates that the bistability domain is shifted to the posterior region of the PSM. Parameter values are as in Figure 3.

PSM cells on branch 1 experience periodic Notch, FGF, and Wnt activation produced by the segmentation clock (Pourquié,2003; Dequeant et al.,2006). Because only one steady state exists in the posterior PSM, cells in this region are unable to respond to the segmentation clock signal. Only when they reach the bistability window (at x2) would PSM cells become able to respond to the clock signal. In this sense, the bistability window defines a competence window for the segmentation clock signal in the PSM. We propose that the segmentation clock triggers the suprathreshold change in RA or FGF signaling causing the abrupt transition from branch 1 to branch 2. The recently described oscillations of FGF signaling are a good candidate to trigger this abrupt transition (Dequeant et al.,2006).

The synchronization of the segmentation clock readouts (such as cyclic gene expression) within a PSM domain suggests that the signal triggering the bistable transition is synchronously delivered to the cohort of cells that entered the bistability window since the previous clock oscillation. All the cells that entered the bistability window during the previous oscillation cycle have already operated their coordinated transition to branch 2 and now define a segmental domain (Fig. 10, hatched black). Thus, even if these cells can still be located in the bistability window (depending on the parameters defining the size of the window), they cannot respond to the clock signal anymore and therefore stay on branch 2.

We propose that the changes caused by the pulsatory signal from the clock induce the cohort of cells that enters the posterior part of the bistability domain during one oscillation cycle (Fig. 10d,f, green) to move simultaneously to branch 2 (Fig. 10g,h), and, therefore, to collectively escape high FGF signaling and become altogether exposed to high RA signaling. This signal only operates in one direction and can only trigger the high FGF (branch 1) to high RA (branch 2) transition. This clock-induced bistable transition would result in the synchronous activation in this cohort of cells of genes such as mesp2 (Fig. 10h, purple), which are repressed by FGF signaling and activated by RA, thereby specifying the next segmental domain (Fig. 10). As a result, the differentiation program becomes synchronously activated in this domain. The posterior border of this domain subsequently defines the future segment boundary (Morimoto et al.,2005).

The question arises as to how the position of the bistability domain along the PSM depends on the gradient in fgf8 mRNA. In Figure 7A and B, we denoted by x1 and x2 the lower and upper limit of the region of bistability along the PSM. Shown in Figure 11 are the positions x1 and x1 as a function of the concentration M0 of fgf8 mRNA at the posterior boundary of the PSM. When the concentration of fgf8 mRNA at the posterior boundary of the PSM decreases, the drift in x1 and x2 indicates that the bistability domain is shifted to the posterior region of the PSM. Numerical simulations thus show that the domain of bistability moves toward the posterior end of the PSM as a result of the progressive flattening of the fgf8 mRNA gradient (Fig. 11). Therefore, the degradation of fgf8 mRNA in the course of time should contribute to the movement of the determination front toward the posterior PSM during somitogenesis. Because of the constant posterior displacement of the bistability window caused by the posterior movement of the FGF gradient resulting from mRNA degradation (Dubrulle and Pourquié,2004), the PSM domains that undergo the bistable transition at each oscillation cycle will appear as a succession of segmental domains (Fig. 10). The periodic nature of the clock signal will ensure that a bistable transition occurs during each oscillation cycle at the same relative PSM level but in sequentially more posterior domains of the paraxial mesoderm. Thus, in this model, the segmental domain is contained within the bistability window and is defined by the distance traveled by x2 during one oscillation of the segmentation clock.

Bistability and the “Clock and Wavefront” Mechanism Controlling Somitogenesis

To account for the progressive formation of somites along the AP axis, Cooke and Zeeman (1976) assumed that a clock, which was later identified experimentally as the segmentation clock (Palmeirim et al.,1997), interacts with a wavefront (created in their abstract model by the sharp transition between two stable steady states) that moves progressively in an anterior-to-posterior direction. Bistability resulting from mutual inhibition of RA and FGF provides a molecular mechanism for the all-or-none transition taking place at the wavefront level and controlled by the clock cycle assumed in the “clock and wavefront” mechanism (Cooke and Zeeman,1976). The postulated displacement of the wavefront can thus be explained by the constant degradation of fgf8 mRNA and FGF8 protein at the determination front level, which triggers the posterior movement of the bistability window (Fig. 11).

Existence of the switch predicted by the model is supported by the sharp segmental expression of genes, such as those of the mesp2 family, which respond to FGF and RA (Moreno and Kintner,2004; Delfini et al.,2005). Also consistent with the model is the observation that in the mouse embryo, the response of a RARE-LacZ transgene, a sensitive reporter for the presence of RA, shows a sharp threshold at the determination front level, whereas RALDH2 protein is expressed in a gradient in the anterior PSM (Vermot et al.,2005). Experimental proof for bistability in other cellular processes, such as the cell cycle, was achieved by demonstrating hysteresis (Pomerening et al.,2003; Sha et al.,2003): the response does not follow the same path when a control parameter is varied back and forth. Here, this would require showing that the spatial position where cells switch on a response gene such as mesp2 will differ depending on whether they go down the FGF gradient (or up the RA gradient) as they exit the posterior region, or in the reverse direction of the gradients if the cells were to move posteriorly. However, for such a demonstration, one would need to reproducibly manipulate the levels of FGF and RA signaling in embryos with a precision that is not within reach using currently available techniques.

Relation to Other Patterning Processes in Embryogenesis

The peculiarity of the present system is that the mutual inhibition of two well-known developmental regulators, such as RA and FGF, can produce a phenomenon of bistability that is localized in a precise position of the embryo when two opposite gradients in these regulators are established along the embryonic axis. A similar situation is also seen in other developmental systems, such as the limb bud, in which opposing gradients of FGF and RA produced by the apical ectodermal ridge (AER) and by the proximal region, control patterning along the proximo-distal (PD) axis (Capdevila and Izpisua Belmonte,2001; Yashiro et al.,2004). In this case, the opposing gradients could lead to the establishment of a sharp developmental threshold when cells reach the end of the bistability domain. Such a threshold could play a role in coordinating cell differentiation during growth along the PD axis. Antagonistic gradients of FGF and RA signaling have also been described in the hindbrain and might be involved in the segmentation of this territory (Gavalas and Krumlauf,2000). The role of two opposed, unspecified gradients has also been considered theoretically in the context of morphogenetic patterns (McHale et al.,2006).

Bistability is a widespread phenomenon in biological systems (Thomas and d'Ari,1990; Ferrell,2002) in which it can readily arise as a result of positive feedback (Lisman,1985; Guidi et al.,1998; Bhalla et al.,2002; Xiong and Ferrell,2003; Ozbudak et al.,2004). It can also originate from mutual inhibition of two regulatory factors (Monod and Jacob,1961), as shown by mathematical models (Meinhardt,1982; Keller,1995; Cherry and Adler,2000) and demonstrated experimentally in a synthetic genetic system (Gardner et al.,2000). Examples of bistable switches operating in development have been reported. Bistability resulting from self-amplification via positive feedback was proposed to play a role in setting morphogen thresholds in some developmental processes (Lewis et al.,1977; Meinhardt,1982), including dorso-ventral and segmental patterning in Drosophila (Lewis et al.,1977; Meinhardt,1982; Von Dassow and Odell,2002; Ingolia,2004; Wang and Ferguson,2005). Bistable switches resulting from mutual inhibition were also proposed to operate in development (Meinhardt,1982), for example, in stabilization of the pair rule stripes in the fly embryo by mutual repression of transcription factors (Edgar et al.,1989) and in delta-notch-mediated lateral inhibition (Collier et al.,1996). The present results based on the coupling of RA and FGF signaling through mutual inhibition in somitogenesis demonstrate that two opposing morphogen gradients mutually inhibiting each other can lead to bistability defining sharp thresholds. Given that such antagonistic morphogen gradients are a common feature of development (Brook and Cohen,1996; DeRobertis and Sasai,1996; Jiang and Struhl,1996; McHale et al.,2006), this mechanism could play a major role in many morphogenetic processes.

MODEL AND SIMULATION PROCEDURES

Model for Bistability Based on Mutual Inhibition of RA and FGF Signaling

We consider the case in which RA and FGF are linked through mutual inhibition, as occurs during somitogenesis in the PSM. As depicted in Figure 2, the synthesis of RA is catalyzed by the RALDH2 enzyme (which is distributed along a gradient, with maximum activity in the anterior part of the PSM), while RA is hydrolyzed by the enzyme CYP26. The inhibitory effect exerted on RA by FGF arises from the induction of cyp26 expression by FGF. Thus, FGF increases RA destruction by increasing the level of the cyp26 mRNA and, consequently, the level of the RA-hydrolyzing enzyme. The negative regulatory effect of RA on FGF is taken into account by subjecting to inhibition by RA the rate of fgf8 mRNA translation into FGF.

To examine the effect of gradients in RALDH2 and fgf mRNA, we kept these as parameters of the model. We, therefore, considered for simplicity and definiteness that the negative effect exerted by FGF on RA signaling is solely mediated by activation of cyp26 expression. Likewise, we considered that the effect of RA on FGF signaling is mediated solely by inhibition of translation of fgf8 mRNA by RA. Other modes of mutual inhibition of FGF and RA signaling should give largely similar results with respect to bistability.

Kinetic Equations

The variables are the concentrations of RA (RA), cyp26 mRNA (MC), CYP26 protein (C), and FGF protein (F). The time evolution of these variables is governed by the following set of kinetic equations:

(1)

(2)

(3)

(4)

In the kinetic eqs. [1]–[4] we assume that the regulatory effects of RA on FGF and of FGF on RA obey cooperative kinetics, described by Hill functions with cooperativity degrees m and n, respectively. In these Hill functions, KA is the constant measuring half-maximum activation of cyp26 expression by FGF, while KI is the constant measuring half-maximum inhibition of fgf8 mRNA translation into FGF. Parameters νs1 and ks2 measure the rate of synthesis of RA (by RALDH2) and of CYP26, while ks3 is the rate constant measuring FGF synthesis; V0 is a basal, FGF-independent rate of expression of cyp26; VsC is the maximum rate of FGF-activated cyp26 expression; kd1, kd2, kd3,kd4, and kd5 represent apparent first-order constants measuring the rates of destruction of RA (by enzyme CYP26), CYP26, cyp26 mRNA, and FGF, and the rate of nonspecific degradation of RA besides that catalyzed by CYP26, respectively.

In addition to the above equations describing the mutual inhibition of RA and FGF, we consider the following opposite, linear gradients in the rate of synthesis of RA by the RALDH2 enzyme, νs1, and in the amount of fgf8 mRNA (MF):

(5)

The first expression indicates that in any given point of the PSM, the rate of synthesis of RA, νs1, is the product of the catalytic constant ks1 times the concentration of RALDH2. The gradients express the assumption that the concentration of the RA-synthesizing enzyme decreases linearly from a maximum value RALDH20 in the anterior extremity (x = 0) down to zero at the posterior extremity (x = L) of the PSM, while the gradient in fgf8 mRNA decreases from a value of M0 down to zero in the opposite direction (Figs. 1, 7, 9, and 10). Nonlinear gradients, such as exponential gradients expected to result from the diffusion of a morphogen from a localized source, have also been considered (Fig. 8). We introduce a coordinate system traveling with the PSM and assume that its length remains constant because incorporation of new cells at the posterior end is balanced by exit of cells at the anterior end due to differentiation, leading to somite formation. We further assume that the gradient conditions are constantly re-established at the boundaries of the PSM. If we assume that the dynamics of the interactions between RA and FGF signaling remain a cell-autonomous process and do not rely on intercellular diffusion of the regulatory factors, we can follow the time evolution of the system by integrating eqs. [1]–[4] in each point in space along the axis of the PSM, i.e., in the range x = 0 (anterior end) to L (posterior end). For a given value of χ, the values of νs1 (rate of RA synthesis) and MF (amount of fgf8 mRNA) are determined by expressions [5] for the two opposite gradients. The time evolution at a particular location in space is obtained by integrating numerically eqs. [1]–[4], using the values of ks1 and MF corresponding to this point.

Extended Model Incorporating Diffusion of RA and FGF

When investigating the effect of diffusion of RA and FGF, we must extend the model not only by including diffusion terms for RA and FGF, but also by distinguishing between the intracellular and extracellular concentrations of the two chemical species. Indeed, we wish to focus on the effect of diffusion of extracellular RA and FGF along the PSM. Because we now distinguish between the intracellular and extracellular concentrations of RA and FGF, the set of four kinetic eqs. [1]–[4] must now be replaced by the set of six kinetic eqs. [6]–[11]:

(6)

(7)

(8)

(9)

(10)

(11)

where RA and F denote the intracellular concentrations of RA and FGF, while RAe and Fe represent their extracellular concentrations. Moreover, kt1 and kr1 represent apparent first-order rate constants for transport of intracellular RA into the extracellular medium and for the reverse process; kt2 and kr2 are the corresponding parameters for FGF transport between the intracellular and extracellular medium. Parameter θ denotes the ratio of extracellular to intracellular volumes. The rate constants ke1 and ke2 measure the degradation of RA and FGF in the extracellular medium. The diffusion coefficients of RA and FGF are denoted by DRA and DF. Diffusion takes place in the spatial dimension x along the AP axis of the PSM.

Measuring the Effects of RA and FGF

To better compare the effects of RA and FGF, it is useful to consider the binding of these factors to their receptors. In the model governed by eqs. [1]–[4], the saturation functions of the RA and FGF receptors, equal to the amount of receptor bound to ligand divided by the total amount of receptor, are given by eqn. [12] where Kr1 and Kr2 represent, respectively, the dissociation constants for RA and FGF binding to their receptors:

(12)

As a relative measure of FGF and RA signaling response, we use the ratio ρ of FGF to RA receptor occupancy, defined by eqn. [13]:

(13)

This allows us to characterize bistability in terms of a single quantity. A value of ρ smaller (larger) than unity will indicate that RA (FGF) signaling predominates.

In the extended model governed by eqs. [6]–[11], which incorporates diffusion of extracellular RA and FGF, we distinguish the intracellular and extracellular concentrations of RA and FGF. Then ρ remains defined by eq. [13] while the expressions for α1 and α2 become:

(14)

to reflect the fact that the regulatory effect of FGF is exerted upon binding of extracellular FGF to a membrane-bound receptor, in contrast to RA, which exerts its effect intracellularly.

Numerical Study of Bistability

We used several methods to determine the steady states admitted by the kinetic eqs. [1]–[4]. The first method consists of integrating numerically the system of ordinary differential equations starting with different initial conditions. When a single, stable steady state exists, the system evolves to this state regardless of initial conditions. In contrast, when bistability occurs, the system evolves to either one of two stable steady states, depending on initial conditions. The set of initial conditions from which the system evolves to one of the two stable steady states represents the basin of attraction of this state. Using such a method, however, we can determine only the stable steady states and not the unstable ones. It is also possible to solve the algebraic equations obtained at steady state, either numerically or graphically, to determine the steady-state solutions as a function of a control parameter such as KI. This approach must be coupled to a linear stability analysis to distinguish between stable and unstable steady states.

An alternative approach rests on the use of the program AUTO (Doedel,1981). This continuation program finds numerically stable and unstable steady states admitted by a system of ordinary differential equations as a function of a control parameter. We used the program AUTO to generate the branches of stable and unstable steady states in Figures 3–6 as a function of parameters KA and KI, RALDH2 and MF, or kd5, and in Figures 7A–E and 9A as a function of position x. The branches of stable steady states in these figures were also determined by the first method described above.

For the determination of bistability in the presence of diffusion (Fig. 9B), we integrated numerically eqs. [6]–[11]. The diffusion terms in the partial differential eqs. [10] and [11] for RA and FGF were approximated in one spatial dimension by finite differences, and space was represented by a mesh of 100 points. The equations were integrated numerically using zero-flux boundary conditions by means of XPPAUT (http://www.math.pitt.edu/∼bard/xpp/xpp.html) and MATLAB. Bistability was determined by integrating the equations starting from different initial conditions after disappearance of (sometimes prolonged) transients.

Parameter Values

The set of parameter values listed in the legend to Figure 3 is physiologically plausible because the degradation rate constants correspond to half-lives of the order of 2 h for the mRNAs and proteins (Hargrove et al.,1991; Lewin,1997; Lewis,2003) in the scheme of Figure 2. It is not straightforward to obtain data from the literature for the other parameter values. The value chosen for the maximum rate of RA synthesis, νs1, is of the order of 1 nM/min. The RA degradation rate by CYP26, equal to the product kd1C, is of the order of 0.2 min−1 (McSorley and Daly,2000). Taking a value of the order of 2 nM min−1 for the maximum rate of RA synthesis by RALDH2, νs1, is compatible with steady-state levels of RA of the order of 10 nM reported experimentally (Maden et al.,1998). The concentration of FGF8 protein in the mouse embryo posterior PSM was estimated to be in the range of 45 nM (Dubrulle and Pourquié,2004). The diffusion coefficient of RA, DRA, was given the experimentally determined value of 6 × 10−6 cm2/min (Eichele and Thaller,1987).

Origin of Bistability

To understand the origin of bistability as a function of position along the opposite gradients of RALDH2 and fgf8 mRNA, it is useful to focus on a particular position along the PSM and to study the behavior of the cross-regulatory scheme of Figure 2 as a function of the strength of the mutual inhibitions. Strength of the inhibition of RA signaling by FGF is measured by KA, which gives the FGF concentration yielding 50% maximum activation of cyp26 expression, while strength of the inhibition of FGF signaling by RA is measured by KI, which gives the concentration of RA yielding 50% inhibition of fgf8 translation. Shown in Figure 3 is the steady-state behavior of the RA/FGF signaling system as a function of KI (Fig. 3A) and KA (Fig. 3B). In both cases we observe a range of values of KA or KI in which bistability occurs. In this range, two stable steady states are separated by an unstable steady state.

As the strength of FGF inhibition of RA signaling decreases when KA progressively increases, we first find a single branch of steady states in which FGF signaling predominates over RA signaling. At large values of KA, a unique branch of steady states again exists in which RA signaling predominates. For intermediate values of KA, the two branches of steady states coexist (Fig. 3B). For KI we find an inverse picture (Fig. 3A), with bistability again occurring at intermediate values of KI. The data shown in Figure 3 thus indicate that in a particular position in space corresponding to a particular pair of values of RALDH2 and fgf8 mRNA defined by the gradients, bistability may arise at comparable strength of the mutual inhibitions of RA signaling by FGF and of FGF signaling by RA.

Robustness of Bistability

The phenomenon of bistability occurs over a wide range of parameter values, generally in a window as a function of a given control parameter (Figs. 3–6). Bistability is, indeed, intimately associated with the pattern of mutual inhibition that links RA with FGF signaling. To stress this point, we used widely different sets of parameter values to obtain bistability (Figs. 3–5 and 7).

Transition From One Steady State to the Other in the Bistability Domain

We have investigated the type of perturbation that can induce the switch from one stable steady state to the other when bistability occurs. Starting from steady state 2 (low FGF, high RA), the switch to steady state 1 (high FGF, low RA) can be triggered rather easily, by increasing above a threshold the initial concentration of FGF or of CYP26 or its mRNA, but not simply by decreasing the level of RA. Starting from steady state 1 (high FGF, low RA), it is more difficult to induce the more physiologically relevant switch to steady state 2 (low FGF, high RA). A combined perturbation of several variables is needed. For example, bringing the levels of CYP26 and its mRNA down to zero and decreasing the level of FGF below a threshold level results in triggering the switch to the (low FGF, high RA) steady state. For the parameter values listed above, this switch occurs in some 100 min, the time window compatible with the time required for the formation of a pair of new somites.

Acknowledgements

We thank T. Lecuit, V. Nanjundiah, R. Li, T. Erneux, P. Kulesa, E. Ozbudak, V. Francois, and R. Baker as well as members of our research groups for their critical reading of the manuscript. Thanks are also due to anonymous referees for their useful suggestions. The work of A.G. and D.G. was supported by grant 3.4636.04 from the Fonds de la Recherche Scientifique Médicale (F.R.S.M., Belgium), by the European Union through the Network of Excellence BioSim, Contract LSHB-CT-2004-005137, and by the Belgian Programme on Interuniversity Attraction Poles, initiated by the Belgian Federal Science Policy Office, project P6/22 (BIOMAGNET). Work in the Pourquié lab is supported by the Stowers Institute for Medical Research and NIH grant R01HD043158 to O.P. During completion of this work, A. Goldbeter held a Blaise Pascal International Research Chair at the Université de Paris Sud-Orsay, supported by the State and the Ile-de-France Region, under the auspices of the Fondation de l'Ecole Normale Supérieure. O. Pourquié is a Howard Hughes Medical Institute Investigator.