Noise drives sharpening of gene expression boundaries in the zebrafish hindbrain.

Morphogens provide positional information for spatial patterns of gene expression during development. However, stochastic effects such as local fluctuations in morphogen concentration and noise in signal transduction make it difficult for cells to respond to their positions accurately enough to generate sharp boundaries between gene expression domains. During development of rhombomeres in the zebrafish hindbrain, the morphogen retinoic acid (RA) induces expression of hoxb1a in rhombomere 4 (r4) and krox20 in r3 and r5. Fluorescent in situ hybridization reveals rough edges around these gene expression domains, in which cells co-express hoxb1a and krox20 on either side of the boundary, and these sharpen within a few hours. Computational analysis of spatial stochastic models shows, surprisingly, that noise in hoxb1a/krox20 expression actually promotes sharpening of boundaries between adjacent segments. In particular, fluctuations in RA initially induce a rough boundary that requires noise in hoxb1a/krox20 expression to sharpen. This finding suggests a novel noise attenuation mechanism that relies on intracellular noise to induce switching and coordinate cellular decisions during developmental patterning.

Introduction

A fundamental feature of developing systems is that cells sense their positions along morphogen gradients and respond collectively to form precise domains of target gene expression (Meinhardt, 2009; Wartlick et al, 2009). How do gene expression domains achieve such sharp boundaries? Cells near a future boundary experience fluctuations or ‘noise' in: (1) morphogen concentration, due to varying synthesis and transport, (2) ability to respond, for example, due to differences in numbers of receptors, (3) transcription and translation rates of target genes, and (4) feedback (Kepler and Elston, 2001; Elowitz et al, 2002; Kaern et al, 2005). Various mechanisms have been proposed to attenuate these sources of noise to generate consistent gene expression domains in every individual. In spatial patterning systems, noise is generally considered as detrimental to the ultimate goal of the system.

However, for systems without spatial constraints, noise can regulate biological switches between high and low gene expression states, and noise can be attenuated by an ultrasensitive signal (Hasty et al, 2000; Thattai and van Oudenaarden, 2002; To and Maheshri, 2010). Could similar switches be operating in spatial patterning systems? Bistability (distinct steady states of a regulatory gene network within a cell) can have a critical role in spatial patterning and lead to sharp borders between gene expression domains in deterministic models (Meinhardt, 1978, 1982; Ferrell, 2002; Lopes et al, 2008). Spatially constrained stochastic models, such as for segmentation of the Drosophila embryo, suggest that noise predominantly depends on transcription and translation dynamics of target gene expression (Holloway et al, 2011), but external fluctuations in signals also have an important role in these downstream responses (He et al, 2012). However, very few studies have addressed mechanisms of noise attenuation in the formation of gene expression boundaries in any system.

Here, we investigate interactions between noise in a morphogen (i.e., retinoic acid—RA) and noise in its downstream, bistable regulatory gene network in boundary sharpening. RA specifies rough boundaries between segments (called rhombomeres) of the zebrafish hindbrain in a concentration-dependent manner, which subsequently become razor sharp (Giudicelli et al, 2001; Cooke and Moens, 2002; White et al, 2007; White and Schilling, 2008). Two genes downstream of RA, hoxb1a (r4) and krox20 (r3 and r5), cross-inhibit one another and auto-activate their own expression to form a bistable switch (Barrow et al, 2000; Giudicelli et al, 2001; Alexander et al, 2009). With a stochastic model that incorporates these interactions we estimate the switching probability between hoxb1a and krox20 expression at different RA concentrations based on an exponential function of Minimum Action Paths (MAPs) between stable and unstable states (Freidlin and Wentzell, 1998). Exploration of the stochastic models reveals that noise in the RA morphogen gradient can lead to rough gene expression boundaries initially, and that sharpening is driven by noise in the expression of hoxb1a and krox20, due to induced switching between expression of one gene and the other. These results reveal an unexpected positive role for noise in boundary sharpening that may be common for many patterning systems.

Results

hoxb1a and krox20 co-expression during rhombomere boundary sharpening

To determine the temporal dynamics of hoxb1a and krox20 expression in the embryonic zebrafish hindbrain, we performed fluorescent in situ hybridization (FISH) analysis. Previous studies showed that initial boundaries of hoxb1a in r4 and krox20 expression in r3 and r5 are rough but become razor sharp between 10 and 14  h post fertilization (h.p.f.) (Figure 1A–F; Cooke and Moens, 2002; Cooke et al, 2005). Cells that find themselves on the wrong side of a boundary (i.e., surrounded by neighbors with a different pattern of gene expression) may go through a transient phase in which they express both genes and subsequently downregulate one or the other to enable sharpening (Schilling et al, 2001; Cooke and Moens, 2002). To quantify sharpness in krox20 expression, confocal stacks were collected for a minimum of 10 embryos at 6 different stages (between 10.7 and 12.7 h.p.f.) (Figure 1A–F) and the fluorescence was measured at different positions along the anterior-posterior (A-P) axis focusing on the r4/5 boundary (Figure 1G–I). This analysis demonstrated quantitatively how krox20 expression sharpens at rhombomere boundaries over time.

