From local resynchronization to global pattern recovery in the zebrafish segmentation clock

  1. Koichiro Uriu  Is a corresponding author
  2. Bo-Kai Liao
  3. Andrew C Oates
  4. Luis G Morelli
  1. Graduate School of Natural Science and Technology, Kanazawa University, Kakuma-machi, Japan
  2. Department of Aquaculture, National Taiwan Ocean University, Taiwan
  3. Department of Cell and Developmental Biology, University College London, Gower Street, United Kingdom
  4. The Francis Crick Institute, United Kingdom
  5. Max Planck Institute of Molecular Cell Biology and Genetics, Germany
  6. Institute of Bioengineering, École polytechnique fédérale de Lausanne (EPFL), Switzerland
  7. Instituto de Investigación en Biomedicina de Buenos Aires (IBioBA) – CONICET – Partner Institute of the Max Planck Society, Polo Científico Tecnológico, Argentina
  8. Departamento de Física, FCEyN UBA, Ciudad Universitaria, Argentina
  9. Max Planck Institute for Molecular Physiology, Department of Systemic Cell Biology, Germany

Abstract

Integrity of rhythmic spatial gene expression patterns in the vertebrate segmentation clock requires local synchronization between neighboring cells by Delta-Notch signaling and its inhibition causes defective segment boundaries. Whether deformation of the oscillating tissue complements local synchronization during patterning and segment formation is not understood. We combine theory and experiment to investigate this question in the zebrafish segmentation clock. We remove a Notch inhibitor, allowing resynchronization, and analyze embryonic segment recovery. We observe unexpected intermingling of normal and defective segments, and capture this with a new model combining coupled oscillators and tissue mechanics. Intermingled segments are explained in the theory by advection of persistent phase vortices of oscillators. Experimentally observed changes in recovery patterns are predicted in the theory by temporal changes in tissue length and cell advection pattern. Thus, segmental pattern recovery occurs at two length and time scales: rapid local synchronization between neighboring cells, and the slower transport of the resulting patterns across the tissue through morphogenesis.

Introduction

Synchronization of genetic oscillations in tissues generates robust biological clocks. To attain synchrony, cells interact with each other locally and adjust their phase of oscillations. How local interactions between oscillators lead to the emergence of collective rhythms has been studied in static tissues and in dynamic tissues with local cell rearrangements, but how collective rhythms are influenced by the more complex deformations of entire tissues typical in embryogenesis remains challenging and is less well understood. A system to explore this is the synchronization of genetic oscillations during the segmentation of the vertebrate embryo’s body axis, a process termed somitogenesis. Cells in the unsegmented tissue, namely the presomitic mesoderm (PSM) and the tailbud, show collective rhythms of gene expression that set the timing of somite boundary formation and are referred to as the segmentation clock (Oates et al., 2012; Pourquié, 2011). In the tailbud, spatially homogeneous oscillations can be observed. In the PSM, kinematic phase waves of gene expression move from posterior to anterior, indicating the presence of a spatial phase gradient along the axis (Delaune et al., 2012; Shih et al., 2015; Soroldoni et al., 2014). Importantly, this unsegmented oscillating tissue undergoes extensive deformations during the time of segment formation, with complex cellular rearrangements, flows and a changing global size and geometry (Gomez et al., 2008; Jörg et al., 2015; Lawton et al., 2013; Mongera et al., 2018; Steventon et al., 2016). However, our current understanding of synchronization in the segmentation clock follows largely from considering a non-deforming tissue with constant size.

In the zebrafish segmentation clock, Her1 and Her7 proteins repress their own transcription, forming negative feedback loops (Hanisch et al., 2013; Schröter et al., 2012; Trofka et al., 2012). These negative feedback loops are considered to generate cell-autonomous rhythms of gene expression (Lewis, 2003; Monk, 2003; Webb et al., 2016). Cells in the PSM interact with their neighbors via Delta-Notch signaling (Horikawa et al., 2006; Jiang et al., 2000; Riedel-Kruse et al., 2007; Soza-Ried et al., 2014). It is thought that Her proteins repress the transcription of deltaC mRNA, causing oscillatory expression of DeltaC protein on the cell surface (Horikawa et al., 2006; Wright et al., 2011). Binding of a Delta ligand to a Notch receptor expressed by neighboring cells leads to the cleavage and release of the Notch intracellular domain (NICD), which is translocated to the nucleus and modulates transcription of her mRNAs.

Several lines of evidence based on the desynchronization of the segmentation clock show that Delta-Notch signaling couples and thereby synchronizes neighboring genetic oscillators in the zebrafish PSM and tailbud. The first collective oscillation of the segmentation clock occurs immediately before the onset of gastrulation at 4.5 hr post fertilization (hpf), independently of Delta-Notch signaling (Riedel-Kruse et al., 2007; Ishimatsu et al., 2010). Thereafter, cells from embryos deficient in Delta-Notch signaling gradually become desynchronized due to the presence of various sources of noise (Horikawa et al., 2006; Keskin et al., 2018). Single-cell imaging of a live Her1 reporter in the Delta-Notch mutant embryos deltaC/beamter, deltaD/after eight and notch1a/deadly seven during posterior trunk formation (~15 hsspf) shows that Her1 protein oscillation is desynchronized across the PSM cells (Delaune et al., 2012). At the tissue level, Delta-Notch mutants form the anterior 4 ~ 6 segments normally, followed by consecutive defective segments (van Eeden et al., 1996). These phenotypes are not caused by a direct failure of segment boundary formation (Ozbudak and Lewis, 2008), but have been explained in terms of the underlying desynchronization of the segmentation clock (Jiang et al., 2000; Riedel-Kruse et al., 2007).

This desynchronization hypothesis has been formalized as a theory based on coupled oscillators (Riedel-Kruse et al., 2007; Liao et al., 2016). The theory postulates a critical value Zc such that if the level of synchrony becomes lower than this critical value, a defective segment boundary is formed, Figure 1A. In the absence of Delta-Notch signaling, the level of synchrony decays over time and eventually becomes lower than Zc. The time point at which the level of synchrony crosses Zc for the first time is considered to set the anterior limit of defects (ALD), that is the anterior-most defective segment along the body axis. Indeed, theory based on this desynchronization hypothesis can quantitatively explain the ALD in Delta-Notch mutants (Riedel-Kruse et al., 2007).

Figure 1 with 1 supplement see all
Segment boundary defects observed in late and early DAPT washout embryos.

(A) Schematic time series of synchrony level during desynchronization and resynchronization. In the presence of DAPT, the synchrony level decreases due to the loss of Delta-Notch signaling (solid line). DAPT is washed out at 14 hr post-fertilization (hpf; ~9 somite stage; ss) in this panel and resynchronization starts from that time point (dotted line). If the synchrony level is higher (lower) than a critical value Zc, normal (defective) segments are formed. (B) Wild-type control embryo treated with DMSO. (C) Embryo with late DAPT washout at 14 hpf (9 ss). Enlargements of (D) broken or fragmented boundaries, (E) incorrect number of boundaries and (F) left-right misaligned boundaries are shown below. (G) Embryo with early DAPT washout at 9.5 hpf (0 ss). Red, blue and green triangles indicate the anterior limit of defect (ALD), first recovered segment (FRS) and posterior limit of defect (PLD), respectively. (H), (I) Histograms of the difference between PLD and FRS (PLD – FRS) for embryos with DAPT washout at (H) late (14 hpf; n = 30) and (I) early (9.5 hpf; n = 28) stages. Numbers of embryos examined in (H) and (I) were 15 and 14, respectively. FRS and PLD were measured separately between left and right sides of embryos. p<0.05 in Kolmogorov-Smirnov test.

Figure 1—source data 1

Segment boundary defects in embryos with different DAPT washout timing.

https://cdn.elifesciences.org/articles/61358/elife-61358-fig1-data1-v2.xlsx

The desynchronization hypothesis implies that the oscillators could be resynchronized by restoring Delta-Notch signaling, with the expectation that resynchronization of the segmentation clock requires several oscillation cycles for the level of synchrony to smoothly increase and surpass the threshold Zc, giving rise to a transition from defective to recovered normal segments. Due to the constitutive lack of coupling, Delta-Notch mutants cannot be used to analyze resynchronization dynamics. A powerful tool for this purpose is the Notch signal inhibitor DAPT, which inhibits the cleavage and release of NICD, blocking cell coupling. Importantly, DAPT can be washed out to recover cell coupling and this allows cells to resynchronize their genetic oscillations, (Figure 1A; Riedel-Kruse et al., 2007; Ozbudak and Lewis, 2008; Liao et al., 2016; Mara et al., 2007). Previous experimental studies showed that after late washout of DAPT at the nine somite stage, embryos start making normal segments again after several oscillation cycles (Riedel-Kruse et al., 2007; Liao et al., 2016; Mara et al., 2007). The position of the first recovered normal segment after DAPT washout represents the time at which the level of synchrony surpasses Zc for the first time (Liao et al., 2016). In previous studies, almost all segments formed normally after the first fully recovered segment (Riedel-Kruse et al., 2007; Liao et al., 2016), consistent with a monotonic increase of the level of synchrony in the vicinity of Zc, Figure 1A, as expected from the desynchronization hypothesis. Despite this success, it remains fundamentally unclear how tissue-scale gene expression patterns underlying segment recovery reorganize from local intercellular interactions and whether they are affected by tissue size and shape changes that occur during development.

Here, we analyze resynchronization dynamics of the zebrafish segmentation clock at different developmental times using both experimental and theoretical techniques. In contrast to late washout mentioned above, we find that washing out DAPT at earlier developmental stages causes a region of scattered segment defects, where normal and defective segments are intermingled. This striking phenotype was not anticipated by previous models (Riedel-Kruse et al., 2007; Uriu et al., 2017). To investigate the processes in the segmentation clock that yield this distinctive recovery behavior, here we develop a new model of the segmentation clock that encompasses two scales, describing resynchronization and pattern recovery in terms of local interactions between cells in the tissue, as well as global properties of the tissue. In concert, we develop observables that allow pattern dynamics to be quantified and compared between simulation and experiment, the vorticity index and local order parameter. Despite its simplicity, this model can describe the intermingling of normal and defective segments. Numerical simulations indicate that persistent phase vortices appear during resynchronization, resulting in scattered and intermingled segment defects along the axis. The length of the PSM and tailbud and advection pattern influence the recovery process via the transport of phase vortices from the posterior to the anterior of the PSM. Moreover, by including temporal changes to tissue length, advection pattern and coupling strength, the model makes predictions about the pattern of resynchronization at both early and late stages that we confirm experimentally in the embryo.

Results

Scattered embryonic defective segments in zebrafish resynchronization assay

To investigate the processes involved in resynchronization of the segmentation clock, we used a resynchronization assay based on washing out the Notch signaling inhibitor DAPT at different developmental stages. In this assay, zebrafish embryos were placed in DAPT for a defined duration, during which time the segmentation clock desynchronized and defective segments were formed, then washed extensively to allow Delta-Notch signaling activity to resume. Subsequently, the segmentation clock gradually resynchronized and normal segments were made.

Throughout this study, we administered DAPT at 4 hpf, a developmental stage before the oscillating genes of the segmentation clock were expressed. This was a treatment duration of at least 5 hr, sufficient to obtain defects on both left and right sides of the embryo for subsequent resynchronization analysis (Ozbudak and Lewis, 2008). A record of the resulting spatiotemporal pattern of somitogenesis was visualized after its completion by whole-mount in situ hybridization for the myotome segment boundary marker gene xirp2a (xin actin-binding repeat containing 2a) in ~36 hpf embryos, Figure 1B,C. This assay was chosen for its high sensitivity for boundary defect detection (Riedel-Kruse et al., 2007), and we analyzed these staining patterns by scoring boundaries as either normal or defective using established criteria (Riedel-Kruse et al., 2007; Liao et al., 2016), with the exception that we scored the left and right embryonic sides separately. Examples of defects observed with DAPT treatment included fragmented, mis-spaced, or mis-aligned boundaries, Figure 1D–F. To prevent incorrect identification of misaligned boundaries in embryos with bent axes, or in tilted samples, we confirmed that the boundaries outside of the defective region were well aligned between left and right sides. To assign the ordinal segment number to defective boundaries when boundaries were severely fragmented, we used the contralateral side or counted either dorsal or ventral boundary ends, which were often clearer, to estimate their axial position.

As described in the introduction, the location of the transition from normal to defective segments resulting from desynchronization is termed the anterior limit of defects (ALD), given by the first segment along the embryo’s axis that shows a defective boundary, Figure 1—figure supplement 1A. After removing DAPT, resynchronization begins and normal segments form eventually. This transition has been recorded by the first recovered segment along the axis (FRS) (Liao et al., 2016) and then the posterior limit of defects (PLD), the most posterior segment along the axis with a boundary defect, Figure 1—figure supplement 1A (Riedel-Kruse et al., 2007). Note that because segments form rhythmically and sequentially along the body axis, FRS and PLD label both an axial position and the developmental time of segment formation.

In late washout experiments, we observed that a normal segment boundary sometimes formed shortly after ALD even when DAPT was still present, possibly due to desynchronization fluctuations. In previous reports, the definition of the FRS avoided counting these early defects because washout was done late, after full desynchronization. However, when DAPT was washed out early, before ALD occurred, we could not discriminate whether a normal segment formed due to desynchronization fluctuations or as a consequence of resynchronization. The frequency of defects kept growing during an early phase in all conditions until reaching a plateau around segment 9, Figure 1—figure supplement 1B, suggesting that the desynchronization phase lasted until segment 9, at least. Hence, in this study we defined FRS as the first recovered segment after segment 9, when the desynchronization phase was over.

We first analyzed the recovery of normal segments when DAPT was washed out at 14 hpf (twash-out = ~9 somite-stage (ss)), as in previous studies, Figure 1C. Several defective segment boundaries were formed after washout, suggesting that the level of synchrony was still lower than the critical value for normal segment formation during that time interval. However, embryos recovered a normal segment boundary after some time, indicating that the level of synchrony surpassed the threshold, Figure 1C. With this late washout time, we often observed contiguous defective segments before FRS, suggesting that cells in the PSM were completely desynchronized by a DAPT pulse of this length. In addition, PLD coincided closely with FRS, with the distribution of PLD – FRS peaking at lowest values, Figure 1H, Figure 1—figure supplement 1A,F, suggesting that once the level of synchrony surpassed the critical value Zc, it remained larger than Zc, as expected. This observation can be interpreted using the desynchronization hypothesis to indicate that the level of synchrony increases monotonically over time, resulting in the formation of consecutive normal segments after the FRS, Figure 1A.

