key: cord-0016242-w5kqv3du
authors: Uriu, Koichiro; Liao, Bo-Kai; Oates, Andrew C; Morelli, Luis G
title: From local resynchronization to global pattern recovery in the zebrafish segmentation clock
date: 2021-02-15
journal: nan
DOI: 10.7554/elife.61358
sha: c6c0a05a5e748c2ddf8a5bec553719c5061dda7a
doc_id: 16242
cord_uid: w5kqv3du

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.

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 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; . 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 , 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 Z c 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 Z c . The time point at which the level of synchrony crosses Z c 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) .

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 Z c , 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 Z c 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 Z c , Figure 1A , as expected from the 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 Z c , 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) Figure 1 continued on next page 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.

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 . 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 Figure 1 continued 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. The online version of this article includes the following source data and figure supplement(s) for figure 1:

Source data 1. Segment boundary defects in embryos with different DAPT washout timing. 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 (t washout =~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 Z c , it remained larger than Z c , 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 t wash-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 Z c in this intermingled region, and persistent fluctuations in synchrony level lead to intermittent defective boundary formation.

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 meanfield 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 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:

is the physical contact force between cells i and j, N is the total cell number and F b x i ð Þ 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:

where i is the phase of cell i, ! x i ð Þ is the autonomous frequency at position x i , k is the coupling strength, n i is the number of neighboring cells interacting with cell i, d c is the cell diameter, D is the phase noise intensity and i t ð Þ 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 x a , Figure 2A ,E. The arrested phase stripes in the region x < x a 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 < x a 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 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 (Z c ) 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 Z c 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 2p 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. 

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 Table 1 . See Materials and methods for the quantification of FRS and PLD in simulations.

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, 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 k 0 influenced both FRS and PLD, Figure 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.

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 k t ð Þ ¼ 0 for t < t wash-out and then switched on coupling at t = t wash-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 = t wash-out and resynchronization begun at that time. We performed 100 realizations of simulations for each washout 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 In the model, we represented the change in the advection pattern in a simplified way such that at earlier somite stages (t < t g ) 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 > t g ), 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, 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.

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 (t wash-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 -t wash-out ) decreased gradually over developmental time, Figure 4E . In contrast, the experimentally observed gradual decrease in time to FRS (FRS -t wash-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- 

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.

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 k 0 between neighboring oscillators. As discussed previously, coupling strength can strongly influence FRS, Table 1 and 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 t wash-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 t wash-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 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 (t wash-out = 9.5 hpf) until 9 ss (t wash-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.

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 t wash-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.

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 timedependent 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 timedelays 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 socalled 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.

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 .

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 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.

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 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 mM 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 mm 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.

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.

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 W l and W r , respectively with radius r. The tailbud is described as a half toroidal domain W 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 ¼ L x . We denote the position of the anterior end of the PSM at time t as x a ðtÞ, so x a ð0Þ ¼ 0. The length of the tissue is L ¼ L x À x a . X c denotes the position of the center of torus core curve. The position of posterior tip of the tailbud can then be written as L x ¼ X c þ R þ r. The y-axis points along the left-right axis and z-axis along the dorsal-ventral axis of an embryo.

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 

PSM length is constant if the velocity of the tailbud is the same as that of the anterior end of the PSM,

In this Lab reference frame, the position of a cell x L ð Þ c t ð Þ may change both due to tissue elongation and due to tissue deformations such as local tissue stretch.

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

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.

For the implementation of the tailbud frame for numerical simulations, we fix the position of the posterior tip of the tailbud at x ¼ L x . 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.

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 x i ðtÞ ¼ x i ðtÞ; y i ðtÞ; z i ð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) :

where v d x i ð Þ is cell advection from posterior to anterior, v 0 x i ð Þ represents the speed of intrinsic cell movement, n i t ð Þ is the cell polarity pointing the direction of intrinsic movement, F x i ; x j À Á is a physical contact force between cells i and j, and F b x i ð Þ is a boundary force that confines cells within the U-shaped domain. We explain each of these terms below.

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 v L ð Þ t . 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 v 

t , 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 x 1 t ð Þ and x 2 t ð Þ 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 Dx t ð Þ ¼ x 2 t ð Þ À x 1 t ð Þ is small. Due to the presence of an advective velocity field, the velocities of neighboring cells may be different, that is v

Þ. Thus, during a short time interval Dt these cells change their relative position due to this velocity difference,

where the velocity field gradient qv=qx ð Þ 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.

Back to the tailbud reference frame, the velocity field is

which we refer to as the cell advection pattern v d 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

x À x a ð Þ=L<x q À Á and the posterior PSM x À x a ð Þ=L>x q À Á , and consider different strain rates qv d =qx for each domain, Figure 2C . The advection field v d x ð Þ in Equation (1) depends on spatial position along the anterior-posterior axis x (x a x L) and is described as ( Figure 2C) :

and v p >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 ¼ L x À x a . The slope changes at the position x q 0<x q <1 À Á . Thus, x q 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 v p in the posterior PSM domain and v a À v p 1 À x q À Á À Á =x q in the anterior PSM domain.

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 v p in Equation (2) at time t ¼ t g . We assume that for t g , advection occurs mostly at the anterior part, v p <v a . For t g , we assume that advection occurs at the posterior part of the tissue v p >v a .

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 v 0 x ð Þ in Equation (1) depends on the spatial position x ¼ x À x a along the anterior-posterior axis ( Figure 2B) :

where v s is the maximum speed at the posterior tip of the tailbud, X v is the lengthscale of the mobility gradient, and the coefficient h determines the steepness of the mobility gradient, Figure 3 

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

where D f is the polarity noise intensity. White Gaussian noise fi t ð Þ and 'i t ð Þ satisfy fi t ð Þ ¼ 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 n i t ð Þ:

where

Þ= m x  n i t ð Þ j j and e z ¼ 0; 0; 1 ð Þ. These unit vectors m x and m y define a plane tangent to the unit sphere at position n i t ð Þ. Random displacement on this tangent plane will move the polarity vector away from the surface of the sphere, so its length ñ i t þ dt ð Þ will be larger than 1,

We normalize the polarity vector to correct for this radial displacement

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

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 n i performs an isotropic and uniform random walk on the surface of a unit sphere.

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

where >0 is the coefficient for intercellular force. (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, 

F b x i ð Þ in Equation

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

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

Then, we define frictionless boundary force in the columns as 

In the tailbud half-toroid

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

where (Figure 2-figure  supplement 1B) . We then define the boundary force in the half-toroidal domain as:

where b is the coefficient and r b is the length scale of the boundary force, dx ¼ R cos p i þ r cos p i cos q i À Dx i j j, dy ¼ R sin p i þ r sin p i cos q i À Dy i j j and dz ¼ r sin q i À Dz i j j.

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) :

where i is the phase of cell i, ! x i ð Þ is the autonomous frequency at position x i , k 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 x j À x i d c . The coupling strength k may be timedependent as described below. i t ð Þ is a white Gaussian noise

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 ! 

where x ¼ x À x a , s 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<x a is arrested (i.e. d i =dt ¼ 0 for x i <x a ). 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.

In this study, we allow the coupling strength in Equation (10) to be a time-dependent function k 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:

where t wash-out represents the time at which DAPT is washed out in simulations. In the presence of DAPT t < t wash-out , there is no interaction between cells. After DAPT washout at t = t wash-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 k s ¼ 0 in Equation (12) describes a constant coupling strength k 0 after DAPT washout. We first analyze resynchronization dynamics with the fixed value of the coupling strength k 0 by setting k s ¼ 0, Figures 3 and 4 

When we model PSM shortening, we consider that the position of the anterior end of the PSM x a ð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 u a ðu a >0Þ, Figure 4 -figure supplement 1A:

Hence, the length of the PSM becomes shorter as LðtÞ ¼ L x À u a t. 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 u a in simulations. For this, we impose the following condition: v a þ u a ¼ c;

where v a is the cell advection speed at the anterior end used in Equation (2) and c is a constant. Hence, v a can be expressed as v a ¼ c À u a . With Equation (14), the segment length S reads S ¼ c  T a where T a is the period of oscillation at the position x a . This anterior period of oscillation T a 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 T a matches the segmentation rate. When the PSM length does not change during simulation (u a ¼ 0), the anterior period T a is equal to the period at the tailbud, T a ¼ 2p=! 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 T a , which becomes shorter. To compensate for this effect here, we set T a ¼ 30 min by tuning ! 0 in Equation (10) if we include PSM shortening.

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 k 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 k t ð Þ ¼ 0 in Equation (10), noise by cell addition alone (D ¼ 0) ruins kinematic phase waves and stripe patterns, 

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 W l , we first compute phase order for cells within a thin slice domain

where n m is the number of cells within the domain W m ðxÞ, Figure 3 -figure supplement 1A. We set Dx equal to the cell diameter Dx ¼ d c so that Equation (16) This definition of Z m 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 Z m over a number M of these domains, to smooth out small finite size fluctuations:

The segment length in the anterior-posterior direction in our parameter sets is~5Dx. 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 W r . In this case, the domain W m ðxÞ in Equation (16) 

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 5Dx would result in low values of the order parameter in the presence of a phase gradient as in the anterior PSM.

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 x a to detect when a segment boundary position is set in simulations, Figure 3 -figure supplement 1A-C. We compute the mean phase É 1 t; x a ð Þ at position x a at time

2r ½ for the right PSM by using Equation (16) (m ¼ 1). As noted in Equation (16), we set Dx as the cell diameter Dx ¼ d c .

We then detect a time t i (i ¼ 1; 2; :::) that satisfies:

where # is a constant that we set # ¼ 3p=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, t i should be the time when the position of the segment boundary i is determined, Figure 3 -figure supplement 1B,C.

For resynchronization simulations starting from random initial conditions, we modify the above procedure as follows. After detecting the time t i when the mean phase of anterior cells becomes #, we check the local phase order parameter Z t; x a ð Þ 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 t i if:

where T a is the period of oscillation at the anterior end of the PSM x a as described in the section "PSM shortening". Since Z t i À T a =2; x a ð Þ is the average local phase order across nearly one segment length at t i À T a =2, it evaluates the integrity of the segment boundary and its neighboring interboundary regions. To suppress a false detection of a normal segment boundary caused by the fluctuation of Z t; x a ð Þ, we monitor Z t; x a ð Þ in a short time interval with a window size h in Equation (19). We set h ¼ 4 min in Equation (19). By visual inspection of stripe patterns in simulations, we set Z c ¼ 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 t i , we then specify the segment boundary number. Note that the subindex i of t 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 t iÀ1 ; x a ð Þ ¼ # at time t iÀ1 (t iÀ1 <t 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 t i and anterior period T a . We assign the current segment boundary with the expected segment number:

where ½ ½ represents rounding of real number to the closest integer value. Note that we assign the number ½ ½t i =T a when the phases of anterior cells become # slightly earlier, t i >T a ½ ½t i =T a À DT a considering the fluctuation after resynchronization, Figure 3 -figure supplement 1F. In this study, we set D ¼ 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.

FRS is defined as j f À 1 where j f is the smallest segment boundary number that was determined to be normal. The subtraction of 1 from j f 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 j p above which all the segment boundaries, including j p , are normal. Then, PLD is j p À 1. In the example simulation shown in Figure 3 -figure supplement 1D,E, FRS is 9 and PLD is 13.

Vortices are regions in space where the phase values circulate from 0 to 2p 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 2p 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 mm and z ¼ 2r À 20 mm (slices 1 and 2 and Figure 3-figure supplement 2A) , and two y-slices located at y ¼ y 0 mm and y ¼ y 0 þ 2r À 20 mm, (slices 3 and 4, Figure 3 -figure supplement 2B), where r is the radius of the PSM. For the left PSM, y 0 ¼ 0 mm, while for the right PSM y 0 ¼ 2R mm where R is the tailbud torus radius. The thickness of these four slices is 20 mm (~2 cell diameters, compare to the 50 mm of PSM diameter). We then project the phase values of cells in each slice to 2D planes, P ðaÞ (a ¼ 1; 2; 3; 4). The two x-y planes P ð1Þ and P ð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 3figure supplement 2A,C). For instance, the x-y plain P ð1Þ for the right PSM contains cells within the