Figure 1

Sharpening of gene expression boundaries in the zebrafish hindbrain. (A–F) Single confocal images of fluorescent in situ hybridization (FISH) for krox20 (red) mRNA, dorsal views, anterior to the left, between 10.7 and 12.7 h post fertilization (h.p.f.). (G–I) Fluorescence measurements at different positions along the anterior-posterior axis (X axis) at 11, 11.7, and 12.7 h.p.f. Lines represent four different samples. (J–L) Single confocal images of two-color FISH for hoxb1a (r4, red) and krox20 (r3 and r5, green). Insets show enlargements of cells co-expressing both (yellow). (M–O) Sample distributions of mis-expressing cells along the r4/5 boundary (black lines) between 10.7 and 12 h.p.f., anterior to the top. Cells mis-expressing krox20—green dots, hoxb1a—red dots and co-expressing cells—orange dots.

To examine this more closely, we used confocal analysis and two-color FISH to colocalize hoxb1a and krox20 near the r3/4 and r4/5 boundaries at 20-min intervals between 10.6 and 12 h.p.f. (Figure 1J–L). hoxb1a expression is initiated broadly in the early gastrula (6.5 h.p.f.; Maves and Kimmel, 2005), and is preceded by its close relative hoxb1b, which is the first gene induced by RA in this system and activates hoxb1a transcription directly (McClintock et al, 2001). By 10.5 h.p.f. hoxb1a expression resolves into a strong r4 stripe 4–6 cells wide along the A-P axis while krox20 is expressed in flanking r3 and r5 stripes that overlap with hoxb1a at its edges (Figure 1J–L). Higher magnification images demonstrated that krox20 and hoxb1a mRNAs colocalize in many of these cells near future boundaries (insets) and occasional colocalization was observed as late as 12.0 h.p.f. This revealed an initial ‘transition zone' containing a mixture of hoxb1a, krox20 and co-expressing cells that was ∼40 μm in length along the A-P axis and later reduced to 5–10 μm (1 cell diameter) by 12 h.p.f. Similar numbers of co-expressing cells were identified at 10.7 h.p.f. (average 7 cells, n=3) and 11.3 h.p.f. (average 7.3 cells, n=3), however, by 12 h.p.f. the number of co-expressing cells decreased (average 3, n=3). Conversely, the percentage of mis-specified cells that were expressing both genes increased from 36% and 34% at 10.7 and 11.3 h.p.f. to 56% at 12 h.p.f. The co-expressing cells were more prevalent in the r5 domain at both 10.7 and 11.3 h.p.f. (Figure 1M and N) while this bias was not observed at 12 h.p.f. (Figure 1O).

Induction of stripes of hoxb1a and krox20 expression by RA requires bistability and initial expression of Hoxb1

RA activates hoxb1a expression in r4 (directly) and krox20 in r3 and r5 (indirectly through Vhnf1 and MafB) in a concentration-dependent manner (Niederreither et al, 2000; Begemann et al, 2001; Hernandez et al, 2004; Labalette et al, 2011). Our deterministic model is based on a previous continuum model of the RA signaling network that consists of diffusive extracellular and intracellular RA, and self-enhanced degradation through the enzyme Cyp26a1 (White et al, 2007), without inclusion of downstream signal responses (see Equation S1.1 in Supplementary information). In the new model, RA activates hoxb1a and krox20 expression, which in turn both positively regulate their own expression and negatively regulate each other (Barrow et al, 2000; Giudicelli et al, 2001; Alexander et al, 2009; Figure 2A). Such positive auto-regulation and mutual inhibition have been modeled and shown to result in only one gene remaining active in a particular cell (Meinhardt, 1978, 1982). Here, the dynamics of both genes are modeled using rate equations along with Hill functions for regulation, with RA as input (see Equation S1.2 in Supplementary information).

Figure 2

Modeling induction of hoxb1a and krox20 expression by a gradient of retinoic acid (RA) in a noise-free system. (A) Diagram illustrating RA movement from extracellular [RA]out to intracellular [RA]in, self-enhanced degradation via Cyp26a1, and induction of hoxb1a and krox20 which undergo auto-activation and cross-inhibition. (B) In the absence of noise, a smooth RA gradient leads to sharp boundaries of gene expression—as long as there is a low initial level of hoxb1a (∼0.1). (C) Three-dimensional graph of krox20 (g, Y axis) and hoxb1a (g, Z axis) expression levels at different points along the RA gradient (X axis). The number of possible gene states is 5 (0<[RA]in<0.22), 3 (0.22<[RA]in<0.85), and 1 ([RA]in>0.85) for a normalized RA concentration. (D) Phase diagram of Hoxb1 (red) and Krox20 (blue) activation illustrating effects of the initial level of Hoxb1 (Y axis) at different segmental positions (X axis). The initial level of krox20 is zero and the RA gradient used to generate the diagram is shown in (B).

