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.
Fluctuations in morphogens often lead to gene expression domains with rough boundaries during development. A mechanism is described whereby intracellular noise can help coordinate cellular decisions during patterning and thereby sharpen expression boundaries.
Fluorescent in situ hybridization reveals rough edges around gene expression domains during development of rhombomeres in the zebrafish hindbrain, in which cells co‐express hoxb1a and krox20 on either side of a future boundary, and these domains sharpen within a few hours.
Computational analysis of stochastic models demonstrates that the initial rough boundaries are caused by stochastic fluctuations in the morphogen, retinoic acid, which patterns the expression domains.
Noise in gene expression induces switching between expression of one gene or the other to narrow the transition zone and enable sharpening of boundaries.
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.
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.
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).
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, gh and gk represent the concentrations of hoxb1a and krox20, respectively, and ω(t)ψ(t) represents white noise with amplitudes ah and ak, 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]in (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.
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)−φ*(Xc)| 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.
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 PH→K and PK→H using the MAP calculation of Equation 2 for a normalized RA concentration. Our models indicate that PH→K increases exponentially when [RA] is high and PK→H 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).
The noise amplitude, corresponding to the size of ah and ak also must be within an appropriate range for noise attenuation and boundary sharpening (Supplementary Figure S6A). When ah and ak 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 ah and ak 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).
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).
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.
We would like to thank Arthur Lander for initial development of the deterministic model and advice throughout this project, as well as Ines Gehring for zebrafish care. This study was supported by NIH P50 GM‐76516 (QN and TS), NIH R01 GM‐67247 (QN), NSF DMS‐0917492 (QN), and NIH R01 NS‐41353 (TS).
Author contributions: LeiZ, LikunZ, AC, and QN performed the computational research; KR and TS performed the experimental research with zebrafish; LeiZ, QN, KR and TS wrote the paper.
Conflict of Interest
The authors declare that they have no conflict of interest.
Supplementary Information [msb201245-sup-0001.doc]
Supplementary Dataset 1 [msb201245-sup-0002.xls]
Supplementary Dataset 2 [msb201245-sup-0003.zip]
This is an open‐access article distributed under the terms of the Creative Commons Attribution License, which permits distribution, and reproduction in any medium, provided the original author and source are credited. This license does not permit commercial exploitation without specific permission.
- Copyright © 2012 EMBO and Macmillan Publishers Limited