(A) Two-color whole mount in situ hybridization of embryonic zebrafish hindbrains for otx2 (purple, anterior region, far left), krox20 (purple segments, center) and aldh1a2 (red, far right) transcripts from 11 to 14 hpf. otx2 marks the midbrain-hindbrain boundary (MHB), krox20 marks r3 and r5 and aldh1a2 marks the RA production region. Embryos are flat-mounted and shown in dorsal view with anterior to the left. Scale bars: 100 μm. (B) Illustration depicting convergent-extension of the hindbrain. The entire rectangular region, including r1-7 and the RA production region, constitutes the morphogen domain. The hindbrain narrows in the L-R direction (width) and elongates in the A-P direction. (C) Gene regulatory network used to model hindbrain patterning in r2-6. Genes and morphogens with black font were previously used for modeling the r3-r5 pattern [17], while additional factors considered in this model are shown in orange font. Pointed arrows depict up-regulation/activation and blunt arrows depict down-regulation/inhibition. Two morphogens, retinoic acid (RA) and Fibroblast Growth Factor (FGF) diffuse and form two distinct gradients to govern downstream gene expression. (D) Illustration depicting r2-6 and distinct identities (i.e. gene expression signatures) underlying selective cell sorting. Cells in r3 and r5 (blue) express krox20 and cells in r4 (red) express hoxb1a, while both krox20 and hoxb1a levels are low in r2 and r6. Cells in r6 (purple) have high vhnf1 expression. Cells with the same segmental identity attract each other and cells with different identities repulse each other.

(A,B) Experimental measurements of hindbrain dimensions along the A-P (A) and L-R (B) axes at 11, 12, 13 and 14 hpf. Error bars represent standard deviation. Cubic interpolation is used to obtain the smooth curves used in the model. (A) A-P hindbrain length was measured from the posterior edge of the mid-hindbrain boundary (MHB) to the anterior edge of the RA production region. A-P length of the RA production region was based on measurements of the aldh1a2 expression domain. (B) L-R hindbrain width was measured at the A-P position of r4. (C,D) Predicted noisy distributions of morphogen signaling at 14 hpf (C) RA ([RA]_in). (D) FGF ([FGF]_signal). (E-G) Time series of gene expression in r2-6 (the hindbrain is represented as a rectangle for simplification): (E)hoxb1a (red) and krox20 (blue), (F)vhnf1 (purple), (G)irx3 (yellow). (H) Quantifications of rhombomere length, number of dislocated cells (DCs) and sharpness indices (SIs) versus time. Rhombomere lengths (r3-5), and SIs for four boundaries (SI(r2/r3), SI(r3/r4), SI(r4/r5), SI(r5/r6)) and DC number in multiple simulations (n = 100): ‘solid line’: quantities for the simulation shown in (E); ‘brown dashed line’ indicates the average and the width of ‘brown shade” indicates standard deviation; ‘black dashed line’ represents rhombomere lengths from experimental measurements and the error bars represent standard deviation.

(A-C) One-dimensional simulations for the two-morphogen model. (A) The upper panel shows spatial distributions of RA, krox20, hoxb1a, vhnf1, irx3 and FGF. The initial hoxb1a level is modeled as a constant 0.21 over the space. In the lower panel, the initial hoxb1a level is randomly sampled over the space independent of the location. The value is randomly uniformly distributed at a level of [0,0.3]. Solid line represents one simulation. Dashed line represents average values and the width of the shading around each line represents the standard deviation (n = 100). Since fluctuations over multiple simulations are small, solid lines overlap with dashed lines and the small standard deviations result in shadings of negligible width around the dashed lines. X-axis, microns; Y-axis, arbitrary units. (B) Phase diagram of hoxb1a and krox20 distributions with different initial hoxb1a levels. (C) Rhombomere lengths with different initial hoxb1a levels. (D-F) Similar one-dimensional simulations for the one-morphogen model. For (D), in the upper panel, the constant initial hoxb1a level is taken as 0.21; in the lower panel, the initial hoxb1a level is randomly sampled with levels in the range [0.19,0.23] with uniform distribution. Corresponding (E) phase diagram and (F) graph of rhombomere lengths with the one-morphogen model.