Exploration of the model reveals that the system robustly resolves into a striped pattern of gene expression with hoxb1a in r4 and krox20 in r3 and r5 (Figure 2B). This demonstrates that by simply including two bistable steady states and an anteriorly declining RA gradient, one can specify alternating gene expression patterns with sharp boundaries. Simulations in two dimensions show a similar striped pattern (Supplementary Figure S1). Auto-activation and mutual inhibition between Hoxb1 and Krox20 allow one to switch from the off to the on state, or vice versa, within a range of RA. In particular, at a low RA concentration (RA<0.22 μM), there are three stable states (hoxb1a-on, krox20-on, or both off) and two unstable critical transition states (Figure 2C). As the RA concentration increases above 0.22 μM, both the off state and one unstable transition state merge and disappear, while other states (hoxb1a-on, krox20-on, and another unstable transition state) remain (Figure 2C). If the RA concentration is high (larger than 0.85 μM in Figure 2C), then the hoxb1a-on and unstable transition states disappear and only krox20 is activated.

Because of the monotonic spatial distribution of RA (and additional influences from Fgf signaling) (Hernandez et al, 2004; White et al, 2007; Labalette et al, 2011), the activation (and auto-activation) of krox20 is likely stronger than hoxb1a, at least in r5. However, our models suggest that enhanced activation alone cannot create this bistability. Rather, the mutual inhibition between hoxb1a and krox20 is essential for generating complementary stripes of gene expression. Lack of either inhibitory interaction eliminates rhombomere-specific gene expression (see Supplementary Figure S2).

Our model suggests that an initial low level of hoxb1a expression is necessary to enable formation of an alternating pattern of hoxb1a and krox20 expression (Figure 2D). krox20 expression can be activated in the absence of initial hoxb1a due to its strong auto-activation. However, as the hoxb1a concentration increases (>0.084 in Figure 2D), it represses krox20 expression and activates its own expression in r4. Interestingly, hoxb1a is first expressed at an intermediate level of RA, not in the anterior hindbrain where RA concentrations are low. This reflects the presence of the initial low level of hoxb1a expression, which compensates for its weaker auto-activation than krox20. This is consistent with previous theoretical findings (Meinhardt, 1978, 1982) and experimental observations of extremely early onset of hoxb1a (and hoxb1b) expression during zebrafish gastrulation (McClintock et al, 2001; Maves and Kimmel, 2005).

Noise in RA can generate rough boundaries of hoxb1a and krox20 expression

Stochastic fluctuations in ligand-receptor binding and morphogen synthesis introduce noise into any morphogen gradient. To study the propagation of such noise and its influence on target gene expression in the hindbrain, we introduced spatial and temporal noise into the deterministic RA model:

where [RA]out and [RA]in represent extracellular and intracellular RA concentrations, and , denote standard white noise in extracellular and intracellular RA concentrations, respectively. [Cyp([RA]in)] represents RA degradation by Cyp26 (see Section 1 in Supplementary information for the other parameters and boundary conditions).

Because of inherent stochasticity in gene expression and other cellular components (Elowitz et al, 2002) in any gene regulatory network, we also include temporal noise in each gene equation:

Here, g and g represent the concentrations of hoxb1a and krox20, respectively, and ω(t)ψ(t) represents white noise with amplitudes a and a, respectively (see Section 1 in Supplementary information for more details).

To investigate the effects of noise on boundary sharpening we have varied each term in Equations (1) and (2) and modeled the outcome. First, if noise is only present in extracellular RA (i.e., ), due to fluctuations in environmental factors (e.g., availability of vitamin A) or in RA synthesis, then the boundaries between Krox20 and Hoxb1 expression domains are sharp from the outset (see Supplementary Figure S3). This is because RA induces Cyp26a1, creating ‘self-enhanced degradation' feedback that makes signaling robust to fluctuations in RA (Eldar et al, 2003; White et al, 2007) resulting in a smooth gradient of intracellular RA concentration along the A-P axis. Consistent with this idea, simulations in which we have varied spatial noise demonstrate that self-enhanced degradation provides excellent noise attenuation for fluctuations in extracellular RA (see Supplementary Figures S3 and S4).