Importantly, however, when we washed DAPT out at an earlier time twash-out = 9.5 hpf (~0 ss), the majority of embryos had an interval along the axis where normal and defective segments were intermingled between FRS and PLD, with the difference between PLD and FRS distributed more uniformly, Figure 1G,I, and Figure 1—figure supplement 1A,C. This result suggests that the segmentation clock has a level of synchrony close to Zc in this intermingled region, and persistent fluctuations in synchrony level lead to intermittent defective boundary formation.

Physical model of the PSM

According to the desynchronization hypothesis, intermingling of normal and defective segment boundaries suggests a fluctuation of the phase order parameter around its critical value for proper segment formation. How could such large and potentially long-lasting fluctuations of the phase order occur? The desynchronization hypothesis (Jiang et al., 2000) was first formalized in a mean-field theory describing synchronization dynamics from global interactions (Riedel-Kruse et al., 2007). Later, synchronization in the tailbud was analyzed with a theory with local interactions and neighbor exchange by cell mobility (Uriu et al., 2017; Uriu et al., 2010). However, a critical prediction of these theories is that once a population of oscillators is synchronized, a large fluctuation of synchrony level is not expected (Hildebrand et al., 2007; Kuramoto and Nishikawa, 1987; Daido, 1987). Instead of such global phase order fluctuations, other hypotheses for the intermingled defects are the emergence of localized disorder, or the existence of local phase order with a mis-orientation to the global pattern. To explore these potential behaviors, following the general lineage of the clock and wavefront (Cooke and Zeeman, 1976), we develop a physical model of the PSM and tailbud that brings together in a novel framework previous descriptions of (i) the local processes of phase coupling (Morelli et al., 2009) and physical forces (Uriu et al., 2017; Uriu and Morelli, 2014) between neighboring oscillators, and (ii) the tissue-level properties of a frequency profile and oscillator arrest front (Jörg et al., 2015; Morelli et al., 2009), changing tissue length (Jörg et al., 2015), and a gradient of cell mixing (Uriu et al., 2017); furthermore, we introduce an advection pattern of the PSM (Jörg et al., 2015) that changes in time, Figure 2, Figure 2—figure supplement 1 and Supplementary file 1.

Figure 2 with 2 supplements see all
Physical model of the PSM and tailbud.

(A) U-shape geometry of the PSM and tailbud (left), and schematics of key ingredients in the model (right). Each sphere represents a PSM cell. The scale bar indicates the mapping of phase θi to color: white is π/2 and blue is 3π/2. R: right. L: left. Scale bar: 50 μm. (B) Intrinsic cell mobility gradient, (C) cell advection speed, and (D) autonomous frequency gradient along the anterior-posterior axis of the PSM and tailbud. In (C), the absolute value of the spatial derivative of advection speed, referred to as strain rate, is indicated by the red line. L is the length of the PSM L = Lxxa. (E) Kinematic phase waves moving from the posterior to anterior PSM in a simulation. Snapshots of the right PSM are shown. See also Figure 2—video 1. t0 = 302 min is a reference time point. T = 30 min is the period of oscillation at the posterior tip of the tailbud. Parameter values for simulations are listed in Supplementary file 1.

The model describes the PSM and tailbud as a U-shaped domain in a 3D space, Figure 2A, Figure 2—figure supplement 1, see Materials and methods for details. We set the posterior tip of the tailbud as a reference point. Cells are represented as particles in the 3D space, subject to physical forces from other particles when they are closer than a typical length scale that we term cell diameter. In addition, tissue boundaries exert confinement forces on cells (Uriu et al., 2017). Although cells are rendered as spheres in simulations, their effective shapes are in fact dynamic and depend on their local physical interactions. In accordance with previous experimental studies, we consider the spatial gradient of intrinsic cell mobility across the PSM and tailbud, with highest mobility in the posterior, Figures 2B (Lawton et al., 2013; Mongera et al., 2018; Uriu et al., 2017). Axis extension as observed in the lab is described here, from the reference point of the posterior tip of the tailbud, as cell advection from posterior to anterior, Figures 2A,C (Jörg et al., 2015; Morelli et al., 2009; Ares et al., 2012). The value of the spatial derivative of the advection velocity at each position effectively represents the local strain rate, Figure 2C. Cellular motion is described by the overdamped equation with four terms that represent the four physical influences listed above:

dxidt=vd(xi)+v0(xi)ni(t)+j=1,jiNF(xi, xj)+Fb(xi),

where xi is the position of cell i, vdxi is cell advection velocity, v0xi is the speed of intrinsic cell movement with direction nit, Fxi,xj is the physical contact force between cells i and j, N is the total cell number and Fbxi is the boundary force that confines cells within the U-shaped domain.

Genetic oscillation in each PSM cell is described as a phase oscillator with noise terms. The phase oscillators are coupled with their neighbors, representing intercellular interaction with Delta-Notch signaling. We define cells to be neighbors when their distance is shorter than the cell diameter. We also consider a left-right symmetric frequency profile along the anterior-posterior axis of the PSM, as observed in zebrafish embryos, to create kinematic phase waves (Figure 2D; Shih et al., 2015; Soroldoni et al., 2014; Jörg et al., 2015). The mobility of cells in the tailbud increases the communication of the phase between left and right halves. The frequency of oscillation is highest at the tip of the tailbud and gradually decreases toward the anterior PSM. The phase equation describing the oscillators has three terms, representing the frequency profile of the cell-autonomous oscillators, the coupling between neighbors, and noise:

dθitdt=ωxi+κtnθitxj-xidcsinθjt-θit+2Dθξθit,

where θi is the phase of cell i, ωxi is the autonomous frequency at position xi, κ is the coupling strength, nθi is the number of neighboring cells interacting with cell i, dc is the cell diameter, Dθ is the phase noise intensity and ξθit is a white Gaussian noise. When a simulation was started from a uniformly synchronized initial condition, the frequency profile generated left-right symmetrical kinematic phase waves due to growing phase differences between the anterior and posterior PSM, Figure 2E and Figure 2—video 1.

To model the formed segments, we arrest the oscillation when cells leave the PSM from its anterior end xa, Figure 2A,E. The arrested phase stripes in the region x < xa are representative of segment boundaries, and the segment length is the wavelength of this arrested phase pattern. Although the determination of the segment boundary in vivo is a complex process (Dahmann et al., 2011; Naganathan and Oates, 2020), for simplicity we consider that these phase stripes correspond to the boundaries of the resulting morphological segments and the chevron-shaped myotomes that are detected in our experiments. Cells in the formed segment region x < xa continue to be advected anteriorly at the same speed as the anterior end of the PSM, due to axial elongation, and their relative positions are fixed. While small heterogeneities in cell density and cell division may exist in the tissue (Mongera et al., 2018; Steventon et al., 2016; Uriu et al., 2017; Zhang et al., 2008), as a simplifying assumption we keep global cell density constant over time. Consequently, we add new cells with random phase values (Horikawa et al., 2006) at random positions in the PSM and tailbud to balance the cells exiting through the anterior end of the PSM.

Simulation results 1. Intermingled defective segments result from spatially heterogeneous resynchronization in the PSM

Using this physical model, we analyzed resynchronization dynamics in simulations. As an initial condition, we described the state of the PSM and tailbud immediately after DAPT washout by assigning random initial phases to cells, top panel in Figure 3A, as an extended treatment of saturating dose of DAPT was expected to cause such complete randomization of oscillator phases (Delaune et al., 2012). In this desynchronized state, normal segment boundaries – as defined by ordered phase stripes – did not form, matching the experimental appearance of embryos with persistent loss of Delta-Notch signaling. For simplicity, we start our analysis with constant tissue parameters. Below, we will introduce temporal changes to the parameters.

Figure 3 with 12 supplements see all
Resynchronization simulations with constant tissue parameters.

(A) Snapshots of a resynchronization simulation. Color scale as in Figure 2A, also in (D) and (F). The black dotted vertical line indicates the position of the anterior end of the PSM xa = 0. Tissue parameters are constant over time. See also Figure 3—video 1. (B-F) Analysis of local phase order and vortex transport in the simulation shown in (A). (B) Kymograph of local phase order parameter of the right PSM shown in (A). (C) Kymographs of phase vorticity for (left) clockwise and (right) counter clockwise rotations. The phase patterns within the black and yellow boxed space-time domains in (B) and (C) are shown in (D). (D) Snapshots of a phase vortex. The yellow arrows indicate the direction of rotation. (E) Local phase order parameters along the anterior-posterior axis of the PSM at different time points. (F) Time series of local phase order parameter at the anterior end of the right PSM xa. The horizontal broken line indicates the threshold Zc = 0.85 for determining normal and defective segments in simulations. The resultant stripe pattern is on top. In (A), (B) and (F), the blue and green triangles mark FRS and PLD, respectively. Red bracket in (F) highlights a segmental defect resulting from vortex in (D). Parameter values for simulations are listed in Supplementary file 1. (G), (H) Defect size distributions for (G) simulation (n = 800) and (H) embryonic experimental data (n = 134). Defect size indicates how many consecutive segment boundaries are defective in between FRS and PLD. In (G) and (H), the data for different washout timing shown in Figure 4 and Figure 1—figure supplement 1 were pooled to make the histograms. The defect size distribution for each washout timing is shown in Figure 3—figure supplement 3.

Figure 3A shows snapshots of a synchrony recovery simulation, see also Figure 3—video 1. To characterize resynchronization dynamics, we computed local phase order at each position along the anterior-posterior axis, Figure 3B and Figure 3—figure supplement 1. Due to local coupling, cells first formed local phase synchronization, Figure 3A,B and Figure 3—video 1. During the first stage of this synchrony recovery, the kymograph of local phase order shows three locally synchronized domains that extended their size due to an increase in number of synchronized cells at the same time as they were advected in an anterior direction, Figure 3B,E. When the size of the most anterior domain with local phase order above a threshold of 0.85 (Zc) exceeded one segment length at its arrival in the anterior end of the PSM, a first recovered segment (FRS) was formed in a simulation, blue triangles in Figure 3A,B,F. However, because local interactions drove resynchronization in a spatially heterogeneous manner, domains where phase order was lower than Zc could exist more posterior to such a well-synchronized domain, Figure 3B,E. Subsequent arrival of the less synchronized domains caused fluctuation of the local order parameter in the anterior end of the PSM, which resulted in defective segment formation after FRS, Figure 3A,F. Note that patterns of synchronized domains in the left and right sides of the PSM were not well correlated during this time interval (Figure 3—video 1) despite a left-right symmetrical frequency profile. This is in contrast to the fully recovered state, where in the presence of well-synchronized oscillators the symmetrical frequency profile and the communication of phase through the tailbud ensures the left-right symmetry of the wave pattern. Thus, the numerical simulation suggests that a sequence of synchronized and less synchronized domains moving along the PSM results in an intermingling of normal and defective segments along the axis. This intermingling matches the experimental observations.

These less-synchronized domains typically formed persistent phase patterns that rotated along an axis as a vortex, as illustrated in the simulation, Figure 3D and Figure 3—video 1. To detect these patterns, we introduced an index referred to as vorticity, Figure 3—figure supplement 2 and Materials and methods. In brief, the vorticity index detects the core of a vortex, where the phase values circulate from 0 to 2π around a point, but does not measure the spatial extent of the vortex. The kymograph of the vorticity indicates that the less synchronized domains in the kymograph of phase order were caused by the phase vortices, Figure 3A–C. When a vortex was brought to the anterior end of the PSM through cell advection, it was converted into a defective segment boundary and delayed the PLD, Figure 3B,C. Although this process is best visualized dynamically in simulations, one example is illustrated for the vortex in Figure 3D and the resulting local phase order and a segment defect in Figure 3F (bracket). The formation of this particular segmental defect corresponds to time points 350-380 min on the right hand side in Figure 3—video 1. Thus, the formation of persistent local phase patterns with a mis-orientation to the global pattern in the posterior PSM caused defective segment boundaries after the FRS, providing an explanation for the early washout experiments.

The passage of each vortex into the anterior PSM in the simulations generated a segment defect with a characteristic length of one or two segments, Figure 3A–D,F. We compared these defect lengths between simulation and embryonic data and found that across all embryonic stages (Figure 1—figure supplement 1) and all simulations (see Figure 4, below), the size distributions were in quantitative agreement, Figure 3—figure supplement 3. The resulting invariant length scale of the defect size is shown for simulation and embryonic data in Figure 3G,H. This finding is consistent with phase vortices as the origins of intermingled segment defects in the embryo.

Simulation results 2. Dependence of FRS and PLD on each tissue parameter

These results show that the model captures the intermingling of normal and defective boundaries frequently observed in the early washout experiments, but can the model also capture the axial distribution of FRS and PLD observed in the late washout experiments, thereby joining these observations into a coherent picture of resynchronization across developmental stages?

The finding that transport of local phase patterns across the tissue can influence segment recovery suggests that in addition to local intercellular interactions, tissue level parameters may also be important. Previous experimental studies showed that the PSM length becomes shorter as segments are added (Soroldoni et al., 2014; Gomez et al., 2008). Convergent extension by cells in the anterior part of the tissue contributes to advection pattern in the PSM at early developmental stages (Yin et al., 2008). At later developmental stages, cellular flows from the tissues dorsal to the tailbud change the advection pattern (Lawton et al., 2013; Mongera et al., 2018; Steventon et al., 2016). These complex rearrangements are represented in our model in a simplified manner by the advection profile. Thus, several lines of experimental data support changes in the PSM length and its advection pattern, and other properties may also vary during development. We therefore studied how the FRS and PLD depend on each of the tissue parameters in the physical model. We begin by shifting a given single parameter to a new constant value, while leaving the others unchanged for the simulation, Table 1 and Figure 3—figure supplements 411, and return to the time-dependent cases in the next section. Simulations were started from complete random initial phases as in Figure 3. We computed FRS and PLD over 100 different realizations of simulations and the results are summarized in Table 1. See Materials and methods for the quantification of FRS and PLD in simulations.

Table 1
Dependence of FRS and PLD on each tissue parameter.

These results were obtained with simulations where all the tissue parameters were constant over time.

InfluenceTissue parameters
Change in only PLDPSM length, cell mixing, advection pattern
Change in only FRSNone
Change in both FRS and PLDCoupling strength, PSM radius (tube radius)
No or weak effectFrequency profile, tailbud size (torus radius), phase noise intensity