The two x-z planes P ð3Þ and P ð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 x s ; y u ð Þ in the planes P ð1Þ and P ð2Þ where x s ¼ x a þ sDx (s ¼ 0; 1; 2; :::) and y u ¼ y 0 þ uDy (u ¼ 0; 1; 2; :::) with the grid size Dx and Dy. Similarly, we set grids x s ; z u ð Þ to the planes P ð3Þ and P ð4Þ where x s ¼ x a þ sDx and z u ¼ uDz. We chose Dx ¼ 5 m and Dy ¼ Dz ¼ 2 m. For each grid point in the plane P ðaÞ , we compute vorticity as follows. Below, we explain the case of P ð2Þ with the grid x s ; y u ð Þ, 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 x s ; y u ð Þ in the plane:

dl ffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffi

x À x s ð Þ 2 þ y À y u ð Þ 2 q <l þ dl

as shown in Figure 3 -figure supplement 2D, E. We set dl ¼ 5:5 m and l ¼ 14 m. The circular ring is subdivided into six domains V i (i ¼ 0; 1; 2; :::; 5) by angles p=3 measured counterclockwise from the x-axis. We then compute average phase i over cells within each subdomain V i , Z i e ii ¼ P k2Vi e ik =n i where n i is the number of cells in V i . If there is no cell within one of the subdomains V i (n i ¼ 0), we do not compute the vorticity for the grid point x s ; y u ð Þ and set ¼ 0.