In contrast, if noise is introduced into the intracellular RA concentration, [RA] (i.e., ), for example, due to fluctuations in RA transport into cells, then boundaries between hoxb1a and krox20 never sharpen (Figure 3A). In this case, the r3 domain of krox20 expression expands as fluctuations in the RA gradient reach the threshold that induces krox20. Monte Carlo simulations indicate that there is a large variation in the distribution of gene expression around the r4/5 boundary over time when [RA]in is noisy (Figure 3B). Two-dimensional simulations also show that hoxb1a and krox20 expression domains initially form a rough r4/5 boundary, which does not sharpen (Figure 3C). An initial noisy distribution of hoxb1a expression at this boundary can also disrupt sharpening (see Supplementary Figure S5). However, if noise is only restricted to later hoxb1a and krox20 expression, and not local RA concentration (i.e., ), then boundaries tend to sharpen from the outset (Figure 3D and F) and Monte Carlo simulations confirm this prediction (Figure 3E). These results suggest that rough boundaries of gene expression between r4 and r5 arise due to noise in [RA]in or initial hoxb1a expression.

Figure 3

Effects of noise either in the RA gradient or in hoxb1a/krox20 expression alone on boundary sharpening. (A–C) With noise in RA alone, boundaries are initially rough and never sharpen. (D–F) With noise in hoxb1a/krox20 expression alone boundaries start out sharp at the outset and remain sharp. (A, D) Single samples at three time points illustrating gene expression levels (Y axis) at different A-P positions in r3-5 (X axis). (B, E) Gene expression distributions (Y axis) at different positions relative to the r4/5 boundary (X axis). Solutions are at the scaled time T=50, which is typically long enough for simulations to reach steady state (1000 samples are taken to calculate the gene distributions). (C, F) 2D simulations at three time points showing the pattern of hoxb1a/krox20 gene expression around the r4/5 boundary (hoxb1a: blue; krox20: red).

Noise in Hoxb1/Krox20 expression enables noise attenuation during boundary sharpening

Rhombomeres form lineally-restricted compartments (Fraser et al, 1990; Jimenez-Guri et al, 2010) and single hindbrain cells can upregulate or downregulate their Hox expression according to their host environment/rhombomere (Trainor and Krumlauf, 2000; Schilling et al, 2001). This suggests that similar gene expression ‘switches' occur in cells on either side of a noisy rhombomere boundary. For example, cells expressing krox20 isolated among neighbors expressing hoxb1a (Figure 1J–L) may downregulate the former and upregulate the latter, thereby attenuating the noise and sharpening the border. To study such switching from one stable gene expression state to another, we employed an MAP analysis based on the Wentzell–Freidlin theory of large deviation (Freidlin and Wentzell, 1998). This theory allows one to estimate the probability of a transition ϕ between two stable states X1, X2 in a stochastic dynamic system (e.g., with a form of Equation (2)). The most probable path ϕ* requires the least action and is called an MAP (see Section 2 in Supplementary information for more details). MAP analysis has been used primarily to model phase transitions between two states in stochastic chemical kinetics (E et al, 2004). Here, we adapt it to estimate the switching probability between two stable gene expression states.

The likelihood that a system switches from X1 to X2 relies on its ability to pass through the unstable critical point Xc that lies between X1 and X2 along the path ϕ*. The distance |ϕ*(X1)−ϕ*(X)| is the minimum barrier to the stochastic transition from one state to the other. For a smooth RA gradient and a simple bistable gene expression state (Figure 2C), we calculate MAPs (E et al, 2004) at different levels of RA. At low RA concentration (e.g., at RA=0.1 μM), three MAPs connect each pair of stable states, with each MAP passing through one unstable critical point (Figure 4A). Based on MAP theory, the activation of Krox20 (Krox20-on) from a ‘both-off' state requires less action (a lower barrier) than activation of Hoxb1, which helps explain why the r3 domain of Krox20 expression expands when noise increases in our models (Figure 3C). In contrast, at intermediate RA concentrations a single MAP connects the two steady states and the action to switching from Hoxb1-on to Krox20-on decreases from RA=0.5 to 0.8 μM (Figure 4A), indicating that it becomes easier to switch in this direction as RA increases. When RA levels are high, Krox20-on is the only stable state.

Figure 4

