A Clock and Wavefront Self-Organizing model explains somitogenesis in vivo and in vitro

During mouse development, presomitic mesoderm cells synchronize oscillations of Wnt and Notch signaling forming sequential waves that pattern somites. Tail explants have shown that these waves can be formed in the absence of global embryonic signals in vitro but the self-organizing mechanism that underlies this process remains unknown. Here, we present a model called Sevilletor that synchronizes two oscillatory reactants via local cell communication to generate a variety of self-organizing patterns. This model provides a unifying minimal framework to study somitogenesis showing that classical models cannot recapitulate the waves formed in explants. Nevertheless, we show that explant patterning can be recapitulated by an extended Clock and Wavefront model where Fgf and Retinoic acid signals modulate a self-organizing oscillator implemented by Wnt and Notch. This model captures the change in relative phase of Wnt and Notch observed in mouse and explains how cells can be excited to oscillate by local cell communication. Introduction Temporal and spatial coordination of cell behaviors is essential for multicellular organization. During embryonic development, this coordination is mediated by cellular signaling pathways that synchronize the dynamic behavior of individual cells to form tissue-level gene expression patterns. A prime example of this dynamical coordination is the segmentation clock, a genetic network that patterns body segments in a rhythmic manner during the elongation of the main embryonic axis. This process, known as somitogenesis, is characterized by genetic oscillations that originate at the posterior tailend of the embryo and arrest sequentially at the anterior side forming somites, the precursors of body segments [3, 32]. A classic hypothesis to explain this sequential patterning is the Clock and Wavefront (CW) model [7], which proposes that the presomitic mesoderm (PSM) possess a smooth oscillator (clock) that freezes when cells exit a differentiation front (wavefront) at the tip of the tail. During embryo elongation, cells leave the front sequentially generating a periodic somite pattern of alternating frozen oscillatory states with a wavelength determined by the relative speed between embryo growth and oscillations. This theoretical model made no assumptions on the molecular nature of the clock but it proposed that the wavefront could be implemented by a posterior morphogen gradient that provides positional information [54]. Experiments in chick, Zebrafish and mouse have provided strong molecular support to the main hypothesis of the CW model. On the one hand, they provided evidence for the clock by showing that cells at the posterior side of the embryo exhibit oscillations of Notch, Wnt and Fgf signaling [35, 33, 47]. On the other hand, they showed that posterior signaling gradients of Wnt and Fgf allow cells to oscillate at the tip of the tail [11, 29, 3] functioning as the wavefront described in the 1 CW model, while Retinoic acid (RA) forms an anterior gradient that localizes within newly formed somites and promotes differentiation [12, 3, 32]. Together, this experimental evidence suggests that PSM cells possess an autonomous oscillator that is modulated by anterior-posterior morphogen gradients that act as global coordinating signals. Other experiments have also shown that the phase of oscillations in the PSM is shifted along the anterior-posterior axis generating kinematic waves of Notch signaling that traverse the embryo from posterior to anterior [23, 5]. This phase shift is lost when different parts of the PSM are dissociated in vitro [23], suggesting that the synchronization that underlies phase waves requires cell communication. In addition, it has been showed that Fgf and Notch oscillate in phase [10] but are out of phase with Wnt signaling [4]. This supported the idea that Fgf and Notch were part of the same oscillator, whereas Wnt oscillations were driven by a redundant genetic clock that might provide robustness [16, 10, 32]. In mouse, however, it has been shown that Wnt and Notch oscillations are coupled and that their relative phase changes at different positions along the embryo, passing from out of phase oscillations at the posterior tip to in phase oscillations at the middle part of the tail [4, 47]. One way to explain how the PSM can form phase waves is by assuming that it behaves as a series of oscillators with a frequency profile that stimulates a gradual slow down of the clock from posterior to anterior [27, 17]. The wave patterns formed by these models are tightly coupled to the specific shape of the frequency profile, whose slope determines the number of waves formed. For example, it has been proposed that in Zebrafish slightly shallower profiles can promote the formation of multiple phase waves in comparison with mouse [34]. Genetic manipulations in mouse, however, have shown that multiple phase waves can form even upon disruption of posterior gradients by ectopic activation of Wnt signaling [5], i.e., the phase of oscillations does not correlate with a monotonically decreasing gradient, but rather it appears to change autonomously in a periodic manner along the anterior-posterior axis. The ability of the PSM to generate phase waves autonomously has also been observed by monitoring tail explants in vitro [21, 49, 18]. These explants can generate periodic circular wave patterns (target patterns) despite local cellular rearrangements and in the absence of anterior-posterior gradients and axial elongation. It has been proposed that target patterns are formed by a system of coupled oscillators that tends to synchronize uniformly (a type IIIo system [9]), but that possesses an amplitude defect in the center of the explant originating from non-oscillating cells [30]. Although explants possess heterogeneities and differences in cell density that could influence the oscillatory behavior of the PSM, there is currently no evidence for a central non-oscillatory cell population [18]. Indeed, traveling waves can be observed even when the center of explants is ablated [18]. Moreover, classical type IIIo systems tend to synchronize uniformly in space and cannot explain the formation of periodic wavefronts observed in explants. Finally, wave formation has also been observed in more heterogeneous cultures obtained by mixing cells from different tailbud explants [18]. Taken together, these data suggest that, in addition to responding to global signals, PSM cells can coordinate the phase of oscillations in a self-organizing manner through local cell-to-cell communication. This idea has been explored in a previous study that proposed a progressive oscillatory reaction-diffusion (PORD) model to explain somitogenesis [8]. In contrast to the CW model, where new somites formed sequentially due to the progressive movement of the wavefront, in the PORD model new somites were formed in a self-organizing manner through a relay mechanism triggered by the last formed somite. In this way, somite patterning became independent of posterior signaling gradients and required instead a pre-patterned somite. In following studies, the PORD model was reformulated by letting a posterior gradient act as a wavefront to promote sequential somite patterning without requiring a pre-patterned somite [36, 20]. Nevertheless, these models could only generate patterns with stripes that had a one cell width, which did not recapitulate the larger phase waves and somite patterns observed in the embryo [5]. In summary, the mechanisms