To detect a vortex with clock-wise rotation, we permutate i based on their values (Figure 3-figure supplement 2D, E): 0 ¼ k where k = arg min{ i : i = 0,...,5}, 1 ¼ kþ1 , 2 ¼ kþ2 ,..., and 5 ¼ kþ5 , (mod 6). For a vortex with counter clock-wise rotation, we permutate i as: 0 ¼ k where k = arg min{ i : i = 0,..., 5}, 1 ¼ kÀ1 , 2 ¼ kÀ2 ,..., and 5 ¼ kÀ5 , where a negative value of k À j (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 x s ; y u ð Þ, i 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 a defined as:

A value of the correlation coefficient a 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 a ! a 0 , 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 x s ; y u ð Þ ¼ 5 À 0 =2p, Figure 3 -figure supplement 2F. If a<a 0 , we set x s ; y u ð Þ ¼ 0 to exclude false positive detection of a vortex by fluctuation of phase values. We used a 0 ¼ 0:75 throughout the article. After calculating vorticity for each grid point, we project x s ; y u ð Þ to x-axis. We use maximum projection, 2 ð Þ x s ð Þ ¼ yu max x s ; y u ð Þ, Figure 3 -figure supplement 2H.

By performing same calculations for the remaining three planes, we obtain 

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 F s of single defects defined below and compare it between simulations and experiment.

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

We defined normal and defective segment boundaries based on the local phase order at the anterior end of the PSM Z t; x a ð Þ 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 N t and the number of single defects N s appeared posterior to the segment 9 as in experimental data. Then, we computed the fraction, F s ¼ N s =N t , Figure 4 -figure supplement 8D.

We solved Equations (1) and (10) with the Euler-Maruyama method with the time step for integration dt ¼ 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).