Noise in hoxb1a/krox20 expression leads to boundary sharpening. (A) Minimum Action Paths (dash lines) at [RA]in=0.1, 0.5, and 0.8 (krox20-on: blue dot, hoxb1a-on: red dot, both-off: black dot, critical point: green dot). (B) Gene switching probability estimated by MAPs reveals that noise in gene expression can drive cells from co-expressing Hoxb1/Krox20 to uniform Krox20 expression when [RA]in is high, and this coincides with the results of Monte Carlo simulations. (C–E) With noise in both [RA]in and hoxb1a/krox20 expression, a transient noisy boundary becomes sharp over time: (C) single sample; (D) gene distribution at the r4/5 boundary (1000 samples are taken to calculate the gene distributions); (E) two-dimensional simulation at the r4/5 boundary (hoxb1a: blue; krox20: red). (F) Sharpness Index versus time. ‘green dashed line': noise only in extracellular RA alone; ‘black dashed-dotted line': noise in both extracellular and intracellular RA; ‘magenta dotted line': noise in gene expression alone; ‘blue solid line': noise interactions between RA and gene expression; ‘red dashed line with green squares': mean value of the Sharpness Index for distributions of Krox20 obtained from the experimental data. The error bar represents the standard error of the mean. The times 11, 11.3, 11.7, 12, and 12.7 h.p.f. correspond 3, 4, 5, 6, and 8 somites, respectively, and are rescaled to 1, 9, 20, 29, and 50.

To quantify such switching capability, we estimate the switching probability from X1 to X2 within a time interval [0, T] through an exponential of the minimal barrier:

The switching probability from X2 to X1 is defined in a similar manner:

We estimate the Hoxb1/Krox20 gene switching probabilities P and P using the MAP calculation of Equation 2 for a normalized RA concentration. Our models indicate that P increases exponentially when [RA] is high and P is low, and cells have a high probability of switching from Hoxb1 to Krox20 expression. On the other hand, Krox20 expression is more stable due to a cell's low switching (to Hoxb1 expression) probability (Figure 4B). Together, this suggests that noise in Hoxb1/Krox20 expression drives cells from occasionally co-expressing Hoxb1/Krox20 expression to a uniform Krox20 expression when RA concentrations are high, leading to a sharp boundary. In support of this analysis, direct Monte Carlo simulations of the gene system (2) of switching probability at the same time intervals provide similar MAP estimates based on Equations (3) and (4) (Figure 4B; see Section 2 in Supplementary information for more details).

Thus, surprisingly, our models suggest that the combination of noise in both [RA]in and Krox20/Hoxb1 expression (i.e., ), synergize to reduce noise during boundary sharpening, at least at the r4/5 boundary (Figure 4C–E). Interestingly, the initial boundary (T=1) is established at 160±10 μm along the A–P axis and, following sharpening, the boundary is located at 144 μm (T=50) (Figure 4E). This suggests that sharpening preferentially drives cells near the initial, rough boundary to krox20 expression due to the irreversibility of gene switching. Similar directional boundary shifts in gene expression have also been observed in Drosophila (Jaeger et al, 2004). This fits well with our in vivo observation of a higher percentage of hoxb1a/krox20 co-expressing cells on the posterior side of the putative r4/5 boundary at 10.7 and 11.3 h.p.f. (Figure 1M and N), which might predict that the forming boundary shifts anteriorly.

To quantify boundary sharpening more systematically, we define a ‘Sharpness Index' (S), which resembles the standard deviation. To define S, we calculate the ‘mean' location of the boundary between Hoxb1 and Krox20 expression domains as the intersection of their distributions at 50% of the normalized value. Using this definition, we can measure the roughness of the boundary, that is, deviation from a sharp boundary (without transition zone) in both one-dimensional and two-dimensional simulations. A decrease in S over time indicates noise attenuation and a sharper boundary while a larger S implies a rougher boundary. Stochastic simulations of four cases that differ in whether or not they include noise in: (a) [RA]out, (b) both [RA]out and [RA]in, (c) Hoxb1/Krox20 expression, and (d) both [RA]in and Hoxb1/Krox20 expression, show that S starts with a large value and decreases significantly over time only when there is both noise in [RA]in and noise in gene expression, which is consistent with S value estimates based on the experimental data for krox20 expression (Figures 1A–F and 4F).

If boundary sharpening depends on gene switching induced by noise in gene expression, then the timing of switching, which is closely related to noise frequency, is likely to be critical. To test this, we have varied noise frequency, both in RA levels and in gene expression, to study effects on boundary sharpening. γ is defined as the ratio of frequency of noise in RA levels over that of gene expression. For a small γ, indicating high frequency noise in gene expression, rough boundaries remain rough (Figure 5A). In this case, gene expression oscillates relatively too fast to reach a critical level of gene switching to enable sharpening. Lower frequency noise in gene expression improves boundary sharpening (Figure 5B). Conversely, high frequency noise in RA (i.e., a larger γ) leads to better noise attenuation (Figure 5C) since the time average of RA signal over a substantial period reduces the influence of noise. As a result, the Sharpness Index, S, decreases as the frequency ratio γ increases, indicating that lower frequency noise in gene expression and higher frequency noise in RA together facilitate boundary sharpening (Figure 5D).

Figure 5