Introduction
Temporal and spatial coordination of cell behaviors is essential for multicellular organization. During embryonic development, this coordination is mediated by cellular signaling pathways that synchronize the dynamic behavior of individual cells to form tissue-level gene expression patterns. A prime example of this dynamical coordination is the segmentation clock, a genetic network that patterns body segments in a rhythmic manner during the elongation of the main embryonic axis. This process, known as somitogenesis, is characterized by genetic oscillations that originate at the posterior tailend of the embryo and arrest sequentially at the anterior side forming somites, the precursors of body segments [3,32].
A classic hypothesis to explain this sequential patterning is the Clock and Wavefront (CW) model [7], which proposes that the presomitic mesoderm (PSM) possess a smooth oscillator (clock) that freezes when cells exit a differentiation front (wavefront) at the tip of the tail. During embryo elongation, cells leave the front sequentially generating a periodic somite pattern of alternating frozen oscillatory states with a wavelength determined by the relative speed between embryo growth and oscillations. This theoretical model made no assumptions on the molecular nature of the clock but it proposed that the wavefront could be implemented by a posterior morphogen gradient that provides positional information [54].
Experiments in chick, Zebrafish and mouse have provided strong molecular support to the main hypothesis of the CW model. On the one hand, they provided evidence for the clock by showing that cells at the posterior side of the embryo exhibit oscillations of Notch, Wnt and Fgf signaling [35,33,47]. On the other hand, they showed that posterior signaling gradients of Wnt and Fgf allow cells to oscillate at the tip of the tail [11,29,3] functioning as the wavefront described in the CW model, while Retinoic acid (RA) forms an anterior gradient that localizes within newly formed somites and promotes differentiation [12,3,32]. Together, this experimental evidence suggests that PSM cells possess an autonomous oscillator that is modulated by anterior-posterior morphogen gradients that act as global coordinating signals.
Other experiments have also shown that the phase of oscillations in the PSM is shifted along the anterior-posterior axis generating kinematic waves of Notch signaling that traverse the embryo from posterior to anterior [23,5]. This phase shift is lost when different parts of the PSM are dissociated in vitro [23], suggesting that the synchronization that underlies phase waves requires cell communication. In addition, it has been showed that Fgf and Notch oscillate in phase [10] but are out of phase with Wnt signaling [4]. This supported the idea that Fgf and Notch were part of the same oscillator, whereas Wnt oscillations were driven by a redundant genetic clock that might provide robustness [16,10,32]. In mouse, however, it has been shown that Wnt and Notch oscillations are coupled and that their relative phase changes at different positions along the embryo, passing from out of phase oscillations at the posterior tip to in phase oscillations at the middle part of the tail [4,47].
One way to explain how the PSM can form phase waves is by assuming that it behaves as a series of oscillators with a frequency profile that stimulates a gradual slow down of the clock from posterior to anterior [27,17]. The wave patterns formed by these models are tightly coupled to the specific shape of the frequency profile, whose slope determines the number of waves formed. For example, it has been proposed that in Zebrafish slightly shallower profiles can promote the formation of multiple phase waves in comparison with mouse [34]. Genetic manipulations in mouse, however, have shown that multiple phase waves can form even upon disruption of posterior gradients by ectopic activation of Wnt signaling [5], i.e., the phase of oscillations does not correlate with a monotonically decreasing gradient, but rather it appears to change autonomously in a periodic manner along the anterior-posterior axis.
The ability of the PSM to generate phase waves autonomously has also been observed by monitoring tail explants in vitro [21,49,18]. These explants can generate periodic circular wave patterns (target patterns) despite local cellular rearrangements and in the absence of anterior-posterior gradients and axial elongation. It has been proposed that target patterns are formed by a system of coupled oscillators that tends to synchronize uniformly (a type IIIo system [9]), but that possesses an amplitude defect in the center of the explant originating from non-oscillating cells [30]. Although explants possess heterogeneities and differences in cell density that could influence the oscillatory behavior of the PSM, there is currently no evidence for a central non-oscillatory cell population [18]. Indeed, traveling waves can be observed even when the center of explants is ablated [18]. Moreover, classical type IIIo systems tend to synchronize uniformly in space and cannot explain the formation of periodic wavefronts observed in explants. Finally, wave formation has also been observed in more heterogeneous cultures obtained by mixing cells from different tailbud explants [18].
Taken together, these data suggest that, in addition to responding to global signals, PSM cells can coordinate the phase of oscillations in a self-organizing manner through local cell-to-cell communication. This idea has been explored in a previous study that proposed a progressive oscillatory reaction-diffusion (PORD) model to explain somitogenesis [8]. In contrast to the CW model, where new somites formed sequentially due to the progressive movement of the wavefront, in the PORD model new somites were formed in a self-organizing manner through a relay mechanism triggered by the last formed somite. In this way, somite patterning became independent of posterior signaling gradients and required instead a pre-patterned somite. In following studies, the PORD model was reformulated by letting a posterior gradient act as a wavefront to promote sequential somite patterning without requiring a pre-patterned somite [36,20]. Nevertheless, these models could only generate patterns with stripes that had a one cell width, which did not recapitulate the larger phase waves and somite patterns observed in the embryo [5]. In summary, the mechanisms that underlies self-organization of periodic phase waves during somitogenesis remains unknown.
The default stability regime of cells during somitogenesis has been investigated by dissociating the PSM of mouse, chick and Zebrafish embryos to study how cells behave in isolation. Early experiments in mouse and chick showed that different PSM parts can undergo Notch signaling oscillations even when separated in vitro, however, their oscillations quickly lose spatial synchronization [22,23]. More recently, similar results have been observed in single cells dissociated from the Zebrafish PSM, which exhibited cell-autonomous oscillations with a noisy phase and amplitude [53]. Although these oscillations were variable across cell populations, it appeared that PSM cells could retain the identity possessed in the embryo. Indeed, while it was initially thought that isolated Zebrafish cells exhibited persistent oscillations [32], a recent study has shown that they stop oscillating after a number of cycles depending on the anterior-posterior position [40]. Similar single cell dissociation experiments in mouse have revealed that PSM cells can oscillate in aggregates but stop oscillating in isolation [18]. Interestingly, this study also showed that oscillations could be recovered by embedding isolated cells on a BSA substrate or by inhibiting actin polymerization with latrunculin A.
In summary, over the past four decades, the genetic characterization of somitogenesis has provided good support for the two main ideas proposed by the CW model: i) That PSM cells possess an autonomous oscillator and ii) That anterior-posterior gradients act as global coordinating signals to influence oscillations. In addition, the phase waves formed in tailbud explants [21,49,18] and upon perturbation of posterior signaling gradients [5] have demonstrated that PSM cells can form periodic waves independently of global signaling gradients in a self-organizing manner. Finally, recent experiments that dissociated the mouse PSM has shown that local cell communications and mechanical stimuli are necessary to excite individual cells to oscillate persistently [18]. Altogether, current evidence suggest that PSM cells undergo excitable self-organizing oscillations that can synchronize together to form periodic phase wave patterns. How Wnt and Notch signaling drive this dynamical behavior in vivo and in explants, however, remains unknown. Moreover, it is unclear how Fgf and RA gradients orchestrate this self-organizing process along the anterior posterior axis to promote changes in phase between Wnt and Notch oscillations and to drive sequential somite patterning.
Here we present a novel self-organizing oscillatory in silico model that is implemented by two self-enhancing reactants coupled by a negative feedback loop. We show that the model can form a wide range of self-organizing patterns by varying just one reaction parameter, which represents the relative strength between self-enhancements. This includes lateral inhibition patterns, rotating spirals, homogeneous oscillatory and static patterns, bi-stable patterns and a novel stripy spiral behavior with a striking visual resemblance to the pattern formed by the Belousov-Zhabotinsky chemical reaction [56]. Discovered in the early 1950s, this autocatalytic oxidative reaction can generate homogeneous sustained oscillations in well-stirred situations and self-organized spiral patterns with heterogeneous initial conditions. Previous theoretical work done in Brussels (Belgium) and Oregon (USA) has shown that this self-organizing dynamics can arise in models named the Brusselator [37,31] and the Oregonator [13]. Continuing the tradition, we named the model presented in this study Sevilletor, as it was developed in Seville (Spain).
By performing complex systems theory analysis on a pair of adjacent cells, we show that the different self-organizing behaviors of the Sevilletor model are based on synchronization and destabilization promoted by diffusion at different stability regimes. A recent study has shown that similar patterns can be formed by a discrete version of a bi-stable FitzHugh-Nagumo model [43]. However, this model generated different patterning dynamics due to changes in the range and strength of the spatial coupling between cells. In contrast, in the Sevilletor model, different self-organizing behaviors can be promoted by changing a single reaction parameter in a continuous manner. This makes our model amenable to study how morphogen gradients modulate self-organizing networks during embryonic development.
Remarkably, the Sevilletor framework can recapitulate the main theoretical models that have been proposed to study somitogenesis, such as the PORD model and the Clock and Wavefront (CW) model. We find that these models cannot recapitulate the patterning behavior observed in tail explants. We use the Sevilletor framework to devise a novel somitogenesis model where an intermediate region between posterior oscillating cells and anterior frozen cells generate self-organizing phase waves. We named this system the Clock and Wavefront Self-Organizing (CWS) model and show that it can recapitulate the formation of the target wave patterns observed in mouse explants. This is done by comparing simulated virtual explants with published experimental data by using a new pattern quantification approach. In agreement with experiments [18], the model predicts that cells in the intermediate regions are excited to oscillate by communicating with neighboring cells and stop oscillating when they are isolated. Finally, we propose that the two self-enhancing reactants of the CWS model correspond to Wnt and Notch signaling. Consistent with experimental data [47], this hypothesis predicts that Wnt and Notch signaling oscillate out of phase at tail tip and move anteriorly in the form of phase waves by oscillating in phase.
The Sevilletor model provides a powerful minimal framework for studying how oscillatory and excitatory behaviors can drive self-organizing spatial patterning. Here we used this framework to develop the CWS model, the first somitogenesis hypothesis that can explain the self-organizing potential of the PSM in vivo and in vitro.

Results
The first subsection of the results introduces a new reaction-diffusion model called Sevilletor, which exhibits a wide range of pattern formation capabilities. The rest of the first half of the result section is dedicated to a theoretical study of the properties of the model, and assumes that the reader is familiar with some aspects of complex systems theory, such as the notion of steady states, stability analysis, bifurcation diagrams and phase portraits [28,46]. Readers primarily interested in somitogenesis might benefit from reading the middle summary section that outlines the theoretical results for a broader audience. The second half of the result section focuses on the application of the Sevilletor model to study somitogenesis.

The Sevilletor model
We developed a minimal reaction-diffusion model to study how oscillatory signaling pathways can self-organize in space trough local cell-to-cell communication ( Figure 1A-B). The model consists of a system of two partial differential equations that represent interactions between two signaling pathways named u and v (equations (1) and (2)). To maintain the model as simple as possible, we consider the case where only the reactant u diffuses while v is immobile. Indeed, we find that a single diffusible reactant is sufficient to promote the emergence of a variety of spatial patterns. The model is implemented by a minimal set of regulatory terms: two positive linear terms associated with the self-enhancement of u and v, which promote instability, two opposite linear regulatory terms that couple u and v with an incoherent negative feedback, and cubic saturation terms for u and v that limit large deviations of concentrations from the steady state without having a significant effect on lower values.
These equations can be viewed as an extension of the first-order formulation of the Van Der Pol oscillator [14] (Section S3 in the supplementary material) with an additional linear self-enhancing feedback and cubic saturation term in equation (2). To characterize the behavior of the Sevilletor model, we focus our analysis on the effect of the two self-enhancement terms associated with parameters k 1 ≥ 0 and k 2 ≥ 0 for u and v respectively. Throughout the study, we set k 2 = 1 unless otherwise mentioned and let k 1 represent the relative strength of the self-enhancement of u and v: k 1 → k 1 k 2 . Finally, the last parameter in equation (1) is the diffusion constant D, which is set to D = 0.3 unless otherwise mentioned and describes how fast u spreads to promote the spatial coupling in the system. The dynamics of the network is centered around the fixed point (u * , v * ) (0,0) = (0, 0), which represents intermediate concentrations. However, the equations can be easily adjusted to form only positive values without affecting the behavior of the system. This can be done by shifting the steady state to (u * , v * ) (0,0) = (c, c) by substituting u → (u − c) and v → (v − c), where c is a positive constant ( Figure S5).

Self-enhancement strength and initial conditions determine Sevilletor patterning dynamics
The Sevilletor model can form a rich variety of two-dimensional wave patterns by varying just one parameter (k 1 ) (supplementary Movie 1). The bifurcation diagram in Figure 1C shows how k 1 affects the number of steady states and their stability (the corresponding bifurcation diagram for v is shown in Figure S2A). For all values of k 1 , there is an unstable steady state at (u * , v * ) (0,0) = (0, 0). However, for k 1 ≥ 2.3 the system undergoes a bifurcation that adds four additional steady states, two on the lower half-plane and two on the upper half-plane, for a total of five steady states. The two steady states furthest away from the center are stable, while the other three steady states are unstable. Importantly, each of the outermost stable states is in the proximity of an unstable saddle point. The higher the value of k 1 , the larger the distance between stable and unstable states. The reason becomes apparent by studying the phase spaces for different values of k 1 ( Figure 1F) where we see that the turning points of the u nullcline (shaped like an upside-down N) are pushed away from (u * , v * ) (0,0) when k 1 increases, generating intersection points between nullclines that are further apart. The calculations of the properties of each steady state is described in the STAR methods, along with further mathematical analyses of the points ( Figure S6).
Next, we investigate which patterning dynamics are promoted by these phase portraits with numerical simulations on a squared spatial domain with representative values of k 1 . We use an Euler method for time integration and a finite difference scheme to integrate space. Finally, we consider zero flux Neumann boundary conditions unless otherwise mentioned (STAR methods).
Starting from homogeneous initial conditions with identical values for (u(t = 0), v(t = 0)) on the whole squared domain ( Figure 1D left), numerical simulations reflect the simple predictions of the phase portraits shown in Figure 1F ( Figure 1E left column). For 0 ≤ k 1 < 2.3 they show homogeneous synchronized oscillations associated with the limit cycle, for k 1 ≥ 2.3 they show homogeneous static patterns associated with the bi-stable regime, where the whole squared domain reaches the closest stable steady state from the initial values (u(t = 0), v(t = 0)). However, starting from random initial conditions ( Figure 1D right), where each point in the squared domain has a random value of u and v chosen from a Gaussian distribution with zero mean and standard deviation 0.1, numerical simulations show a variety of complex oscillatory and static patterns ( Figure 1E right column). These include lateral inhibition patterns for k 1 = 0, which are characterized by cells with alternating opposite concentrations, rotating spirals when k 1 = 1, stripy spiral patterns that resemble the one formed by the Belousov-Zhabotinsky reaction [56] for k 1 = 2.3, propagating bi-stable fronts that generate homogeneous static patterns when k 1 = 3 and bi-stable frozen states for k 1 = 4 (supplementary Movie 1).
In summary, the phase portraits of the Sevilletor model correlate with the patterning dynamics observed in simulations started from homogeneous initial conditions, nevertheless they cannot explain the patterns formed with random initial conditions. For example, lateral inhibition patterns (k 1 = 0) and rotating spirals (k 1 = 1) are formed with similar phase portraits. This is because the spatial heterogeneity of random initial conditions introduces diffusion, which can have a significant impact on patterning dynamics. To overcome this limitation, we extend our complex system analysis to study how phase portraits change with diffusion between two cells.