draft, Writing -review and editing; Andrew C Oates, Conceptualization, Supervision, Funding acquisition, Visualization, Writing -original draft, Writing -review and editing; Luis G Morelli, Conceptualization, Formal analysis, Funding acquisition, Visualization, Methodology, Writing -original draft, Writing -review and editing . Transparent reporting form

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

The total number of defects in embryo i is N it ¼ N is þ N id . Then, we take the population average of these quantities. The average number of double defects in the population is

using left/right independence. Similarly, the average number of single defects N is h i is

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

and its population average

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

Using the left/right symmetry of defect distributions p l x ð Þ ¼ p r x ð Þ ¼ p x ð Þ, we can further simplify this expression

The good agreement observed between this probabilistic calculation and direct measurement of F s 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 

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 F s . 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.

Collective modes of coupled phase oscillators with delayed coupling

Short-lived her proteins drive robust synchronized oscillations in the zebrafish segmentation clock

A random cell motility gradient downstream of FGF controls elongation of an amniote embryo

Order and stochastic dynamics in Drosophila planar cell polarity

A clock and wavefront model for control of the number of repeated structures during animal morphogenesis

Boundary formation and maintenance in tissue development

Scaling behaviour at the onset of mutual entrainment in a population of interacting oscillators

Single-cell-resolution imaging of the impact of notch signaling and mitosis on segmentation clock dynamics