A larger noise frequency ratio improves boundary sharpening. (A–C) Two-dimensional simulations of hoxb1a/krox20 gene expression at T=50 for the ratio of the frequency of noise in RA over the frequency of noise in gene expression at three different values: (A) γ=0.01; (B) γ=1; (C) γ=100 (hoxb1a: blue; krox20: red). (D) Sharpness Index versus frequency ratio.

The noise amplitude, corresponding to the size of a and a also must be within an appropriate range for noise attenuation and boundary sharpening (Supplementary Figure S6A). When a and a are relatively small, switching between expression of Hoxb1 and Krox20 rarely occurs and oscillations in gene expression caused by noise in [RA]in are barely affected by noise in gene regulatory networks (Supplementary Figure S6B). Conversely, when a and a are relatively large, cells do not commit to expressing one gene or the other, leading to oscillations that make the boundary even less sharp (Supplementary Figure S6C) or isolated krox20 expressing cells within the hoxb1a expression domain (Supplementary Figure S13).

Because tissue growth potentially changes the RA gradient and affects the pattern, we also explore a two-dimensional stochastic model incorporating growth along the A-P axis. As measured in our fluorescent in situs (Figure 1), the length of the r3-r5 field increases from 106±4 μm at 10.7 h.p.f. to 127±11 μm at 12.7 h.p.f., corresponding to approximately a 20% increase in A-P length. When incorporated into the stochastic PDE model (see Section 1.3 in Supplementary information), simulations reveal that as the RA gradient fades due to growth (Figure 6A) sharpening of the r4/5 boundary is delayed and less accurate (Figure 6B). This reflects dependence of gene switching on RA concentration—lower RA levels reduce the switching probability for cells around the future boundary. Furthermore, without noise in the gene expression, a fading RA due to growth becomes unable to induce sharp boundaries of gene expression, indicating that this noise is essential (Supplementary Figure S9).

Figure 6

Growth of the embryo delays boundary sharpening. (A) Mean values of the intracellular RA gradient along the X axis at different times, showing a growth-dependent fading of RA levels (see one sample of two-dimensional RA gradients in Supplementary Figure S9A). (B) The pattern of hoxb1a/krox20 gene expression during growth (hoxb1a: blue; krox20: red).

In addition to the stochastic continuum model studied above, we also use a discrete cell model containing a one-dimensional array of cells within which the reactions and regulatory interactions (Figure 2A) are calculated using a Stochastic Simulation Algorithm (Gillespie, 1977). The number of RA molecules is obtained either by the continuum approach shown above or through a stochastic reaction-diffusion process (Kang et al, 2012; see Section 3 in Supplementary information for more details). We found that when the number of RA molecules is large, noise in gene expression can drive boundary sharpening (Supplementary Figures S7 and S8), while when the number of RA molecules is very small, leading to large fluctuations in the RA gradient, the rough boundary remains. Both results are consistent with the stochastic continuum approach.

Time delays, such as those that occur during transcriptional regulation, often affect the dynamics of gene switching and sharpening (Kepler and Elston, 2001; Bratsun et al, 2005). By incorporating a constant time delay in the expression of hoxb1a and krox20 in our model (see Section 1.4 in Supplementary information; Kuang, 1993), we find additional stochastic oscillations that reduce the speed and efficiency of boundary sharpening (Supplementary Figure S10). One possible explanation is that the time delay averages (or filters) noise such that gene switching becomes difficult, particularly when the time delay is long relative to the dynamics of gene expression.

Boundaries may also sharpen through cell movements and differential adhesion that can lead to cell sorting (Xu et al, 1999; Firtel and Meili, 2000; Dormann and Weijer, 2003, 2006). Previous studies in the zebrafish hindbrain have demonstrated some limited sorting at rhombomere boundaries (Cooke et al, 2005; Kemp et al, 2009). Simulations with two-dimensional stochastic discrete cell models, which incorporate both directional cell movements and stochastic reaction diffusion (see Section 4 in Supplementary information), suggest that the speed of cell sorting must be strongly regulated to facilitate sharpening (Supplementary Figure S11). Moreover, directed cell movements without noise in gene expression lead to variability in the location of the r4/5 boundary (Supplementary Figure S12).

Discussion