The effect of cell-to-cell communication on patterning behaviors
We consider a 1D version of the Sevilletor model implemented in a system with only two cells. To study how diffusion affects the patterning behavior, we simulate this simplified 1D model starting from a heterogeneous initial state with (u 1 (t = 0), v 1 (t = 0)) = (0.1, 0) and (u 2 (t = 0), v 2 (t = 0)) = (−0.1, 0) for the two cells respectively and run two simulations for each representative value of k 1 : one without diffusion (D = 0 in equation (1)), and one with diffusion (D = 0.3 in equation (1)) ( Figure 2A-B). Without diffusion, each cell acts as an individual unit and its behavior depends solely on changes driven by reaction (black arrows in the phase spaces of Figure 2A). This case corresponds to the behavior seen in two-dimensional simulations with a homogeneous initial state ( Figure 1E left column), since when every cell has the exact same amounts of u and v there is no active contribution from diffusion, i.e. ∇ 2 u = 0. By including cell communication in the form of diffusion of u, the combined effect of reaction and diffusion coordinates the behavior of the two cells giving rise to large scale patterns (black and green arrows in the phase spaces in Figure 2B). To characterize these behaviors we also quantify phase changes over time (right columns in Figure 2A-B, details of the phase calculation are provided in the STAR methods. With k 1 = 0, in the absence of diffusion, the two cells undergo continuous oscillations with independent trajectories in the limit cycle (first row in Figure 2A). In the presence of diffusion, however, these oscillatory trajectories are counterbalanced by the effect of diffusion, which stabilizes the two cells in a frozen out-of-phase configuration (first row in Figure 2B) [45] (supplementary Movie 2). Our two-dimensional simulations show that this type of behavior generates a chessboard pattern where each cell has a frozen phase opposite to that of its neighbors, as in lateral inhibition patterns (right column of the first row in Figure 1E and first column in supplementary Movie 1).
When k 1 = 1, in the absence of diffusion, the two cells oscillate independently (second row in Figure 2A). In this case, however, when diffusion is added to the model, the cells continue to oscillate by following a smaller limit cycle that in the long run synchronizes the two cells together as in a type IIIo system [9] (second row in Figure 2B). Using random initial concentrations that lay on one of the half planes, this behavior is associated with homogeneous oscillations ( Figure S8). However, starting from random initial concentrations spread on the two half planes, the system generates rotating spirals similar to those formed by a diffusive Van der Pol oscillator (right column of the second row Figure 1E, Figure S7 and second column in supplementary Movie 1).
At the bifurcation point k 1 = 2.3, in the absence of diffusion, the two cells do not oscillate and are trapped at the nearest stable state on the upper or lower half plane (third row in Figure 2A). When diffusion is added to the model, however, the equilibrating effect of diffusion pushes cells out of stability towards the trajectory of the closest saddle point (green arrows in the third row of Figure  2B). Following this trajectory, each cell goes to the opposite half-plane towards the stable steady state, and it is again destabilized by diffusion. The repetition of this process generates a novel limit cycle that keeps cells oscillating (third row of Figure 2B and supplementary Movie 3). This new limit cycle generates in-phase oscillations of u and v, because the permanence time around the stable states is greater than the time it takes to follow the trajectory to the opposite half-plane (see trajectories in the phase space in supplementary Movie 3). In two-dimensional simulations with random initial conditions, this dynamic gives rise to a new type of spiral patterning behavior that to the best of our knowledge has not been described previously (right column of the third row Figure 1E and third column in supplementary Movie 1). The pattern looks very similar to those formed by classic models of the the Belousov-Zhabotinsky reaction [56] such as the Brusselator [37] and Oregonator [13], however it emerges from a different dynamical behavior.
Classic BZ models [37,13] generate spiral patterns based on a limit cycle around a single unstable steady state ( Figure S7), and the stripy spiral patterns formed by the Sevilletor model with k 1 = 2.3 arise from a diffusion-driven excitation of a bi-stable regime. This difference is reflected by different patterning dynamics. Unlike the spirals generated by the Brusselator model [37], which are formed at the beginning and disappear progressively over time, the Sevilletor model produces spirals that appear and disappear dynamically at random locations over time while moving in space. Numerical simulations show that these spirals are formed only with zero flux boundary conditions, while with periodic boundary conditions, the long term behavior of the system is to generate periodic phase waves that are straight ( Figure S1). In contrast, in classic BZ models the spiral formation is independent of boundary conditions. These results suggest that in the Sevilletor model, the appearance of spirals is due to crashing of phase waves promoted by reflective zero flux boundary conditions, as it has been observed in a previous study for a discrete version of the FitzHugh-Nagumo model [43] (Section S3 in the supplementary material).
The limit cycle that underlies stripy spiral patterns is possible only for values of k 1 in proximity of the bifurcation point ( Figure 1C). Indeed, larger values of k 1 increase the distance between the stable steady states and saddle points, thereby requiring a stronger effect of diffusion to destabilize the stable states and induce oscillations. For example, for a slightly larger self-enhancement strength with k 1 = 3, the effect of diffusion is still strong enough to stimulate a few oscillations but the two cells eventually freeze together at the same stable state (fourth row of Figure 2B). In two-dimensional simulations, this behavior generates a propagating front of one stable state that covers the whole domain, giving rise to a homogeneous static pattern (right column of the fourth row Figure 1E and fourth column of supplementary Movie 1). For even stronger self-enhancement with k 1 = 4, the distance between stable states and saddle points is too large for diffusion to impact the trajectory of the cells. Therefore, the two cells are frozen on the closest steady state (fifth row of Figure 2B). In two-dimensional simulations this behavior generates a static bi-stable frozen pattern that amplifies the noise present in random initial conditions (right column of the fifth row Figure 1E and fifth column of supplementary Movie 1).
A major benefit of the Sevilletor model over previous models is that its patterning dynamics can be changed smoothly by varying just the parameter k 1 . This property can be exploited to switch between different patterning behaviors over time, see for example the simulation in the supplementary Movie 4, where the model switches between a lateral inhibition pattern to a rotating spiral pattern, a stripy spiral pattern, a bi-stable frozen spiral pattern and finally back to lateral inhibition pattern k 1 = 0.0 → 1.0 → 2.3 → 4.0 → 0.0. In addition, it allows morphogen gradients to modulate the patterning behaviour of the model, as presented in the following section. Gradients can coordinate the behavior of the Sevilletor model to create complex global patterns During embryonic development, morphogen gradients can orchestrate tissue patterning in space by inhibiting or promoting specific parts of gene regulatory networks in a concentration-dependent manner. Our complex systems analysis revealed that the Sevilletor model can switch between different patterning behaviors upon progressive changes in the relative strength of the self-enhancements of u and v. Therefore, to characterize how the Sevilletor model behaves upon morphogen gradient modulation, we performed a two-dimensional simulation where the strength of k 1 increased linearly along the y-axis ( Figure 2C). The simulation recapitulated the results presented in Figure 1E and supplementary Movie 1 by forming lateral inhibition patterns, spirals, stripy spirals, homogeneous static patterns and bi-stable frozen patterns as k 1 increased. Moreover, it revealed that for intermediate k 1 values, the gradient promotes the formation of straight waves that move towards increasing values of k 1 (white arrows in Figure 2C and supplementary Movie 5).
To analyze more in-depth how self-enhancements promote different patterning behaviors, we performed a two-dimensional simulation where k 1 and k 2 increased linearly along the y-and x-axis respectively ( Figure 2D), and calculated the corresponding two-dimensional bifurcation diagrams of k 1 and k 2 ( Figure 2E). Finally, to take into account the effect of diffusion, we performed a series of two-cell simulations for each combination of parameters (k 1 , k 2 ) in a 100x100 grid with intervals of (k 1 , k 2 ) = ([0, 4], [0, 4]), for a total of 10000 two-cell simulations. In each case, we quantified the number of oscillations exhibited by the two cells and colored the corresponding point in the bifurcation diagram accordingly, see Figure 2E and the STAR methods for details on the calculation. This allowed us to further refine the bifurcation diagram into regions associated with different patterning behaviors depending on the number of oscillations ( Figure 2F).
This extended bifurcation diagram correlated well with the graded two-dimensional simulations presented in Figure 2D and reveled that the effect of diffusion on patterning depends on the relative self-enhancement strengths of u and v in each stability region. Inside the region with one steady state, diffusion can stop oscillations to form lateral inhibition patterns only when k 1 k 2 < 1, see the cyan and pink regions in Figure 2F. On the other hand, inside the five steady state region with k1 > k 2 , diffusion can promote stripy spiral patterns and homogeneous static patterns when k 1 k 2 has low and intermediate values respectively, see the green and purple regions in Figure 2F.

Summary of the theoretical analysis
We developed a minimal model called Sevilletor that can generate self-organizing spatial patterns by coupling oscillations with diffusion. To characterize the capabilities of the model, we performed numerical simulations by varying the relative self-enhancement strength of u and v, which is a key parameter (k 1 ) that promotes different patterning behaviors.
Our analysis reveled that, in the absence of diffusion and/or when every cell has the same initial values of u and v, the model can be either in an unstable state where the concentrations of u and v oscillate, or in a bi-stable state where the concentrations of u and v are frozen on high and low values (left column in Figure 1E and Figure 2A). Starting from random initial conditions, in the presence of diffusion, the model gives rise to five different behaviors depending on the relative self-enhancement strength of u and v: lateral inhibition patterns, spiral patterns, stripy spirals, homogeneous patterns and bi-stable frozen patterns ( Figure 1E on the right and supplementary Movie 1). An analysis with two-cell simulations revealed that diffusion plays a different role in each of these scenarios, as it can be seen by comparing the trajectories with and without diffusion in phase space (Figure 2A-B and supplementary Movies 2 and 3). It can be either stabilizing (lateral inhibition), synchronizing (spirals) or destabilizing (stripy spirals, homogeneous patterns). Finally, two-dimensional simulations revealed that spatial gradients of the self-enhancement strengths of u and v can promote a smooth transition between different patterning behaviors in space and can simultaneously orient phase waves towards regions with stronger self-enhancement ( Figure 2C and supplementary Movie 5). In the following section, we show that this feature makes the model amenable to simulate the different hypotheses that have been proposed to explain somitogenesis.

Sevilletor model as a framework to study somitogenesis
The variety of spatial patterning behaviors exhibited by the Sevilletor model and its ability to change behaviors upon spatial modulations makes it a suitable framework for examining different somitogenesis hypotheses. In this context, the two oscillatory reactants of the Sevilletor model can be interpreted as Wnt and Notch signaling. Indeed, these two signaling pathways play an important role in somitogenesis and exhibit coupled oscillations as the two reactants of the Sevilletor model [47]. Wnt is a long range paracrine signaling, which can be represented by the diffusible reactant u, while Notch is a juxtacrine signaling [32], which can be interpreted as the non-diffusive reactant v ( Figure 3B). We suggest that the coupling between Wnt and Notch oscillations supports the idea that these two signaling pathways interact to generate oscillation. According to this hypothesis, the two equations of the Sevilletor model can be rewritten as equations (3) and (4), where the cross-regulatory terms represent a coupling by a feedback between Wnt and Notch [47]. The parameter k 1 in the model is associated with the strength of the positive self-enhancement of Wnt, while the parameter D corresponds to the diffusion coefficient of Wnt. Negative self-regulatory cubic terms in the equations should be interpreted as phenomenological saturation around the central steady state of the system, rather than specific molecular reactions or gene regulations. Indeed, the goal of this model is not to build a precise molecular-level description of the gene regulatory network that underlies somitogenesis, but rather to provide a general framework to represent somitogenesis hypotheses as dynamical systems. These descriptions can be used to investigate how well alternative hypotheses fit the qualitative patterns observed in different experiments.
We analyze different somitogenesis hypotheses by simulating equations (3) and (4) in a 2D growing rectangular grid of virtual cells that elongates with constant speed along the y-axis, which represents the anterior-posterior axis of the developing tail (see the STAR methods for details). In this model, tail growth is implemented by proliferation of the posterior-most cells, which generate a new line of cells that inherit the concentrations of Wnt and Notch (details in the STAR methods). Other important ingredients in the simulation are two opposing morphogens of Fgf and RA signaling that have high values at the posterior and anterior side respectively ( Figure 3A). The two morphogens are introduced in the model with the terms Fgf and RA in equation (3), which correspond to sigmoid functions with formula ∆k 1+e −y , where y represents an anterior-posterior spatial coordinate and ∆k the amplitude of the gradient that promotes changes in k 1 . Over time, the morphogens move in coordination with tail elongation: Fgf is maintained at the tip of the tail and RA follows tail growth on the anterior side, giving rise to regions with different values of RA and Fgf. Similar to the simulations presented in Figures 1 and 2, we use an Euler method for time integration and a finite difference scheme to integrate space, and we consider zero flux Neumann boundary conditions (STAR methods). In the following sections, we show that this theoretical framework can recapitulate the PORD and CW somitogenesis hypotheses simply by using different anterior-posterior modulations of k 1 . Finally this theoretical framework suggests a novel somitogenesis model called the Clock and Wavefront Self-Organizing model (CSW).  Figure S9). Right: average Wnt (lime) and Notch (cyan) values over time at the tip, middle and top part of the tail. Check-marks indicate whether the measurements match quantifications presented in [4,47].
C) The PORD model is implemented with k 1 = 0 and relies on lateral inhibition to generate a periodic pattern of stripes with a one cell width via a rely mechanism triggered at the top by a stripe of high Notch signaling in the initial conditions. This happens in the absence of anterior-posterior signals. A Sevilletor implementation of the PORD model The first somitogenesis hypothesis that we examine is the Progressive Oscillatory Reaction-Diffusion (PORD) model originally presented in [8]. This model is characterized by the sequential appearance of somites independent of global gradients. It also predicts homogeneous oscillations at the posterior side of the tail. The model consists of two reaction-diffusion equations associated with an immobile reactant called "A" and a diffusive reactant called "R" (Section S7 in the supplementary material). Although in its original implementation the model produces stripes separated by several cells, a distinctive feature of all PORD models is that they form stripes with a one cell width. To investigate the reason behind this patterning behavior, we analyzed the original PORD model introduced in [8] by performing two-cell simulations similar to ones presented in Figure 2B, and found that the stripes are formed by freezing oscillations of neighboring cells ( Figure S11) with the same diffusion-driven mechanisms of lateral inhibition patterns formed in the Sevilletor ( Figure 1E and supplementary Movie 2). To form stripes rather than chessboard patterns, the original PORD model relies on specific initial conditions that consist of an anterior line of cells with high R concentration. This initial state promotes a sequential freezing of neighboring cells in a stripe-like conformation.
To implement the PORD model in the Sevilletor framework, we consider that the diffusible reactant "R" corresponds to Wnt and the immobile reactant "A" to Notch ( Figure 3C). We use equations (3) and (4) with parameters k 1 = 0 and D = 1 for all the cells, and set the gradient terms Fgf and RA to zero. Both the original PORD implementation [8] and our model have a unique unstable steady state that promotes oscillations. In contrast to the original PORD model, however, in our implementation the orbit around this central steady state is symmetric ( Figure S11). In the original PORD model and in our implementation, the two reactants oscillate out of phase at the tip of the tail but, in contrast to the original PORD model, in our implementation the oscillations arrest out of phase on the anterior side. In phase freezing can be promoted by switching the signs of the cross-regulatory linear terms that couple Wnt and Notch, however, this changes also the phase of oscillations at the tip ( Figure S2).
In agreement with the original PORD model, the model generates periodic somite patterns trough a relay mechanism that is triggered by initial conditions with a line of high Notch concentration at the anterior boundary. This initial condition promotes a sequential freezing of neighboring cells based on lateral inhibition that promotes horizontal stripes with a one cell width ( Figure 3C and supplementary Movie 6). In the original model, these stripes are separated by several cells while in our model by a single line of cells. This difference depends on the simpler equations of our model that give rise to a symmetric phase space and lack of a heavy side function to cut low reactant values.
Importantly, as discussed in [36], our simulations reveal that this pattern is fragile to spatial noise. Indeed, noise added during the simulation breaks up the stripe configuration promoting a chessboard-like pattern instead ( Figure S13). Robustness improves when the PORD model is combined with global modulations by posterior gradients [36]. Nevertheless, the overall fragility of the model against noise and its tendency to form somites with a one-cell width fails to recapitulate the two dimensional somite patterns observed in vivo. Finally, even under the influence of a posterior gradient, the PORD model cannot form a central region where Wnt and Notch oscillate in phase as observed in the mouse tail [4,47].