In vitro characterization of the human segmentation clock

Dictyostelium: the mathematician's Organism

Interplay of cell dynamics and epithelial tension during morphogenesis of the Drosophila pupal wing

Stochastic Methods: A Handbook for the Natural and Social Sciences

Control of segment number in vertebrate embryos

Multiple Delta genes and lateral inhibition in zebrafish primary neurogenesis

The elongation rate of RNA polymerase II in zebrafish and its significance in the somite segmentation clock

Intercellular coupling regulates the period of the segmentation clock

A multi-cell, multi-scale model of vertebrate segmentation and somite formation

Kinetic theory of coupled oscillators

Instability of Hes7 protein is crucial for the somite segmentation clock

Noise-resistant and synchronized oscillation of the segmentation clock

Excitable dynamics and Yap-Dependent mechanical cues drive the segmentation clock

Emergence of traveling waves in the zebrafish segmentation clock

Size-reduced embryos reveal a gradient scaling-based mechanism for zebrafish somite formation

Stochastic regulation of her1/7 gene expression is the source of noise in the zebrafish somite clock counteracted by notch signalling

Sustained oscillations and time delays in gene expression of protein Hes1

Time-delayed spatial patterns in a two-dimensional array of coupled oscillators

Notch signalling and the synchronization of the somite segmentation clock

Continuum theory of gene expression waves during vertebrate segmentation

Noise in the vertebrate segmentation clock is boosted by time delays but tamed by notch signaling

Stages of embryonic development of the zebrafish

Kosterlitz-Thouless physics: a review of key issues

Ordering, metastability and phase transitions in two-dimensional systems

Adjoint method provides phase response functions for delay-induced oscillations

Nonlinear phase-amplitude reduction of delay-induced oscillations

Chemical Oscillations, Waves, and Turbulence

Statistical macrodynamics of large dynamical systems. Case of a phase transition in oscillator communities

Regulated tissue fluidity steers zebrafish body elongation

Autoinhibition with transcriptional delay: a simple mechanism for the zebrafish somitogenesis oscillator

Faster embryonic segmentation through elevated Delta-Notch signalling

Segmentation of the zebrafish axial skeleton relies on notochord sheath cells and not on the segmentation clock

Priming, initiation and synchronization of the segmentation clock by deltaD and deltaC

Recapitulating the human segmentation clock with pluripotent stem cells

Foundations of Synergetics I