Morphogen gradients activate target genes in a concentration-dependent manner to generate distinct spatial domains of expression in developing tissues. Extrinsic noise in morphogen concentration and tissue geometry, as well as intrinsic noise in signal transduction and gene expression, reduces precision of patterning (Bollenbach et al, 2008; Balazsi et al, 2011; Kang et al, 2012). Here we show, surprisingly, that noise in intracellular signal transduction actually improves precision and robustness of patterning. In particular, we find that the initially rough expression boundaries of hoxb1a and krox20 in the developing zebrafish hindbrain are due to extrinsic noise in the morphogen (RA) that induces them, but that intracellular fluctuations in their expression lead to gene switching that enables sharpening. The resulting noise attenuation progressively narrows the transition zone between the two gene expression domains, similar to the sharpening that occurs in vivo. Within this transition zone cells transiently co-express both genes before adopting one segmental fate or the other. This result is consistent with experimental evidence showing that cells at these stages can upregulate or downregulate Hox expression, rather than having to migrate to sort themselves out at the future boundary (Trainor and Krumlauf, 2000; Schilling et al, 2001) and that rhombomeres form lineally restricted compartments at early embryonic stages (Fraser et al, 1990; Jimenez-Guri et al, 2010). Similar rules may also apply for other RA target genes (e.g., Vhnf1) and other gene expression boundaries where cross-inhibition and/or auto-activation between target genes occur in morphogen systems (Rivera-Pomar and Jackle, 1996; Perkins et al, 2006; Balaskas et al, 2012). Notably, this property of the system does not depend on the genes being direct transcriptional targets of the morphogen signal—Hoxb1 is a direct target while Krox20 is induced indirectly through Vhnf1, Mafb and other transcription factors (Barrow et al, 2000; Giudicelli et al, 2001; Hernandez et al, 2004; Alexander et al, 2009).

Our deterministic model shows that an initial low level of Hoxb1 is required to generate the alternating striped pattern of Hoxb1 and Krox20 expression in response to a spatially monotonic RA gradient (see Figure 2). During gastrulation in zebrafish, hoxb1a expression is initiated several hours before krox20, and induced by the even earlier expression of hoxb1b in response to RA (McClintock et al, 2001; Maves and Kimmel, 2005). This early onset is important genetically (it fine-tunes target genes within its expression domain; Labalette et al, 2011) and computationally (i.e., pre-steady-state decoding) during embryonic patterning (Bergmann et al, 2007; Saunders and Howard, 2009). This helps explain how alternating, mutually exclusive domains of target gene expression can be induced by the same signal, a common feature of many boundary-forming morphogen systems (Lander, 2011). Another major assumption of the model is the irreversibility of cell switching such that gene expression remains stable when the morphogen gradient decreases or disappears (Gould et al, 1998; Grapin-Botton et al, 1998). This irreversibility has also been pointed out in previous modeling studies (Meinhardt, 1978, 1982).

Intracellular noise may arise from two sources: an intrinsic one due to small numbers of molecules or stochastic fluctuations inherent in biochemical reactions and an extrinsic one driven by fluctuations in cellular environment (Swain et al, 2002). Gene expression noise in some morphogen systems depends predominantly on fluctuations in transcription and translation (Holloway et al, 2011). From our combination of spatial SSA simulations, which take into consideration the number of molecules, and stochastic continuum PDE models, both types of intracellular noise can be utilized to induce switching between gene expression states (e.g., the switching from Hoxb1 to Krox20 is dominant in the transition region between r4 and r5), leading to boundary sharpening. Of course, the level of noise in the morphogen signal must be within a range that allows switching in the transition region through this system of intracellular noise.

Our models also suggest that the cells in the transition region near the r3/r4 boundary utilize somewhat different noise attenuation mechanisms despite undergoing similar boundary sharpening. Unlike the case at higher levels of RA (i.e., at the r4/r5 boundary) where the switching probability from Hoxb1 to Krox20 is significantly higher than the one from Krox20 to Hoxb1 (Figure 4B), near the r3/r4 boundary where RA levels are lower, switching probabilities are also lower. Since we can detect cells co-expressing hoxb1a and krox20 at the r3/4 boundary during sharpening (Figure 1J–L), it seems likely that differences lie upstream, at the level of the inductive signals. Cells at this boundary may have intrinsic differences in their RA responses. Alternatively, these cells integrate responses to additional morphogens such as Fgfs (e.g., Fgf3 and Fgf8), which are produced in the anterior hindbrain during sharpening and influence both RA degradation and expression of hoxb1a and krox20 (Hernandez et al, 2004; Labalette et al, 2011).

In particular, our simulations suggest that both time delays and cell movements can affect rhombomere boundary sharpening. We find that time delay in the expression of hoxb1a and krox20 may introduce additional stochastic effects, leading to reduced speed and efficiency in boundary sharpening. Likewise, cell movements and differential adhesion can lead to cell sorting (Xu et al, 1999; Firtel and Meili, 2000; Dormann and Weijer, 2003, 2006) and some limited sorting has been observed at rhombomere boundaries in the hindbrain (Fraser et al, 1990; Cooke et al, 2005; Kemp et al, 2009). However, clones derived from single progenitors in the neural plate are for the most part lineally restricted to individual rhombomeres (Fraser et al, 1990; Jimenez-Guri et al, 2010) and do not move across boundaries. Hindbrain cells at these stages are also capable of upregulating or downregulating Hox expression if they find themselves on the wrong side of a boundary (Trainor and Krumlauf, 2000; Schilling et al, 2001). Our simulations also reveal that: (1) the speed of cell sorting is critical for sharpening and (2) directed cell movements without noise in gene expression disrupt the location of the r4/5 boundary. Our simulation data suggest that cell movements alone are unlikely to account for rhombomere boundary sharpening, and that noise in gene expression is critical for this process. However, more comprehensive experiments are needed to quantify the amount of sorting that occurs and modeling to understand the roles of this sorting in the establishment and maintenance of precise boundaries.