A Sevilletor implementation of the Clock and Wavefront model
The second somitogenesis hypothesis that we consider is the Clock and Wavefront (CW) model [7]. In its simplest formulation, this model is characterized by homogeneous oscillations at the posterior side of the tail (clock) and by the sequential arrest of oscillations on the anterior due to tail outgrow (exit of the wavefront). Mathematical descriptions of this hypothesis as a system of coupled oscillators [27,17] have captured the main qualitative aspects of somitogenesis. However, they did not provide a mechanistic explanation for the emergence and the arrest of oscillations. Recent studies have described this patterning event by using one-dimensional dynamical systems, where the freezing of oscillations was modeled as a bifurcation from an unstable state to a bi-stable regime along the anterior posterior axis [15,19]. Our two-dimensional implementation of the Clock and Wavefront model is based on the same principle and has a phase space that resemble the ones in [19].
To implement the CW hypothesis in the Sevilletor framework, we developed a two-dimensional dynamical model with a posterior gradient implemented by the sigmoid function Fgf in equation (3) and set the term RA = 0 ( Figure 3D and supplementary Movie 7). In this model, Fgf acts as a wavefront by reducing the parameter k 1 associated with Wnt self-enhancement, to promote a transition from an oscillatory state to a bi-stable regime. Indeed, the theoretical analysis presented in Figure 2B reveals that with low values of k 1 the Sevilletor model has a single unstable steady state, while for higher values of k 1 the system is in a bi-stable regime. In the model, cells have a default strong Wnt self-enhancement with k 1 = 4.0, which is associated with bi-stability and is decreased by Fgf to k 1 = 1 promoting oscillations at the posterior side (see network diagram in Figure 3D). Starting from random initial conditions that lie on the lower half plane, this model gives rise to homogeneous oscillations at the posterior side and does not form spirals ( Figure S8). From a biological point of view, this initial conditions correspond to low Wnt and Notch signaling, which is a plausible initial state of the PSM that can explain why spiral waves have never been observed in vivo. As shown in Figure S2, we find that the same results can be obtained if Fgf modulates the self-amplification of Notch instead of Wnt, when the interactions connecting Wnt and Notch are inverted.
In contrast to the PORD model, our dynamical implementation of the CW model generates robust somite patterns that span across multiple cells and do not depend on pre-patterns in the initial conditions. Similarly to the PORD case, the CW model give rise to out-of-phase oscillations of Wnt and Notch at the posterior side of the tail but freezes Wnt and Notch in phase at the anterior side, as observed experimentally [4,47] ( Figure 3D and supplementary Movie 7). Nevertheless, the model cannot recapitulate the formation of a central region in the tail where Wnt and Notch oscillate in phase, and it is not able to form multiple phase waves.
The Clock and Wavefront Self-Organizing model  Figure 1E on the right) giving the name to the model: the Clock and Wavefront Self-Organizing model (CWS). This model generates robust periodic somite patterns that span multiple cells. The initial state of the CWS model is the same as the one used for the CW model. As shown in Figure  S2, we find that the same results can be obtained by letting the gradients affect the self-amplification of Notch instead of Wnt, when the signs of the interactions connecting Wnt and Notch are switched. Nevertheless, in this case we observe that the precision of the somite boundary decreases slightly.
The pattern generated by the CWS model appear similar to the one formed by the CW model. However, a more careful examination of the patterning dynamics of the two models reveals that, in the middle tail region, the CWS model forms phase waves that travel anteriorly, as it can be seen by comparing the simulations in Figures 3D and 3E and supplementary Movies 7 and 8. Quantification of the oscillation phases along the anterior-posterior axis reveals that in the middle region Wnt and Notch oscillate in phase, in agreement with experimental data [4,47] (graphs in Figure 3E). The switch from out of phase to in phase oscillations of Wnt and Notch is driven by the excitable dynamics of the model. Indeed, as cells move from the posterior to the middle section of the tail, they enter a bi-stable region where the destabilization driven by diffusion ( Figure 2B and supplementary Movie 3) generates a limit cycle with a long permanence time around the two steady states where Wnt and Notch are in phase. Then, as the cells transit to the anterior side, the oscillations will remain frozen at the two stable states forming stripes where Wnt and Notch signaling are in phase.
The self-organizing capabilities of the CWS model guarantee robustness against noise during the patterning processes ( Figure S13). In addition, the model has the ability to regenerate patterns upon drastic perturbations. For example, we find that the model can regenerate somite patterns upon tissue ablations in the middle tail region as well as on the anterior part where oscillations are frozen ( Figure S10). Further confirming the validity of our extension of the CW hypothesis, we observed that the model can recapitulate the formation of wider somites patterns observed upon Fgf inhibition and faster growth of the tail [42] ( Figure S10).
In summary, the CWS model can recapitulate the formation of realistic somite patterns by combining modulations from global opposite gradients of Fgf and RA with an intrinsic self-organization potential based on local cell communications. The different self-organizing regimes of the model along the anterior-posterior axis can explain for the first time the change of relative phase between Wnt and Notch observed in the tail middle region during mouse somitogenesis [4,47]. We next wish to explore to which extent the different somitogenesis models can recapitulate the formation of self-organizing somitogenesis waves observed in embryonic explants.