We found that the speed of cell mixing at the tailbud, PSM length L and the cell advection pattern in the PSM did not affect FRS, whereas they strongly affected PLD, Figure 3—figure supplements 46. Faster cell mixing in the tailbud reduced PLD, Figure 3—figure supplement 4, because it prevented the formation of persistent phase vortices at the tailbud. The PSM length and cell advection pattern determined the time taken for a phase vortex to reach the anterior end of the PSM. Longer PSM increased PLD because a phase vortex needed more time to arrive at the anterior end by the cell advection, Figure 3—figure supplement 5. Cell advection underlies the transport of phase vortices, Figure 2C and Figure 3—figure supplement 6. If cell advection was slower at the posterior part of the PSM, PLD became larger because vortices stayed longer at the posterior part, Figure 3—figure supplement 6B. In contrast, if cell advection was faster at the posterior part, PLD became smaller due to faster transport of phase vortices to the anterior part, Figure 3—figure supplement 6B. Thus, these parameters are important for understanding the new phenotype, because each can increase the difference between FRS and PLD, producing the intermingled defects observed experimentally.

In contrast, the PSM radius r and the coupling strength between neighboring oscillators κ0 influenced both FRS and PLD, Figure 3—figure supplements 7 and 8. Smaller PSM radius r decreased FRS and PLD, Figure 3—figure supplement 7. Smaller cell number at each position in the PSM with a smaller radius r allowed local coupling to more rapidly generate a synchronized domain as large as the tissue diameter, leading to a normal segment. Larger values of κ0 reduced FRS and PLD, Figure 3—figure supplement 8. A larger coupling strength reduced the time for local order to form, including vortex patterns. As a result, the last-formed vortex departed the posterior PSM earlier, decreasing PLD.

Coupling keeps phase differences between neighboring oscillators in check. There are two sources of local phase fluctuations in the model: (i) the noise term in individual phase dynamics and (ii) the addition of new cells with random phase values. Desynchronization simulations, where the coupling between cells was absent, demonstrated that the addition of new cells alone was enough to disrupt the wave pattern and compromise the integrity of segment boundaries, Figure 3—figure supplement 9. Additionally, a larger phase noise intensity further contributed to a faster decay of the pattern. However, in resynchronization simulations with coupling, both FRS and PLD only weakly depended on the noise intensity for biologically realistic values, Figure 3—figure supplement 9.

Finally, we found that the shape of frequency profile and the torus size for the tailbud R did not influence either FRS or PLD, Figure 3—figure supplements 10 and 11. Note that there was no parameter that influenced only FRS, Table 1. In summary, FRS was determined by parameters that influence local synchronization of oscillators. PLD, on the other hand, was influenced by parameters that control local synchronization, formation of phase vortices, and their arrival at the anterior end of the PSM.

Simulation results 3. Prediction of PLD from DAPT washout timing, PSM shortening and changing advection pattern

In the previous section, we considered constant values of parameters defining coupling and tissue properties. However, as noted above, some features like PSM length vary during development on timescales that may be relevant for resynchronization (Soroldoni et al., 2014; Gomez et al., 2008; Jörg et al., 2015). To further investigate whether the early and late segmentation phenotypes shown in Figure 1 could result from a common underlying set of processes, we introduced the washout process into the model, and examined the effect of different washout times in simulations in which tissue properties changed over the course of the simulation.

To model differences in timing of DAPT washout, we started with coupling strength κt=0 for t < twash-out and then switched on coupling at t = twash-out. We assigned random phases to the oscillators in the model as an initial condition, assuming that all DAPT treatments completely desynchronized oscillators, as above. Hence, the phase disordered state lasted until t = twash-out and resynchronization begun at that time. We performed 100 realizations of simulations for each washout time twash-out, and recorded the developmental time taken from washout to observation of FRS and PLD, termed the time to FRS (FRS – twash-out) and time to PLD (PLD – twash-out).

In the absence of tissue shortening or a changing cell advection pattern, the times to FRS and PLD were not affected by washout time, as expected, Figure 4A. We analyzed the consequence of PSM shortening on PLD while keeping all the other parameters constant over time, Figure 4B, Figure 4—figure supplement 1 and Figure 4—video 1. For simplicity, we assumed that the PSM length decreased linearly over time in the simulation, Figure 4—figure supplement 1A. PSM shortening decreased the time to PLD for later washout times, Figure 4B, because in a shorter PSM at later somite stages phase vortices reached the anterior end of the PSM more quickly, Figure 4—figure supplement 1B. With higher speed of PSM shortening, the time to PLD after washout became shorter, Figure 4—figure supplement 1C–F. As expected from the independence of FRS on constant PSM length, PSM shortening did not affect FRS, Figure 4B and Figure 3—figure supplement 5.

Figure 4 with 12 supplements see all
Gradual transition from early to late washout boundary phenotypes captured by the physical model.

(A-C) Dependence of times to FRS and PLD on DAPT washout time for different conditions in simulations. (A) Constant tissue where all the tissue parameters remain unchanged during a simulation. (B) PSM length becomes shorter with time. All the other parameters are constant. See also Figure 4—video 1. (C) Cell advection pattern changes at 9 somite stage (ss). Before 9 ss, the strain rate is larger in the anterior than posterior PSM. After 9 ss, the strain rate becomes larger in the posterior PSM. See also Figure 4—video 2. All the other parameters are constant. The box-whisker plots indicate 5, 25, 75, and 95 percentiles. The white bars mark the median. In (B) and (C), the gray dotted lines mark the medians of FRS and PLD in the constant tissue shown in (A). (D) Whole-mount in situ hybridization for the myotome segment boundary marker gene xirp2a in ~36 hr post-fertilization (hpf) embryos. DAPT washout time is 9.5 hpf (0 ss; n = 28), 11 hpf (3 ss; n = 22), 12.5 hpf (6 ss; n = 28), 14 hpf (9 ss; n = 30), and 15.5 hpf (12 ss; n = 26) from top to bottom. Red, blue and green triangles indicate the ALD, FRS, and PLD, respectively. (E) Dependence of times to FRS and PLD on DAPT washout time. Light blue and green box-whisker plots indicate 5, 25, 75, and 95 percentiles for embryonic experimental data (exp.). Dark blue and green box-whisker plots indicate those for simulation data (sim.). The white bars mark the median. The PSM shortening, change in cell advection pattern and increase in the coupling strength are combined in the model, see also Figure 4—videos 3 and 4. The lack of information about the formation of final segments in embryos precludes simulations for the latest washout (12 ss), see the text. (F) Spatial distribution of segment boundary defects. Symbols indicate embryonic experimental data and lines indicate simulation data. Grey dashed vertical line across panels is a guide to the eye. In (A-C), (E), results of 100 realizations of simulations with each washout timing are plotted. Parameter values for numerical simulations are listed in Supplementary files 1 and 2. Source data for (D-F) is available in Figure 1—source data 1.

We next analyzed the effect of a change in cell advection pattern on PLD, keeping all the other parameters constant over time, Figure 4C, Figure 4—figure supplement 2 and Figure 4—video 2. In the model, we represented the change in the advection pattern in a simplified way such that at earlier somite stages (t < tg) the local strain rate was larger in the anterior region of the PSM, whereas the strain rate became larger in the posterior region at later stages (t > tg), Figure 4—figure supplement 2A. We found that such a change in advection pattern increased time to PLD for earlier DAPT washout, Figure 4C. If the change in advection pattern occurred at later developmental stages, time to PLD for later washouts was also increased, Figure 4—figure supplement 2D–G. As described above, the cell advection pattern in the PSM underlies the transport of phase vortices. When advection did not occur in the posterior PSM at earlier somite-stages, the movement of phase vortices relative to the tailbud was slowed in that region, delaying PLD for earlier washout timing, Figure 4—figure supplement 2A–C. As expected from the independence of FRS on constant advection patterns, a change in pattern did not affect FRS, Figure 4C and Figure 4—figure supplement 2D–G.

Taken together, these results predict that the changes in tissue length and cell advection pattern observed in the embryo over developmental time may have an impact on resynchronization dynamics, reflecting in the decrease in difference between FRS and PLD. Thus, changing tissue properties may explain the difference between early and late washout phenotypes.

Embryonic segment recovery in zebrafish depends on timing of DAPT washout

To test the theoretical predictions about the influence of tissue-level changes on the time from washout to FRS and PLD over developmental stages, we next performed DAPT washout experiments, as illustrated in Figure 4D, in which DAPT was removed at different times (twash-out) between 9.5 and 15.5 hpf. We visualized the resulting distribution of segment defects along the axis on both left and right-hand sides of the embryo and identified the FRS and PLD, arrowheads in Figure 4D.

We found that FRS and PLD both increased with later washout times, Figure 4—figure supplement 3. In accordance with the prediction of the model, the time to PLD (PLD – twash-out) decreased gradually over developmental time, Figure 4E. In contrast, the experimentally observed gradual decrease in time to FRS (FRS – twash-out) in Figure 4E was not expected from the simulations, Figure 4B,C. After washout, it took ~13 segments to observe FRS for an earlier washout time, whereas it took eight segments for a later washout time. Embryos yielded more scattered, non-continuous defects with earlier than with later washouts, Figure 4D and Figure 1—figure supplement 1C–G.

Combined, these results revealed the transition between early and late washout segmentation phenotypes. The observed decrease in the time to PLD was qualitatively consistent with the theoretical predictions. However, the expectation that FRS would be independent of washout timing, based on its insensitivity to global tissue properties of PSM shortening and changing cell advection patterns in the model, was not observed experimentally. We therefore hypothesized that FRS might instead be affected by mechanisms that determine the level of local synchronization.

Prediction of embryonic FRS from simulated DAPT washout timing and increasing coupling strength

Local synchronization is thought to be driven by local intercellular interactions (Delaune et al., 2012; Horikawa et al., 2006; Jiang et al., 2000; Riedel-Kruse et al., 2007). The intensity of such local interactions is described in the theory by the coupling strength κ0 between neighboring oscillators. As discussed previously, coupling strength can strongly influence FRS, Table 1 and Figure 3—figure supplement 8. Although FRS can also be influenced by the PSM radius, which becomes smaller with developmental stage, the effect of change in the PSM radius was weaker than the coupling strength within the biologically plausible range, Figure 3—figure supplement 7 and Figure 4—figure supplement 4. Therefore, we tested whether a changing coupling strength could describe the dependence of time to FRS on twash-out observed in experiments.

For simplicity, we assumed that, in the absence of any perturbation, the coupling strength increased as a linear function of time in the simulation, Figure 4—figure supplement 5A. An increase in coupling strength over developmental stages in the embryo could be caused by an increase in the abundance or activity of Delta and Notch proteins in cells (Wright et al., 2011; Liao et al., 2016; Haddon et al., 1998; Westin and Lardelli, 1997), an increase in the contact surface area between neighboring cells, or some other mechanism.

A temporal increase in coupling strength in the simulation, allowing it to double by ~15 ss, led to a decrease in the time to FRS that reproduced the experimental results, Figure 4—figure supplement 5B. The time to PLD also decreased with twash-out. However, the magnitude of reduction in time to PLD was similar to that in time to FRS, meaning that this effect would not be expected to contribute to the observed experimental reduction in the difference between PLD and FRS. These results indicate that the increase in the coupling strength over somite stages alone is sufficient to realize the dependence of time to FRS on twash-out, whereas the global tissue parameters contribute to the behavior of PLD observed in experiments.

Embryonic segmentation defect patterns are captured in simulations by PSM shortening, change in cell advection pattern and an increase in coupling strength

Finally, we simulated a resynchronizing PSM with all the three effects described above combined: the changes in PSM length and cell advection pattern over time were imposed by existing experimental data, the change in coupling strength was motivated by the results in the previous section, and all other parameters remained unchanged, Figure 4E, Figure 4—figure supplement 6 and Figure 4—videos 3 and 4.

The beginning and the end of somitogenesis are special. The segmentation clock becomes active and rhythmic before somitogenesis starts, during epiboly (Riedel-Kruse et al., 2007). At this stage, the embryo undergoes dramatic morphological changes that we do not describe with the current model, which instead describes the PSM shape from 0 ss onwards. For the end of somitogenesis, there is a lack of quantitative information about the formation of the final segments that precludes constraining the theory at this late stage. Therefore, we simulated DAPT washout times from the beginning of somitogenesis 0 ss (twash-out = 9.5 hpf) until 9 ss (twash-out = 14 hpf), where the model could be well parametrized and provided a fair description of tissue shape changes.

The requirement for many realizations to compute FRS and PLD precluded application of standard fitting procedures for determining parameter values in the model. Instead, we used parameter values close to those observed in embryos. We found a decrease in time to PLD with the magnitude of the decrease greater than time to FRS, as observed in the experimental data, Figure 4E. Inclusion of the PSM shortening and change in cell advection pattern in simulations recapitulated the experimental observation that the difference between PLD and FRS became smaller with later washout time. PSM shortening decreased time to PLD, without affecting FRS, thereby reducing the difference between PLD and FRS for later washout time, Figure 4B. The slower cell advection in the posterior PSM at earlier time in simulations delayed PLD without affecting FRS, in principle enlarging the difference between PLD and FRS for earlier washout time, Figure 4C.

In summary, a change in the coupling strength was sufficient to reproduce the behavior of FRS. Combined effects of the PSM shortening and cell advection pattern were the dominant factors that generated the behaviors of PLD in simulations. Thus, the physical model quantitatively reproduced the behaviors of FRS and PLD, suggesting that the physiologically plausible changes in these tissue parameters may underlie behaviors observed in the experiment.

Predicted segment defect distribution confirmed in zebrafish resynchronization assay

We showed that the model could capture the onset of segment boundary recovery and its completion, quantified by FRS and PLD, respectively. However, segment recovery is a complex gradual process reflected in intermingled segment defects. Therefore, we further tested whether the model captured this gradual recovery process between its onset and completion with data that were not used to develop it: the spatial distribution of segment boundary defects along the embryonic body axis.

We used the same parameters listed in Supplementary file 2 that we established to describe FRS and PLD, Figure 4E. Since simulations started from completely random initial phases, the initial fraction of defective segments was one. The fraction of defective segments decreased from one to zero along the body axis after DAPT washout, shifting posteriorly for later twash-out, Figure 4F. We then compared the simulated axial distribution of defective segments with embryonic DAPT washout experiments. We restricted the comparison to the washout phase of the experiment. We counted the number of defective segments along the axis in embryos and defined the fraction of defective boundaries at each segment position, Figure 4F. The distributions of defective segments were similar between left and right sides of embryos, Figure 4—figure supplement 7A. After washout, the fraction of defective segment boundaries gradually decreased, and eventually it became zero, suggesting that synchronization was fully recovered at that time. As DAPT was washed out at increasingly later times, defective segment boundaries continued to more posterior locations, in agreement with simulations, Figure 4F.

From this distribution in embryos, we could compute ALD, FRS, and PLD using a probabilistic theory that assumed left-right independence, see Appendix and Figure 4—figure supplement 7. This distribution also explained the ratio of single defects, where either the left or right segment was defective at a segment locus, to double defects, where both left and right segments were defective, Appendix and Figure 4—figure supplement 8. This agreement between experimental data and probability theory for the fraction of single defects indicates that recovery occurred independently between left and right PSMs, Figure 4—figure supplement 8C. In summary, the physical model predicted the segment defect distribution, providing a thorough description of synchrony recovery.