Our model is limited to two spatial dimensions. Interestingly, in other systems the precision of responses to morphogen gradients rapidly increases when considered in three dimensions (Bollenbach et al, 2008). The zebrafish hindbrain is several cell diameters thick during the sharpening period considered here. Therefore, incorporating more accurate tissue geometries in our stochastic model will undoubtedly reveal new features of the system—most likely including more robust patterning and boundary sharpening that can resist larger amplitude fluctuations. There exist other potential mechanisms that may facilitate boundary sharpening and noise attenuation, including cell-to-cell communication (e.g., Notch signaling) (Louvi and Artavanis-Tsakonas, 2006; Ozbudak and Lewis, 2008; Koseska et al, 2009; VanHook, 2011) and averaging (Tanouchi et al, 2008).

While molecular noise often introduces fluctuations in biochemical reactions and increases uncertainty in cellular decisions, it also can improve performance objectives of biological systems in surprising ways. For example, noise can: (1) help create synchronous oscillations in cell–cell signaling systems (Zhou et al, 2005; Springer and Paulsson, 2006); (2) enhance sensitivity in intracellular regulation (e.g., stochastic focusing; Paulsson et al, 2000), and (3) through reversible progression, help reliable cellular decision making (Kuchina et al, 2011). Our results add to the growing body of evidence that points to important roles for molecular noise in cell fate decisions, and reveals a novel mechanism by which intracellular molecular noise reduces uncertainty in the ability of a fluctuating morphogen to induce precise domains of target genes.

Materials and methods

Numerical methods of stochastic PDE model

The stochastic PDEs (1) were solved with a finite-difference approximation in both one- and two-dimensional spaces. The stochastic system (2) was solved with Milstein's method (Higham, 2001). Noise was added at each grid point in the space at a specified time interval. We used spatial resolutions of 100 grid points in the one-dimensional model or 100 × 20 grid points in the two-dimensional model and temporal resolution of h=0.1s. Numerical tests were conducted to ensure sufficient resolution. The solutions were typically observed at the scale time T=1.25 and 50, which is the steady state for all of the variables at each spatial point to reach an approximately invariant distribution.

To estimate the stationary distribution from one realization, we performed Monte Carlo simulations, which are repeated computations of the stochastic models. We took 1000 samplings in the Monte Carlo simulations to calculate the Sharpness Index. We explored a range of sample numbers from 50 to 5000 to ensure consistent results. Please see Section 1 in Supplementary information for more details and Dataset 2 in Supplementary information for the simulation codes.

Numerical methods for finding MAPs

We followed the algorithm in E et al (2004) to find the MAPs for the gene switch. First, a discrete time interval was used to form a mesh and the path ϕ is approximated on the mesh. Next, the action S[ϕ] of this path was approximated according to the midpoint rule. The steepest descent method was then applied to minimize the discrete action S[ϕ]. Please see Section 2 in Supplementary information for the equations and Dataset 2 in Supplementary information for the simulation codes.

Numerical methods for spatial SSA

We partitioned the one-dimensional or the two-dimensional space into identical compartments based on a computational strategy previously developed (Kang et al, 2012), and applied the Gillespie algorithm to the stochastic reaction-diffusion simulations (Gillespie, 1976). Please see Sections 3 and 4 in Supplementary information for more details and Dataset 2 in Supplementary information for the simulation codes.

Gene expression analysis in zebrafish

Wild-type zebrafish embryos (TL) were collected from natural matings and staged as previously described (Kimmel et al, 1995). Fluorescent whole mount in situ hybridization was performed as previously described for Fast Red alone (Thisse et al, 2004) or with two colors using tyramide amplification (Zuniga et al, 2010). Probes were synthesized from cDNA clones of krox20 (Oxtoby and Jowett, 1993) and hoxb1a (McClintock et al, 2001). Embryos were imaged on an Olympus Fluoview FV1000 confocal microscope, processed in Image J, and post-processed on Adobe Photoshop CS3. Fast Red in situ images for krox20 were processed and analyzed by the Matlab Image Processing Toolbox (The MathWorks, Natick, MA, USA). Graphs of cell location were made in Microsoft Excel based on confocal stacks analyzed in Image J. Please see Dataset 1 in Supplementary information for the processed cell location data. Full image data are available from the authors upon request.