Self-Organizing somitogenesis waves in ex-vivo tail explants
In recent years, tail explants have been established as an in vitro assay to study how PSM cells can synchronize oscillations when decoupled from global embryonic signals [18,21,49]. A typical explant is made by dissecting and culturing a portion of the embryonic tail during somitogenesis with the posterior end of the tail pointing upwards. Over time, these explants spread in the culture, forming a thin circular quasi-monolayer [2] ( Figure 4A). Other types of explants are made by mixing together cells from different individual tails and by allowing the mixture to form a circular colony. Depending on the anterior-posterior position of the dissected region, explants generate different types of Notch signaling patterns including multiple phase waves emitted from the center of the explant that are called "target patterns" [18,21]. The fact that explants form coherent Notch signaling patterns despite the heterogeneities induced by mixing cells from different tails suggests that PSM cells have the capacity to generate phase waves in a self-organizing manner [49,18].

Virtual explants to test different somitogenesis hypotheses
To test the self-organizing capabilities of different somitogenesis hypotheses, we perform virtual explant simulations of the PORD, CW and CWS models. A virtual explant simulation is constructed to reflect experiments ( Figure 4A). This is done by first performing a y-projection of a simulated tail to calculate the average intensity of Wnt and Notch along the anterior-posterior axis. Secondly, by projecting a portion of the anterior-posterior axis radially to obtain a circular grid of cells. The posterior part of the tail is in the center of the circle and the anterior part is on the outer edge. These circular simulations begin with the concentrations possessed in the tail, which means that each cell retains a memory of the state and the gradients experienced in vivo [21]. Finally, the positions of each cell in the virtual explant is switched randomly with its neighbors to take into account cell mixing ( Figure 4A). Using this method, we generate virtual explants of the posterior tip, the middle, and the whole posterior part (tip+middle) of simulated tails.

Quantification of explant patterns
In order to compare the virtual and experimental explants, we devise a novel quantification method to characterize two-dimensional wave patterns. Space-time plots of circular profiles defined around the center of the patterns are used to distinguish between homogeneous oscillations, target patterns, spirals and bi-stable patterns ( Figure STAR2 in the STAR methods). The space-time plot of bi-stable patterns is characterized by vertical straight lines, which represent concentration profiles that do not change over time. In contrast, oscillatory patterns are characterized by alternating horizontal lines of high and low values that reflect the oscillation. Finally, spiral patterns are characterized by space-time plots with diagonal stripes that reflect the movement of spiral arms along the circle. A single circular space-time plot, however, cannot distinguish between homogeneous oscillations and target patterns. Instead we use two circular profiles with different radius. In the case of homogeneous oscillations, the two space-time plots have horizontal lines that are in the same phase, while target patterns have horizontal lines in opposite phase ( Figure STAR2A-B in the STAR methods).
We compare the behavior of PSM cells in in vitro experiments obtained from Matsuda et al. [24], Lauschke et al. [21], Hubaud et al. [18] and Tsiairis and Aulehla [49] with simulated explants of the tip, middle and the whole posterior part of the tail in the PORD, CW and CWS model. The simulations are performed by generating multiple virtual explants in each situation starting from different development stages, which are randomly distributed between the three and four somite stage. This provides an overview of the different behaviors that can be obtained by varying initial conditions ( Figure S3).

Explants of the posterior tail tip
The first type of explants that we generate are obtained from the posterior tip of the tail in the PORD, CW and CWS models (first rows in Figure 4B-D). In all models, this portion of the tail contains cells with a unique unstable state that oscillate synchronously. Explants predict that in all models cells will continue global synchronized oscillations despite cell mixing (first row in Figure 4B-D). Unfortunately, this prediction cannot be easily tested experimentally because it is difficult to completely isolate the tip of the tail from the anterior part. As an alternative, cell culture experiments can be used to differentiate a group of cells uniformly into tip-like PSM cells that exhibit synchronized oscillations of Notch [24] in agreement with all three models.
Explants of the middle and the whole posterior part of the tail The second type of explants that we consider includes explants obtained from the middle portion, and the whole posterior part of the tail. This type of explant experiments have been performed experimentally by Hubaud et al. in [18] and Lauschke et al. in [21] by first dissecting a posterior portion of the tail that included the tip and the middle part the tail [18] and by further cutting away the tip to retain only the middle part [21]. In both cases, the explants gave rise to concentric waves of Notch signaling that propagated from the central core towards the periphery.
The behavior of virtual explants that contain only the middle portion of the tail varies depending on the somitogenesis model that is considered. In the PORD model, both middle and whole posterior parts of the tail explants generate chessboard patterns, due to the intrinsic propensity of cells to generate lateral inhibition patterns starting from heterogeneous initial conditions (second and third row in Figure 4B). In the CW model, the middle and the whole posterior part of the tail explants give rise to homogeneous oscillations in the inner core of the explant and a frozen bi-stable pattern on the outer part (second and third row in Figure 4C). This due to the fact that cells in the CW model can undergo only two types of behaviors: homogeneous oscillations and bi-stable frozen behavior. Therefore, neither PORD nor CW explants can recapitulate the formation of target patterns observed in the experiments.
The virtual explant of the middle region in the CWS model can produce self-organizing phase waves originating from spiral centers (second row in Figure 4D). These centers arise due to the mixing of cells, which creates regions where neighboring cells have phases that lie on different half-planes. Depending on the initial conditions, these phase waves can dissolve into a homogeneous pattern or be maintained in the presence of spiral centers, see the possible outcomes in Figure S3. In the latter case, waves originating from different spiral centers can merge together to give rise to circular phase waves that resemble target patterns (second row in Figure 4D). As expected, quantification with two radial profiles of different radius shows space-time plots with out of phase horizontal lines that are characteristic of target patterns. The lines contain small dipping points due to the irregular shape of the circular pattern created by spirals. The same behavior can be observed in the experiment with the mid-tail section by Lauschke et al. [21] where a target pattern of signaling waves propagates out from the center. A virtual CWS explant of the middle tail part is shown in the supplementary Movie 9.
Explants of the CWS model obtained from the whole posterior tail give rise to regular target patterns consistently in all simulations ( Figure S3). This is due to the fact that cells derived from the tip, which are positioned at the center of the explant, undergo homogeneous oscillations, while cells on the outer part of the explant, which corresponds to the middle part of the tail, have the ability to generate self-organizing phase waves (third row in Figure 4D). The homogeneous oscillations at the center act as a guide to promote radially symmetric waves propagating towards the periphery. Indeed, quantification with circular profiles shows straight horizontal lines with opposite phases in the two space-time plots. The same behavior can be observed in explants of the posterior part of the tail that include the middle and the tip in the experiment by Hubaud et al. in [18]. A virtual CWS explant of whole posterior tail is shown in the supplementary Movie 10.