Control of waves, patterns and turbulence in chemical systems

A fluid-to-solid jamming transition underlies vertebrate body Axis elongation

Oscillatory expression of Hes1, p53, and NF-kappaB driven by transcriptional time delays

Delayed coupling theory of vertebrate segmentation

Patterning and mechanics of somite boundaries in zebrafish embryos

Patterning embryos with oscillations: structure, function and dynamics of the vertebrate segmentation clock

Waiting on the fringe: cell autonomy and signaling delays in segmentation clocks

Notch signalling synchronizes the zebrafish segmentation clock but is not needed to create somite boundaries

Vertebrate segmentation: from cyclic gene networks to scoliosis

Synchrony dynamics during initiation, failure, and rescue of the segmentation clock

Dynamics of zebrafish somitogenesis

Topology and dynamics of the zebrafish segmentation clock core circuit

Dynamics of the slowing segmentation clock reveal alternating two-segment periodicity

Oscillatory control of Delta-like1 in Somitogenesis and neurogenesis: a unified model for different oscillatory dynamics

Wave propagation in nonlocally coupled oscillators with noise

Spatial fold change of FGF signaling encodes positional information for segmental determination in zebrafish

Genetic oscillations. A doppler effect in embryonic pattern formation

Pulses of notch activation synchronise oscillating somite cells and entrain the zebrafish segmentation clock

Species-specific contribution of volumetric growth and tissue convergence to posterior body elongation in vertebrates

High-resolution in situ hybridization to whole-mount zebrafish embryos

The Her7 node modulates the network topology of the zebrafish segmentation clock via sequestration of the Hes6 hub

Self-Organization of embryonic genetic oscillators into spatiotemporal wave patterns

Random cell movement promotes synchronization of the segmentation clock

A framework for quantification and physical modeling of cell mixing applied to oscillator synchronization in vertebrate somitogenesis

Nü sslein-Volhard C. 1996. Mutations affecting somite formation and patterning in the zebrafish

Retinoic acid controls the bilateral symmetry of somite formation in the mouse embryo

Single-Cell reconstruction of emerging population activity in an entire developing circuit

Persistence, period and precision of autonomous cellular oscillators from the zebrafish segmentation clock

Timing by rhythms: daily clocks and developmental rulers

Three novel notch genes in zebrafish: implications for vertebrate notch gene evolution and function

Spiral waves of chemical activity

The Geometry of Biological Time

DeltaC and DeltaD interact as Notch ligands in the zebrafish segmentation clock

Cooperation of polarized cell intercalations drives convergence and extension of presomitic mesoderm during zebrafish gastrulation

Coupling delay controls synchronized oscillation in the segmentation clock

Cell cycle progression is required for zebrafish somite morphogenesis but not segmentation clock function

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 