Discussion

The segmentation clock produces dynamic patterns that determine the formation of vertebrate segments. This multicellular clock, consisting of thousands of cells that make the PSM and tailbud, produces a kinematic wave pattern. The integrity of this wave pattern relies on local synchronization of the oscillators mediated by Delta-Notch cell-cell signaling. Our current view of synchrony is largely informed by desynchronization experiments in which cells were uncoupled by interfering with Delta-Notch signaling, resulting in a loss of wave patterns. In contrast, resynchronization, where oscillators re-establish coherent rhythms from a desynchronized state, can be used to probe how tissue-scale collective patterns arise from local interactions during morphogenesis.

In this study, we applied a combination of experiments and theory to explore the recovery of normal body segments during the resynchronization of oscillators. Our surprising experimental discovery was regions of normal and defective segments intermingled along the body axis following DAPT washout. Since we seek to capture pattern recovery at multiple scales, our new physical model draws from previous work describing cellular oscillations with an effective phase variable together with local intercellular interactions (Uriu et al., 2017; Morelli et al., 2009), as well as larger scale mechanics such as cell movements (Uriu et al., 2017; Uriu and Morelli, 2014) and tissue shape changes (Jörg et al., 2015). One of the novelties of this work is to combine these descriptions in a coherent framework for the first time. The second key novelty is the framework itself, which uses tissue geometry linked with changing mechanical and biochemical properties, such as cellular advection profile and coupling strength.

The choice of a phase description for the cellular oscillator is motivated by its generality and by prior success in analyzing coupling in the segmentation clock (Riedel-Kruse et al., 2007; Uriu et al., 2017; Herrgen et al., 2010). Core genetic components of the cellular oscillator have been identified and their dynamics can be described by detailed delay models (Schröter et al., 2012; Lewis, 2003; Monk, 2003; Horikawa et al., 2006; Ay et al., 2013; Hirata et al., 2004; Jensen et al., 2003), but since we do not measure any of the components in these networks, the choice of phase oscillators captures the core oscillatory behavior without additional underconstrained parameters (Kotani et al., 2012; Kotani et al., 2020). Thus, although we do not anticipate any qualitative differences between these modeling approaches, future work could include such a detailed description of oscillatory genes, potentially allowing a more direct connection to mutant conditions or imaging experiments.

The model qualitatively explains the formation of these intermingled defective segments by the emergence of persistent phase vortices in the posterior PSM and their advection through the tissue to the anterior. The intermingled defects span 1 or 2 segments in the embryonic data, but are not multiples of the local segment length. In simulations with the parameter set we determined here, the defects have the same characteristics. The vortices arise from the local coupling of desynchronized oscillators and the advection is a consequence of axis elongation. As vortices arrive at the anterior end, their mis-oriented local phase patterns result in segment defects, before global pattern recovery is achieved. Thus, the intermingled regions are explained by the intermittent arrival of vortices with a size of approximately one segment. Note that vortices only form during resynchronization, being seeded by the random phases of the desynchronized oscillators and requiring local coupling, and not during desynchronization, which starts by the loss of coupling in a population with a locally smooth distribution of phase.

Importantly, for quantitative comparison between simulations and experiments, we needed to introduce observables such as an index of off-lattice vorticity, and a local order parameter to distinguish normal and defective segments. It was also necessary to incorporate global features of the PSM and tailbud that are observed in the embryo, such as changing tissue length, a change in the advection pattern and a gradient of cell mixing. To achieve a description of faster recovery at later developmental stages, as indicated by shorter time to embryonic FRS, our model includes a time-dependent increase in effective coupling strength. This is plausible from existing data of Delta and Notch gene expression, for example, but remains an expectation of the current work. Nevertheless, simulations of the model confirm that stronger local coupling leads to faster resynchronization and recovery of the gene expression pattern, as expected from previous experimental work (Liao et al., 2016). To reduce computational complexity and time, in our model we neglected the coupling time-delays that are thought to occur in the segmentation clock (Herrgen et al., 2010; Oates, 2020; Shimojo and Kageyama, 2016; Yoshioka-Kobayashi et al., 2020). We do not anticipate their inclusion would affect our main results, since other systems with time-delayed coupling exhibit vortex formation (Jeong et al., 2002), although a quantitative description of some phenotypes may need to account for coupling delays in Delta-Notch signaling.

Because our data was captured from snapshots at a late developmental stage, for simplicity we assumed a constant segment length. However, it has been long recognized that segment length at the time of formation gradually decreases in the posterior of all extending vertebrate embryos (Gomez et al., 2008; Schröter et al., 2008). These early length differences are gradually reduced and eventually eliminated by subsequent growth in the zebrafish (Lleras Forero et al., 2018). It is possible that a quantitative description should also take the phenomenon of changing segment length into account. This would require a timelapse analysis of the forming segments during perturbations and a greater understanding of the mechanism underlying the determination of and change in segment lengths during development (Ishimatsu et al., 2018; Simsek and Özbudak, 2018), and could be accommodated in the framework of the current model. In summary, the formulation of the model allows a quantitative comparison to experimental data and a phenomenological understanding of pattern recovery.

Vortices could arise in other models of the segmentation clock with similar local interactions (Uriu et al., 2017; Morelli et al., 2009; Uriu and Morelli, 2014; Hester et al., 2011). Indeed, vortex formation is a common feature of systems that can be described with locally coupled polar variables, in both excitable and oscillatory media (Mikhailov and Showalter, 2006). These include biological systems such as cAMP patterns in aggregating populations of Dictyostelium cells (Durston, 2013) and spiral patterns in heart tissue (Winfree, 1980), planar cell polarity (Burak and Shraiman, 2009), chemical systems like the BZ reaction (Winfree, 1980; Kuramoto, 1984; Winfree, 1972) and in physical systems generally (Kosterlitz, 2016; Kosterlitz and Thouless, 1973). In unconstrained and homogeneous oscillatory media, vortices can grow to system-size length scales and are very stable on their own (Mikhailov, 1990). However, the effect on size and stability of concurrent features in our model such as the particular geometrical confinement, as well as the existence of a frequency profile and non-uniform advection has not been examined. Thus, the relationships between vortex size, frequency of arrival, orientation and the underlying processes such as coupling strength and tissue geometry in our model are interesting topics for future exploration.

Recent experiments using cells isolated from mouse tailbuds and reaggregated to form so-called emergent PSM have shown the emergence of striking wave patterns that depend on Notch signaling, suggesting that local interactions can drive large-scale dynamical features (Hubaud et al., 2017; Tsiairis and Aulehla, 2016). The relationship between such features and normal versus defective segment boundary formation remains to be explored. Since phase vortices arise naturally in systems of locally coupled oscillators starting from disordered initial conditions, Video 1, we predict that these structures will form also in mammalian PSM tissue culture systems (Hubaud et al., 2017; Tsiairis and Aulehla, 2016; Diaz-Cuadros et al., 2020; Matsuda et al., 2020). The framework of our model with different geometries will facilitate analysis of dynamics in these and other collective cellular systems with both local interactions and tissue-level deformations.

Video 1
Formation of phase vortices in a resynchronization simulation in a cuboid domain, 110 × 110 × 55 μm3.

The color indicates 1+sinθi/2. The number of oscillators is N = 998. Frequencies of all the oscillators are identical, ω = 0.2094 min−1. Values of the other relevant parameters in Equation (1) in the supporting information are: κ0 = 0.07 min–1, κs = 0 min–1, Dθ=0.0013 min–1, μ = 8.71 μm min−1, dc = 11 μm, μb = 20 μm min−1, and rb = 1 μm.

In simulations, the kinematics of phase vortices across the PSM is determined by global tissue properties, including the PSM length and its cell advection pattern. In this way, the recovery of a gene expression wave pattern across the PSM and subsequent normal boundary formation is set both by the timescales of local synchronization through intercellular interactions, and also by those of morphological processes. In addition, this result implies that the quantification of global tissue parameters will be necessary to obtain quantitative agreement between theory and experiment. Temporal changes in PSM length have been measured in various species (Soroldoni et al., 2014; Gomez et al., 2008; Ishimatsu et al., 2018), and cell advection patterns across the PSM have been investigated (Lawton et al., 2013; Mongera et al., 2018; Steventon et al., 2016; Bénazéraf et al., 2010). Involvement of these global tissue properties discriminates the resynchronization of the segmentation clock from its desynchronization, which is dominated by local cellular properties such as noise in clock gene expression (Jiang et al., 2000; Riedel-Kruse et al., 2007; Keskin et al., 2018; Jenkins et al., 2015). Furthermore, morphological tissue development distinguishes the synchronization of the segmentation clock from that of other biological oscillator systems in spatially static tissues, such as the suprachiasmatic nucleus in the mammalian circadian clock (Webb and Oates, 2016).

Numerical simulations of the physical model showed that the first recovered segment FRS and the posterior limit of defects PLD contain different information about resynchronization. From FRS, we estimated properties of local synchronization, such as the coupling strength between genetic oscillators. In contrast, PLD is influenced by global tissue properties, such as tissue length and cell advection pattern, that determine the transport of persistent mis-oriented phase patterns. Hence, it is important to choose these measures appropriately depending on the question addressed. For example, recent studies suggested that the level of cell mixing observed in the tailbud promotes synchronization of genetic oscillation (Uriu et al., 2017). This effect appears in PLD, but not in FRS, because the elevated mixing in the tailbud prevents formation of the local phase patterns at the posterior PSM. On the other hand, FRS can be a better measure for the coupling strength (Liao et al., 2016), because it is less affected by the other tissue parameters than PLD.

Based on the new theoretical framework developed in this study, we have proposed that the distinctive intermingling of well-formed and defective segments seen during recovery arises from persistent phase vortices. Furthermore, we expect the kinetics of phase vortex formation and transport to change throughout development. We have presented quantitative experimental evidence on the size of defects, the shape of the spatial defect distribution and the left-right independence of defects that supports these theoretical predictions. The occurrence of vortices and the role that these transient patterns have in causing intermingled defects remain to be directly observed. Such transient patterns would be difficult to recognize in snapshots of the segmentation clock due to the restricted system geometry, the discrete character of cellular tissue, and the difficulty of determining phase from a single time point, and consequently may have been overlooked in previous experiments. However, they should be within reach of techniques for perturbation of cell coupling combined with live imaging over the long durations required to observe them in the embryo (Yoshioka-Kobayashi et al., 2020; Wan et al., 2019). Such experiments will be challenging nevertheless, as quantitatively testing the model will require combining (i) imaging, segmentation and tracking of single cells within the embryo across appropriate developmental stages; (ii) an approach to reliably estimate phase from segmentation clock transgenic reporter signals with amplitude fluctuations; (iii) data that allows the calculation of vorticity; and (iv) a method of correlating the passage of vortices with the resulting intermittent segment defects in the same embryo. The use of mutants with reduced coupling strength (Liao et al., 2016) may increase the time interval over which vortices are produced, and may thus facilitate their observation.

In conclusion, our study of resynchronization has revealed how the segmentation clock can be influenced by global developmental processes such as tissue length change, cell advection pattern and cell movement, as well as local coupling strength between cells. Our findings suggest that segmental pattern recovery occurs at two scales: local pattern formation and transport of these patterns through tissue deformation. Other developing systems such as Dictyostelium colony aggregation show similar dynamics, with individual motile cells interacting via local signaling to generate spiral waves that guide the formation of large multicellular fruiting bodies (Durston, 2013). In the case of the elongation of the Drosophila pupal wing, local cell rearrangements within the epithelium combine with tissue-level pulling forces applied by a neighboring tissue to form the final pattern (Etournay et al., 2015). The hallmark of these systems is an interplay between locally driven interactions and global morphological changes, pointing to a common principle of pattern dynamics within developing tissues.

Materials and methods

Animals and embryos

Request a detailed protocol

Zebrafish (Danio rerio) adult stocks were kept in 28°C fresh water under a 14:10 hr light:dark photoperiod. Embryos were collected within 30 min following fertilization and incubated in petri dishes with E3 media. Until the desired developmental stages (Kimmel et al., 1995), embryos were incubated at 28.5°C. For whole mount in situ hybridization, PTU (1-phenyl 2-thiourea) at a final concentration of 0.003% was added before 12 hr post fertilization (hpf) to prevent melanogenesis. All wildtypes were AB strain.

DAPT treatment and washout

Request a detailed protocol

DAPT treatment was carried out as previously described (Riedel-Kruse et al., 2007). 50 mM DAPT stock solution (Merck) was prepared in 100% DMSO (Sigma) and stored in a small volume at −20°C. Embryos in their chorions were transferred to 12-well plates at 2 hpf in 1.4 ml E3 medium with 20 embryos per well. 50 μM DAPT in E3 medium was prepared immediately before the treatment. To prevent precipitation, the DAPT stock solution was added into E3 medium while vortexing, and then filtered by 0.22 μm PES syringe filter (Millipore). DAPT treatment was initiated by replacing E3 medium with E3/DAPT medium at desired stage. At 9.5 (0 somite stage: ss), 11 (3 ss), 12.5 (6 ss), 14 (9 ss), or 15.5 hpf (12 ss), DAPT was washed out at least twice with fresh E3 medium + 0.03% PTU. Embryos were dechorionated and fixed at 36 hpf. All experimental steps were incubated at 28.5°C, except for short operations, for example, washing out, which were at room temperature.

Whole-mount in situ hybridization and segmental defect scoring

Request a detailed protocol

In situ hybridization was performed according to previously published protocols (Thisse and Thisse, 2008). Digoxigenin-labeled xirp2a (clone: cb1045) riboprobe was as previously described (Riedel-Kruse et al., 2007). Stained embryos were visually scored under an Olympus SZX-12 stereomicroscope and images were acquired with a QImaging Micropublisher 5.0 RTV camera. Defective segment boundaries were scored as previously described (Riedel-Kruse et al., 2007), with the addition that left and right (LR) sides of the embryo were scored and assessed separately. Since boundary formation in unperturbed embryos is extremely reliable, with errors occurring in less than 1 in 1000 embryos, any interruption or fragmentation to the boundary, and/or alterations in spacing, or alignment was recorded as a defect. The following observables were collected for each LR side: an anterior limit of defects (ALD), that is, the position of the first defective boundary; a posterior limit of defects (PLD), that is, the position of the last defective boundary; and the first recovered segment (FRS), that is, the position of the first normal segment after the segment 9, as described below.