Explants obtained by mixing cells
A way to further test the ability of the PSM to self-organize is by randomly mixing cells to disrupt any pre-pattern in the tail. This has been done in two studies [18,49] that generated explants by centrifuging together cells taken from tails of different individual embryos. In the experiment by Hubaud et al. [18], these mixed explants gave rise to continuous rotatory oscillatory movements that synchronized in space. Virtual explants obtained by mixing cells from the tail of the PORD and CW model cannot recapitulate this self-organizing behavior ( Figure S4B-C). In contrast, explants generated by mixing cells from the whole posterior part of the tail of the CWS model induce spiral rotatory patterns that are characterized by diagonal lines in the space-time plots (first row in Figure  S4D).
Experiments presented in [49], however, showed that explants derived by mixing cells from different tails can also give rise to static Notch signaling patterns after undergoing few global oscillations. Virtual explants derived by mixing cells of the CWS model suggest that this type of behavior arises when explants are obtained by mixing large portions of the tail that contain anterior cells in the frozen bi-stable regime. The presence of these cells promotes patterns that remain frozen after a small initial rearrangement (second row in Figure S4D).  Figure S12). Right: two space-time plots showing Notch signaling over time along the two circular profiles in the last time point of virtual explant simulations (cyan and green circles). The graph on the right shows the average projection of the two space-time plots along the x-axis with correspondent confidence interval (details on radial profiles in the STAR methods). Check-marks indicate whether the simulations fit experimental data. B) Virtual explants of the PORD model. Tip explants make homogeneous oscillations marked by in phase oscillations in the small and large circular profiles, in agreement with cell culture experiment with tip-like cells showed in [24] (green mark). Middle explants and whole posterior tail explants, however, create a lateral inhibition pattern that propagates in from the outer edge and do not resemble the Notch signaling patterns observed in experiments showed in [21] and [18] (red cross). C) Virtual explants of CW model. Tip explants generate homogeneous oscillations in agreement with the cell culture experiment in [24] (green mark). Middle explants and whole posterior tail explants, however, do not reproduce the circular wave patterns observed in experiments in [21] and [18] (red cross) generating instead homogeneous oscillations in the tip cells and frozen patterns in cells coming from the middle. D) Virtual explants of CWS model. Tip explants generate homogeneous oscillations that resemble cell culture experiments in [24] (green mark) with in phase oscillations in small and large circular profiles. In addition, middle and whole posterior tail explants generate circular wave patterns that resemble the experiments showed in [21] and [18] (green marks). The presence of circular wave patterns is confirmed by out of phase oscillations between small and large circular profiles on the space-time plots and graph on the right.

Discussion
During embryonic development, cells need to coordinate their behaviors to form coherent spatial patterns that drive tissue specification. One way to achieve this coordination is by responding to global signals, such as morphogen gradients that provide positional information to the cells [54]. Alternatively, self-organizing spatial patterns can be formed by coupling cell-autonomous behaviors at the tissue level through local cell communication [50]. In addition, increasing evidence is showing that these two patterning strategies are not exclusive and that embryonic development is often controlled by self-organizing processes guided by external signals [26,39,1].
The waves of oscillating gene expression observed during somitogenesis are an example of this coordination. Models based on the Clock and Wavefront hypothesis [7] suggested that these waves arise when posterior morphogen gradients modulate the frequency of oscillations along the anterior-posterior axis [27,17]. However, multiple waves can be formed despite ectopic anterior activation of posterior morphogens [5]. In addition, explant experiments showed that presomitic mesoderm cells can form waves that resemble somitogenesis in vitro [21,49,18] despite cellular re-arrangements and in the absence of global embryonic signals. This evidence suggests that, during somitogenesis, presomitic mesoderm cells can coordinate genetic oscillations in a self-organizing manner via local cell communication.
To investigate this hypothesis, we developed a minimal self-organizing model called Sevilletor that couples the oscillations of two self-enhancing reactants through local cell communication. Our analysis showed that, starting from random initial conditions, this model can form a variety of spatial patterns depending on the relative strength of self-enhancements. Without self-enhancement of the diffusible reactant, the model formed chessboard lateral inhibition patterns. With a weak self-enhancement it formed spiral patterns and with a strong self-enhancement it formed bi-stable salt and pepper patterns. In addition, we found that with an intermediate self-enhancement the model generated stripy spiral patterns that resemble the Belousov-Zhabotinsky reaction [56,38] that were formed due to a diffusion-driven destabilization of a bi-stable regime. We used this model as a backbone to investigate if the patterning dynamics observed during somitogenesis could be recapitulated by self-organizing behaviors modulated by Fgf and RA signaling. To this end, we performed simulations of the developing tail as well as of virtual tail explants.
Although the main goal of this study was not to derive a detailed molecular model of somitogenesis, we proposed that the diffusible reactant of the model is related to Wnt, which is a long range paracrine signaling pathway, while the immobile reactant is related to Notch, which is a juxtacrine signaling. Indeed, these two signaling pathways play a central role in somitogenesis and oscillate in coordination at the posterior side of the embryo [47]. Previous somitogenesis models have coupled these two pathways in larger biochemical networks that produced oscillations at the single cell level [51]. Here, for the first time, we considered an oscillator that couples both signaling pathways to recapitulate the two-dimensional patterning dynamics of somitogenesis in vivo and in vitro. The crucial part of the model is a negative feedback between Wnt and Notch, which is a known regulatory topology that promotes coupled oscillations [6]. Throughout the study, we considered a negative feedback implemented by a positive regulation from Wnt to Notch and a negative regulation from Notch to Wnt. However, equivalent results can be obtained by switching the signs of these two interactions ( Figure S2). In addition, we assumed that both signaling pathways have self-enhancing feedback that helps to promote an oscillatory instability. To keep the complexity of the model to a minimum, we considered equations with simple linear terms and added negative cubic saturation terms for both reactants. These non-linear terms do not represent specific molecular interactions and are introduced only to impose symmetric saturation around the central steady state of the system. Finally, to maintain the model as simple as possible, we considered Wnt as an oscillatory signal and ignored its role as a posterior gradient, and for the same reason, we consider Fgf only as a posterior signaling gradient and neglect its oscillating behavior. Remarkably, this framework could recapitulate the somitogenesis dynamics of both the PORD self-organizing model [8] and the Clock and Wavefront hypothesis [7].
The PORD model could be implemented by setting Wnt self-enhancement to zero. Indeed, in agreement with [8], these parameters produced periodic somite patterns via a rely mechanism triggered at the anterior boundary and independent of global morphogen gradients. Numerical simulations, however, showed that these periodic patterns were fragile against noise and required initial conditions with an anterior straight line of high concentrations to avoid the formation of chessboard patterns ( Figure S13), similar to the reaction-diffusion model in [52]. Two recent studies have shown that the robustness of PORD-like patterns can be improved with additional modulation by posterior morphogenic gradients [36,20], which can also make somite patterning independent of specific initial conditions. However, our analysis shows that the underlying self-organizing behavior of PORD models corresponds to one of the six pattern-forming regimes presented by Alan Turing in [50] as "stationary waves of extreme short wave-length", where cells tend to form chessboard lateral inhibition patterns due to a diffusion-driven stabilization of cells in opposite phase, as proposed in Singh et al. [45]. This tendency was highlighted by simulations of virtual tail explants ( Figure 4B), where the mixing between cells stimulated the formation of chessboard patterns, which cannot be reconciled with the circular signaling waves observed in tail explants in [21] and [18].
The Clock and Wavefront hypothesis could be modeled by considering an Fgf-driven inhibition of Wnt self-enhancement or Notch self-enchantment depending on the type of negative feedback between Wnt and Notch ( Figure S2). In both cases, the modulation promoted homogeneous oscillations at the tip of the tail and the arrest of oscillations when cells exited the influence of Fgf. However, in contrast to several previous CW models [27], the freezing of oscillations was modeled as a bifurcation from an unstable oscillatory state to a bi-stable frozen state. This bifurcation is similar to the one proposed in a recent model of metazoan body segmentation [19] but in our model is induced by modulating only one reaction parameter in a simpler network. The Fgf modulation could also explain the switch from out of phase oscillations of Wnt and Notch at the tail tip to an in phase freezing at the anterior side [47]. Nevertheless, it could not recapitulate the formation of waves with in phase oscillations of Notch and Wnt observed in the middle part of the mouse tail [47]. Finally, virtual explants revealed that the model was also unable to form the circular signaling waves observed in tail explants in vitro [21,18].
To overcome these limitations, we developed an extended CW model where Fgf decreased Wnt self-enhancement to promote posterior oscillations and RA increased Wnt self-enhancement to promote anterior bi-stable patterns. This defined a novel central region in the tail with low levels of both Fgf and RA signaling, where cells had intermediate Wnt self-enhancement strength. We named this model the Clock and Wavefront Self-Organizing model (CWS), to highlight that in this central region cells possess the self-organizing potential to form periodic phase waves. Indeed, in contrast to previous CW models [27,17], in the CWS model the signaling waves that characterize somitogenesis were formed independent of global gradients due to local cell-communication, as showed for the Sevilletor model in Figure 1E and Figure S1. Strikingly, we found that the model could also explain the switch from out of phase to in phase oscillations of Wnt and Notch observed in the middle part of the tail during mouse somitogenesis [4,47]. The presence of this central self-organizing region was also supported by virtual explant simulations, which in agreement with experimental data generated circular signaling waves that propagated from the center of the explants [21,18]. Finally, the CWS model exhibited regeneration pattern capabilities across the whole tail and, in contrast to the PORD model, it is robust to noise ( Figure S10 and Figure S13).
The default self-organizing regime of the CWS model fitted also the idea that the segmentation clock behaves as an excitable system [18]. Indeed, in the absence of diffusion, isolated cells with intermediate Wnt self-enhancement were trapped into one of the two stable states of the system (Figure 2A for k 1 = 2.3). Nevertheless, if the same cells were coupled together via diffusion they could excite each other to maintain sustained oscillations ( Figure 2B for k 1 = 2.3 and supplementary Movie 3). A similar behavior has been observed in mouse where isolated PSM cells stop oscillating in low density cultures but could be excited to oscillate in the presence of neighboring cells [18]. It has been proposed that the cells sensed neighbors with Yap signaling, which together with Notch signaling determine the excitability threshold of the cells. The authors modeled this idea by using the FitzHugh-Nagumo equations in a single cell simulation, where the behavior of isolated cells and high density cell cultures was modeled by manually changing the excitability threshold [18]. In the CWS model, in contrast, the switching from a quiescent state to an excited oscillatory state happens with the same excitability threshold and is due to a diffusion-driven destabilization of cells in a bi-stable regime.
A prediction of the two-dimensional Sevilletor simulations presented in Figure 1E (supplementary Movie 1) is that, starting from random initial conditions, the CWS model should form spirals and stripy spiral patterns with weak and intermediate Wnt self-enhancement respectively. So why does the CWS model not form spirals at the tip and in the middle of the tail? At the tail tip, initial conditions consist of biologically plausible low levels of Wnt and Notch signaling that are perturbed by noise. This initial conditions lie in the same half-plane and therefore synchronize oscillations in space ( Figure S8). On the other hand, in the central region of the tail cells are in a stripy spiral regime but form straight phase waves without spirals. The reason for the lack of spirals in this region can be understood by considering that spirals emerge when propagating wavefronts break up into opposite waves that curl into rotatory movements [57]. The breaking of these wavefronts is initiated by different refractory times in space, which can be promoted by heterogeneities in excitable media or obstacles [57]. In the central part of the CWS tail model, this phenomenon is prevented by the homogeneous oscillations at the neighboring tip, which act as a peacemaker population to remove heterogeneity in refractory time. Moreover, as showed by the two dimensional simulation in Figure 2C, the graded self-enhancement promoted by Fgf and RA guides the alignment and movement of waves anteriorly helping to prevent breaking of wavefronts.
In summary, cells in the CWS model do not form spirals because they are constrained by initial conditions and interactions with the neighboring cells at the tip. Therefore, the CWS model is an example of a guided self-organizing process [26], where the intrinsic capacity of the PSM to form periodic phase waves is guided by global signals. One way to uncover the full self-organizing potential of the PSM is by removing the influence of these global embryonic signals in tail explants. Our simulations, however, suggest that if cells inherit the identity and concentrations possessed in the tail, explants can retain guiding cues. For example, virtual explants generated from the whole posterior tail can give rise to circular wave patterns due to a central cell population from the tip that oscillates homogeneously and guides the outer edge to propagate circular waves (third row in Figure 4D and supplementary Movie 10). On the other hand, explants generated from the middle part of the tail posses less guiding cues, because all the cells share the same identity. Nevertheless, they can still form circular wave patterns [21] (second row in Figure 4D and supplementary Movie 9). In this case, virtual explant simulations suggest that the circular patterns are promoted by spiral centers induced by cell mixing that arise primarily at the center of the explant due to the concentrations inherited from the tail. Finally, our simulations predicted that the tendency to form spiral patterns can be further uncovered by increasing the degree of cell mixing in the explant to reveal the self-organizing dynamics of the model (first row in Figure  S4D). In an experiment by Hubaud et al. [18], explants obtained by mixing cells from different individual tails gave rise to continuous oscillations that synchronized in space with an anti-clockwise rotation, suggesting that PSM cells have the self-organizing capacity to coordinate oscillations trough local cell communications. However, the lack of a clearly defined spiral signal could be associated with remaining global guides that can be induced by heterogeneity in the adherent cultures. This observation helps to address the concern that the PSM lacks hallmarks of excitable media such as spiral formation [32] by showing that the self-organizing behavior of the PSM can be hidden by globally guiding cues. For this reason, we propose that future experiments aimed to study the self-organizing behavior of the PSM should be based on non-adherent cultures and should maximize the mixing and randomization of cells.
Recently, non-adherent stem cell cultures have been used to generate human organoids that give rise to somites in vitro [55,41,25]. The main goal of these studies, however, was to improve somite reproducibility by reducing cell heterogeneities and by inducing an anterior-posterior axis. One of the studies showed that when organoids were placed in adherent cultures, they first generated inner-outer signaling waves and then formed periodic somite rosettes via cell sorting [25]. This demonstrates that adhesive culture conditions can indeed modulate somitogenesis and suggests that additional self-organizing behaviors drive the periodic morphogenesis of somites independent of the clock, as recently showed in Zebrafish [44].
In conclusion, the Sevilletor model provides a minimal theoretical framework to study multi-cellular pattern formation by synchronized oscillations. We applied this framework to investigate mouse somitogenesis and proposed the Clock and Wavefront Self-Organizing model (CWS) to recapitulate the behavior of the PSM both in vivo and in vitro. The model forms somitogenesis phase waves in a self-organizing manner and recapitulates for the first time the formation of circular waves observed in tail explants. In addition, it provides the first mechanistic explanation for the switch from out phase to in phase oscillations of Wnt and Notch observed in the mouse tail. Finally, it explains how PSM cells can be excited to oscillate from a quiescent state via local cell communication.