(A) Three convergence rates during the 11–14 hpf period are considered in the model, rapid (from experimental measurements, Fig 2B), medium and slow. All start and terminate with the same L-R width. The curve of medium convergence is depicted as a linear function. The curve of slow initial convergence is symmetric to the curve of rapid initial convergence with respect to the curve of linear function. (B-E) Time series of cell distributions with different convergence rates from 11–14 hpf. hoxb1a (red) and krox20 (blue) expression. Dislocated cells (DCs) are highlighted by yellow edges. (B) Three simulations start with the same initial cell distribution (11 hpf) generated by the gene expression model (see Methods). Cell distributions with (C) rapid, (D) medium and (E) slow initial convergence rates from 12–14 hpf. (F) The boundary sharpness index (SI) for four boundaries (SI(r2/r3), SI(r3/r4), SI(r4/r5) and SI(r5/r6)) and DC number versus time. (G-I) Histograms depicting three convergence rates analyzed for (G) rhombomere lengths of r3-5, (H) SI and (I) DC number. Each represents 100 independent stochastic simulations for each convergence rate based on the same parameters. Error bars represent standard deviation.

(A-D) Time series of cell distributions with different convergence rates from 11 to 14 hpf. hoxb1a (red) and krox20 (blue) expression. Dislocated cells (DCs) are highlighted by yellow edges. (A) Three simulations start with the same initial cell distribution (11 hpf) generated by the Gaussian mixture distribution. Cell distributions with (B) rapid, (C) medium and (D) slow initial convergence from 12–14 hpf. (E) The boundary sharpness index (SI) for four boundaries (SI(r2/r3), SI(r3/r4), SI(r4/r5) and SI(r5/r6)) and number of DCs versus time. (F-H) Histograms depicted three convergence rates analyzed for their (F) rhombomere lengths of r3-5, (G) SI and (H) the DC number. Each represents 100 independent stochastic simulations for each convergence rate are based on the same parameters. Error bars represent standard deviation.

The statistics of the dynamics of (A) intracellular RA [RA]_in, (B) vhnf1, and (C) FGF signaling [FGF]_signal at different A-P lengths of the tissue domain: 50 μm, 100 μm, 150 μm and 200 μm. Lines represent average values and the width of the shading around each line represents the standard deviation. (D) The temporal dynamics of total percentages of cells that have committed in each rhombomere (r2-r6). Each panel shows the dynamics of the percentages of cells that have committed in each rhombomere. Data are collected from the full models, see Fig 4.

Parameters for gene regulation were randomly perturbed and a total of n = 1000 simulations are displayed for each convergence rate. There are 513, 563 and 452 simulations for rapid, medium and slow initial convergence, respectively, which successfully generate the r2-6 pattern with four boundaries. (A-F) Dot plots showing the relationship between rhombomere length and boundary sharpness. Each point represents the corresponding quantities for each simulation. (A-C) Length of r4 versus sharpness index (SI) of the r4/r5 boundary with: (A) rapid initial convergence, (B) slow initial convergence and (C) a comparison between rapid and slow initial convergence. (D-F) Length of r5 versus SI of r5/r6 boundary with: (D) rapid initial convergence, (E) slow initial convergence and (F) a comparison between rapid and slow initial convergence. (G) Fractions of simulations achieving roughly equal rhombomere lengths versus the deviation d%. With a deviation d, a simulation has roughly equal rhombomere lengths if the length of each rhombomere is deviated at most d% from its average experimental length (i.e. r3, r4 and r5 are in the ranges of 42*(100%±d%) μm, 34*(100%±d%) μm and 37*(100%±d%) μm, respectively).

Schematic illustration for rapid initial convergence improves pattern robustness by comparing with slow initial convergence (A) Gene regulation/cell fate: rapid initial convergence produces a steeper RA distribution to induce more cells switching from r4 (red) to r5 (blue) identity than slow initial convergence. Consequently, r5 has similar size with r3, consistent with the experimental measurements. (B) Cell sorting: rapid initial convergence increases cell-cell contacts to enhance efficiency of cell sorting, leading to sharper boundaries comparing to slow initial convergence. Number of green lines represent efficiency of cell sorting. (C-C’) Synergy between cell sorting and gene regulation: rapid initial convergence induces an early peak of morphogens for both RA and FGF, leading to an early commitment of cell fates. Cell sorting mechanisms fully function to sharpen boundaries with sufficient time without disrupting cell fate switching.