Due to cell advection, cells reach the anterior end of the PSM x ¼ x a . These cells exit from the domains W l and W r at the advection speed jv d x a ð Þj ¼ v a . 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 ðx a ðtÞ<x<X c Þ, respectively, and that of the tailbud % t ðX c <xÞ. Then, we compute the difference between the measured density and the target density % 0 , d l ¼ % l À % 0 , d r ¼ % r À % 0 and d t ¼ % t À % 0 . 2. If all d l , d r , and d t are positive, we do not add a cell to the PSM and tailbud when a cell leaves the PSM from its anterior end x a ðtÞ due to advection. 3. In contrast, if some of d l , d r , and d t are negative, we add one cell to the region of which density is smallest. For example, if d 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 2p. 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 x x a þ z for the left and right PSM. In this study, we set z = 100 mm. 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 2p for ', and between 0 and p for f. The autonomous frequency, speed of intrinsic cell movement and advection speed for the added cell are determined depending on its added position.

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:where s r (s r >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 x i and the radius size rðtÞ. When updating rðtÞ into rðt þ DtÞ by Equation (15) above, some cells may be left outside the U-shaped domain, r i ðt þ DtÞ>rðt þ DtÞ where r i ðt þ DtÞ is the position of cell i in the radial direction defined in Equation (6) or (Equation 8) at time t þ Dt. For such cells, we correct their positions after the update of PSM radius as: r i ðt þ DtÞ ! rðt þ DtÞ À 2r b where r b is the lengthscale of boundary force used in Equations (7) and (9).

We set the initial position of the anterior end of the PSM x a ð0Þ as x a ð0Þ ¼ 0. We choose the cell number Nð0Þ to satisfy the cell density % ¼ % 0 with a given tissue geometry. Initial positions of cells x i ð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 ¼ x a ¼ 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Þ ¼ 3p=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 2p. The values of initial cell polarity angles in Equation (4) are also chosen randomly from the uniform distributions between 0 and p for f i 0 ð Þ, and between 0 and 2p for ' i 0 ð Þ. The autonomous frequency, speed of intrinsic cell movement and advection speed are determined by the initial position x i ð0Þ. The funders had no role in study design, data collection and interpretation, or the decision to submit the work for publication. 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.

We introduce a state variable that accounts for the presence or absence of a defective segment S ik x ð Þ, where i labels the embryo, k ¼ fl; rg 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 S ik x ð Þ takes the values zero and one depending on whether the segment is normal or defective.

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,andThese 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 p l x ð Þ ¼ p r x ð Þ ¼ p x ð Þ. In the following sections, we use probability theory to compute ALD, PLD and FRS from this spatial distribution of defective segments.

Let q a x a ð Þ be the probability to find an ALD at position x a . q a x a ð Þ can be expressed in terms of p x ð Þ as:The first factor Q xaÀ1 ¼0 1 À p ð Þ ð Þ in the second line represents the probability of normal segment boundaries from position x ¼ 0 to x ¼ x a À 1. The second factor p x a ð Þ is the probability of a defective segment boundary at position x ¼ x a . The resulting ALD distribution q a x a ð Þ 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,x a Á q a x a ð Þ:

Let q p x p À Á be the probability of PLD at position x p . It can be written as:The first factor p x p À Á in the first line is the probability of a defective boundary at x p . The second factor Q M ¼xpþ1 1 À p ð Þ ð Þ represents the probability that all the remaining segment boundaries posterior to x p are normal. The resulting PLD distribution q p x p À Á 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,

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 s d . We define FRS as the first normal segment after s d À 1. This is the definition we use to measure the embryonic FRS.With this definition of FRS, the probability q f x f À Á of the first normal segment boundary at locus x f is:The first factor of the second line represents the probability for all the segment boundaries between s d and x f À 1 to be defective. The second factor is the probability for a normal segment boundary to form at x f .To compute q f x f À Á , we set s d ¼ 9 as in the main text, Figure 1 -figure supplement 1B. The resulting distribution q f x f À Á 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 asIn Equation (30), we subtract one because FRS is defined with the anterior boundary of the first normal segment after s d À 1, see definition of FRS for experimental data in Materials and methods. q a , q p and q f 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.

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 distributionthe single defects distributionand the normal segments distributionThese 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) ,The two terms in this covariance can be written using the segment defect distributions of the two frameworks,Similarly, we obtain for single defects p s x ð Þ ¼ p l x ð Þ 1 À p r x ð Þ ð Þ þ 1 À p l x ð Þ ð Þ p r x ð Þ;and for normal segments p n x ð Þ ¼ 1 À p l x ð Þ ð Þ 1 À p r x ð Þ ð Þ:We can further simplify these expressions with the assumption p l x ð Þ ¼ p r x ð Þ ¼ p x ð Þ. The axial distribution of double defects p d x ð Þ gradually grows to a plateau at the onset of the defective region and then decays at its end, while the distribution of single defects p s x ð Þ 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 p d x ð Þ and p s 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.

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 iand