Acknowledgments
This work is supported by grant PGC2018-095170-A-100 from the Spanish Ministry of Science and Innovation and J.K. is supported by a doctoral fellowship PRE2019-087911, from the Spanish Ministry of Science and Innovation, Spain. We thank Dr. Miki Ebisuya for providing feedback.

Declaration of interests
The authors declare no competing interests.

Data and code availability
The data and codes, and any additional information required to reanalyze the data reported in this paper is available from the lead contact upon request.

Method details
Simulations are run with a finite difference solver written in Julia 1.6.4. The complex systems analysis and two cell simulations have been run in Mathematica 12, the reverse-LUT method has been written in Python 3.8.10, and Fiji (ImageJ) has been used to analyze the patterns seen in the simulated explants. Python 3.8.10 is used to create all plots, excluding the phase spaces and bifurcation diagrams which have been generated using Mathematica 12.

Space discretization
The diffusion term ∇ 2 u(x, y) in the system is calculated using a first order finite difference scheme to calculate the discrete Laplace. This is done by using a Taylor expansion of the functions u(x, y, t) and v(x, y, t) around the steady state point (0,0). ǫ = dx = dy in the following. For u the calculations are: u(x, y + ǫ) = u(x, y) + ǫ du(x, y) dx + ǫ du(x, y) dy We isolate the first order differential terms in the equations and combine the equations to get: Similar calculations are also performed for v.

Boundary conditions
In all simulations, except the one in Figure S1B, we consider zero flux Neumann boundary conditions. This means that nothing diffuses between the outside and inside of the system, i.e. the Laplace is equal to zero at the boundaries. The Laplace at the boundary is derived as: du(x, y) This leads to: u(x, y + ǫ) = u(x, y − ǫ).
Inserting the appropriate ones for the edge/corner into the discrete Laplace equation (13) gives the correct formula, for example for the bottom left corner (x, y) = (1, 1) in equation (18): The periodic boundary conditions used for the simulation in Figure S1B are also calculated using the discrete Laplace, treating cells on the boundaries as direct neighbors with cells on the corresponding boundary. For example for the cell in the bottom left corner (x, y) = (1, 1) in a system where the lengths of the axes are L: -Clock and Wavefront Self-Organizing model: every 100th time step (for the perturbed version with faster growth in Figure S10 it is 50th time step)

Stability and types of steady states
The steady states of the Sevilletor model and their properties are calculated with a phase plane analysis [28]. The steady states of a system are all pairs (u * , v * ) for which f (u * , v * ) = du dt (u * ,v * ) = 0 and g(u * , v * ) = dv dt (u * ,v * ) = 0. The stability of the steady states are found by using linear stability analysis from the determinant and trace of the matrix A: For Det(A (u * ,v * ) ) > 0 and Trace(A (u * ,v * ) ) < 0 the steady state is stable, otherwise it is unstable. The type of steady state is identified as discussed in Appendix A in [28] and describes the shape of the vector field around the steady state. Figure 2 The phase is calculated, at every time point, as the angle between the initial position and the position of the cell in phase space around the central unstable steady state in (u, v) = (0, 0).

Calculation of the phase and the number of oscillations performed in
The number of oscillations is determined by counting the local minima that are close to zero along the phase profile, excluding the first minima for t = 0.  Figure STAR1: Details of the calculation of phase and number of oscillations in Figure  2 A) The timeline shows the path of a cell in the phase space. The phase is calculated for each time point, starting from 0, increasing to π for half a loop, and decreasing back to 0 for the second half of the loop.
B) The phase is plotted as a function of time and the number of oscillations is measured as the number of local minima in phase = 0 excluding the initial for t = 0, marked by black points.

Explant quantification with circular profiles
We developed a novel method to quantify the patterns formed in explants by deriving space time plots from circular profiles. These profiles are generated by a macro written in Fiji, which defines a ring of pixels with a determined radius and width that is placed around a manually defined point corresponding to the center of the pattern. Average values along the ring are measured with the Plot Profile function of Fiji at each time point to generate space-time plots. In the analysis in Figure  4, we consider two concentric profiles with a smaller and a larger radius to generate two space-time plots. In Figure 4 on the right we also calculated phase profiles by averaging the values of the space-time plots along the x-axis. The resulting graphs showed the average values and standard deviation of the concentration over time for smaller and larger circular profiles. Comparison of the phase between the two space-time plots allowed us to distinguish between homogeneous oscillations and circular waves (target wave patterns).  [4,47].
A) The standard implementation of the models presented in Figure 3 with a positive interaction from Wnt to Notch and an inhibition from Notch to Wnt. In this case Fgf and RA modulate the self-enhancement of Wnt represented by the parameter k 1 . Only the CWS model can generate realistic patterns with the correct relative phase between Notch and Wnt at the posterior, middle and anterior part of the tail (three green marks). B) When the incoherent negative feedback between Notch and Wnt is inverted and Fgf and RA modulate k 1 as in the standard implementation, none of the models can reproduce the relative phase between Notch and Wnt observed in experiments, i.e. all contain at least two red cross marks. C) When the incoherent negative feedback between Notch and Wnt is inverted and Fgf and RA modulate k 2 , which represents the self-amplification of Notch, similar results to the standard implementation can be obtained. Also in this case, only the CWS model can generate patterns with a relative phase between Notch and Wnt that fit experiments in the posterior, middle and anterior part of the tail (three green marks). The variability across simulations depends on the tail simulation stage used to generate the explant, which is randomly selected around the time point showed on the left, and on the stochastic shuffling between cells performed in the making of the explants. A) Explants obtained from the tip give rise to homogeneous oscillations in 100% of the simulations (22/22). B) Explants obtained from the middle section of the tail (without the tip) give rise to three possible outcomes: circular wave patterns with a single center made from two spirals in 27% of simulations (6/22), a single concentric wave that dissipates into a homogeneous pattern in 59% of the cases (13/22), or two concentric waves that merge together and are promoted by four spirals (3/22). In real middle tail explants, the second and the third case may be prevented by modulations provided by adherent cultures that may promote a behavior more similar to the tip+middle explants showed in panel C. C) Explants obtained from tip+middle part of the tail generate circular wave patterns in 100% of the simulations (22/22). In this case, a small portion of the explant center contains cells coming from the tip that prevent spiral formation by oscillating homogeneously. The space-time plots on the right show that in both explants cells form frozen patterns that do not recapitulate the rotatory waves of the experiment in [18]. D) In the tip+middle case (first row), mixed explants of the CWS model give rise to rotatory spiral waves that resemble the rotatory patterns observed in [18]. In the simulation, this rotatory pattern arises due to the self-organizing properties of the middle section of the tail and the heterogeneity induced by cell randomization. In the tip+middle+top case (second row), mixed explants of the CWS model give rise to patterns with frozen regions that resemble the frozen patterns observed in [49]. This suggests that the mixed tail explants generated in [49] contain cells from the anterior part of the tail that promote the freezing of oscillations.