Each segment has anterior and posterior boundaries. For the segment i (i = 1, 2, 3,. .,), the anterior and posterior boundaries were numbered as i – 1 and i, respectively, Figure 1—figure supplement 1A. Both ALD and PLD were numbered using the posterior boundary of the defective segment. For example, if the jth segment boundary was the last defective boundary, PLD was numbered as j. FRS was numbered using the anterior boundary of the recovered segment. For example, if the kth segment boundary was the first normal boundary after washout, the first normal segment was segment k, and FRS was numbered as k – 1, Figure 1—figure supplement 1A. With this definition of FRS, if the first normal segment boundary after the washout was located immediately after the last defective boundary, the difference between PLD and FRS, PLD – FRS, was 0, Figure 1—figure supplement 1C–G. Definitions of FRS and PLD for simulation date were the same as those for experimental data written here, and described in a later section.

Physical model of PSM and tailbud cells

3D tissue geometry

Request a detailed protocol

We model the PSM and tailbud as a U-shaped domain in a 3D space, Figure 2—figure supplement 1A. The left and right PSMs are represented as two tubular domains Ωl and Ωr, respectively with radius r. The tailbud is described as a half toroidal domain Ωt with a larger radius R and smaller radius r.

We implement this U-shaped domain in the spatial coordinate system as follows. The x-axis is along the anterior-posterior axis of the PSM. We set the initial position of the anterior end of the PSM at x=0 and the posterior tip of the tailbud at x=Lx. We denote the position of the anterior end of the PSM at time t as xa(t), so xa(0)=0. The length of the tissue is L=Lxxa. Xc denotes the position of the center of torus core curve. The position of posterior tip of the tailbud can then be written as Lx=Xc+R+r. The y-axis points along the left-right axis and z-axis along the dorsal-ventral axis of an embryo.

Reference frames: Lab reference frame

Request a detailed protocol

To describe cell movements and tissue deformations it is important to define a reference frame. A natural choice may be a reference frame which is at rest in the Lab, termed a Lab reference frame. For instance, the origin of x-axis can be set at the initial position of the anterior end of the PSM at t=0. We write the position of the posterior tip of the tailbud in this Lab reference frame as xt(L)t, where superscript (L) indicates variables in the Lab reference frame. Because an embryo elongates posteriorly, dxt(L)t/dt>0. The position of anterior end of the PSM in the Lab reference frame xa(L)t also changes over time due to the formation of new segments. In this reference frame, PSM length is Lt=xt(L)t-xa(L)t. The rate of change of PSM length is dLt/dt=dxt(L)t/dt-dxa(L)t/dt. PSM length is constant if the velocity of the tailbud is the same as that of the anterior end of the PSM, that is dxt(L)t/dt=dxa(L)t/dt, but it will be changing over time whenever dxt(L)t/dtdxa(L)t/dt. In this Lab reference frame, the position of a cell xc(L)t may change both due to tissue elongation and due to tissue deformations such as local tissue stretch.

Reference frames: tailbud reference frame

Request a detailed protocol

In this study however, we mostly use the tailbud reference frame (t) unless noted otherwise. In this reference frame we measure the position of cells and tissue landmarks from the tailbud. The position of the tailbud is fixed at xt(t)t=0 for all time t, so dxt(t)t/dt=0. The position of a tissue landmark l is xl(t)t=xl(L)t-xt(L)t. In this tailbud reference frame, the velocity of a landmark is zero if it moves with the same velocity as the tailbud in the Lab frame. In other words, a tissue landmark moves relative to the tailbud when their velocities in the Lab reference frame are different. The positions of cells in the tailbud reference frame may change over time, and this will enter as cell advection in the cell equations of motion as we describe below. In the following, we omit the superscript (t) for the tailbud reference frame to alleviate the notation.

Reference frames: computational implementation

Request a detailed protocol

For the implementation of the tailbud frame for numerical simulations, we fix the position of the posterior tip of the tailbud at x=Lx. In contrast, cells positions and other tissue landmarks such as the anterior end of the PSM change over time. We model cell displacement relative to the tailbud as cell advection in the anterior direction as explained below.

Cell mechanics and equation of motion

Request a detailed protocol

We model PSM and tailbud cells as particles in the 3D domain. We denote the number of cells in the PSM and tailbud at time t as N(t). State variables for these cells are their position x in the domain, the phase θ of their genetic oscillators and the unit vector n representing cell polarity for intrinsic cell movement. Position of cell i at time t is denoted as xi(t)=(xi(t),yi(t),zi(t))(i=1,2,3,...,N(t)). We describe cellular motion by the following over-damped equation based on the previous study (Uriu et al., 2017):

(1) dxidt=vd(xi)+v0(xi)ni(t)+j=1,jiNF(xi,xj)+Fb(xi),

where vdxi is cell advection from posterior to anterior, v0xi represents the speed of intrinsic cell movement, nit is the cell polarity pointing the direction of intrinsic movement, Fxi,xj is a physical contact force between cells i and j, and Fbxi is a boundary force that confines cells within the U-shaped domain. We explain each of these terms below.

Elongation, strain rate, and cell advection

Request a detailed protocol

To introduce tissue elongation and deformations we first discuss the movement of tissue landmarks and cells in the Lab reference frame. Then, we switch to the tailbud reference frame where the tissue elongation results in cell advection from the posterior to the anterior directions.

In the Lab reference frame, segments and cells within do not move. The tailbud, together with cells at that point, moves away from segments with a velocity vt(L). We call elongation to this outgrowth of the tissue. As cells differentiate into a segment, the anterior end of the PSM also moves after the tailbud, with a velocity va(L). If this velocity matches that of the tailbud va(L)=vt(L), the length of the PSM remains constant. If this velocity differs from that of the tailbud va(L)vt(L), the length of the PSM changes over time, causing a global tissue deformation. At these PSM ends we have boundary conditions for the cells velocities: (i) a cell at the tailbud moves with the same velocity as the tailbud and (ii) a cell at the anterior end of the PSM is at rest, like neighboring cells within a segment, even though the anterior end of the PSM itself is moving in the Lab frame. Within the PSM, cells are subject to an advection velocity field v(L)(x) that satisfies these boundary conditions. This velocity field producing internal local deformations is caused by internal strains described below. The shape of this velocity field determines the kind of local deformations. For example, a linear velocity field produces uniform deformations, while a non-linear velocity field produces non-uniform deformations, as the piecewise linear function described below in the implementation.

To see the relation between the velocity field and underlying internal strains, let x1t and x2t be the x positions of cells 1 and 2 at time t in the Lab reference frame, where we drop the (L) superscript for notational convenience. We assume that their positions are close to each other, so the distance between them Δx(t)=x2(t)x1(t) is small. Due to the presence of an advective velocity field, the velocities of neighboring cells may be different, that is v1v2, with v1=dx1/dt=v(x1(t)) and v2=dx2/dt=v(x2(t))=v(x1(t)+Δx(t)). Thus, during a short time interval Δt these cells change their relative position due to this velocity difference, Δx(t+Δt)=x2(t+Δt)x1(t+Δt)Δx(t)+(v2v1)Δt. Thus, there is a local strain (Δx(t+Δt)Δx(t))/Δx(t)(v/x)Δt, where the velocity field gradient v/x is a local strain rate. Thus, the advection velocity field effectively describes the presence of local strains and determines the local tissue deformation. Although such local deformations make cells move apart from each other along the x-axis, intercellular and boundary forces in Equation (1) constrain cell distances and we do not observe large density fluctuations due to the velocity field gradient, as can be seen in simulations. Furthermore, where local cell density fluctuations do happen, cell addition described below will bring back the density to its average value.

Cell advection patterns relative to the tailbud reference frame

Request a detailed protocol

Back to the tailbud reference frame, the velocity field is v(t)(x)= v(L)(x)-vt(L), which we refer to as the cell advection pattern vd(x) in Equation (1) and in the main text. We model different advection patterns of the PSM effectively by changing the spatial derivative of the advection speed in subdomains of the tissue. The advection pattern of the PSM may change at a certain developmental stage as we discussed in the main text. Below, we model the spatial profile of cell advection speed and its temporal change. For simplicity, we divide the PSM into two subdomains, namely the anterior PSM ((xxa)/L<xq) and the posterior PSM ((xxa)/L>xq), and consider different strain rates vd/x for each domain, Figure 2C. The advection field vd(x) in Equation (1) depends on spatial position along the anterior-posterior axis x (xaxL) and is described as (Figure 2C):