S1 Simulation with only positive values of u and v
To obtain simulations with only positive concentration of u and v centered around a positive value c, equations (1) and (2) can be rewritten by u → (u − c) and v → (v − c). The result is the equations (S22) and (S23). An example is shown in Figure S5B, where a stripy spiral pattern is formed centered around the value c = 2 with k 1 = 2.3 and k 2 = 1. Real part Imaginary part Figure S6: (Caption on next page.) Figure S6: Dispersion relation diagrams for the steady states of the example systems in Figure 1E. A) Phase spaces (as shown in Figure 1F). B) Dispersion relation diagrams for the steady states of the patterns shown in Figure 1E. The dispersion relation diagrams show the real (solid line) and imaginary (dashed line) parts of the eigenvalues as a function of the wavenumber (wavenumber = 2Π/wavelength). Where the real part of the eigenvalue is positive the wavelength (corresponding to the wavenumber) is unstable. If the real part of the eigenvalue is positive for wavenumber=0, the system is unstable without diffusion.
If it is negative, the system is stable without diffusion.
Brusselator equations [38]: We have simulated spirals with the Brusselator equations with the parameters found in the paper by Torabi and Davidsen [48]: A = 1.9, B = 4.8, D X = 1.0 and D Y = 0.7, L = 256. We set the initial values around (X, Y ) = (1.9, 2.52632) and used zero flux boundaries.
Difference between the Sevilletor model and classic models of the Belousov-Zhabotinsky reaction: The pattern generated by the Sevilletor model near the bifurcation point (k 1 = 2.3) is a new type of spiral patterning behavior that to the best of our knowledge has not been described previously ( Figure 1E third row, right column and supplementary Movie 1). Although the pattern looks similar to those formed by classical models of the the Belousov-Zhabotinsky reaction [56] such as the Brusselator [37] and Oregonator [13], it emerges from a different dynamical behavior. The classic models [37,13] generate spiral patterns based on a limit cycle around a single unstable steady state ( Figure S7), as shown for the spiral pattern presented in this paper with k 1 = 1 ( Figure 2B second row). In contrast, the stripy spiral patterns formed by the Sevilletor model with k 1 = 2.3 arises from a diffusion-driven excitation of a bistable regime (third row of Figure 2B and supplementary Movie 3). Simulations of spiral formation (top row) and corresponding phase spaces (bottom row). The Van der Pol [14] (A), Brusselator [38,48] (B) and Sevilletor "smooth" spirals (C) have a single unstable steady state and a limit cycle in the phase space. The Sevilletor stripy spirals (D) has five steady states, including two stable ones. However, visually the Van der Pol is similar to the smooth spirals, and the Brusselator similar to the stripy spirals.
This difference is also reflected by different patterning dynamics. Our analysis shows that in the Brusselator model, spirals are formed out of local heterogeneities in the initial conditions. Analysis of the phase space shows that these neighboring cells follow a spiral trajectory towards the central steady state, temporarily transforming the unstable steady state of the system into an effectively stable one. Following this trajectory, neighboring cells reach the steady state and eventually synchronize their phases following the original limit cycle. This event is associated with the disappearance of spiral centers. These initial spiral centers are formed independently of boundary conditions and are required to produce a stripy phase pattern.
In contrast, with k 1 = 2.3, the Sevilletor model generates spirals that appear and disappear dynamically over time while moving in space (third column in supplementary Movie 1). Numerical simulations show that these spirals are formed only with zero flux boundary conditions, suggesting that the stripy spiral appearance is due to crashing of phase waves reflected by zero flux boundary conditions. With periodic boundary conditions, the long term behavior of the system is to generate a periodic phase wave pattern of straight stripes ( Figure S1).

FitzHugh-Nagumo model:
Like the Sevilletor model, another model that was also developed from the Van der Pol equations is the famous FitzHugh-Nagumo model, proposed by FitzHugh in 1961 [14] to study neuronal firing. In this model, an additional excitation term in equation (1) and a negative linear term in equation (2) was added to investigate how a stable system could be excited to become unstable. Unlike the FitzHugh-Nagumo model, the Sevilletor model always has an instability in the steady state (u * , v * ) (0,0) = (0, 0) that is driven by the two positive self-enhancing feedbacks of u and v. Moreover, while the FitzHugh-Nagumo model can have only up to three steady states, due to the linear shape of the v nullcline, in our model both v and u have nullclines with non-linear shapes that can intersect on more than three points.

S4
The initial values of u and v affect the outcome of the patterning process

S7
The original PORD Model by Cotterell et al. [8] The PORD model by Cotterell et al. [8] freezes cells in opposite ends of the phase space as observed in the lateral inhibition patterning regime of the Sevilletor model, as shown by the two cell simulations in Figure S11A. The cell starting in (0,0) makes one single loop around the steady state, before freezing at the bottom left corner while the other cell freezes at the top right corner. This is a non-symmetric phase space, while in our PORD implementation ( Figure S11B) the phase space and the path of cells is symmetric. This difference together with use of a heavy side function to prevent negative values may explain why the periodic somite patterns formed by the original PORD model are separated by several cells, while in our model and in other implementations of the PORD hypothesis [20] are separated by only one cell. Despite this difference, all PORD implementation exhibit the same underlying behaviour based on the freezing of oscillations in opposite phase due to diffusion between neighboring cells.
The Sevilletor equations for the lateral inhibition behavior used to make an implementation of the PORD model are (equations (1) and (2) repeated for convenience): With k 1 = 0.0 and k 2 = 1.0. D = 1.0.

S10 Supplementary movies
Movies of two-dimensional square simulations of the Sevilletor model, two-cell simulations in phase space and somitogenesis simulations of the tail and tail explants for the PORD, CW and CWS model. The number of time steps in each movie varies to illustrate the dynamics of the patterning processes.
Movie 1: Self-organizing patterning behaviors of the Sevilletor model From left to right, two dimensional simulations of the Sevilletor model started with random initial concentrations for increasing values of k 1 , as shown in Figure 1E. With k 1 = 0.0 the model forms a lateral inhibition chessboard pattern, with k 1 = 1.0 spirals, with k 1 = 2.3 stripy spirals, with k 1 = 3.0 a homogeneous pattern and with k 1 = 4.0 bi-stable salt and pepper pattern. Without diffusion (graph on the left) the two cells oscillate by following a limit cycle around the central steady state. In the presence of diffusion (graph on the right) the two cells stop at opposite phase where the changes promoted by reaction (black arrows) are balanced by the changes promoted by diffusion (green arrows). This behavior generates chessboard patterns as shown in Figure 1E.

Movie 3: Phase space trajectories of two cells in stripy spiral patterning
Two-cell simulations in phase space for the stripy spiral regime with k 1 = 2.3 in the absence or presence of diffusion, as showed in Figure 2A-B. The trajectories of the two cells with initial values (u, v) = (±0.1, 0.0) are shown in teal colored lines, the gray points are unstable steady states and white points are stable steady states. Without diffusion (graph on the left) each cell is pushed by reaction (black arrows) into one of the stable steady states. In the presence of diffusion (graph on the right) the influence of reaction is counterbalanced by diffusion (black and green arrows) that pushes each cell onto the trajectory of the other unstable gray point creating coordinated loops. This behavior generates stripy spiral patterns as shown in Figure 1E. The behavior of the Sevilletor model can be easily modulated to transition between patterning regimes over time by varying the parameter k 1 in all cells. In this example, the model transitions from lateral inhibition pattern, to rotating spirals, to stripy spirals, to a bi-stable frozen spiral pattern, and back to lateral inhibition pattern: k 1 = 0.0 → 1.0 → 2.3 → 4.0 → 0.0.

Movie 5: A gradient of k 1 guides self-organizing patterning
Square two-dimensional simulation with random initial conditions and a gradient of the parameter k 1 along the y-axis, as showed in Figure 2C. The gradient starts with low values at the bottom (k 1 = 0) that increase linearly towards the top (k 1 = 4). From bottom to top, the Sevilletor model smoothly transitions between lateral inhibition pattern, spirals, stripy spiral patterns, homogeneous patterns and bi-stable frozen patterns. In addition, over time, the gradient modulates the stripy spiral region to form straight phase waves that move towards the top following increasing values of k 1 .

Movie 6: Sevilletor implementation of the PORD model
Two dimensional simulation of the PORD model implemented in the Sevilletor framework, showed in Figure 3C. Wnt and Notch signaling are shown side-by-side. All cells have k 1 = 0.0 and the simulations start with high signaling at the top boundary (anterior region). The virtual tail grows with a constant speed by addition of a line of cells at the bottom (posterior). The model generates a periodic somite pattern with a one-cell width independent of gradients due to a rely mechanisms based on lateral inhibition. Simulated explant obtained from the middle part of the tail in the Clock and Wavefront Self-Organizing model. The simulation recapitulates the circular Notch signaling waves that propagate from the center of the colony towards the periphery observed in experiments, simulated Wnt and Notch signaling are shown side-by-side. The virtual explant is generated by making a radial projection of average values along the y-axis for the middle part of the tail, with cells having k 1 = 2.3 that generate circular wave patterns due to two spiral centers that rotate in opposite directions. This type of pattern depend on initial condition and arise in 27% (6/22) of simulations ( Figure S3).