(2a) vdx=-vdx-/L,0, 0T, (2a)
(2b) vd(χ)={vavp(1xq)xqχ+vaχxq,vpχ+vpχ>xq,

where x-=x-xa, va>0 and vp>0 are the parameters that determine the advection speed, the superscript T denotes transposition of the vector, and L is the length of the PSM L=Lxxa. The slope changes at the position xq(0<xq<1). Thus, xq divides the PSM into anterior and posterior subdomains. Within each domain, the strain rate is uniform, whereas it may be different between these two domains. The strain rate, given by the magnitude of the spatial gradient of advection speed, is vp in the posterior PSM domain and va-vp1-xq/xq in the anterior PSM domain.

Temporal change in advection pattern

Request a detailed protocol

The advection pattern of the PSM in embryos may change at a certain developmental stage. To model the change in the advection pattern of the PSM, we change the value of vp in Equation (2) at time t=tg. We assume that for t<tg, advection occurs mostly at the anterior part, vp<va. For t>tg, we assume that advection occurs at the posterior part of the tissue vp>va.

Gradient of intrinsic cell movement speed

Request a detailed protocol

A cell mixing gradient is observed along the anterior-posterior axis of the PSM (Lawton et al., 2013; Mongera et al., 2018; Mara et al., 2007; Uriu et al., 2017). The degree of cell mixing is higher in the tailbud and posterior region of the PSM than anterior region. To model the cell mixing gradient, we assume that the speed of intrinsic cell movement v0(x) in Equation (1) depends on the spatial position x-=x-xa along the anterior-posterior axis (Figure 2B):

(3) v0x=vs1+1-x-/LXvh-1, (3)

where vs is the maximum speed at the posterior tip of the tailbud, Xv is the lengthscale of the mobility gradient, and the coefficient h determines the steepness of the mobility gradient, Figure 3—figure supplement 4A.

Cell polarity

Request a detailed protocol

The unit vector ni(t) in Equation (1) represents the polarity of cell i and determines the direction of intrinsic cell movement. In spherical coordinates, ni(t)=(sinϕi(t)cosφi(t),sinϕi(t)sinφi(t),cosϕi(t))T with 0ϕi(t)π and 0φi(t)<2π. We assume random change of cell polarity by letting the polarity angles ϕi and φi perform a random walk. The time evolution of ni is described by Langevin equations for these two polarity angles (Uriu et al., 2017):

(4a) dϕitdt=Dϕtanϕit+2Dϕξϕit, (4a)
(4b) dφitdt=2Dϕξφitsinϕit, (4b)

where Dϕ is the polarity noise intensity. White Gaussian noise ξϕit and ξφit satisfy ξϕit=0, ξϕitξϕjt'=δijδt-t', ξφit=0, ξφitξφjt'=δijδt-t' and ξϕitξφjt'=0.

To obtain Equations (4a) and (4b), we first consider isotropic diffusion in a plane that is locally tangent to the sphere. Next, we normalize the polarity vector to keep its unit length. Finally, we transform to global spherical coordinates. Formally, we first update the vector nit:

(4c) n~it+dt=nit+2Dξφitmxdt+2Dξϕitmydt, (4c)

where mx=(ni(t)×ez)/ni(t)×ez, my=(mx×ni(t))/mx×ni(t) and ez=(0, 0, 1). These unit vectors mx and my define a plane tangent to the unit sphere at position ni(t). Random displacement on this tangent plane will move the polarity vector away from the surface of the sphere, so its length n~it+dt will be larger than 1,

(4d) n~it+dt=1+2Ddt+Odt2. (4d)

We normalize the polarity vector to correct for this radial displacement

(4e) nit+dt=n~it+dtn~it+dt, (4e)

and substitute Equations (4c) and (4d) into Equation (4e) to obtain the differential equation for the unit vector ni(t)

dnitdt=2Dξφitmx+2Dξϕitmy-2Dnit.

Finally, we transform variables from cartesian to spherical coordinates using Ito calculus, and obtain Equations (4a) and (4b). With these Langevin equations, the polarity vector ni performs an isotropic and uniform random walk on the surface of a unit sphere.

Intercellular force

Request a detailed protocol

For simplicity, we consider a linear elastic force F(xi,xj) to model volume exclusion of cells in Equation (1). If the distance between two cells becomes shorter than a threshold dc which we term cell diameter, they repel each other (Uriu et al., 2017; Uriu and Morelli, 2014):

(5a) Fxi, xj=Fxi, xjxj-xixj-xi, (5a)

where F(xi,xj) is the modulus of intercellular force

(5b) F(xi,xj)={μ(|xixj|/dc1)|xixj|/dc1 0|xixj|/dc>1,

where μ>0 is the coefficient for intercellular force.

Boundary force

Request a detailed protocol

Fb(xi) in Equation (1) is a confinement force that a cell receives from the boundaries of the domains. Below, we specify this confinement force depending on which tissue domains cells are located in, Figure 2—figure supplement 1B.

In the PSM tubes

Request a detailed protocol

In the PSM regions Ωl and Ωr, we consider boundary force in y and z directions. The position of cell i in the right PSM region Ωr can be written as

(6a) xi=xi, r+2R+ricosqi(c), Zc+risinqi(c), (6a)

see Figure 2—figure supplement 1B. If cell i is in the left PSM region Ωl, its position can be described as

(6b) xi=xi, r+ricosqi(c), Zc+risinqi(c). (6b)

Then, we define frictionless boundary force in the columns as

(7) Fb=FbxFbyFbz=0-μbe-δy(c)rbcosqi(c)-μbe-δz(c)rbsinqi(c), (7)

where δy(c)=r-ricosqi(c) and δz(c)=r-risinqi(c).

In the tailbud half-toroid

Request a detailed protocol

Position of cell i within the half-toroidal domain Ωt can be expressed as:

(8) xi=Xc+Δxi=(Xc+ΔxiYc+ΔyiZc+Δzi)=(Xc+Rcospi+ricosqicospiYc+Rsinpi+ricosqisinpiZc+risinqi),

where pi=tan1(Δyi/Δxi), qi=tan1(Δzicospi/(ΔxiRcospi)) and ri=Δzi/sinqi (Figure 2—figure supplement 1B). We then define the boundary force in the half-toroidal domain as:

(9) Fb=FbxFbyFbz=-μbe-δxrbcospicosqi-μbe-δyrbsinpicosqi-μbe-δzrbsinqi, (9)

where μb is the coefficient and rb is the length scale of the boundary force, δx=|Rcospi+rcospicosqiΔxi|δy=|Rsinpi+rsinpicosqiΔyi| and δz=|rsinqiΔzi|.

Phase equation

Request a detailed protocol

The phase dynamics of genetic oscillators in single PSM cells is described by a phase oscillator model (Riedel-Kruse et al., 2007; Uriu et al., 2017; Morelli et al., 2009; Kuramoto, 1984):

(10) dθitdt=ωxi+κtnθitxj-xidcsinθjt-θit+2Dθξθit, 10

where θi is the phase of cell i, ωxi is the autonomous frequency at position xi, κ is the coupling strength, nθi is the number of neighboring cells interacting with cell i and Dθ is the phase noise intensity. Interactions occur between touching cells |xjxi|dc. The coupling strength κ may be time-dependent as described below. ξθit is a white Gaussian noise satisfying ξθit=0 and ξθitξθjt'=δijδt-t'.

We assume the frequency profile ωx along the anterior-posterior axis of the PSM to generate traveling phase waves (Jörg et al., 2015; Morelli et al., 2009; Ares et al., 2012). It is described as ωx=ω0Ux where ω0 is the frequency at the posterior tip of the tailbud. For simplicity, we scale the frequency profile with the PSM length L. The function U(x) reads, (Figure 2D; Jörg et al., 2015):

(11) Ux=σ+1-σ1-e-kx-/L1-e-k, (11)

where x-=x-xa, σ denotes the difference in the frequency between anterior and posterior ends of the tissue, and k determines the shape of the frequency profile, Figure 2D and Figure 3—figure supplement 10.

In Equation (10), we introduced some simplifications. For example, time delays in intercellular interactions with Delta-Notch signaling play important roles in setting the period of collective rhythms and synchronization (Herrgen et al., 2010; Yoshioka-Kobayashi et al., 2020). However, Equation (10) does not include it to reduce computational time. In addition, a recent study suggested the presence of a spatial gradient of noise intensity along the PSM (Keskin et al., 2018), but we assume that Dθ is constant across the tissue. Wildtype embryos with Delta-Notch signaling do not spontaneously form defective segments, suggesting that coupling strength is sufficiently strong to overcome phase noise. Therefore, we assume that the phase noise intensity is sufficiently lower than coupling strength and approximate its spatial gradient by the zeroth order term.

The phase of cells anterior to the PSM x<xa is arrested (i.e. dθi/dt=0 for xi<xa). Then, we obtain advecting stripes for normal traveling waves in the PSM, Figure 2E. We consider that the stripes represent segment boundaries and the interval between two consecutive stripes represents segment length as described below.

DAPT washout at different time points

Request a detailed protocol

In this study, we allow the coupling strength in Equation (10) to be a time-dependent function κt. The value of the coupling strength is varied in the presence or absence of DAPT in embryos. Besides the influence of DAPT, the coupling strength in embryonic cells may change intrinsically with developmental stages due to, for example, gradual changes in Delta and Notch protein levels on cell membrane, and/or changes in contact surface areas between neighboring cells. To model such changes in the coupling strength, we assume the following time dependence:

(12) κ(t)={0t<twashoutκst+κ0ttwashout,

where twash-out represents the time at which DAPT is washed out in simulations. In the presence of DAPT t < twash-out, there is no interaction between cells. After DAPT washout at t = twash-out, cells restore coupling immediately and interact with each other at finite rates. For simplicity, we assume that the coupling strength changes linearly with time in Equation (12) to consider its intrinsic dependence on developmental stages, Figure 4—figure supplement 5A. Setting κs=0 in Equation (12) describes a constant coupling strength κ0  after DAPT washout. We first analyze resynchronization dynamics with the fixed value of the coupling strength κ0 by setting κs=0, Figures 3 and 4, Figure 3—figure supplements 411 and Figure 4—figure supplements 1, 2 and 4. Then, we consider a positive value of κs>0 and let the value of the coupling strength gradually increase to study its effect on time to FRS, Figure 4E and Figure 4—figure supplements 3, 5, 6 and 8.

PSM shortening

Request a detailed protocol

When we model PSM shortening, we consider that the position of the anterior end of the PSM xa(t) changes over time. For simplicity, we assume that the anterior end of the PSM moves in the posterior direction (x>0) at a constant velocity ua(ua>0), Figure 4—figure supplement 1A:

(13) xat=uat. (13)

Hence, the length of the PSM becomes shorter as L(t)=Lxuat. The speed of the anterior end of the PSM may influence the segment length. For better comparison, we fix the segment length S constant for different values of ua in simulations. For this, we impose the following condition:

(14) va+ua=c, (14)

where va is the cell advection speed at the anterior end used in Equation (2) and c is a constant. Hence, va can be expressed as va=cua. With Equation (14), the segment length S reads S=c×Ta where Ta is the period of oscillation at the position xa. This anterior period of oscillation Ta is the period that one would measure by recording the phase at the tissue boundary at the anterior end as waves – and cells – pass across this boundary, not to be confused with the period of a single cell there. With this definition, this anterior period Ta matches the segmentation rate. When the PSM length does not change during simulation (ua=0), the anterior period Ta is equal to the period at the tailbud, Ta=2π/ω0. When the PSM becomes shorter over time, there is a Doppler effect because waves are traveling toward an approaching anterior boundary (Soroldoni et al., 2014; Jörg et al., 2015). This Doppler effect changes the readout of anterior period Ta, which becomes shorter. To compensate for this effect here, we set Ta=30 min by tuning ω0 in Equation (10) if we include PSM shortening.

Cell influx and outflux

Request a detailed protocol

Due to cell advection, cells reach the anterior end of the PSM x=xa. These cells exit from the domains Ωl and Ωr at the advection speed |vd(xa)|=va. We make the simplifying assumption of constant cell density ϱϱ0. We implement this assumption in simulations by local density dependent cell addition to the PSM and tailbud, as described below:

  1. We measure cell density of the left and right PSM ϱl and ϱr (xa(t)<x<Xc), respectively, and that of the tailbud ϱt (Xc<x). Then, we compute the difference between the measured density and the target density ϱ0, δl=ϱlϱ0, δr=ϱrϱ0 and δt=ϱtϱ0.

  2. If all δl, δr, and δt are positive, we do not add a cell to the PSM and tailbud when a cell leaves the PSM from its anterior end xa(t) due to advection.

  3. In contrast, if some of δl, δr, and δt are negative, we add one cell to the region of which density is smallest. For example, if δl is negative and the smallest, we add a cell to the left PSM. The phase θ of the added cell is assigned randomly from the uniform distribution between 0 and 2π. The position of the added cell is randomly assigned within the chosen domain. In the anterior end of the embryonic PSM, cells divide less frequently (unpublished observations). For this reason, we do not add cells in the region xxa+ζ for the left and right PSM. In this study, we set ζ = 100 μm. We note that adding cells with random phases in this region would be detrimental for the segmented pattern, given that phase disturbances would not have time to synchronize to their neighbors. The cell polarity angles of the added cell are assigned randomly from the uniform distribution between 0 and 2π for φ, and between 0 and π for ϕ. The autonomous frequency, speed of intrinsic cell movement and advection speed for the added cell are determined depending on its added position.

Change in the PSM radius

Request a detailed protocol

We examine how change in the PSM radius over time influences resynchronization, Figure 4—figure supplement 4. For simplicity, we assume that the PSM radius r decreases uniformly in the U-shaped domain and linearly over time:

(15) rt=-srt+r0,

where sr (sr>0) is the magnitude of change in the PSM radius. Because we fix cell density as described in the subsection ‘cell influx and outflux’, the number of cells in the U-shaped domain decreases as the PSM radius becomes narrower. Note that although the thinning of PSM in one direction may cause a length extension in other directions to preserve a volume, we do not consider such extension for simplicity.

To simplify the implementation of the model, we separately update cell positions xi and the radius size r(t). When updating r(t) into r(t+Δt) by Equation (15) above, some cells may be left outside the U-shaped domain, ri(t+Δt)>r(t+Δt) where ri(t+Δt) is the position of cell i in the radial direction defined in Equation (6) or (Equation 8) at time t+Δt. For such cells, we correct their positions after the update of PSM radius as: ri(t+Δt)r(t+Δt)2rb where rb is the lengthscale of boundary force used in Equations (7) and (9).

Initial condition of simulations

Request a detailed protocol

We set the initial position of the anterior end of the PSM xa(0) as xa(0)=0. We choose the cell number N(0) to satisfy the cell density ϱ=ϱ0 with a given tissue geometry. Initial positions of cells xi(0)(i=1,2,,N) are within the PSM and tailbud domains. To set initial positions of cells, we first randomly locate cells within the U-shaped domain. Then, local stresses caused by the random cell positioning are relaxed by simulating Equation (1) without the advection term for 10 min with the integration time step of 0.01 min. In this relaxation process, we include a boundary force at x=xa=0 to constrain cells within the U-shaped domain. For simulations of embryos without DAPT treatment (i.e. synchronized initial condition), we use a synchronized initial condition θi(0)=3π/2 for i=1,,N. For resynchronization simulations, the initial phase value of cell i, θi(0), is chosen randomly from a uniform distribution between 0 and 2π. The values of initial cell polarity angles in Equation (4) are also chosen randomly from the uniform distributions between 0 and π for ϕi(0), and between 0 and 2π for φi(0). The autonomous frequency, speed of intrinsic cell movement and advection speed are determined by the initial position xi(0).

Parameter values

Request a detailed protocol

The values of parameters used in this study are listed in Supplementary files 1 and 2. Cell density and diameter of cells are based on the estimation by Uriu et al., 2017. The values of parameters that determine intrinsic cell movement and physical forces between cells are set to reproduce experimental data, Figure 3—figure supplement 4B (also refer to Uriu et al., 2017 for details). We use the PSM length within the range reported in Soroldoni et al., 2014. The parameters for the frequency profile in Equation (10) are based on Soroldoni et al., 2014; Jörg et al., 2015. Also, the value of the coupling strength κ in Equation (10) is in the range estimated by Riedel-Kruse et al., 2007; Herrgen et al., 2010.

Previous studies estimated noise intensity in the segmentation clock employing different theoretical formalisms (Riedel-Kruse et al., 2007; Keskin et al., 2018; Jenkins et al., 2015). The present physical model has two noise sources for phase of oscillation. One is the white Gaussian noise with intensity Dθ in the phase Equation (10). The other is cell addition with a random phase value described in the section 'Cell influx and outflux'. In desynchronization simulations with κt=0 in Equation (10), noise by cell addition alone (Dθ=0) ruins kinematic phase waves and stripe patterns, Figure 3—figure supplement 9A,D. First five stripes are typically recognizable with Dθ=0, Figure 3—figure supplement 9A. With the increase in Dθ, the local phase order decays more quickly and stripes are lost earlier, Figure 3—figure supplement 9B–D. Because ALD with a saturated dose of DAPT is around five in current experiments, we set Dθ=0.0013 throughout this study, Figure 3—figure supplement 9B. In resynchronization simulations with κt=κ0=0.07, both FRS and PLD do not depend on Dθ within its examined rage, Figure 3—figure supplement 9E. Hence, even if a different value of Dθ is used in resynchronization simulations, only a slight modification of the parameter values would be enough to better fit experimental FRS and PLD. Qualitative aspects of FRS and PLD in simulations do not depend on the value of Dθ.

Calculation of local phase order

Request a detailed protocol

To measure the level of synchrony at local domains along the anterior-posterior axis of the PSM, we introduce a local phase order parameter (Shiogai and Kuramoto, 2003). It is based on the Kuramoto phase order parameter (Kuramoto, 1984) with some modifications in computing average over cells to capture the presence of spatial phase waves in the PSM. For the left PSM Ωl, we first compute phase order for cells within a thin slice domain Ωmx=x+m-1Δx, x+mΔx×0, 2r×0, 2r for m=1,2,,M by:

(16) Zmt,xeiΨmt,x=1nmxjΩmxeiθjt, (16)

where nm is the number of cells within the domain Ωm(x), Figure 3—figure supplement 1A. We set Δx equal to the cell diameter Δx=dc so that Equation (16) measures the phase order for cells in same position in the anterior-posterior axis. Zm indicates the level of local synchrony of this slice domain. If the phases of cells in the domain are synchronized, Zm is close to one. In contrast, if they are not synchronized, Zm is close to zero. Ψm indicates the mean phase of the cells in the slice domain.

This definition of Zm measures local order with a high spatial resolution along x-axis. However, only a few cells are present in each slice, introducing finite size fluctuations. Thus, we define the local phase order parameter at position x by taking the average of Zm over a number M of these domains, to smooth out small finite size fluctuations:

(17) Zt,x=1Mm=1MZmt,x. (17)

The segment length in the anterior-posterior direction in our parameter sets is 5Δx. Therefore, we set M=5 in the calculation of the local phase order Equation (17). We calculate the local phase order parameter in a similar manner for the right PSM Ωr. In this case, the domain Ωm(x) in Equation (16) is Ωmx=x+m-1Δx, x+mΔx×2R, 2R+r×0, 2r.

Note that computing local order in thin slices first as in Equation (16), and then averaging as in Equation (17), we can capture high local order values even in the presence of a phase gradient. Computing a local order parameter directly in a thicker slab domain of width 5Δx would result in low values of the order parameter in the presence of a phase gradient as in the anterior PSM.

Definition of a normal segment boundary in simulations

In computational simulations, we use the phase order parameter Equation (17) to define normal and defective segment boundaries. We first define a segment boundary by considering simulations for wild-type embryos untreated with DAPT. Such simulations are started from a synchronized initial condition. Then, we describe how to detect normal segment boundaries in resynchronization simulations started from random initial conditions.

In a untreated embryo where cells in the PSM are locally synchronized, kinematic phase waves can be observed across the tissue. In such embryos, the position of a new segment boundary is specified when a wave of gene expression arrives in the anterior end of the PSM. Namely, a position of a segment boundary is determined when the phase of cells near the anterior end of the PSM attains a certain value ϑ, Figure 3—figure supplement 1B,C.

Based on this, we monitor the mean phase at the anterior end of the PSM xa to detect when a segment boundary position is set in simulations, Figure 3—figure supplement 1A–C. We compute the mean phase Ψ1(t,xa) at position xa at time t for Ω1xa=xa, xa+Δx×0, 2r×0, 2r for the left PSM and Ω1xa=xa, xa+Δx×2R, 2R+r×0, 2r for the right PSM by using Equation (16) (m=1). As noted in Equation (16), we set Δx as the cell diameter Δx=dc.

We then detect a time τi (i=1,2,...) that satisfies:

(18) Ψ1τi,xa=ϑ, (18)

where ϑ is a constant that we set ϑ=3π/2 without loss of generality in this study, Figure 3—figure supplement 1C. For simulations for control embryos where DAPT is not added and, therefore, the level of synchrony is high, τi should be the time when the position of the segment boundary i is determined, Figure 3—figure supplement 1B,C.

Identification and numbering of defective segment boundaries

Request a detailed protocol

For resynchronization simulations starting from random initial conditions, we modify the above procedure as follows. After detecting the time τi when the mean phase of anterior cells becomes ϑ, we check the local phase order parameter Z(t,xa) defined in Equation (17) to determine whether these cells can form a normal boundary, Figure 3—figure supplement 1D. We define that the anterior cells can form a normal boundary at time τi if:

(19) Z(t,xa)Zc for τiTa/2tτiTa/2+η,

where Ta is the period of oscillation at the anterior end of the PSM xa as described in the section “PSM shortening”. Since Z(τiTa/2,xa) is the average local phase order across nearly one segment length at τiTa/2, it evaluates the integrity of the segment boundary and its neighboring inter-boundary regions. To suppress a false detection of a normal segment boundary caused by the fluctuation of Z(t,xa), we monitor Z(t,xa) in a short time interval with a window size η in Equation (19). We set η=4 min in Equation (19). By visual inspection of stripe patterns in simulations, we set Zc=0.85 for the recovery simulations throughout this paper, see Figure 3F and Figure 3—figure supplement 1D,E. Note that this threshold value is simply for detecting a normal segment boundary in simulations. It may be different from the critical value of the order parameter for normal segment boundary formation in actual embryonic tissues.

If Equation (19) is satisfied for τi, we then specify the segment boundary number. Note that the subindex i of τi does not specify the segment number in resynchronization simulations due to the fluctuation of the average phase Ψ1 for earlier time when cells are not synchronized yet, Figure 3—figure supplement 1D. If the previous anterior cell population that satisfied Ψ1τi-1,xa=ϑ at time τi1 (τi1<τi) was also satisfied Equation (19) and numbered as segment boundary j, the current one is numbered as j+1. If not, we infer the segment boundary number based on τi and anterior period Ta. We assign the current segment boundary with the expected segment number:

(20) {τi/Ta,ifτi/Taτi/Ta<Δ,τi/Ta1,otherwise,

where χ represents rounding of real number χ to the closest integer value. Note that we assign the number τi/Ta when the phases of anterior cells become ϑ slightly earlier, τi>Taτi/TaΔTa considering the fluctuation after resynchronization, Figure 3—figure supplement 1F. In this study, we set Δ=0.3 in Equation (20). After detecting the normal segment boundaries and identifying their numbers in this way, we assign all the remaining segment boundaries to be defective.

Definition of FRS and PLD in simulations

Request a detailed protocol

FRS is defined as jf1 where jf is the smallest segment boundary number that was determined to be normal. The subtraction of 1 from jf is to match the definition of FRS for experimental data, see the section ‘Whole-mount in situ hybridization and segmental defect scoring’. PLD is defined as follows. We find the minimum segment boundary number jp above which all the segment boundaries, including jp, are normal. Then, PLD is jp1. In the example simulation shown in Figure 3—figure supplement 1D,E, FRS is 9 and PLD is 13.

Calculation of phase vorticity

Request a detailed protocol

Vortices are regions in space where the phase values circulate from 0 to 2π around some point. Vorticity could be detected taking a closed path around cells and computing the accumulating change in the phase of neighboring cells around it. However, it is challenging to detect vorticity in a tissue where phase is not continuous in space, but only defined at points where cells are. Besides, there are phase fluctuations that can introduce local variations of phase change. Therefore, here we discretize the closed path in angular steps and average the phase over the resulting domains, Figure 3—figure supplement 2. The phase within these domains may grow linearly from 0 to 2π when one turns around a position close to the center of a vortex. Below, we describe the definition of vorticity that is shown in Figure 3C, Figure 3—figure supplements 5 and 6 and Figure 4—figure supplements 1, 2 and 6.

Vortex axis can have different spatial orientations. To detect vortices with different rotating axis, we set several planes in the space and compute phase vorticity at each plane. Then, we project phase vorticity on x-axis to obtain its trajectory along the anterior-posterior axis of the PSM. We first choose either the left or right PSM for the calculation of vorticity. Subsequently, we consider the four slices in the PSM, Figure 3—figure supplement 2A,B. These slices are two z-slices located at z=0 μm and z=2r20 μm (slices 1 and 2 and Figure 3—figure supplement 2A), and two y-slices located at y=y0 μm and y=y0+2r20 μm, (slices 3 and 4, Figure 3—figure supplement 2B), where r is the radius of the PSM. For the left PSM, y0=0 μm, while for the right PSM y0=2R μm where R is the tailbud torus radius. The thickness of these four slices is 20 μm (~2 cell diameters, compare to the 50 μm of PSM diameter). We then project the phase values of cells in each slice to 2D planes, Π(α) (α=1,2,3,4). The two x-y planes Π(1) and Π(2) obtained by the projection of the two z-slices (slices 1 and 2) are used to detect vortices with a rotating axis parallel to z-axis (Figure 3—figure supplement 2A,C). For instance, the x-y plain Π(1) for the right PSM contains cells within the z-slice 0,Lx×2R,2R+2r×0, 20. The two x-z planes Π(3) and Π(4) obtained by projection of the y-slices (slices 3 and 4) are used to detect vortices with a rotating axis parallel to y-axis, Figure 3—figure supplement 2B. We hardly observed phase vortices with the rotating axis parallel to x-axis in simulations. Therefore, we do not consider y-z planes in the calculation of vorticity.

Next, we set grids (xs,yu) in the planes Π(1) and Π(2) where xs=xa+sΔx (s=0,1,2,...) and yu=y0+uΔy (u=0,1,2,...) with the grid size Δx and Δy. Similarly, we set grids (xs,zu) to the planes Π(3) and Π(4) where xs=xa+sΔx and zu=uΔz. We chose Δx=5 μm and Δy=Δz=2 μm. For each grid point in the plane Π(α), we compute vorticity ψ as follows. Below, we explain the case of Π(2) with the grid (xs,yu), Figure 3—figure supplement 2A,C. Same calculations were performed in the other planes as well.

We set a circular ring for the grid point (xs,yu) in the plane:

(21) δl (xxs)2+(yyu)2<l+δl

as shown in Figure 3—figure supplement 2D, E. We set δl x-xs2+y-yu2<l+δl and δl=5.5 μm. The circular ring is subdivided into six domains Vi (i=0,1,2,...,5) by angles π/3 measured counterclockwise from the x-axis. We then compute average phase l=14 μm over cells within each subdomain Vi, θ-i where ni is the number of cells in Vi. If there is no cell within one of the subdomains Vi (ni=0), we do not compute the vorticity for the grid point (xs,yu) and set ψ=0.

To detect a vortex with clock-wise rotation, we permutate θ¯i based on their values (Figure 3—figure supplement 2D, E): θ-i where k = arg min{θ¯i: i = 0,...,5}, θ-i, θ^1=θ-k+1,..., and θ^2=θ-k+2, (mod 6). For a vortex with counter clock-wise rotation, we permutate θ^5=θ-k+5 as: θ-i where k = arg min{θ^0=θ-k: i = 0,..., 5}, θ-i, θ^1=θ-k-1,..., and θ^2=θ-k-2, where a negative value of kj (j=1,2,...,5) should be replaced as −1 → 5, −2 → 4, …, and −5 → 1.

We assume that when a phase vortex is present near the grid point (xs,yu), θ^5=θ-k-5 for the grid point increases linearly with i, Figure 3—figure supplement 2D, F. To detect this linear increase of θ^i, we compute the correlation coefficient θ^i defined as:

(22) α

where 5/2=i=05i/6, 5/2=i=05i/6, θ^-=i=05θ^i/6, and σi=i=05i-5/22/6=35/12. A value of the correlation coefficient σθ^=i=05θ^i-θ^-2/6 close to one means that the phase increases linearly along a perimeter of a circle, indicating the existence of a phase vortex, Figure 3—figure supplement 2D,F. If the correlation coefficient is larger than a threshold α, we consider that the phase value consistently increases along the circumference of the ring and rotates along the z-axis. In this case we define vorticity for the grid point as ψ(xs,yu)=(θ^5θ^0)/2π, Figure 3—figure supplement 2F. If xs,yu=θ^5-θ^0/2π, we set ψ(xs,yu)=0 to exclude false positive detection of a vortex by fluctuation of phase values. We used xs,yu=0 throughout the article. After calculating vorticity for each grid point, we project ψ(xs,yu) to x-axis. We use maximum projection, ψ(2)(xs)=maxyu ψ(xs,yu), Figure 3—figure supplement 2H.

By performing same calculations for the remaining three planes, we obtain {ψ(1)(xs),ψ(2)(xs),ψ(3)(xs),ψ(4)(xs)} for position xs, top panel of Figure 3—figure supplement 2I. Finally, we take their maximum value and adopt it as the vorticity at the position xs: ψ(max)(xs)=maxα ψ(α)(xs), bottom panel of Figure 3—figure supplement 2I. For visualization, we make a density plot as a kymograph by using the data (xs,t,ψ(max)).

Quantification of single and double defects

In both experiments and simulations, we sometimes observe that a defective segment boundary appears either only left or right side of an embryo, which we refer to as a single defect. We also observe another instance where left and right boundaries are both defective, which we call double defect. To examine whether the theory can account for the emergence of single defects in embryonic experiments, we compute the fraction Fs of single defects defined below and compare it between simulations and experiment.

Experimental data

Request a detailed protocol

We first counted the total number Nt of defective segment boundary loci for an embryo. We then counted the number of single defects Ns and computed the fraction Fs=Ns/Nt. When we compared the experimental data with simulation data, Figure 4—figure supplement 8D, we measured Nt and Ns after segment 9 that marks the onset of resynchronization, Figure 1—figure supplement 1B.

Simulation data

Request a detailed protocol

We defined normal and defective segment boundaries based on the local phase order at the anterior end of the PSM Z(t,xa) as described in the previous section ‘Definition of a normal segment boundary in simulations’. For single realizations of simulation, we counted the total number of defective segment boundary loci Nt and the number of single defects Ns appeared posterior to the segment 9 as in experimental data. Then, we computed the fraction, Fs=Ns/Nt, Figure 4—figure supplement 8D.

Implementation for numerical simulations

Request a detailed protocol

We solved Equations (1) and (10) with the Euler-Maruyama method with the time step for integration δt=0.01 min. Custom simulation codes were written in C language (Source code 1). Videos of numerical simulations, calculations of local phase order and vorticity, and analysis of left-right segment boundary defects were done with custom Mathematica (Wolfram) codes (Source code 1).

Appendix 1

Segment statistics from the spatial distribution of defective segments

The distribution of defects along the embryonic axis can be parametrized in different ways. Here, we introduce complementary pictures, one relies on the fraction of left and right defects, and the other on fractions of single and double defects. Using a random defect hypothesis we show how to relate these pictures with probability theory. Taking into account the spatial distribution of defects we compute ALD, PLD and FRS from these statistics, and we predict and test the values of the fraction of single defects per embryo.

Segment state variables

We introduce a state variable that accounts for the presence or absence of a defective segment Sik(x), where i labels the embryo, k={l,r} labels the side of the embryo and x is the segment locus along the axis. We consider N embryos with M+1 segment boundaries, so i=1,2,...,N and x=0,1,...,M. The state variable Sik(x) takes the values zero and one depending on whether the segment is normal or defective.

Segment defect distribution on both sides of the embryo

When a segment locus is defective on the left (right) side of the embryo we call it a left (right) defect independently of the state of the other side. We define left and right defect distributions taking the population average of segment state variables,

(23) plx=Sil(x)=1Ni=1NSil(x), (23)

and

(24) prx=Sirx=1Ni=1NSirx. (24)

These spatial distributions of defective segments on the left and right sides of embryos are very similar, for different DAPT washout timing, Figure 4—figure supplement 7A. This agreement between left and right distributions supports the assumption plx=prx=px. In the following sections, we use probability theory to compute ALD, PLD and FRS from this spatial distribution of defective segments.

Probabilistic calculation of ALD

Let qaxa be the probability to find an ALD at position xa. qaxa can be expressed in terms of p(x) as:

(25) qa(xa)={p(xa)xa=0ξ=0xa1(1p(ξ))×p(xa)xa>0. 

The first factor ξ=0xa-11-pξ in the second line represents the probability of normal segment boundaries from position x=0 to x=xa1. The second factor pxa is the probability of a defective segment boundary at position x=xa. The resulting ALD distribution qaxa presents a clear peak at the onset of the defective region, Figure 4—figure supplement 7B. The ALD is then calculated as the mean value for this probability distribution,

(26) ALD=xa=0Mxaqaxa. (26)

Probabilistic calculation of PLD

Let qpxp be the probability of PLD at position xp. It can be written as:

(27) qp(xp)={p(xp)×ξ=xp+1M(1p(ξ))xp<Mp(xp),xp=M.

The first factor pxp in the first line is the probability of a defective boundary at xp. The second factor ξ=xp+1M1-pξ represents the probability that all the remaining segment boundaries posterior to xp are normal. The resulting PLD distribution qpxp peaks at the end of the defective region, Figure 4—figure supplement 7B. The PLD can be written as the mean value for this distribution,

(28) PLD=xp=0Mxpqpxp. (28)

Probabilistic calculation of FRS

Occurrence of a recovered segment is conditioned to the previous occurrence of defective segments. In this study, we measure the FRS after the desynchronization phase. The desynchronization phase is determined based on the distribution of defective segments, Figure 1—figure supplement 1B. Suppose that the desynchronization phase ends by the formation of segment sd. We define FRS as the first normal segment after sd-1. This is the definition we use to measure the embryonic FRS. With this definition of FRS, the probability qf(xf) of the first normal segment boundary at locus xf is:

(29) qf(xf)={1p(xf),xf=sd ξ=sdxf1p(ξ)×(1p(xf)),xf>sd.

The first factor of the second line represents the probability for all the segment boundaries between sd and xf-1 to be defective. The second factor is the probability for a normal segment boundary to form at xf.

To compute qfxf, we set sd=9 as in the main text, Figure 1—figure supplement 1B. The resulting distribution qfxf has a peak that precedes that of PLD and partly overlaps with it, Figure 4—figure supplement 7B. From these results, the FRS then can be expressed as

(30) FRS=xf=sdMxfqfxf-1. (30)

In Equation (30), we subtract one because FRS is defined with the anterior boundary of the first normal segment after sd-1, see definition of FRS for experimental data in Materials and methods. qa, qp and qf calculated with Equations (25), (27) and (29) agree well with the direct measurement of these distributions, Figure 4—figure supplement 7B. Furthermore, expressions obtained for ALD, PLD and FRS from the spatial distribution of defects p(x), Equations (26), (28) and (30), are in very good agreement with direct measurements of these quantities, Figure 4—figure supplement 7C.

The spatial distribution of single and double defects

In a complementary framework, we introduce double defects, single defects and normal segments. Double defects occur when both sides of the embryo at a given locus x are defective. Single defects occur when at a given locus there is a defect either left or right, but not on the other side. Normal segments have no defects on either side. Taking population averages, we can define the double defect distribution

(31) pdx=Sil(x) Sir(x)=1Ni=1NSil(x) Sir(x), (31)

the single defects distribution

(32) psx=1Ni=1NSilx 1-Sirx+1-Silx Sirx, (32)

and the normal segments distribution

(33) pn(x)=1Ni=1N1-Silx1-Sirx. (33)

These distributions can be related to the left/right framework introduced above. Key to this is that the two sides of the anterior PSM are physically unconnected, separated by another tissue –the notochord. It has been shown that interfering with Retinoic Acid, which controls somitogenesis bilateral symmetry by shielding asymmetric cues, results in asymmetric left/right segmentation and clock waves (Vermot et al., 2005). This indicates that segmentation clock oscillations are independent in the left and right sides of the PSM. In the physical model, we tacitly assume this independence since there is no coupling between oscillators on one side and the other. A consequence of this left/right independence should be a vanishing covariance of the segment defect variables at opposite sides of the embryo (Gardiner, 2009),

(34) Sil(x) Sir(x) -SilxSirx=0. (34)

The two terms in this covariance can be written using the segment defect distributions of the two frameworks,

(35) pdx=plxprx. (35)

Similarly, we obtain for single defects

(36) psx=plx 1-prx+1-plx prx, (36)

and for normal segments

(37) pnx=1-plx1-prx. (37)

We can further simplify these expressions with the assumption plx=prx=px. The axial distribution of double defects pdx gradually grows to a plateau at the onset of the defective region and then decays at its end, while the distribution of single defects psx peaks both at the onset of the defective region and at its end, Figure 4—figure supplement 8A. The good agreement observed between direct measurement of pd(x) and ps(x) and results obtained from p(x) together with probabilistic arguments, provides a test for the vanishing of the covariance and left/right independence, Figure 4—figure supplement 8A.

Fraction of single defects from the axial distribution of defects

As a further test of the reach of the defect distribution p(x), we use it to compute the fraction of single defects in a population of embryos. We first count the number of double and single defects in a given embryo i

(38) Nid=x=0MSilxSirx, (38)

and

(39) Nis=x=0MSilx 1-Sirx+1-Silx Sirx. (39)

The total number of defects in embryo i is Nit=Nis+Nid. Then, we take the population average of these quantities. The average number of double defects in the population is

(40a) Nid=1Ni=1NNid=1Ni=1Nx=0MSil(x) Sir(x)=x=0M1Ni=1NSil(x) Sir(x), (40a)

so

(40b) Nid=x=0MSil(x) Sir(x)=x=0Mpd(x)=x=0Mplxprx, (40b)

using left/right independence. Similarly, the average number of single defects Nis is

(41) Nis=x=0Mpsx=x=0Mplx 1-prx+1-plx prx. (41)

We define the fraction of single to total defects in an embryo

(42) Fis=Nis/Nit, (42)

and its population average

(43) Fs=Fis=Nis/Nit. (43)

If the coefficient of variation of the total number of defects in the population is small, we can approximate

(44a) Fs=Nis/NitNis/Nit, (44a)

so

(44b) FsNisNit=x=0Mplx 1-prx+1-plx prxx=0M1 -1-plx 1-prx . (44b)

Using the left/right symmetry of defect distributions plx=prx=px, we can further simplify this expression

(45) FsNisNit=2x=0Mpx 1-pxx=0M1 -1-px2 . (45)

The good agreement observed between this probabilistic calculation and direct measurement of Fs counting single defects in individual embryos, Figure 4—figure supplement 8C, suggests that the recovery of segments after DAPT washout occurs independently between left and right sides of the PSM.

Finally, the physical model of the PSM reproduces both the single and double defect distributions, Figure 4—figure supplement 8B, and the dependence of Fs on DAPT washout timing, Figure 4—figure supplement 8D.

Conclusion

Taken together, these results suggest that knowledge of p(x) is enough to compute some of the key observables that we use to quantify segment recovery, such as ALD, PLD, FRS and Fs. Together with the result that the physical model of the PSM reproduces p(x) at different washout timings, Figure 4F, this is evidence for the breadth of the physical model results and predictions.

Data availability

All data generated or analyzed during this study are included in the manuscript and supporting files.

References

  1. Book
    1. Gardiner C
    (2009)
    Stochastic Methods: A Handbook for the Natural and Social Sciences
    Springer-Verlag.
    1. Haddon C
    2. Smithers L
    3. Schneider-Maunoury S
    4. Coche T
    5. Henrique D
    6. Lewis J
    (1998)
    Multiple Delta genes and lateral inhibition in zebrafish primary neurogenesis
    Development 125:359–370.
    1. van Eeden FJ
    2. Granato M
    3. Schach U
    4. Brand M
    5. Furutani-Seiki M
    6. Haffter P
    7. Hammerschmidt M
    8. Heisenberg CP
    9. Jiang YJ
    10. Kane DA
    11. Kelsh RN
    12. Mullins MC
    13. Odenthal J
    14. Warga RM
    15. Allende ML
    16. Weinberg ES
    17. Nüsslein-Volhard C
    (1996)
    Mutations affecting somite formation and patterning in the zebrafish, Danio rerio
    Development 123:153–164.

Article and author information

Author details

  1. Koichiro Uriu

    Graduate School of Natural Science and Technology, Kanazawa University, Kakuma-machi, Kanazawa, Japan
    Contribution
    Conceptualization, Software, Formal analysis, Funding acquisition, Investigation, Visualization, Methodology, Writing - original draft, Writing - review and editing
    Contributed equally with
    Bo-Kai Liao
    For correspondence
    uriu@staff.kanazawa-u.ac.jp
    Competing interests
    No competing interests declared
    ORCID icon "This ORCID iD identifies the author of this article:" 0000-0003-1802-2470
  2. Bo-Kai Liao

    1. Department of Aquaculture, National Taiwan Ocean University, Keelung, Taiwan
    2. Department of Cell and Developmental Biology, University College London, Gower Street, London, United Kingdom
    3. The Francis Crick Institute, London, United Kingdom
    4. Max Planck Institute of Molecular Cell Biology and Genetics, Dresden, Germany
    Contribution
    Conceptualization, Resources, Funding acquisition, Validation, Investigation, Visualization, Writing - original draft, Writing - review and editing
    Contributed equally with
    Koichiro Uriu
    Competing interests
    No competing interests declared
    ORCID icon "This ORCID iD identifies the author of this article:" 0000-0003-4742-4487
  3. Andrew C Oates

    1. Department of Cell and Developmental Biology, University College London, Gower Street, London, United Kingdom
    2. The Francis Crick Institute, London, United Kingdom
    3. Max Planck Institute of Molecular Cell Biology and Genetics, Dresden, Germany
    4. Institute of Bioengineering, École polytechnique fédérale de Lausanne (EPFL), Lausanne, Switzerland
    Contribution
    Conceptualization, Supervision, Funding acquisition, Visualization, Writing - original draft, Writing - review and editing
    Competing interests
    No competing interests declared
    ORCID icon "This ORCID iD identifies the author of this article:" 0000-0002-3015-3978
  4. Luis G Morelli

    1. Instituto de Investigación en Biomedicina de Buenos Aires (IBioBA) – CONICET – Partner Institute of the Max Planck Society, Polo Científico Tecnológico, Buenos Aires, Argentina
    2. Departamento de Física, FCEyN UBA, Ciudad Universitaria, Buenos Aires, Argentina
    3. Max Planck Institute for Molecular Physiology, Department of Systemic Cell Biology, Dortmund, Germany
    Contribution
    Conceptualization, Formal analysis, Funding acquisition, Visualization, Methodology, Writing - original draft, Writing - review and editing
    Competing interests
    No competing interests declared
    ORCID icon "This ORCID iD identifies the author of this article:" 0000-0001-5614-073X

Funding

Japan Society for the Promotion of Science (KAKENHI grant number 17H05762)

  • Koichiro Uriu

Japan Society for the Promotion of Science (KAKENHI grant number 19H04955)

  • Koichiro Uriu

Japan Society for the Promotion of Science (KAKENHI grant number 19H04772)

  • Koichiro Uriu

Ministry of Science and Technology, Taiwan (MOST 108-2311-B-019-001-MY3)

  • Bo-Kai Liao

SNSF (31003A_176037)

  • Andrew C Oates

Wellcome Trust (WT098025MA)

  • Andrew C Oates

European Research Council (ERC-2007-StG: 207634)

  • Bo-Kai Liao
  • Andrew C Oates

Francis Crick Institute

  • Bo-Kai Liao
  • Andrew C Oates

Agencia Nacional de Promoción Científica y Tecnológica (PICT 2012 1954)

  • Luis G Morelli

Agencia Nacional de Promoción Científica y Tecnológica (PICT 2013 1301)

  • Luis G Morelli

Agencia Nacional de Promoción Científica y Tecnológica (PICT 2017 3753)

  • Luis G Morelli

FOCEM-Mercosur (COF 03/11)

  • Luis G Morelli

Japan Society for the Promotion of Science (Short Term Grant S17064)

  • Koichiro Uriu
  • Luis G Morelli

The funders had no role in study design, data collection and interpretation, or the decision to submit the work for publication.

Acknowledgements

We thank the fish facility staff of the MPI-CBG, Dresden and Arianne Bercowsky, Gabriela Petrungaro, Sundar Naganathan, and Olivier Venzin for comments on the manuscript. Research was sponsored by JSPS KAKENHI grant number 17H05762, 19H04955, 19H04772 to KU; Ministry of Science and Technology, Taiwan (MOST 108–2311-B-019–001-MY3) to B-KL; SNSF Project funding division III (31003A_176037) and Wellcome Trust Senior Research Fellowship in Basic Biomedical Science (WT098025MA) to AO; European Research Council Starting Independent Research Grant (ERC-2007-StG: 207634) and the Francis Crick Institute to B-KL and AO; ANPCyT PICT 2012 1954, PICT 2013 1301, PICT 2017 3753 to LGM, FOCEM-Mercosur (COF 03/11) to IBioBA; and JSPS Short Term Grant S17064 to KU and LGM.

Ethics

Animal experimentation: Zebrafish experimentation was carried out in strict accordance with the ethics and regulations of the Saxonian Ministry of the Environment and Agriculture in Germany under licence Az. 74-9165.40-9-2001, and the Home Office in the United Kingdom under project licence PPL No. 70/7675.

Version history

  1. Received: July 23, 2020
  2. Accepted: January 27, 2021
  3. Accepted Manuscript published: February 15, 2021 (version 1)
  4. Version of Record published: March 22, 2021 (version 2)

Copyright

© 2021, Uriu et al.

This article is distributed under the terms of the Creative Commons Attribution License, which permits unrestricted use and redistribution provided that the original author and source are credited.

Metrics

  • 1,395
    Page views
  • 237
    Downloads
  • 6
    Citations

Article citation count generated by polling the highest count across the following sources: Crossref, PubMed Central, Scopus.

Download links

A two-part list of links to download the article, or parts of the article, in various formats.

Downloads (link to download the article as PDF)

Open citations (links to open the citations from this article in various online reference manager services)

Cite this article (links to download the citations from this article in formats compatible with various reference manager tools)

  1. Koichiro Uriu
  2. Bo-Kai Liao
  3. Andrew C Oates
  4. Luis G Morelli
(2021)
From local resynchronization to global pattern recovery in the zebrafish segmentation clock
eLife 10:e61358.
https://doi.org/10.7554/eLife.61358

Share this article

https://doi.org/10.7554/eLife.61358

Further reading

    1. Developmental Biology
    2. Neuroscience
    Guangqin Wang, Yunpeng Gu, Zhiyong Liu
    Research Article

    Mammals harbor a limited number of sound-receptor hair cells (HCs) that cannot be regenerated after damage. Thus, investigating the underlying molecular mechanisms that maintain HC survival is crucial for preventing hearing impairment. Intriguingly, Pou4f3-/- or Gfi1-/- HCs form initially but then rapidly degenerate, whereas Rbm24-/- HCs degenerate considerably later. However, the transcriptional cascades involving Pou4f3, Gfi1, and Rbm24 remain undescribed. Here, we demonstrate that Rbm24 expression is completely repressed in Pou4f3-/- HCs but unaltered in Gfi1-/- HCs, and further that the expression of both POU4F3 and GFI1 is intact in Rbm24-/- HCs. Moreover, by using in vivo mouse transgenic reporter assays, we identify three Rbm24 enhancers to which POU4F3 binds. Lastly, through in vivo genetic testing of whether Rbm24 restoration alleviates the degeneration of Pou4f3-/- HCs, we show that ectopic Rbm24 alone cannot prevent Pou4f3-/- HCs from degenerating. Collectively, our findings provide new molecular and genetic insights into how HC survival is regulated.

    1. Developmental Biology
    2. Evolutionary Biology
    Jaaved Mohammed, Neha Arora ... Joanna Wysocka
    Research Article

    Genome-wide association studies (GWAS) identified thousands of genetic variants linked to phenotypic traits and disease risk. However, mechanistic understanding of how GWAS variants influence complex morphological traits and can, in certain cases, simultaneously confer normal-range phenotypic variation and disease predisposition, is still largely lacking. Here, we focus on rs6740960, a single nucleotide polymorphism (SNP) at the 2p21 locus, which in GWAS studies has been associated both with normal-range variation in jaw shape and with an increased risk of non-syndromic orofacial clefting. Using in vitro derived embryonic cell types relevant for human facial morphogenesis, we show that this SNP resides in an enhancer that regulates chondrocytic expression of PKDCC - a gene encoding a tyrosine kinase involved in chondrogenesis and skeletal development. In agreement, we demonstrate that the rs6740960 SNP is sufficient to confer chondrocyte-specific differences in PKDCC expression. By deploying dense landmark morphometric analysis of skull elements in mice, we show that changes in Pkdcc dosage are associated with quantitative changes in the maxilla, mandible, and palatine bone shape that are concordant with the facial phenotypes and disease predisposition seen in humans. We further demonstrate that the frequency of the rs6740960 variant strongly deviated among different human populations, and that the activity of its cognate enhancer diverged in hominids. Our study provides a mechanistic explanation of how a common SNP can mediate normal-range and disease-associated morphological variation, with implications for the evolution of human facial features.