Arnold tongue entrainment reveals dynamical principles of the embryonic segmentation clock

  1. Paul Gerald Layague Sanchez
  2. Victoria Mochulska
  3. Christian Mauffette Denis
  4. Gregor Mönke
  5. Takehito Tomita
  6. Nobuko Tsuchida-Straeten
  7. Yvonne Petersen
  8. Katharina Sonnen
  9. Paul François  Is a corresponding author
  10. Alexander Aulehla  Is a corresponding author
  1. European Molecular Biology Laboratory (EMBL), Developmental Biology Unit, Germany
  2. McGill University, Canada
  3. European Molecular Biology Laboratory (EMBL), Transgenic Service, Germany

Abstract

Living systems exhibit an unmatched complexity, due to countless, entangled interactions across scales. Here, we aim to understand a complex system, that is, segmentation timing in mouse embryos, without a reference to these detailed interactions. To this end, we develop a coarse-grained approach, in which theory guides the experimental identification of the segmentation clock entrainment responses. We demonstrate period- and phase-locking of the segmentation clock across a wide range of entrainment parameters, including higher-order coupling. These quantifications allow to derive the phase response curve (PRC) and Arnold tongues of the segmentation clock, revealing its essential dynamical properties. Our results indicate that the somite segmentation clock has characteristics reminiscent of a highly non-linear oscillator close to an infinite period bifurcation and suggests the presence of long-term feedbacks. Combined, this coarse-grained theoretical-experimental approach reveals how we can derive simple, essential features of a highly complex dynamical system, providing precise experimental control over the pace and rhythm of the somite segmentation clock.

Editor's evaluation

This manuscript undertakes a fundamental study of the entrainment of the mouse segmentation clock to an external oscillating signal. This provides insight into both the nonlinear properties of the clock as well as how it might integrate global signals. By combining experimental measurement of the finite pulse phase response curve with theoretical analysis using tools from nonlinear dynamics, this work produces a compelling coarse-grained model of the segmentation clock, albeit in an in vitro setting. In future, it will be of great interest to see how the predictions of this coarse-grained view, in particular the consequences of operating near a bifurcation, hold up in an in vivo system. This work is likely to be of interest to developmental biologists, synthetic biology experts, theoreticians and others interested in the emergence of complex patterning in living systems.

https://doi.org/10.7554/eLife.79575.sa0

Introduction

How do we gain insight into a complex system, which exhibits emergent properties that reflect the integration of entangled interactions and feedback regulation? As pointed out in the late 1970s by David Marr and Tomaso Poggio in their seminal paper (Marr and Poggio, 1976), understanding the complexity encountered when studying the ‘nervous system or a developing embryo’ requires the analysis at multiple levels of organization. Their core tenet is that also in biological systems, different levels of organization, while obviously causally linked, exhibit only a loose connection and importantly, can be studied and understood independently from each other.

Such observations are not specific to biology and have been made more quantitative in other fields. In physics, renormalization techniques coarse-grain degrees of freedom to obtain scale-free theories, allowing to define universality classes independent of the precise details of interactions (Cardy, 1996; Bialek, 2015). Another recent example is the parameter space compression theory, showing how complex systems (in biology or physics) can be typically reduced to simpler descriptions with few parameters (Machta et al., 2013; Proulx-Giraldeau et al., 2017; Hsu et al., 2020).

Going one step further, this suggests that one might be able to study – and control – complex systems provided we identify the essential, macro-level behavior. This is possible because only a limited number of universal descriptions exist, with defining behaviors and properties, that do not depend on the detailed implementation. A central challenge that remains is to implement these theoretical ideas to the experimental study of biological complexity.

Here, we develop a coarse-grained approach combining theory and experiments to study a cellular oscillator ensemble that constitutes the embryonic somite segmentation clock. Functionally, this clock controls the periodic formation of somites, the precursors of the vertebral column, and other tissues (Palmeirim et al., 1997; Krol et al., 2011). Molecularly, the segmentation clock comprises the oscillatory activity of several major signaling pathways, such as the Notch, Wnt, and Fgf signaling pathways, which show oscillatory dynamics with a period matching somite formation, that is, ∼2 hr in mouse embryos (Dequéant and Pourquié, 2008; Hubaud and Pourquié, 2014; Aulehla et al., 2003; Dequéant et al., 2006; Aulehla et al., 2008). More recently, segmentation clock oscillations with a period of ∼5 hr have been identified in human induced pluripotent stem cells differentiated into paraxial mesoderm, identifying a set of ∼200 oscillating genes, including targets of Notch and Wnt signaling (Chu et al., 2019; Diaz-Cuadros et al., 2020; Matsuda et al., 2020).

Strikingly, as individual oscillating cells are coupled to their neighbors via Notch-Delta signaling, the oscillations occur synchronized and wave-like activity patterns appear to periodically sweep along the embryonic anterior-posterior axis (Maroto et al., 2005; Masamizu et al., 2006; Aulehla et al., 2008; Herrgen et al., 2010; Webb et al., 2016; Yoshioka-Kobayashi et al., 2020; Rohde et al., 2021).

Adding to the complexity, these periodic spatiotemporal wave patterns are linked to an underlying spatial period gradient along the embryonic axis, that is, signaling dynamics in cells close to the posterior of the embryo oscillate faster compared to those in cells located more anteriorly. Such a period gradient linked to the segmentation clock has been identified in several species (Gomez et al., 2008; Morelli et al., 2009; Oates et al., 2012; Lauschke et al., 2013; Soroldoni et al., 2014; Tsiairis and Aulehla, 2016; Falk et al., 2022) and also in in vitro assays culturing intact or even dissociated presomitic mesoderm (PSM) (Lauschke et al., 2013; Tsiairis and Aulehla, 2016).

Of note, an analogous oscillatory system was also described during segmentation in arthropods (Clark et al., 2019) and while distinct at molecular level, it also exhibits spatiotemporal wave patterns traversing the embryo axis, again with indication of a period gradient (Sarrazin et al., 2012; El-Sherif et al., 2012; Brena and Akam, 2013; Schönauer et al., 2016).

In this work, we coarse-grain these underlying complexities and take a dynamical systems, macro-perspective on the segmentation clock, studying it as a single phase oscillator (Figure 1A). We build on the theory of synchronization and entrainment (see below) to first perform a systematic experimental characterization of its response to perturbation. We compare the outcome to qualitative and quantitative theoretical predictions. In turn, these experimental quantifications allow to derive a phase response curve (PRC) that uniquely characterizes the dynamical properties of the segmentation clock. This new insight provides the means to understand – and control – the timing of a complex embryological patterning process.

Figure 1 with 1 supplement see all
Entrainment of an embryonic oscillator with a zeitgeber: setup and theory.

(A) Schematic of the experimental entrainment setup using a microfluidics device previously described in Sonnen et al., 2018, and an overview of the analysis approach used in this present study. (1) The posterior presomitic mesoderm (PSM) is recovered from E10.5 mouse embryo and used in a quasi-two-dimensional segmentation assay (‘2D-assay’) (Lauschke et al., 2013). (2) The embryonic tissue is loaded into a microfluidics device bonded to fibronectin-coated cover glass, and is imaged using a ×20 objective at 37°C and 5% CO2. (3) Simultaneously, the 2D-assay is subjected to periodic pulses of a drug (e.g. DAPT for Notch signaling), which serves as zeitgeber. (4) Using dynamic fluorescent reporter for segmentation clock genes (e.g. LuVeLu for Notch signaling), we quantify endogenous signaling oscillations during the experiment. (5) We generate a single timeseries using reporter signal from the entire sample (‘global region of interest [ROI]’) that allows to quantify the segmentation clock rate and rhythm (see validation in Figure 2 and Figure 7—figure supplement 3). Illustration by Stefano Vianello. A photo of the actual microfluidics device and its design are shown in Figure 1—figure supplement 1. (B) Abstract definition of phase: two different points in the plane (x,y) have the same phase if they converge on the same point on the limit cycle (indicated in gray). (C) Perturbations change the phase of the cycle by an increment Δϕ (here phase is defined by the angle in the plane). (D) Time courses with similar perturbation as C showing the oscillations of x and phase difference Δϕ. (E) Illustration of generic Arnold tongues, plotted as a function of zeitgeber strength (ϵ) and zeitgeber period (Tzeit), mapping n:m entrainment where the entrained oscillator (with natural period of Tosc) goes through n cycle/s for every m cycle/s of the zeitgeber. Three different points in the 1:1 and 2:1 Arnold tongues are specified with corresponding graphical illustration of an autonomous oscillator as it is subjected to zeitgeber with different periods (Tzeit): when Tzeit is less than Tosc (F), when Tzeit is greater than Tosc (G), and when Tzeit is much much greater than Tosc but is close to twice of Tosc (H). Free-running rhythm of the oscillator (i.e. before perturbation) is marked by a dashed line, while solid line illustrates its rhythm during perturbation with zeitgeber. Magenta bars represent the zeitgeber pulses. Illustration by Stefano Vianello.

Theory of synchronization guides the experimental study of segmentation clock entrainment

Our experimental study is based on and guided by the theory of entrainment of oscillators by an external periodic signal – a subset of the general theory of synchronization (Pikovsky et al., 2001).

Entrainment is observed when an autonomous oscillator adapts its behavior to lock to an external periodic signal (called zeitgeber in the circadian rhythm literature). The general theoretical framework to understand entrainment requires the definition of oscillator phases (Figure 1B–D), and their response to perturbation (Kuramoto, 2003). Assuming the zeitgeber consists in periodic pulses, entrainment is observed when the phase of the oscillator ϕent at the time of the zeitgeber is constant (technically, a fixed point of the Poincaré return map [Cross and Siggia, 2005]). This defines period-locking (also termed mode-locking) (Pikovsky et al., 2001). Entrainment is not always manifested and conditions for its existence can be derived. Quantitatively, when entrainment occurs, the zeitgeber induces a periodic phase perturbation (or response) of the entrained oscillator, which exactly compensates the detuning (or period mismatch, Tzeit-Tosc) between the zeitgeber and the free-running oscillator. For this reason, when the detuning is very small, a weak external perturbation is enough to entrain an oscillator. Conversely, if the detuning is big, a strong signal and associated response is required for entrainment. One can then plot the minimal strength of the zeitgeber (ϵ) versus the corresponding detuning (or simply Tzeit, if Tosc is constant): these maps are more commonly known as Arnold tongues (Figure 1E). Arnold tongues predict the period- and phase-locking behavior in oscillatory systems as different as electrical circuits (Balanov, 2008), oscillatory chemical reactions (Ito, 1979; Kai and Tomita, 1979; Buchholz et al., 1985; Makki et al., 2014), or living systems such as circadian rhythms (Granada et al., 2013).

Lastly, more complex patterns of entrainment can be observed: for instance, stable phase relationships can be established where the entrained oscillator goes through n cycles for every m cycles of the external signal, defining n:m period-locking. In that case, the instantaneous period of the oscillator matches m/n the period of the zeitgeber (Tzeit). Corresponding Arnold tongues can be obtained, leading to a rich structure for entrainment in parameter space (Figure 1F–H).

Previously, periodic activation of Notch signaling via heat shock-driven expression of the ligand Delta was used to modulate the segmentation of PSM in zebrafish (Soza-Ried et al., 2014). In the cited study, the readout to assess the effect of the periodic perturbation relied on somite size and morphology. While this gave important insight into how morphological segmentation of the PSM and somite length are affected by a zeitgeber, experimental investigation on the underlying signaling dynamics upon periodic perturbation remains limited.

In this work, to experimentally apply the theory of synchronization to the segmentation clock, we make use of a microfluidics-based entrainment setup, which we had established previously in the lab (Figure 1—figure supplement 1; Sonnen et al., 2018).

We showed before that using a quasi-two-dimensional in vitro segmentation assay (hereafter referred to as a 2D-assay), which recapitulates segmentation clock dynamics and PSM patterning (Lauschke et al., 2013), the microfluidics-entrainment approach allowed us to take control of Notch and Wnt signaling oscillations, providing direct functional evidence that the oscillation phase shift between Wnt and Notch signaling is critical for PSM patterning (Sonnen et al., 2018).

Results

A coarse-grained, single oscillator description of the segmentation clock

To perform a systematic analysis of entrainment dynamics, we first introduced a single oscillator description of the segmentation clock. We used the segmentation clock reporter LuVeLu, which shows highest signal levels in regions where segments form (Aulehla et al., 2008). Hence, we reasoned that a global region of interest (ROI) quantification, averaging LuVeLu intensities over the entire sample, should faithfully report on the segmentation rate and rhythm, essentially quantifying ‘wave arrival’ and segment formation in the periphery of the sample. Indeed, our validation confirmed that the oscillation period based on the global ROI analysis using the LuVeLu reporter closely matched the rate of morphological segment boundary formation (Figure 2D, Figure 2—figure supplement 1). In addition, global ROI quantifications using additional reporters to measure Wnt signaling oscillations (i.e. Axin2, Figure 2C–D, Figure 2—figure supplement 1) and the segmentation marker Mesp2 (Figure 2D, Figure 2—figure supplement 1) also showed close correspondence with LuVeLu global ROI measurements of segmentation clock period.

Figure 2 with 4 supplements see all
Quantifying the segmentation clock rhythm.

(A) Comparison of measurements obtained from a region of interest (ROI) spanning the entire field (‘global ROI’) with those obtained from ROIs in which central regions of increasing size were excluded. The excluded area is represented in terms of the diameter (in pixels, px) of a circular region at the center. The corresponding timeseries are shown in the top panel and marked with different colors. Bottom panel shows the global ROI and excluded region in a snapshot of signal of LuVeLu, a dynamic reporter of Notch signaling driven from the Lfng promoter (Aulehla et al., 2008), in a 2D-assay at 700 min from the start of the experiment. Henceforth, timeseries is obtained using ‘global ROI’, unless otherwise specified. (B) Detrended timeseries of the segmentation clock (obtained using global ROI) in 2D-assays, which express the LuVeLu reporter, subjected to 170 min periodic pulses (magenta bars) of either 2 µM DAPT (in solid line) or DMSO (for control, in dashed line). Timelapse movie and corresponding timeseries from global ROI are available at https://youtu.be/fRHsHYU_H2Q. (C) Detrended timeseries of Axin2-linker-Achilles reporter in 2D-assays, subjected to 170 min periodic pulses (magenta bars) of either 2 µM DAPT (in solid line) or DMSO (for control, in dashed line). (D) Period of morphological segment boundary formation, LuVeLu oscillation, Axin2-linker-Achilles oscillation, and rhythm of Mesp2-GFP in 2D-assays subjected to 170 min periodic pulses of either 2 µM DAPT (or DMSO for control). Each sample is represented as a dot, the median is denoted as a solid horizontal line. For morphological segment boundary formation, period was determined by taking the time difference between two consecutive segmentation events (for CTRL: 17 segmentation events in six samples from four independent experiments, and for DAPT: 24 segmentation events in seven samples from five independent experiments, with p-value < 0.001) in the brightfield channel. For period quantifications based the reporters, the mean period per sample from 650 to 850 min after start of the experiment was plotted. LuVeLu: (CTRL: n=24 and N=7) and (DAPT: n=34 and N=8) with p-value < 0.001, Axin2-linker-Achilles: (CTRL: n=3 and N=1) and (DAPT: n=5 and N=1) with p-value = 0.107, Mesp2-GFP: (CTRL: n=8 and N=3) and (DAPT: n=9 and N=3) with p-value = 0.01. Data were visualized using PlotsOfData (Postma and Goedhart, 2019). To calculate the p-value, two-tailed test for absolute difference between medians was done via a randomization method using PlotsOfDifferences (Goedhart, 2019). The timeseries and corresponding period evolution during entrainment, obtained from wavelet analysis, are in Figure 2—figure supplement 1A and Figure 2—figure supplement 1B, respectively. Timelapse movies are available at https://youtu.be/edFczx_-9hM and https://youtu.be/tQeBk0_U_Qo, respectively.

Hence, we conclude that LuVeLu global ROI timeseries analysis provides a valid coarse-grained quantification of the pace (period ∼ rate of segment formation) and rhythm (phase) of the segmentation clock.

The pace of the segmentation clock can be locked to a wide range of entrainment periods

Having established a quantitative, coarse-grained readout for segmentation clock pace and rhythm, we next analyzed whether pace and rhythm can be experimentally tuned using microfluidics-based entrainment (Figure 1A, Figure 1—figure supplement 1).

First, we tested whether the segmentation clock can be entrained to periods different and far from the endogenous period of ∼140 min, which we refer to as its free-running, natural period or Tosc. To address this question, we modified the entrainment period from 120 to 180 min, while keeping the drug concentration and the pulse duration (i.e. 30 min/cycle) constant. Our results show that, while controls cycled close to Tosc (Figure 3A–C, Figure 3—figure supplement 1), the segmentation clock rhythm in DAPT-entrained samples closely adjusted to Tzeit over the specified range (Figure 3D, Figure 3—figure supplement 1, Figure 3—figure supplement 2). Hence, we were able to speed up and slow down the pace of the segmentation clock system using entrainment.

Figure 3 with 2 supplements see all
The segmentation clock can be locked to a wide range of entrainment periods.

(A) Detrended (via sinc filter detrending) timeseries of the segmentation clock in 2D-assays subjected to 170 min periodic pulses of 2 µM DAPT (or DMSO for controls). Periodic pulses are indicated as magenta bars and the timeseries of each sample (for CTRL: n=24 and N=7, for DAPT: n=34 and N=8) is marked with a dashed line. (B) Period evolution during entrainment, obtained from wavelet analysis. The period evolution for each sample and the median of the periods are represented here as a dashed line and a solid line, respectively. The gray shaded area corresponds to the interquartile range. Magenta dashed line marks Tzeit. (C) Mean period from 650 to 850 min after start of the experiment of samples subjected to 170 min periodic pulses of 2 µM DAPT (or DMSO for controls), with p-value <0.001. Each sample is represented as a dot, while the median of all samples is denoted as a solid horizontal line. This plot is the same as the plot for the LuVeLu condition in Figure 2D. (D) Mean period from 650 to 850 min after start of the experiment of samples entrained to periodic pulses of 2 µM DAPT. Each sample is represented as a dot, while the median of all samples is denoted as a solid horizontal line. The period of the DAPT pulses is specified (for 120 min: n=14 and N=3, for 130 min: n=39 and N=10, for 140 min: n=15 and N=3, for 150 min: n=17 and N=4, for 160 min: n=15 and N=3, for 170 min: n=34 and N=8, for 180 min: n=6 and N=1). Data were visualized using PlotsOfData (Postma and Goedhart, 2019), and a summary is provided in Table 1. A similar plot including each condition’s respective control is in Figure 3—figure supplement 1. The analysis of period and wavelet power across time is summarized in Figure 3—figure supplement 2. To calculate the p-value, two-tailed test for absolute difference between medians was done via a randomization method using PlotsOfDifferences (Goedhart, 2019).

Notably, period-locking was less precise (i.e. higher standard deviation as shown in Table 1) with 120 and 180 min zeitgeber periods, a possible indication that the limit of entrainment range is approached at these conditions.

Table 1
Summary statistics on period-locking of the segmentation clock in E10.5 2D-assays to periodic pulses of 2 µM DAPT.

This table summarizes the median, 95% confidence interval (CI) of the median, mean, standard deviation (SD), and standard error of the mean (SEM) of the segmentation clock in 2D-assays subjected to periodic pulses of 2 µM DAPT. These summary statistics were determined using PlotsOfData (Postma and Goedhart, 2019). A plot of these data is shown in Figure 3D.

Pulse periodnNMedian, min95% CI of median, minMean, minSD, minSEM, min
120 min143127.36123.16–148.03139.3922.286.18
130 min3910133.44132.21–135.97134.754.560.74
140 min153147.68146.10–153.02150.007.421.98
150 min174155.14154.33–160.74158.187.771.94
160 min153167.09166.16–169.84167.335.281.41
170 min348176.14173.06–181.20175.568.521.48
180 min61184.43164.54–194.21181.0617.888.00

We also tested the effect of varying DAPT concentration, as means of changing zeitgeber strength (ϵ in Figure 1), on entrainment dynamics. Synchronization theory predicts that zeitgeber strength correlates with the time it takes to reach period-locking (Pikovsky et al., 2001; Granada and Herzel, 2009). To test this prediction, we entrained samples with periodic DAPT pulses at fixed intervals of 170 min and varied drug concentration. We indeed found that the time needed to show period-locking was shortened in samples using higher DAPT concentrations (Figure 4A–B, Figure 4—figure supplement 1A). Additionally, as expected, higher drug concentrations also resulted in more robust entrainment, indicated by the quantifications of the first Kuramoto order parameter R, a measure for in-phase synchrony between samples (Figure 4C–D, Figure 4—figure supplement 1B, Figure 4—figure supplement 2) (see definition in Appendix 1, notice that R is the modulus of the mean field variable used in Riedel-Kruse et al., 2007).

Figure 4 with 2 supplements see all
Effect of varying DAPT concentrations on entrainment dynamics.

(A) Left: Detrended (via sinc filter detrending) timeseries of the segmentation clock in 2D-assays entrained to 170 min periodic pulses of either 1 or 2 µM DAPT, zoomed in from 500 to 950 min. Periodic pulses are indicated as magenta bars and the timeseries of each sample (for 1 µM: n=18 and N=5, for 2 µM: n=34 and N=8) is marked with a dashed line. The median of the oscillations is represented here as a solid line, while the gray shaded area denotes the interquartile range. Right: Detrended (via sinc filter detrending) timeseries of the segmentation clock in 2D-assays entrained to 170 min periodic pulses of either 1 µM (n=18 and N=5) or 2 µM (n=18 and N=8) DAPT represented as heatmaps, generated using PlotTwist (Goedhart, 2020). Periodic pulses are indicated as magenta bars. Each row corresponds to a sample. (B) Period evolution during entrainment, obtained from wavelet analysis. The period evolution for each sample and the median of the periods are represented here as a dashed line and a solid line, respectively. The gray shaded area corresponds to the interquartile range. The plot for the 2 µM condition is the same as the plot for the DAPT condition in Figure 3B. Magenta dashed line marks Tzeit. (C) Evolution of first Kuramoto order parameter (R) in magenta and circular standard deviation (circSD) in green over time, showing change in coherence of multiple samples during the experiment. An R equal to 1.0 means that samples are in-phase. circSD is equal to -2lnR. (D) Polar plots at different timepoints showing phase of each sample and their first Kuramoto order parameter, represented as a magenta dot along the circumference of a circle and a magenta line segment at the circle’s center, respectively. A longer line segment corresponds to a higher first Kuramoto order parameter, and thus to more coherent samples. The direction of the line denotes the vectorial average of the sample phases. Time is indicated as minutes elapsed from the start of the experiment.

Higher-order entrainment of the segmentation clock

In theory, a non-linear oscillator should be amenable not only to 1:1 entrainment, but also to higher-order n:m entrainment, in which n cycles of the endogenous oscillation lock to m cycles of the zeitgeber (Pikovsky et al., 2001; Mitarai et al., 2013). In practice, demonstration of higher-order entrainment is challenging due to the narrow permissive parameter region, that is, Arnold tongues are progressively narrower away from the 1:1 regime. Strikingly, we found experimental evidence for higher-order 2:1 entrainment (Figure 5). Samples entrained with either 300 or 350 min pulses showed evidence of 2:1 entrainment (Figure 5A–C), that is, the segmentation clock oscillated twice per each zeitgeber pulse and hence the segmentation clock rhythm adjusted to a period very close to 150 min for Tzeit=300 min and approaching 175 min for Tzeit=350 min (Figure 5C). We also found evidence for phase-locking during 2:1 entrainment, that is, a narrowing of the phase distribution among the individual samples (Figure 5D).

The segmentation clock can be entrained to a higher order.

(A) Left: Detrended (via sinc filter detrending) timeseries of the segmentation clock in 2D-assays subjected to 300 min periodic pulses of 2 µM DAPT (or DMSO for controls). Periodic pulses are indicated as magenta bars and the timeseries of each sample (for CTRL: n=5 and N=1, for DAPT: n=5 and N=1) is marked with a dashed line. The gray shaded area denotes the interquartile range. Hypothetical pulses at half the zeitgeber period are indicated as light pink bars. Right: Period evolution during entrainment, obtained from wavelet analysis. The period evolution for each sample and the median of the periods are represented here as a dashed line and a solid line, respectively. The gray shaded area corresponds to the interquartile range. Magenta dashed line marks 0.5Tzeit. (B) Left: Detrended (via sinc filter detrending) timeseries of the segmentation clock in 2D-assays subjected to 350 min periodic pulses of either 5 µM DAPT or 10 µM DAPT. Periodic pulses are indicated as magenta bars and the timeseries of each sample (for 5 µM DAPT: n=8 and N=2, for 10 µM DAPT: n=10 and N=2) is marked with a dashed line. The gray shaded area denotes the interquartile range. Hypothetical pulses at half the zeitgeber period are indicated as light pink bars. Right: Period evolution during entrainment, obtained from wavelet analysis. The period evolution for each sample and the median of the periods are represented here as a dashed line and a solid line, respectively. The gray shaded area corresponds to the interquartile range. Magenta dashed line marks 0.5Tzeit. (C) Mean period either from 1000 to 1150 min (for Tzeit=300 min) or from 800 to 950 min (for Tzeit=350 min) after start of the experiment of samples subjected to periodic pulses of DAPT (or DMSO for controls). Each sample is represented as a dot, while the median of all samples is denoted as a solid vertical line. For 300 min 2 µM DAPT: (CTRL: n=5 and N=1) and (DAPT: n=5 and N=1) with p-value = 0.191, for 350 min 5 µM DAPT: (CTRL: n=7 and N=2) and (DAPT: n=8 and N=2) with p-value = 0.049, for 350 min 10 µM DAPT: (CTRL: n=10 and N=2) and (DAPT: n=10 and N=2) with p-value = 0.016. The period and concentration of the DAPT pulses are specified. Data were visualized using PlotsOfData (Postma and Goedhart, 2019). To calculate p-values, two-tailed test for absolute difference between medians was done via a randomization method using PlotsOfDifferences (Goedhart, 2019). (D) Polar plots at different timepoints showing the phase of each sample in (A) and the first Kuramoto order parameter, represented as a magenta dot along the circumference of a circle and a magenta line segment at the circle’s center, respectively. A longer line segment corresponds to a higher first Kuramoto order parameter, and thus to more coherent samples. The direction of the line denotes the vectorial average of the sample phases. Time is indicated as minutes elapsed from the start of the experiment.

The segmentation clock establishes a stable entrainment phase relative to the zeitgeber

We next analyzed the phase-locking behavior between the segmentation clock and the zeitgeber. According to dynamical systems theory, phase-locking, that is, the entrainment phase (ϕent), can be characterized as an attractor, that is, it is a stable fixed point (Cross and Siggia, 2005; Izhikevich, 2006). To quantify the entrainment phase, we plotted the data as stroboscopic maps (Figure 6A), which take a snapshot of the segmentation clock phase at regular intervals based on zeitgeber pulses (Isomura et al., 2017; Cross and Siggia, 2005; Izhikevich, 2006). Stroboscopic maps enable determination of ϕent as a stable fixed point that lies on the diagonal, where there is phase-locking. Plotting the stroboscopic maps shows that in samples entrained to periodic pulses of DAPT the segmentation clock phase dynamics converge to a region on the diagonal, which marks ϕent (Figure 6B). Such convergence toward a stable entrainment phase is further highlighted by looking at consecutive pulses and the corresponding phase trajectories of individual samples. As exemplified in Figure 6C and as previously theoretically predicted (Granada and Herzel, 2009; Granada et al., 2009), at the same zeitgeber strength and zeitgeber period, faster (or slower) convergence toward this fixed point (i.e. entrainment) was achieved when the initial phase of the endogenous oscillation (ϕinit) was closer (or farther) to ϕent.

The segmentation clock establishes a stable phase relationship with the zeitgeber.

(A) Schematic of how to generate a stroboscopic map, where the phase of the segmentation clock just before a DAPT pulse (old phase, ϕn) is iteratively plotted against its phase just before the next pulse (new phase, ϕn+1).The position of each point in a stroboscopic map thus denotes a stepwise change in phase of the segmentation clock as it undergoes entrainment to the zeitgeber. Upon entrainment and phase-locking, the new phase is equal to the old phase, thus marking a point that lies on the diagonal of the stroboscopic map. This point on the diagonal is the entrainment phase (ϕent). Illustration by Stefano Vianello. (B) Stroboscopic map of samples subjected to 170 min periodic pulses of 2 µM DAPT (or DMSO for controls). Colors mark progression in time, from purple (early) to yellow (late). Note that while in control samples, points remain above the diagonal (reflecting that endogenous oscillations run faster than Tzeit=170 min as shown in Figure 3B–C), in entrained samples, the measurements converge toward a point on the diagonal (i.e. the entrainment phase, ϕent), showing phase-locking. This localized region marks the entrainment phase (ϕent). This is highlighted with a magenta star, which corresponds to the centroid (xc,yc). The centroid was calculated from the vectorial average of the phases of the samples at the end of the experiment, where xc = vectorial average of old phase, yc = vectorial average of new phase. The spread of the points in the region is reported in terms of the circular standard deviation (-2lnR, where R is the first Kuramoto order parameter). (C) Stroboscopic maps of the segmentation clock entrained to 170 min periodic pulses of 2 µM DAPT (n=34 and N=8), for all samples (ALL) and for three individual samples (SAMPLE A, SAMPLE B, SAMPLE C). The numbers and colors (from purple to yellow) denote progression in time. The initial phase (old phase at point 2) and the entrainment phase (new phase at point 5) of each sample are specified.

The phase of entrainment varies according to zeitgeber period

One fundamental property of entrained oscillatory systems is that its phase of entrainment varies as a function of the detuning, that is, the period mismatch between the endogenous oscillator and the entrainment period (Pittendrigh and Bruce, 2015; Aschoff and Pohl, 1978; Rémi et al., 2010; Bordyugov et al., 2015; Granada et al., 2013). In addition, theoretical studies have supported a ‘180 degree’ rule (Granada et al., 2013), stating that the phase of entrainment varies by half a cycle (180° or π) within the range of permissible zeitgeber entrainment periods, Tzeit (Aschoff and Pohl, 1978; Granada et al., 2013).

To test these theoretical predictions, we quantified ϕent across a wide range of detuning, that is, from 120 to 180 min zeitgeber period. We found that ϕent, based on the centroid localization close to the diagonal in the stroboscopic maps, gradually changed its position as zeitgeber period, that is, detuning, was varied (Figure 7A, Figure 7—figure supplement 1). From 120 to 180 min entrainment periods, ϕent systematically shifted from π/2 to 3π/2 (Figure 7A–B), spanning a range of almost π (Figure 7B–C). Hence, as predicted by theory, our results show that the segmentation clock phase of entrainment varies as a function of the detuning relative to zeitgeber pulses.

Figure 7 with 3 supplements see all
The entrainment phase varies according to zeitgeber period within a range of π.

(A) Stroboscopic maps for different periods of DAPT pulses (i.e. zeitgeber periods) placed next to each other. In these maps, only samples that were phase-locked by the end of the experiment are considered (for 120 min: n=13/14 and N=3/3, for 130 min: n=38/39 and N=10/10, for 140 min: n=10/15 and N=3/3, for 150 min: n=16/17 and N=4/4, for 160 min: n=11/15 and N=3/3, for 170 min: n=28/34 and N=8/8, for 180 min: n=4/6 and N=1/1). Here, a sample was considered phase-locked if the difference between its phase at the time of the final drug pulse considered and its phase one drug pulse before is less than π/8. The localized region close to the diagonal in each map marks the entrainment phase (ϕent) for that zeitgeber period. This is highlighted with a magenta star, which corresponds to the centroid (xc,yc). The centroid was calculated from the vectorial average of the phases of the samples at the end of the experiment, where xc = vectorial average of old phase, yc = vectorial average of new phase. The spread of the points in the region is reported in terms of the circular standard deviation (-2lnR, where R is the first Kuramoto order parameter). The zeitgeber period is indicated above the maps. Colors mark progression in time, from purple to yellow. Stroboscopic maps of all samples and their respective controls are shown in Figure 7—figure supplement 1. Drug concentration and drug pulse duration were kept constant between experiments at 2 µM and 30 min/cycle, respectively. Phase 0 is defined as the peak of the oscillation. (B) Entrainment phase (ϕent) at different zeitgeber periods, each calculated from the vectorial average of the phases of phase-locked samples at the time corresponding to last considered DAPT pulse. The spread of ϕent between samples is reported in terms of the circular standard deviation (-2lnR, where R is the first Kuramoto order parameter). (C) Detrended (via sinc filter detrending) timeseries of the segmentation clock in 2D-assays entrained to either 130 or 170 min periodic pulses of 2 µM DAPT, zoomed in from 500 to 950 min. Periodic pulses are indicated as magenta bars and the timeseries of each sample (for 130 min: n=39 and N=10, for 170 min: n=34 and N=8) is marked with a dashed line. The median of the oscillations is represented here as a solid line, while the gray shaded area denotes the interquartile range. The detrended timeseries for the 130 and 170 min conditions are the same as the detrended timeseries for the 2 µM condition in Figure 4—figure supplement 1A and Figure 4A, respectively. The full detrended timeseries for the 170 min condition can be seen in Figure 3A.

It was previously documented that the period of tissue segmentation matches the period of oscillations localized in the posterior PSM, that is, center of of the 2D-assays (Tsiairis and Aulehla, 2016; Soza-Ried et al., 2014). We were then curious whether or not the dynamics from our coarse-grained analysis match that for oscillations in the center of the 2D-assays, and if the behavior of oscillations in this localized region also exhibits dependence to zeitgeber period. To study this, we compared the timeseries acquired from (a) global ROI spanning the entire field of view and from (b) a circular ROI restricted at the center of the 2D-assays (‘center ROI’, radius: 25 pixels, 1 pixel = 1.38 µm). Upon comparison, we noted that the oscillations at the center of 2D-assays adjusted its rhythm to the zeitgeber with strong correspondence to the coarse-grained oscillator (Figure 7—figure supplement 2), regardless of whether the segmentation clock was sped up (Figure 7—figure supplement 2A–B) or slowed down (Figure 7—figure supplement 2C–D). Moreover, the phase-locking between the zeitgeber and either the segmentation clock or the oscillations at the center of the samples are similar (Figure 7—figure supplement 2B1–B2, D1–D2, Figure 7—figure supplement 3). Concomitantly, for oscillations in the center of 2D-assays (i.e. corresponding to the posterior PSM), the entrainment phase at Tzeit=130 min is different to the entrainment phase at Tzeit=170 min (Figure 7—figure supplement 3B), here too exemplifying the effect of zeitgeber period on the phase of entrainment.

PRCs derivation and period change

Building on our finding that the phase of entrainment varies as a function of detuning (Figure 7B), our goal was to extract the quantitative information embedded in this dynamic behavior. This is possible since the phase of entrainment dynamics reflect, at a quantitative level, the fundamental properties of a dynamical system that responds to external perturbations. This behavior, in turn, can usually be captured with a single function, the PRC (Izhikevich, 2006; Granada et al., 2013). The PRC describes the change of phase induced by a perturbation, and a priori depends on both the nature of the perturbation received and the phase of the cycle. Because of the direct dependence of the phase of entrainment dynamics on the PRC, our goal was to use the experimental data to gain insight into the segmentation clock PRC.

We follow Izhikevich, 2006, to define a generalized PRC not tied to the shape or duration of the perturbation. Without loss of generality, we then model perturbations in form of pulses of amplitude ϵ (Pikovsky et al., 2001) and write:

(1) PRC(ϕpulse,ϵ)=ϕafterpulse(ϕpulse,ϵ)ϕpulse

where ϕpulse is the phase of the segmentation clock on the cycle at the moment of the pulse, and ϕafter pulse the phase after the pulse (which can be defined even if the system transiently moves outside of the limit cycle, via isochrons, see Winfree, 1967). Amplitude ϵ is a dimensionless parameter, quantifying zeitgeber strength in the mathematical formalism we develop below (see Appendix 1, section Arnold tongue computations to see how we relate ϵ to DAPT concentrations). If, following a perturbation, the oscillator relaxes quickly toward the limit cycle, one can use the PRC to compute response to periodic pulses with period Tzeit. The sequence of phases at each pulse is then given by the stroboscopic map, introduced in Figure 6A:

(2) ϕn+1=(ϕn+PRC(ϕn,ϵ)+2πTzeitTosc)mod2π

and 1:1 entrainment occurs when this stroboscopic map converges toward a single fixed point ϕent for a given Tzeit. When Tzeit is varied, different phases of entrainment ϕent(Tzeit) are observed, here plotted in Figure 7B.

To derive the PRC directly from the experimental stroboscopic maps, we invert Equation 2 into:

(3) PRC(ϕn,ϵ)=(ϕn+1-ϕn-2πTzeitTosc)mod2π

By estimating ϕn,ϕn+1 for a given Tzeit, one can estimate the PRC. The advantage of this approach is that it allows to estimate the PRC at any observed phase at the pulse ϕn , even far from the entrainment phase ϕent. Notice, however, that phases far from ϕent will be only sampled over the first few pulses (so with possibly much noise) while conversely, ϕent will be quickly oversampled (and as such better defined statistically).

Figure 8A1 shows the PRC computed from the data as well as Fourier series fits for different entrainment periods. PRCs computed for different entrainment periods have a similar shape, with similar locations for maxima and minima. Strikingly, those PRCs are not sinusoidal but essentially 0 or strongly negative, an unusual situation from a dynamical systems theory standpoint associated with special classes of oscillators (see more details below and in Discussion). However, contrary to theoretical predictions, the inferred PRCs at different entrainment periods do not overlap and appear shifted vertically as Tzeit is changed.

Figure 8 with 4 supplements see all
Modeling the segmentation clock entrainment response.

(A) Phase response curve (PRC) from the data for different zeitgeber periods. (1) PRCs calculated at different Tzeit (points) and Fourier series fitted to them (lines). (2) Original PRCs are shifted vertically to collapse the data points on one curve. (B) Oscillator model, optimized by fitting the vertically shifted PRC. The limit cycle is an ellipse (blue) with eccentricity λ=0.5, the region with speeding up s*=5.6 is shaded. (C) Optimized model PRC (line) overlaid to the vertically shifted data points. (D) Modeled intrinsic period Tosc as a function of entrainment period Tzeit. Points were inferred from the PRCs data in A1 by matching the detuning and entrainment phase ϕent, the function was interpolated with cubic splines. (E) Entrainment phase ϕent from the model with changing intrinsic period (line) agrees with the experimental data (points and error bars). Experimental data are the same as those in Figure 7B. (F) 1:1 Arnold tongue and isophases calculated with the model. Stars correspond to experimental data (Tzeit,ϵ) with observed entrainment, in agreement with ϕent data for different DAPT concentrations (Figure 8—figure supplement 4A). Black stars represent the experiments used for optimizing the model. Crosses correspond to experimental conditions with no entrainment. (G) Fit of the model (solid lines) overlaid with experimental phase differences (dots) at different periods Tzeit.

Such vertical shifts in PRCs as a function of entrainment period have been previously observed in so-called ‘overdrive suppression’ for cardiac cell oscillators (Kunysz et al., 1995). Shifts occur when the entrainment signal impacts the biochemical control of the system, leading to a change of intrinsic period (here Tosc). Here, such change of intrinsic period can not be directly measured experimentally, as the system is entrained to another period, however, we found several lines of experimental and theoretical evidence, in addition to the vertical PRC shift, in full agreement with such an effect.

First, the slope of the ϕent(Tzeit) curve in Figure 7B is unusual: PRC theory predicts that for high and low detuning this curve should have vertical slopes. Here, we observed plateaus in the data, which can be simply explained if the intrinsic period changes (see mathematical explanation in Appendix 1). In addition, we performed entrainment release experiments, where the segmentation clock is first entrained, then released (Figure 7—figure supplement 1): we observed a slow recovery over several cycles to a period matching control samples, compatible with a transient change of the intrinsic period. This effect was confirmed by a more precise study of the phase return map after release (Figure 8—figure supplement 3A). Combined, the inference of the PRC based on entrainment quantification data thus reveals two properties: a highly asymmetrical, mainly negative PRC and second, a change of the intrinsic period during entrainment.

ERICA model and Arnold tongue

To gain understanding of these findings from a theory perspective, we next build a minimal model of our data.

First, we take the PRC computed at Tzeit=140 min to be the ‘reference’ curve, since this Tzeit coincides with the natural period of the process. We can then estimate for each entrainment period the period shift most consistent with the data, to obtain a single, period-independent PRC in Figure 8A2 (see Appendix 1, section Theoretical methods). From there, we build upon the simplest non-linear phase-oscillator, that is, the classical Radial Isochron Cycle (RIC). To modulate its sine-form PRC, which is incompatible with our data, we perturb it into an Elliptic Radial Isochron Cycle with Acceleration or ERICA.

ERICA is designed to have radial isochrons, meaning that the phase (and the PRC) can be analytically computed in the entire plane. ERICA then allows for the speeding up of the cycle for some angle range (parameter s*) or for changes of the limit cycle shape into an ellipse of increasing eccentricity (parameter λ). Both high values of λ and s* generate high negative asymmetry in the PRC (see Appendix 1, section Theoretical methods; the full justification of the ERICA construction, the study of its multiple properties and its connection to biological oscillators will be described elsewhere). We then used Monte Carlo (MC) optimization to find parameters best fitting the experimental PRC. The results from this MC optimization put the oscillator far from the standard RIC oscillator (see cycle and corresponding flow in Figure 8B), consistent with strongly negative PRCs, with a moderate value of perturbation ϵ=0.4, high value of λ=0.5 (indicating a strong elliptical shape), and high value of s*=5.6 over half a cycle (indicative of a relaxation type dynamics in x with a slow decay phase followed by a fast reset). Figure 8C compares the PRC of this optimized model with the multiple data points, showing excellent agreement. In addition, we combined the ERICA framework with a simple fit for intrinsic period change (Figure 8D) to account for the experimental phase transition curves, as illustrated in Figure 8G (see also Figure 8—figure supplement 2).

With the ERICA model at hand, we derived numerically the Arnold tongues of the system and the phase/detuning curves for all entrainment parameters (Figure 8E–F). We notice that the Arnold tongues are heavily skewed toward the right, meaning that the system can be more easily slowed down than sped up, consistent with the negative PRC shape. Remarkably, while we build the model using only one entrainment drug concentration, we can also largely explain data obtained at other concentrations (Figure 8F). In particular, the Arnold tongues/our model predict a specific change of the entrainment phase as the entrainment strength is varied (Figure 8—figure supplement 3G). The comparison to experimental data (Figure 8—figure supplement 4) shows that this prediction is, qualitatively, verified. More generally, having the PRC and a minimal, coarse-grained model that captures the essential dynamical features of the segmentation clock during entrainment, including the change in intrinsic period, enables predictable control over the pace and the rhythm of the segmentation clock using the entrainment strategy.

Discussion

In this work, we used a coarse-graining, entrainment approach to gain new insights into the properties of the segmentation clock from dynamical systems’ perspective. We mode-lock the segmentation clock to various entrainment periods and use the information about the dynamic phase-locking behavior to derive the somite segmentation clock PRC from the experimental data.

A coarse-graining approach captures essential dynamical features using a simple ‘one oscillator’ phase description

Given the complexity underlying the somite segmentation clock, comprised of several, interconnected, signaling pathways with countless molecular interactions, it was a priori unclear whether a simple ‘one oscillator’ phase description and a perturbation approach would capture the essential dynamical characteristics at the systems level. Another potential difficulty arises from the fact that cellular oscillators define a phase gradient, controlling the segmentation pattern (Lauschke et al., 2013). For all those reasons, one could imagine that no single phase could describe the system’s behavior, a general concern for systems of interacting oscillators already mentioned by Winfree, 1967.

Remarkably, one key finding of this work is that indeed, using a systems-level single oscillator phase description, we observe a consistent entrainment behavior, that is, period-locking with convergence toward a well-defined entrainment phase, as predicted from the oscillator phase response theory. It is important to point out that theoretical predictions are non-trivial and quantitative, such as the higher order 2:1 entrainment (Figure 5) and the dependence of the entrainment phase on detuning (Figure 7A–B). Given these experimental findings we conclude that the coarse-graining approach and the description of the entire system using a single oscillator phase is justified and enables to extract the essential dynamical properties, which are captured by a mathematical model including the system’s PRC.

Asymmetrical segmentation clock PRC compatible with SNIC bifurcation

Insight into the system’s PRC allows to infer about the nature and characteristics of the oscillatory system, independent of its specific molecular realization (Izhikevich, 2006). For instance, excitable systems can be functionally categorized with their PRC, that is, in neuroscience, Class I excitable systems have PRCs with constant sign (Izhikevich, 2006). The canonical example of (Class I) oscillators with constant-sign PRCs is the quadratic fire and integrate model (Izhikevich, 2006), which can be used to model neuronal spiking.

In contrast, Class II excitable systems exhibit sinusoidal PRCs, examples include many biological systems such as the circadian clock (Granada et al., 2009).

This difference in PRC reflects fundamental properties of the underlying dynamics, including the associated bifurcations: Class I systems are associated with infinite period bifurcations, a typical example being the saddle-node on invariant cycle (SNIC) bifurcation. For such bifurcations, the period can be tuned and can become arbitrarily long depending on the proximity to the bifurcation, and oscillations present relaxation-type dynamics (i.e. combination of slow dynamics with fast resetting). The simplest example, the quadratic fire and integrate model, is a one variable x oscillator, where x goes from - to + in a finite time before resetting to - (Izhikevich, 2006). The PRC for this oscillator is ‘naturally’ of constant sign because there is only one variable, which can only advance (or delay) in response to a constant signal. On the other hand, Class II systems are associated with Hopf bifurcations (Izhikevich, 2006), have a fixed period, and more sinusoidal dynamics close to the bifurcation.

Remarkably, the segmentation clock PRC that we computed here shares features of systems close to an infinite period bifurcation, that is, a highly asymmetrical, constant-sign PRC with an extended flat region. Consistent with this, the analytical ERICA model that we developed to capture the experimental PRCs in a wide range of conditions exhibits relaxation-type dynamics (Figure 8B, slow phase in white and fast resetting phase in gray). It is nevertheless important to point out that a given PRC can result from different oscillators, that is, a constant sign does not prove the existence of a SNIC oscillator. For instance, delayed oscillators close to a Hopf bifurcation can show PRCs with an (almost) constant sign (Kotani et al., 2012). Infinite period can also be obtained with saddle homoclinic bifurcations (Izhikevich, 2006). However, such models are mathematically more complex, for instance delays include long-term memory in the system, and homoclinic bifurcations require the existence of other attractors outside of the cycle.

In agreement with SNIC oscillators, the segmentation clock does show a tunable period and a slowing down behavior at multiple levels: it is known that the overall rate of segmentation slows down over developmental time, a feature described in several species. In addition, the oscillations also slow down along the embryonic axis, that is, a period gradient is present within the PSM, reflecting tunability of oscillations at cellular level, as PSM cells progress toward differentiation and eventually stop oscillations. How the tunability at cellular level is linked to the system-level change in segmentation clock rate is not understood (see also below for macro- versus micro-scale comparison). However, our findings of mostly negative PRCs might reflect that the segmentation clock has a natural bias of slowing down the oscillations.

From a network standpoint, biological networks with tunable periods have been associated with specific structural features, combining negative and positive feedback loops (Tsai et al., 2008). Positive feedback in particular puts oscillators close to a multistable regime, leading to relaxation oscillations (Ferrell et al., 2011). SNICs appear naturally when the regulatory logic of a system moves from a negative feedback oscillator to multistability (Jutras-Dubé et al., 2020), reflecting network modularity. Such modularity and robustness with changes in regulatory logic is a hallmark of developmental plasticity (West-Eberhard, 2003) and hence one would anticipate to find SNIC oscillators to be abundant in developing systems. First examples are indeed emerging, such as during Caenorhabditis elegans development (Meeuse et al., 2020).

Findings not predicted by the theory of PRC

We also made several unexpected findings, not predicted by the theory of PRC.

First, we found evidence that during entrainment, not only does the segmentation clock adjust its observed period to the zeitgeber pulses, but also changes its intrinsic period in the direction of the zeitgeber rhythm. Hence, during entrainment that slows down the clock (i.e. 170 min), we find evidence that the intrinsic period lengthens (not just the observed rhythm), while during entrainment that speeds up the clock (i.e. 120 min), we find evidence that the intrinsic period decreased. Again, such a change of intrinsic period is not predicted in entrained systems. The simplest explanation could be that the drug pulses change the period of the oscillator by changing some biochemical parameter in the system, similar to overdrive suppression in cardiac cells (Kunysz et al., 1995). However, as stated above, we do not find evidence for a consistent slowing down or speeding up effect: the intrinsic period is decreased for short entrainment cycles and increased for longer entrainment cycles. This rather suggests the existence of feedbacks of the clock on itself, leading to higher-order adaptations beyond the rapid Notch phase response. One can only speculate on the mechanisms underlying such adaptation, but this is compatible with the idea that two interacting oscillators control the intrinsic period. Here, it is possible that entraining Notch with a zeitgeber modifies the Wnt oscillation period, which in turn feeds back on Notch on a longer time scale. So the internal period change that we see might in fact come from the induced change of Wnt period. We have shown previously that Wnt and Notch oscillators are coupled but are not phase-locked, with functional impact on tissue patterning (Sonnen et al., 2018). Alternatively, a similar role could also be played by the long inter- or intracellular delays in the system, postulated in multiple theoretical works (Lewis, 2003; Morelli et al., 2009; Ares et al., 2012). Such delays could effectively couple multiple cycles, changing clock parameters (such as the period) beyond the instantaneous phase response. More experimental and theoretical work is needed to explore these ideas.

Second, a striking outcome we obtained was that even after entrainment, a spatial period gradient and phase waves emerged over time (Appendix 1—figure 1). Put differently, while the overall system is entrained, as evidenced by a control of the timing of morphological segmentation and of segmentation clock rhythm, the underlying cellular oscillators show a divergent, yet, coherent response, that is, a spatial period gradient emerges within the tissue. A recent entrainment model of coupled PSM oscillators (Juul et al., 2018) does not account for the emergence of a period gradient during entrainment. This unexpected behavior hence awaits further experimental and theoretical analysis. In particular, the role of intra- and intercellular coupling needs to be addressed to reveal how the macro-scale behavior relates to the underlying cellular scale oscillations.

Conclusions

Our work demonstrates how, despite all the molecular and functional complexities, coarse-graining and theory can be used to effectively take control of complex biological processes. A molecular mechanism is not needed to exert control, as long as we have a mathematical one, one that captures the essential features of a system at a meaningful coarse-grained level.

We also aim to illustrate the potential of an integrated, theoretical-experimental approach to complex biological systems: while from theoretical viewpoint, the fact that entrainment phase varies as detuning is altered is a mathematical necessity and hence ‘obvious’, this outcome is not at all intuitive, a priori, from an experimental-observational viewpoint. Theory guides experimentation, that is, Figure 1 Theory (Phillips, 2015), leading to new insight and understanding of a complex biological phenomenon.

Future studies are needed to gain further insight into the response to other perturbation regimes. In addition, it will be essential to perform the analysis both at the macro-scale, as done in this study, but also extending it to the cellular scale, that is, taking spatial differences, such as the period gradient, fully into account. We consider that such a dynamical systems’ theoretical-experimental approach is a promising and, as we show here, feasible way forward with the goal in categorizing the segmentation clock in its universality class.

Materials and methods

Please refer to the Appendix 1 for detailed materials and methods.

Appendix 1

Appendix 1—key resources table
Reagent type (species) or resourceDesignationSource or referenceIdentifiersAdditional information
Genetic reagent (Mus musculus)LuVeLuAulehla et al.DOI:10.1038/ncb1679
Genetic reagent (Mus musculus)Mesp2-GFPMorimoto et al.DOI:10.1016/j.ydbio.2006.08.043
Genetic reagent (Mus musculus)Axin2-linker-
Achilles
This worksee Materials and methods >
Mouse lines
Chemical compound, drugBovine serum albumin (BSA)Equitech-BioBAC62
Chemical compound, drugDMEM/F12Cell Culture TechnologiesCustomizedNo glucose, glutamate,
L-glutamine, and phenol red
Chemical compound, drug1 M HEPESGibco15630-106
Chemical compound, drug10,000 U/mL PenStrepGibco15140-122
Chemical compound, drug45% glucoseSigmaG8769
Chemical compound, drug200 mM L-glutamineGibco25030-081
Chemical compound, drugFibronectinSigma-AldrichF1141
Chemical compound, drugCascade BlueInvitrogenC-3239
Chemical compound, drugDMSOSigma-AldrichD8418
Chemical compound, drugDAPTSigma-AldrichD5942-5MGStock concentration
in DMSO: 10 mM
Commercial assay, kitSYLGARD
184 silicone kit
Dow16739219:1 (w/w) elastomer
base: curing agent
Software, algorithmMicroscopy Pipeline
Constructor (MyPiC)
Politi et al., 2018DOI:10.1038/nprot.2018.040https://github.com/manerotoni/mypic/
Software, algorithmFijiSchindelin et al.DOI:10.1038/nmeth.2019
Software, algorithmMatplotlibHunterDOI:10.1109/MCSE.2007.55
Software, algorithmNumPyvan der Walt et al.DOI:10.1109/MCSE.2011.37
Software, algorithmpandasMcKinneyDOI:10.25080/Majora-92bf1922-00a
Software, algorithmscikit-imagevan der Walt et al.DOI:10.7717/peerj.45310.7717/peerj.453
Software, algorithmSciPyVirtanen et al.DOI:10.1038/s41592-019-0686-2
Software, algorithmseabornWaskom et al.DOI:10.5281/zenodo.12710
Software, algorithmPlotTwistGoedhartDOI:10.1371/journal.pbio.3000581https://huygens.science.uva.nl/PlotTwist/
Software, algorithmPlotsOfDataPostma and GoedhartDOI:10.1371/journal.pbio.3000202https://huygens.science.uva.nl/PlotsOfData/
Software, algorithmPlotsOfDifferencesGoedhartDOI:10.1101/578575https://huygens.science.uva.nl/PlotsOfDifferences/
Software, algorithmpyBOATMönke, 2020aDOI:10.1101/2020.04.29.067744https://github.com/tensionhead/pyBOAT/
Software, algorithmSpatial pyBOAT (SpyBOAT)Mönke, 2020bhttps://github.com/tensionhead/SpyBOAT/
Software, algorithmEntrainment
Analysis Jupyter notebook
This work, Sanchez, 2022EntrainmentAnalysis in https://github.com/PGLSanchez/EMBL_OscillationsAnalysis/
OtherCover glassMarienfeld0107999 09870 mm × 70 mm, 1.5H, see
Materials and methods >
Preparations for microfluidics-
based entrainment of 2D-assays
OtherPTFE tubingAPTAWG24TInner diameter: 0.6 mm, see
Materials and methods >
Preparations for microfluidics-
based entrainment of 2D-assays
OtherSyringe needleBD Microlance 330090022G×1.25" – Nr. 12, see
Materials and methods >
Preparations for microfluidics-
based entrainment of 2D-assays
Other3 mL syringeBD Luer-Lok309658Diameter: 8.66 mm, see
Materials and methods >
Preparations for microfluidics-
based entrainment of 2D-assays
Other10 mL syringeBD Luer-Lok300912Diameter: 14.5 mm, see
Materials and methods >
Preparations for microfluidics-
based entrainment of 2D-assays
OtherProgrammable pumpWorld Precision InstrumentsAL-400see Materials and methods >
Setting up automated pumping
OtherLSM 780 laser-scanning microscopeCarl Zeiss Microscopysee Materials and methods >
Confocal microscopy
Appendix 1—figure 1 with 5 supplements see all
A spatial period gradient emerges in the tissue even upon entrainment of the segmentation clock.

Left: Schematic of pipeline to generate period and phase movies using pixel-by-pixel wavelet analysis (A). Illustration by Stefano Vianello. Right: Snapshot of period wavelet movie of control subjected to 170 min periodic pulses of DMSO (B) or entrained sample subjected to 170 min periodic pulses of 2 µM DAPT (C), taken at 840 min after the start of experiment. Period is either shown using heatmap (B1,C1) or as a surface plot (B2,C2). Bottom: Snapshots of period wavelet movie of samples in (B) and (C) at different time points. Time is indicated as minutes elapsed from the start of the experiment. Sample is rotated so that the dorsal side is up. Timelapse of period wavelet movies and corresponding surface plots are available at https://youtu.be/3Y53OhXacKI.

Materials and methods

Mouse lines

Request a detailed protocol

For most of the experiments in this study, we used a transgenic mouse line expressing a dynamic Notch signaling reporter driven from the Lfng promoter, more commonly known as LuVeLu. The generation of LuVeLu was previously described (Aulehla et al., 2008). Briefly, the expression of Venus, an improved version of YFP (Nagai et al., 2002), was driven from a 2 kb fragment of the Lfng promoter (Morales et al., 2002; Cole et al., 2002). The locus was flanked by Lfng 3’-UTR and a modified PEST domain (Rogers et al., 1986) to destabilize the reporter mRNA and protein, respectively. For experiments visualizing Mesp2 during entrainment, we used the Mesp2-GFP mouse line generated by Morimoto et al., 2006. Axin2-linker-Achilles, a mouse line expressing Axin2 fused with a GSAGS linker to Achilles, a fast-maturing variant of YFP (Yoshioka-Kobayashi et al., 2020), was generated in-house. To generate the knock-in alleles, we targeted the stop codon of endogenous Axin2 locus with vector containing the reporter sequence coding for Achilles and a selection cassette. This targeting vector was constructed as follows: linker-Achilles-loxP-PGK Neo-loxP. The selection cassette was flanked by loxP- sites for eventual Cre-mediated excision. Axin2-linker-Achilles knock-in reporter line was generated by standard gene targeting techniques using R1 embryonic stem cells. Briefly, chimeric mice were obtained by C57BL/6 blastocyst injection and then outbred to establish the line through germline transmission. The Achilles/pRSETB plasmid was a gift from the lab of Atsushi Miyawaki at RIKEN Center for Brain Science (RIKEN-CBS) in Japan. Mice were kept in an outbred background and were housed in the EMBL Laboratory Animal Resources (LAR). All animal experiments were conducted under veterinarian supervision and after project approval by European Molecular Biology Laboratory, following the guidelines of the European Commission, Directive 2010/63/EU and AVMA Guidelines 2007.

Media preparation

Request a detailed protocol

On the day of the experiment, dissection medium and culture medium were freshly prepared as indicated in Appendix 1—table 2. Culture medium was filter sterilized using a PVDF filter (pore size: 0.22 µm, Merck). Both dissection medium and culture medium were equilibriated to 37°C for at least 15 min, and were kept in a 37°C incubator under 5% CO2 until use.

Appendix 1—table 1
Recipe for media preparation.

Formulations specified here are for the preparation of approximately 50 mL of medium. Special DMEM/F12* used in these media does not contain glucose, L-glutamine, sodium pyruvate, and phenol red. Culture medium is filter sterilized after preparation.

ComponentDissection mediumCulture medium
BSA (Equitech-Bio, BAC62)0.5 g0.02 g
DMEM/F12* (Cell Culture Technologies)50 mL50 mL
1 M HEPES (Gibco, 15630-106)1 mL
10,000 U/mL PenStrep (Gibco, 15140-122)500 µL
45% Glucose (Sigma, G8769)44.4 µL44.4 µL
200 mM L-Glutamine (Gibco, 25030-081)500 µL500 µL

Mouse dissection and embryo recovery

Request a detailed protocol

Female mice were sacrificed on 10.5 dpc (days post coitum) via cervical dislocation. The skin on their ventral side (belly area) was wiped with 70% ethanol, and an incision was made using a clean pair of surgical scissors. The uterine horns were harvested and were washed once with dissection medium. In dissection medium, under a stereo microscope (Leica M80), the deciduae were cut open using clean forceps and the embryos were recovered. The embryos were again washed with fresh dissection medium and their tails were clipped using forceps. Clipped tails were then transferred to a new dish with fresh dissection medium. The volume of dissection medium in the dish was kept to minimum to lessen autofluorescence, which could interfere with subsequent screening. The tails were then screened for presence of the reporter-of-interest (e.g. LuVeLu) using a stereo fluorescence microscope. The tailbud was cut from the rest of the tail and was immediately transferred to pre-equilibriated culture medium that was supplemented with HEPES (170 µL of 1 M HEPES in 10 mL culture medium). These isolated embryonic tissues were then directly loaded in the microfluidics device.

Preparations for microfluidics-based entrainment of 2D-assays

We used in this study a microfluidics-based experimental entrainment platform with a general protocol elaborated in a recent publication (van Oostrom et al., 2021). The PDMS microfluidics device was made using standard soft lithography (El Debs et al., 2012) and as previously described (Sonnen et al., 2018). The ratio of Sylgard 184 silicone elastomer base (Dow) to curing agent (Dow) was 9:1 (w/w). PDMS chip was attached to cover glass (70 mm × 70 mm, 1.5H, Marienfeld 0107999 098) via plasma bonding. UV irradiation of PTFE tubing and PDMS device PTFE tubing (inner diameter: 0.6 mm, APT AWG24T) and syringe needles (22 G 1 1/4 – Nr. 12), as summarized in Appendix 1—table 3, were prepared a day prior to actual microfluidics-based entrainment experiment. In addition to one PDMS microfluidics device (Figure 1—figure supplement 1), each experiment required four 3 m PTFE tubing (each with syringe needle inserted in one end) for the drug/medium inlets, two 1 m PTFE tubing (each with syringe needle inserted in one end) for the outlets, and twenty-four 1 cm plugs made from cut PDMS-filled PTFE tubing for the sample inlets and unused drug/medium inlets. These were all sterilized under UV for at least 20 min.

Appendix 1—table 2
List of PTFE tubing needed for a microfluidics-based entrainment experiment.

For the first two items, a syringe needle is to be inserted inside one end of each tubing. Plugs are made from cut PDMS-filled PTFE tubing, and are used to seal inlets after sample loading. Controls are already taken into account in the specified quantities.

ItemWith needle?QuantityUse
3 m PTFE tubingY4Drug/medium inlet
1 m PTFE tubingY2Outlet
1 cm plugN24Sample inlet+unused drug/medium inlet
Fibronectin-coating of PDMS device and overnight soaking in buffer
Request a detailed protocol

While waiting, 5 mL of PenStrep (Gibco, 15140-122) was added to 495 mL of ×1 PBS (PBS + PenStrep); 5.6 mL of the buffer was set aside to prepare fibronectin solution, while the rest was poured in a glass dish. After UV irradiation, the sterilized PDMS device was immersed in the PBS + PenStrep and bubbles were removed by flushing the channels with buffer. PDMS-filled PTFE plugs were also immersed in PBS + PenStrep. To prepare the fibronectin solution, 280 µL of fibronectin (Sigma-Aldrich, F1141) was added to the set aside 5.6 mL PBS + PenStrep. At least 2.5 mL of fibronectin solution was loaded into a 3 mL syringe (diameter: 8.66 mm, BD Luer-Lok REF 309658). A UV-irradiated needle, which was earlier inserted into a 1-meter PTFE tubing, was attached to the filled syringe. The tubing was then inserted into an outlet in the PDMS device, carefully avoiding introduction of bubbles. Fibronectin was flowed (flow rate: 50 µL/hr) into the PDMS device at room temperature overnight.

Preparation of syringes containing drug/medium
Request a detailed protocol

For a microfluidics-based experiment with periodic pulses of drug, four 10 mL syringes (diameter: 14.5 mm, BD Luer-Lok REF 300912) were filled with either the drug, DMSO control, or culture medium (see Appendix 1—table 2). Components of solution in each of these syringes are specified in Appendix 1—table 4. Drug used in this study was DAPT (Sigma-Aldrich, D5942-5MG). To prepare 10 mM stock of DAPT, 5 mg DAPT (MW = 432.46 g/mol) was dissolved in 1156.2 µL DMSO (Sigma-Aldrich, D8418). This was aliquoted and stored at –20°C until use.

Appendix 1—table 3
List of syringe containing drug/medium for a microfluidics-based entrainment experiment.

Formulations specified here are for any drug with final concentration X µM. For recipe to prepare culture medium, please refer to Appendix 1—table 2.

ComponentSyringe ASyringe BSyringe CSyringe D
DRUGCONTROLMEDIUMMEDIUM
Drug (10 mM in DMSO)X µL
DMSO (Sigma-Aldrich, D8418)X µL
Cascade Blue (Invitrogen, C-3239)2 µL2 µL
Culture medium (see Appendix 1—table 2)Dilute to 10 mLDilute to 10 mL10 mL10 mL
Degassing drug/medium and PDMS device
Request a detailed protocol

After coating and overnight soaking, the PTFE tubing was cut away from the needle and was immediately immersed in the buffer. The dish containing immersed PDMS device, with attached tubing for the two outlets, and plugs were placed inside a vacuum desiccator chamber. The plunger of each syringe containing the drug/medium was pulled to maximum. The syringes were then also placed in the vacuum chamber, almost vertically, with the plunger resting on the desiccator. These were degassed under vacuum for at least 1.5 hr.

Loading embryonic tissues in microfluidics device and mounting for live imaging

Request a detailed protocol

Before mouse dissection and recovery of mouse embryos, syringes containing the degassed drug/medium were each connected to a UV-irradiated 3 m PTFE tubing via attached syringe needle, and were carefully mounted on programmable pumps (World Precision Instruments, AL-400) next to the microscope. Gas in the tubing was displaced with drug/medium by careful pushing of the syringes’ plunger. Pumps were turned on and flow rate was set to 900 μL/hr. The microscope was then equilibriated to 37°C and 5% CO2. After recovery of embryonic tissues, using a pipette (i.e. P200 for 2D-assays and P1000 for intact PSM), each sample was carefully loaded into the microfluidics device, which was already coated with fibronectin, immersed in a buffer of PBS and PenStrep, and degassed. Each sample inlet was plugged with a PDMS-filled PTFE tubing immediately after sample loading. Unused inlets, if any, were also plugged. Flow rate of drug/medium was set to 20 μL/hr and the tubings were carefully inserted into the drug/medium inlets in the PDMS device while it is immersed in buffer. The tubings connected to the syringes with medium were inserted first, and the tubing connected to the syringe with the drug was inserted last. A 15 min timer was started after insertion of the drug tubing. PDMS device was then removed from the buffer and excess liquid on the cover glass was removed carefully with lint-free wipes (Kimberly-Clark). The PDMS device, with more than 1 m of attached tubings, was carefully placed inside a pre-equilibriated microfluidics holder (EMBL Mechanical Workshop, Appendix 1—figure 2) customized to fit the 70 mm × 70 mm cover glass (Marienfeld 0107999 098) and some PTFE tubing. Putting some of the tubing inside the microfluidics holder was necessary to equilibriate the drug/medium to desired environmental conditions before they are perfused into the microfluidics device. The cover glass was secured in place with grease and a U-shaped metal clamp. The microfluidics holder was then carefully mounted on the stage of the microscope, and the end of the outlet tubings were placed in a beaker. Fifteen minutes after insertion of the drug tubing in the PDMS device, each pump was tilted on its side and was equilibriated for another 15 min. Afterward, the pump mounting the syringes with the drug and DMSO control was turned off, and the pump mounting the syringes with the medium was set to flow rate of 60 µL/hr. Samples were equilibrated at these conditions for at least 30 min before start of imaging/entrainment.

Appendix 1—figure 2
Customized box to mount microfluidics device on microscope for simultaneous culture, entrainment, and live imaging.

Computer-aided design (CAD) of metal box customized to hold a microfluidics device bonded to a 70 mm × 70 mm cover glass. The box fits the stage of an LSM 780 laser-scanning microscope (Carl Zeiss Microscopy). Design by Katharina Sonnen, the EMBL Mechanical Design Office, and the EMBL Mechanical Workshop.

Setting up automated pumping

Request a detailed protocol

Entrainment via periodic pulses of drug was performed through alternate perfusion of medium and drug into the microfluidics device. Perfusion of drug/medium was done using syringes mounted on programmable syringe pumps (World Precision Instruments, AL-400). Diameter of syringe (14.5 mm for 10 mL syringe, BD Luer-Lok REF 300912) was accordingly set to match a defined flow rate and a specified volume of solution to be perfused with the duration of perfusion. Standard pumping programs of medium and drug/control are summarized in Appendix 1—table 5, respectively, considering an entrainment experiment with Tzeit = 170 min.

Confocal microscopy

Request a detailed protocol

For most of the experiments here, samples were imaged in an LSM 780 laser-scanning microscope (Carl Zeiss Microscopy) fitted with an incubation chamber (EMBL Mechanical Workshop). A Plan-Apochromat ×20 air objective with a numerical aperture (NA) of 0.8 (Carl Zeiss Microscopy) was used for imaging, and the zoom was set to 0.6. Three z-stacks (spacing: 8 µm) were scanned for each sample every 10 min to acquire timelapse movies. Imaging of multiple samples in multiple locations was done with a motorized stage, controlled using Zen Black software (Carl Zeiss Microscopy), and automated using a VBA macro developed by Antonio Politi (Politi et al., 2018), which is available at https://github.com/manerotoni/mypic. The dimension of the images was 512 pixels × 512 pixels, with a pixel size of 1.38 µm and bit depth of 16 bit. Detection of drug pulses, using Cascade Blue (excited with 405 nm) added to the solution, was also done every 10 min with lower image resolution: 1 z-stack, 32 pixels × 32 pixels (pixel size = 22.14 µm). Imaging and automated pumping of drug/medium through the microfluidics device were started simultaneously.

Appendix 1—table 4
Pumping program of medium for entrainment to 170 min periodic pulses of drug.
PhaseFunctionRateVolumeRemark
PH:01LP:ST
PH:02RATE60 µL/hr140 µLPerfuse medium for 140 min
PH:03LP:ST
PH:04PS:60Pause for 60 s
PH:05LP:30Loop PH:04 30 times (pause for 30 min=60 s × 30)
PH:06LP:30Loop PH:02 to PH:05 30 times (total: 30 periodic pulses)
PH:07STOP
Appendix 1—table 5
Pumping program of drug/control for entrainment to 170 min periodic pulses of drug.
PhaseFunctionRateVolumeRemark
PH:01LP:ST
PH:02LP:ST
PH:03LP:ST
PH:04PS:60Pause for 60 s
PH:05LP:70Loop PH:04 70 times (pause for 70 min=60 s × 70)
PH:06LP:02Loop PH:04 to PH:05 2 times (pause 140 min=70 mins × 2)
PH:07RATE60 µL/hr30 µLPerfuse drug/control for 30 min
PH:08LP:30Loop PH:02 to PH:07 30 times (total: 30 periodic pulses)
PH:09STOP

Data analysis

Details on the coarse-graining strategy and analyses of entrainment dynamics are specified here.

Extracting timeseries from global intensity of 2D-assays for subsequent analyses
Request a detailed protocol

To extract the timeseries corresponding to the segmentation clock (Figure 2, Figure 1—figure supplement 1A), global intensity analysis of timelapse fluorescence imaging was done using Fiji (Schindelin et al., 2012). Z-stacks were first projected based on maximum intensity (in Fiji: Image > Stacks > Z Project). Subsequently, the timeseries was obtained by plotting the z-axis profile (in Fiji: Image > Stacks > Plot Z-axis Profile). Timeseries of replicate samples were compiled in a .txt file for subsequent analyses. To extract timeseries corresponding to signaling oscillations localized at the center of 2D-assays, samples that were centered in the field of view were considered. A 50 pixel × 50 pixel oval ROI was specified at the center (x=256 and y=256) of a re-oriented (and registered, if necessary) 2D-assay timelapse file (in Fiji: Edit > Selection > Specify). Timeseries was extracted after specifying the center ROI (in Fiji: Image > Stacks > Plot Z-axis Profile). Timeseries of replicate samples were compiled in a .txt file for subsequent analyses. Definition of Kuramoto order parameter We follow the definition from Pikovsky et al., 2001: introducing the mean field variable

(4) Z=1Nk=1Neiϕk=ReiΘ

(see also Riedel-Kruse et al., 2007), the first Kuramoto order parameter is defined by

(5) R=|Z|
Monitoring period-locking and phase-locking
Request a detailed protocol

Entrainment was evaluated based on period-locking and phase-locking of the signaling oscillations to the periodic drug pulses. Oscillatory components were extracted from timeseries using a wavelet analysis workflow that was developed by Gregor Mönke, which was recently implemented as a Python-based standalone software (Mönke et al., 2020c, available at https://github.com/tensionhead/pyBOAT; Mönke, 2020a). In this workflow, timeseries was first detrended using a sinc filter and then subjected to continuous wavelet transform. The Morlet wavelet is used as the mother wavelet (Mönke et al., 2020c). Time-resolved frequency analysis was done by cross-correlating the signal to wavelet functions of known frequencies, generating a power spectrum. A high power score was assigned to wavelets that correlated well with the signal relative to white noise. Instantaneous period and phase were extracted upon evaluation of the power spectrum along a ridge tracing wavelet with maximum power for every timepoint. Phase dynamics of signaling oscillations, upon subjecting them to periodic perturbation, were analyzed using stroboscopic maps (Balanov, 2008; Cross and Siggia, 2005; Isomura et al., 2017). Briefly, the phase difference (Δϕ) was defined as

Δϕ=ϕ(t+Tzeit)-ϕ(t)

where t is the time of perturbation and Tzeit is the zeitgeber period (i.e. one cycle after time = t). ϕ(t) and ϕ(t+Tzeit) denote the old_phases and their corresponding new_phases, respectively. The stroboscopic maps were then plotted as new_phases versus old_phases (for scheme, please refer to Figure 6A). The centroid, marking the entrainment phase (ϕent), was determined considering only phase-locked samples (where the difference between a sample’s phase at the time of the final drug pulse considered and its phase one drug pulse before is less than π/8). This was quantified from the average phases of the final (old_phase,new_phase) pairs considered, and the circular standard deviation (circSD) was calculated using the formula:

circSD=-2lnR

where R is the first Kuramoto order parameter. As wavelets only partially overlap the signal at the edges of the timeseries, resulting in deviations from true phase values (Mönke et al., 2020c), the first and last pulse pairs were not considered in the generation of stroboscopic maps. Polar plots were also generated summarizing the instantaneous phase of replicate samples and their first Kuramoto order parameter as shown in Figure 4D, Figure 5D, and Figure 4—figure supplement 1B. The Python code is available as a Jupyter notebook (.ipynb) at https://github.com/PGLSanchez/EMBL_OscillationsAnalysis/tree/master/EntrainmentAnalysis (Sanchez, 2022; copy archived at swh:1:rev:bd4a258570e7b7b5f0c179bc976f0d45988efa1a). This code uses Matplotlib (Hunter, 2007), NumPy (van der Walt et al., 2011; Harris et al., 2020), pandas (McKinney, 2010), scikit-image (van der Walt et al., 2014), SciPy (Virtanen et al., 2020), and seaborn (Waskom et al., 2014).

Detrended timeseries of replicate samples were in some part of the study represented as a heatmap using PlotTwist (Goedhart, 2020), as shown in Figure 4A. Average periods (from 650 to 850 min after start of experiment) were meanwhile plotted using PlotsOfData (Postma and Goedhart, 2019), as shown in Figure 3C–D and Figure 3—figure supplement 1. These apps are available at https://huygens.science.uva.nl/PlotTwist/ and at https://huygens.science.uva.nl/PlotsOfData/, respectively. To calculate p-values, two-tailed test for absolute difference between medians was done via a randomization method using PlotsOfDifferences (Goedhart, 2019). While not recommended, p-values were also calculated similarly for conditions with n<10. PlotsOfDifferences is available at https://huygens.science.uva.nl/PlotsOfDifferences/.

Generating wavelet movies

Period and phase wavelet movies were generated using the Wavelet Processing and Export workflow developed by Gregor Mönke, which runs on the EMBL cluster and is implemented in Galaxy with technical assistance from Jelle Scholtalbers (EMBL Genome Biology Computational Support). This workflow extracts timeseries of every pixel in a timelapse movie and subjects them to sinc filter-based detrending and subsequent continuous wavelet transform (Mönke et al., 2020c). This results in extraction of instantaneous period and phase of each pixel, recovering period and phase wavelet movies corresponding to the input timelapse movie (for scheme, please refer to Appendix 1—figure 1A). The workflow is available at https://github.com/tensionhead/SpyBOAT, (Mönke, 2020b).and can be used via a public Galaxy server at https://usegalaxy.eu/root?tool_id=/spyboat. Settings used to generate wavelet movies in this study were: sigma of 8.0, sample interval of 10 min, period range from 100 to 250 min, and number of periods analyzed of 151.

Examining period gradient in 2D-assays

Re-oriented 2D-assays (dorsal side up) and their corresponding period and phase wavelet movies were used for the analyses. Temporal evolution of the period gradient during the course of entrainment experiments was evaluated from the period wavelet movies of the 2D-assays. As the period wavelet movies also contained wavelet transformations for pixels in the background, a binary mask was first created to differentiate pixels corresponding to signal. To create the mask, re-oriented (and registered, if necessary) timelapse movies of 2D-assays were blurred using a Gaussian blur (in Fiji: Process > Filters > Gaussian Blur [sigma radius: 8, scaled units in µm]). Then, signal was specified by thresholding (in Fiji: Image > Adjust > Threshold [Default method, dark bakground]). After thresholding, pixels corresponding to the signal were assigned a value of 255, while those corresponding to the background were assigned a value of 0. If opposite, the values were inverted (in Fiji: Edit > Invert). Then, the binary mask (signal = 1 and background = 0) was created by dividing all values by 255 (in Fiji: Process > Math > Divide [value: 255]), and was used to mask the period wavelet movie.

Theoretical methods

PRC correction and computation

PRCs at different entrainment periods appear shifted with respect to one another. To correct for this, we fitted the PRCs at entrainment periods Tzeit=(120,130,140,170) min (those for which we have enough experimental data points to build a continuous curve) using three Fourier modes (Figure 8A1) and manually chose the phase shift optimizing their overlap for Figure 8A2. The shift values used were (-0.45,-0.25,0,0.45) for Tzeit=(120,130,140,170) min, respectively. The data points shown in Figure 8C are used for MC optimization (see next section). Our choice was later validated by recomputation of PRCs and comparison with data in Figure 8G. We also show the corresponding PTC and stroboscopic maps in Figure 8—figure supplement 3D-E.

Model

Our model (Appendix 1—figure 3A) is based on simple modifications of the classical Poincaré oscillator, otherwise called RIC model. The main motivation is to find the simplest possible model able to reproduce the general shape of the experimental PRC, with both a flat and a negative region, while being analytical with simple isochrons. The full description and analytical calculations for this model will be described elsewhere, here we outline the basic equations and some general properties of the model. For PRC computations, we will assume that a perturbation is a horizontal shift of magnitude ϵ (Appendix 1—figure 3B).

Appendix 1—figure 3
ERICA model.

(A) Our model consists of an elliptic limit cycle (blue) with increased angular velocity in one sector (shaded region). (B) The effect of such acceleration is to change the spacing of isochrons (radial lines). We compute the phase response curve (PRC) of the model by introducing a perturbation (arrow) at points on the limit cycle and looking at the starting and ending isochrons.

To flatten the PRC for half the cycle, we first modify the limit cycle of the standard RIC model into an ellipse, keeping the origin as one of the foci. This is done through the introduction of a parameter λ, so that the corresponding equation of the ellipse is

(6) rcycle=11+λsinθ

combined with radial equidistant isochrons, θ˙=1, where θ is the polar angle. This is however not enough to define the entire flow in the plane since one needs to specify how r converges on this cycle. We thus impose a radially uniform convergence rate equal to 1 for the radius, which leads to the following differential equation in polar coordinates:

(7) r˙=r˙cycle+r(rcycle-r)

where r˙cycle=-λrx(r+λy)2 (again assuming θ˙=1). This leads to the following equations in cartesian coordinates:

(8) x˙=(rr+λy-r-λx(r+λy)2).x-y=dxλ(x,y)

y˙=x+(rr+λy-r-λx(r+λy)2).y=dyλ(x,y)(9)

We name this model Elliptic Radial Isochron Cycle, or ERIC. An important feature of ERIC is that the angle in the plane is the phase of the oscillator, in particular since the phase is defined by the planar angle, the PRC following a horizontal perturbation of size ϵ toward the right can be computed in a straightforward way and is:

(10) PRCϵ,λ(θ)=cot-1(ϵ(csc(θ)+λ)+cot(θ))-θ

As said above, the effect of introducing λ is to smoothly flatten the PRC over half of the cycle, as illustrated in Appendix 1—figure 4A.

Appendix 1—figure 4
Effects of model parameters on the phase response curve (PRC).

(A) Effect of eccentricity λ: (A1) Increasing λ makes the limit cycle more elongated, while keeping the upper focus at the origin. (A2) The positive part of the PRC is then flattened; for all PRCs here s*=1 (no acceleration). (B) Effect of speeding up s*: (B1) When s* is increased, the spacing between isochrons in the sped up region increases, so this region contains less and less isochrons. Shown here is the case when the acceleration happens in the lower half-plane (α=1.5π, β=π), which is shaded in gray. (B2) The effect on the PRC is to shrink the portion where the oscillator is sped up, thus emphasizing the other part of the curve. Here, λ=0. (C, D) Effects of the sped up sector parameters α,β: the location and width of the sector determine which parts of the PRC get rescaled. (C1) The shaded regions are sectors with different values of α and β=π. (C2) The corresponding PRCs with,λ=0.s*=10. (D1) Sped up sectors located at α=15/8π with different widths β. (D2) The corresponding PRCs for the case λ=0.s*=10. All PRCs in this figure were computed with perturbation amplitude A=0.6.

We then introduce a second modification, allowing us to rescale the portion of the cycle where the PRC is flattened. We modify the equations by introducing a ‘speeding factor’ s so that

(11) θ˙=s(θ)

This will keep isochrons radial but will change their spacing. With the elliptic limit cycle, this leads to the differential equation in cartesian coordinates

(12) x˙=s(θ)dxλ(x,y)
(13) y˙=s(θ)dyλ(x,y)

We name this class of model ERICA. For simplicity, and to keep the system analytical, we restricted ourselves first to s functions linear by piece, that is, s=s>1 for one sector and s=1 otherwise. The sped up sector is centered at angle α and has width β, so that the modified period of the cycle is Ts*=β/s*+(2π-β). It is also convenient to define the rescaled angular velocity ωs*=2π/Ts*=2πs*β+(2π-β)s*. From there, the phase of the cycle as a function of the angle θ in the plane is the simple linear transformation (defining θ0=α-β/2):

(14) ϕs(θ)=ωs*θ=θ2πs*β+(2π-β)s*

for 0<θ<θ0,

(15) ϕs(θ)=ωs*(θ0+θ-θ0s*)=(θ0+θ-θ0s*)2πs*β+(2π-β)s*

for θ0<θ<θ0+β, and

(16) ϕs(θ)=ωs*(θ0+βs*+θ-(θ0+β))=(θ0+βs*+θ-(θ0+β))2πs*β+(2π-β)s*

for θ0+β<θ<2π. Those functions ensure that the rate of phase evolution in sector θ0<θ<θ0+β is 1/s* times the rate in the other sectors (compare θ coefficients in Equations 14-16), that angle 0 in θ is phase ϕs=0, and that phase is continuous so that ϕs(θ0)=ωs*θ0 and ϕs(θ0+β)=(θ0+β/s*)ωs*. Notice also that for θ=2π we get from Equation 16 ϕs(2π)=(βs*+2π-β)2πs*β+(2π-β)s*=2π as expected after one full cycle.

A full study of the possible behaviors of the ERICA model will be published elsewhere. Appendix 1—figure 4 illustrates different shapes of PRC obtained by varying λ,s*, α, and β independently, for various choices of sped up sectors.

For optimization purpose, it is easier to numerically compute the PRC in the following way. We consider an ensemble of angles θi linearly spaced on the interval [0,2π]. We then compute the corresponding position on the ERIC limit cycle, and compute numerically the angle θiϵ of the point at distance ϵ on the right. By definition of ERICA, the PRC for each index i is then

(17) PRCmodel(ϕs(θi))=ϕs(θiϵ)-ϕs(θi)

We checked that this PRC coincides both with the one computed from the integration of the ODEs (simulating the full stroboscopic map procedure) and with the analytical expression. To fit the experimental PRC, we run MC simulations to optimize parameters ϵ, λ, s*, sped up sector location α and width β, and the location of the zero-phase reference point on the cycle ϕ0. We minimize χ2:

(18) χ2=i(PRCicorr-PRCimodel)2σ2,

where PRCicorr are data points of the corrected experimental PRC and σ is the noise, estimated as the root mean square error of a Fourier fit to data (see Figure 8A). We use Markov chains to generate distributions of parameters and take the average values to plot the model limit cycle and PRC. The optimized parameter values are: ϵ=0.43, λ=0.53, s*=5.64, α=1.51π, β=1.15π, ϕ0=1.17π.

Evidence for intrinsic period changes

The PRC shifts are interpreted as change of period of the intrinsic oscillator because of entrainment.

We first checked that such PRC shifts are really needed to explain data. To do this, we computed an average PRC by only taking into account entrainment periods between 130 and 150 min. If we assume there is no period change, we can compute the maximum entrainment period compatible with such PRC, and found it is equal to 170 min. However it is clear from our data that we can go beyond this entrainment period, confirming the shift of PRC, in agreement with a change of intrinsic period.

We further have two evidences for such changes from release experiments. In those experiments, we entrain the oscillators for three cycles at 170 min, then let it go, and measure both the phase and instantaneous periods of the oscillator (Figure 7—figure supplement 1). When we measure the period, we first find it takes several cycles to come back to the intrinsic period, showing a long lived effect, unlikely coming from some artifact of period computation (Figure 7—figure supplement 1). A more direct way to estimate period change is to compute the return maps after the last pulse of DAPT. In Figure 8—figure supplement 3A, we consider all release experiments, which were done with zeitgeber period of 170 min. We take the phase one full cycle after the last pulse, and compute the return map from this phase (so for the first full cycle with no pulses). From this, we can compute the length of the cycle to have an average phase difference of 0 (corresponding to the period two full cycles after the last pulse). We found a cycle length of 150 min, compatible with the instantaneous period measurement and incompatible with a fast return to a 140 min intrinsic period.

We further show the Arnold tongue computed numerically with a similar PRC (1:1 entrainment region) and corresponding isophases assuming a fixed intrinsic period of 140 min in Figure 8—figure supplement 3B. We see much narrower width clearly incompatible with the data, again confirming the need for an intrinsic period change.

The last evidence is both more indirect and more mathematically grounded. It comes from the sigmoidal shape of the experimental entrainment phase as a function of zeitgeber period (Figure 7B, Figure 8E). In a nutshell, flatness of this curve is incompatible with the classical PRC theory, but can be easily explained if the internal period Tosc changes linearly with Tzeit.

To see this, we start with the classical condition for stability of the entrainment phase (see, e.g., Izhikevich, 2006),

(19) 1<dPTCdϕ|ϕ=ϕent<1

or equivalently:

(20) 2<dPRCdϕ|ϕ=ϕent<0

Now we also have by definition of the entrainment phase ϕent:

(21) PRC(ϕent)=2πTosc(Tosc-Tzeit),

where the factor 2πTosc is the unit conversion from time to phase. This equation defines in an implicit way a function ϕent(Tzeit), and the corresponding experimental curve is plotted in Figure 8E for our data. Taking the derivative on both sides with respect to the zeitgeber period, we get:

(22) dPRCdϕ|ϕ=ϕentdϕentdTzeit=-2πTosc

This allows to implicitly define the curve ϕent(Tzeit) with the help of the PRC:

(23) dϕentdTzeit=-2πTosc(dPRCdϕ|ϕ=ϕent)-1

There is a geometric interpretation here: the curve ϕent(Tzeit) is in fact a (rescaled) +π/2 rotation of the PRC! This can be seen here because Equation 23 combines a mirror image of the PRC along the y axis (minus sign) with a mirror image along the first diagonal (power −1 which corresponds to the inversion function exchanging ϕ and PRC(ϕ,PRC(ϕ))(PRC(ϕ),ϕ)). We know from classical group theory that two mirror images give one rotation with an angle equal to twice the angle between the axis (so here 2π/4=π/2). For instance, if the PRC is sinusoidal, ϕent(Tzeit) looks itself like (half a period) a vertical sinusoidal, as can be clearly seen, for example, in Granada et al., 2009. In Figure 8—figure supplement 3C, we further illustrate this rotation using the PRC computed from the data: the PRC rotated by π/2 is compared with the curve ϕent(Tzeit) computed numerically from it, showing perfect agreement.

Now, since the geometry of the PRC is constrained by Equation 20, this in turn imposes geometric constraints on ϕent(Tzeit) , for example, combining Equations 20-23 and counting time in units of the intrinsic period Tosc, 2π/Tosc=1, we get

(24) dϕentdTzeit>12

This defines an absolute, minimum slope of the curve ϕent(Tzeit). However, we see experimentally in Figure 8E that the entrainment phase is becoming almost constant for low and high values of the entrainment period, indicating a zero slope, which is thus theoretically impossible from Equation 24.

Now the simplest way to reconcile this observation with this calculation is to assume that the intrinsic period Tosc in fact depends on the zeitgeber. We then get a generalized version of Equation 22 with changing intrinsic period

(25) dPRCdϕ|ϕ=ϕentdϕentdTzeit=2πTosc(dToscdTzeitTzeitTosc-1)

giving the new implicit definition

(26) dϕentdTzeit=2πTosc(dToscdTzeitTzeitTosc-1)(dPRCdϕ|ϕ=ϕent)-1

Now we see that when dToscdTzeitTzeitTosc1, the right-hand side of Equation 26 can be very small, so that one can get dϕentdTzeit0, meaning that the entrainment phase barely changes as observed experimentally in Figure 8E. This effect will also come with a considerable enlargement of the Arnold tongue as discussed below. Also we notice that dToscdTzeitTzeitTosc=dlnToscdlnTzeit1 could be indicative of a mechanism where the oscillator adapts its intrinsic period (here to the zeitgeber period).

Practically, we obtain the Tosc(Tzeit) dependency by fitting the stroboscopic maps and entrainment phase from all experiments with 2.0 µM DAPT, using the PTC of our optimized model (Figure 8—figure supplement 3D-E). The values of Tosc are shown as data points in Figure 8D and Figure 8—figure supplement 3F.

Arnold tongue computations

To build the Arnold tongues, we first need to interpolate the period changes for period values not used in entrainment experiments. We used cubic spline interpolation to draw Figure 8D. For periods outside of the range of entrainment, in the absence of data some arbitrary choices have to be made. We know experimentally that entrainment does not occur below 120 min and above 200 min which indicates that dToscdTzeit is becoming smaller around those zeitgeber periods. More biologically, this indicates that the system does not adjust to any entrainment period. Again some arbitrary choices have to be made, but based on those constraints and our own interpolation, a cut-off for Tosc of ±20% of the natural period of 140 min was assumed, and is further consistent with the relative period changes observed experimentally in zebrafish segmentation clock mutants (Herrgen et al., 2010) and theoretically derived from delayed coupled models (Moreno-Risueno et al., 2010). With such cut-off, the obtained entrainment range of zeitgeber periods is from 118 to 200 min.

Computation of Arnold tongue and isophases in Figure 8F was done by iterations of the return map Equation 2 and detection of fixed points, using as the intrinsic period the interpolation and extrapolation shown in Figure 8—figure supplement 3F.

Lastly, we can use this model to compute the entrainment phase for different DAPT concentrations. For the main concentration of 2.0 µM DAPT, we have the optimized amplitude of perturbation ϵ=0.43. Assuming that the DAPT concentration is reflected in the strength of perturbation, we can find cross-sections of the isophases at different values of ϵ to have an excellent agreement with experiments, as illustrated in Figure 8—figure supplement 3G. We take ϵ=0.55 for 3.0 µM, ϵ=0.31 for 1.0 µM, and ϵ=0.13 for 0.5 µM DAPT, as shown on the Arnold tongue in Figure 8F. We capture both the presence/absence of entrainment and the ϕent(Tzeit) relation. In particular, we see that for lower DAPT concentrations, the range of entrainment is smaller but we still get plateaus in entrainment phase consistent with the change of period we computed.

Data availability

Source data and Python codes are made available via GitHub (https://github.com/PGLSanchez/EMBL_OscillationsAnalysis [copy archived at swh:1:rev:bd4a258570e7b7b5f0c179bc976f0d45988efa1a] and https://github.com/PGLSanchez/EMBL-files).

References

  1. Book
    1. Balanov A
    (2008)
    Synchronization: From Simple to Complex
    Springer Berlin Heidelberg. Google-Books-ID y6_UDAEACAAJ.
    1. Cross FR
    2. Siggia ED
    (2005) Mode locking the cell cycle
    Physical Review. E, Statistical, Nonlinear, and Soft Matter Physics 72:021910.
    https://doi.org/10.1103/PhysRevE.72.021910
    1. Granada A
    2. Hennig RM
    3. Ronache B
    4. Kramer A
    5. Herzel H
    (2009)
    Methods in Enzymology
    1–27, Chapter 1 Phase Response Curves: Elucidating the Dynamics of Coupled Oscillators, Methods in Enzymology, Academic Press, 10.1016/S0076-6879(08)03801-9.
  2. Book
    1. Kuramoto Y
    (2003)
    Chemical Oscillations, Waves, and Turbulence
    Courier Corporation. Google-Books-ID.
    1. Marr D
    2. Poggio T
    (1976)
    From understanding computation to understanding neural circuitry
    AI Memos 1:AIM-357.
  3. Book
    1. Pittendrigh CS
    2. Bruce VG
    (2015) V. An oscillator model for biological clocks
    In: Rudnick D, editors. Rhythmic and Synthetic Processes in Growth. Princeton University Press. pp. 75–110.
    https://doi.org/10.1515/9781400876167-006
  4. Conference
    1. van der Walt S
    2. Colbert SC
    3. Varoquaux G
    (2011)
    The NumPy Array: A Structure for Efficient Numerical Computation
    Computing in Science & Engineering. pp. 22–30.

Decision letter

  1. Sandeep Krishna
    Reviewing Editor; National Centre for Biological Sciences­‐Tata Institute of Fundamental Research, India
  2. Aleksandra M Walczak
    Senior Editor; CNRS LPENS, France

In the interests of transparency, eLife publishes the most substantive revision requests and the accompanying author responses.

[Editors' note: this paper was reviewed by Review Commons.]

https://doi.org/10.7554/eLife.79575.sa1

Author response

Reviewer #1 (Evidence, reproducibility and clarity (Required)):

Major comments:

Are the key conclusions convincing?

We discuss 4 key conclusions.

1. A PRC of the segmentation clock was constructed.

Although the authors have produced an interesting phase map, the regulation function F(\phi) of the circle map does not give the phase response curve (PRC) (Hoppensteadt & Keener 1982, Guevara & Glass 1982). This holds only when the system is stimulated with very short pulses (ideally Dirac δ), but the experimental pulses here are a quarter of the intrinsic period.

There are several definitions of the PRC (Dirac pulses PRCs, linear PRCs, etc.). We use the general definition from Izhikievich, 2007: “In contrast to the common folklore, the function PRC (θ) can be measured for an arbitrary stimulus, not necessarily weak or brief. The only caveat is that to measure the new phase of oscillation perturbed by a stimulus, we must wait long enough for transients to subside“.

The corresponding equation from Izhikievich (section 10.1.3) is

PRC(θ)= θnewθ

which is equivalent to our Equation 1.

Hence, the key assumption we make is that after perturbing the system, we are back on the limit cycle as pointed out by Izhikievich. We think this is a reasonable assumption, because the perturbation we impose is relatively weak, despite pulsing for almost one quarter of the intrinsic period. The concentrations of DAPT we used in this current study are just enough to elicit a measurable response, and further lowering the concentration does not result in entrainment within our experiment time (0.5uM, Figure S7B in submitted version of the manuscript). Additionally, we previously reported that periodic pulsing with 2uM DAPT did not result in change of the Notch signaling activity with respect to control samples (Sonnen et al., 2018). Along similar lines, the DAPT drug concentrations we used are much lower compared to what has been used in previous studies aiming to perturb signaling levels, e.g. 100uM and 50uM used in study of segmentation clock in zebrafish embryos (Özbudak and Lewis, 2008 and Liao et al., 2016, respectively), and 25uM used in study of the segmentation clock in mouse PSM cells (Hubaud et al., 2017). Combined, we reason that we apply weak perturbations that allow to extract the PRC of the segmentation clock during entrainment. Additional evidence that indeed we have revealed a meaningful PRC is provided below, please see our response to point #3.

2. Furthermore, in eq. 1 T_ext must be the winding number, and the modulus must be in units of phase, either one or two pi, for the circle map to be correct. Thus, calling the measured response of the system a PRC is not convincing.

We thank the reviewer for pointing this out. We indeed rescaled everything to express the PRC in units of phase. We made this more explicit and updated equations throughout the text.

3. The system is being entrained. Technically, It would also be easier to get the stroboscopic maps in the quasi-periodic regime since all the points in the circle will be sampled. Since no quasi-periodic response was demonstrated, the claim of entrainment is not convincing.

While, in principle, PRC can be indeed obtained from responses in the “quasi-periodic” regime, such an approach is, in practice, challenging due to the intrinsic noise. The closest approximation to this is the phase response after the first pulse, that we reproduce in Author response image 1 and compare to our inferred PRC, where we indeed clearly see a high noise level. Nevertheless, also the PRC based on the first pulse is in agreement with the PRC we derived from the entrainment data.

Author response image 1

In the entrained regime, one can get a much more reliable estimate of the phase response despite the noise. The level of noise in the stroboscopic map lowers as the samples approach entrainment (Figure S12), and the entrainment phase itself is a reliable statistical quantity that can be used to infer regions of the PRC as the detuning is varied.In addition, and maybe even more importantly, we identify several key features characteristics of entrainment, such as the change of entrainment phase as a function of detuning (Figure 7, Figure S6-S7 in submitted version of the manuscript) and the dependency of the time to entrainment as a function of initial phase (Figure 6). While additional features can be linked, in theory, with entrainment, i.e. period-doubling, higher harmonics (Figure 5), quasi-periodicity, we do not agree with the reviewer that all of these need, or in fact, can be found in the experimental data, in particular because of the influence of the noise. Conversely the positive experimental evidence that we provide for the presence of entrainment, combined with the theoretical framework we develop, justifies, in our view, the conclusions we make.

4. The response of the system to external pulses is compatible with a SNIC. This is compatible, but it is equally compatible with other explanations. Assuming that the PRC is the same as the regulation function F(\phi), the PRC in Kotani 2012 (PRL 2012 Figure 3C) would be a similar shape as that shown by the authors. Similar models to that in Kotani et al., have been studied, but a SNIC has not been found (an der Heiden & Mackey 1982). It is relatively straightforward to construct a phenomenological model with a SNIC, but having underlying biological insight is not guaranteed. No argument for choosing a SNIC is given, so this emphasis of the paper is not convincing.

It is true that the mapping of PRCs to oscillators is undetermined, in the sense that many systems could potentially give rise to similar PRCs. That said, there is value in parsimonious models, which often generalize very well despite their simplicity. This explains why in neuroscience, constant sign PRCs are generally associated with SNIC. There is a mathematical reason for this: 1-D oscillators with resetting (such as the quadratic fire-and-integrate model) are the simplest models displaying constant sign PRCs, and are the “normal” form for SNICs. In other words, SNIC bifurcations are among the simplest ones compatible with constant sign PRCs, and we think it is informative to point this out. In our manuscript, we go one step further by actually fitting the experimental PRC with a simple, analytical model that allows us to compute Arnold tongue for any values of the perturbation (contrary to more complex models).

Other models such as Kotani 2012 can display similar PRC shapes, but they are of mathematically higher complexity, and furthermore it is not clear how such systems might behave when entrained. For instance that model in particular uses delayed differential equations, and as such contains long term couplings, so that a perturbation might have effects over many cycles, which is not consistent with the hypothesis we here make of a relatively rapid return to the limit cycle. Furthermore, for more complex models, PRCs are analytical only in the linear regime, while our model is analytical for all perturbations. That said, we agree that other types of oscillators can be associated with constant sign PRCs, and we have given more details in this part, in particular we better emphasize the Class I vs Class II oscillators as a way to broaden our discussion on PRC, and emphasize the “infinite period” bifurcation category which is more intuitive and further includes saddle node homoclinic bifurcations.

5. The work demonstrates coarse graining of complex systems.

This conclusion is correct, but coarse graining theory-driven analysis and control of dynamical systems has been established for many years. What is new here is that it is applied specifically to the in vitro culture system of the mouse segmentation clock.

We agree it is new to successfully apply coarse-graining analysis and, importantly, control, to the in vitro culture system of the mouse segmentation clock. We also agree that such an approach has been pioneered and established for many years, especially in (theoretical) physics, but indeed, the key question is whether and how this can be applied to complex biological systems. Insights coming from theoretical considerations on idealized physical systems might not necessarily apply to biology, as already pointed out by Winfree.

There are still very few examples in biology with coarse graining similar to what we do here. We think there is immense value in demonstrating that quantitative insights, and control of the biological systems, can be obtained without precise knowledge of molecular details, which is still counter-intuitive to many biologists. In this sense, we think our report will be of interest to both colleagues within the field of the segmentation clock and also to anyone interested to in the question, how theory and physics guided approaches can enable novel insight into biological complexity.

Should the authors qualify some of their claims as preliminary or speculative, or remove them altogether?

Following on the points above, each of these needs to be corrected or re-done, and/or the conclusions need to be modified accordingly.

We have modified the manuscript in response to all those points.

6. Would additional experiments be essential to support the claims of the paper? Request additional experiments only where necessary for the paper as it is, and do not ask authors to open new lines of experimentation.

If the authors wish to make the strong claim of determining a true PRC, Dirac δ-like perturbation needs to be applied, or approximated by short time duration pulses compared to the intrinsic period.

Please refer to our response to point #1 and #3.

7. Are the suggested experiments realistic in terms of time and resources? It would help if you could add an estimated cost and time investment for substantial experiments.

It's not clear to this reviewer if it is feasible to deliver a very short pulse and record a response. But this may not be relevant, see above.

Please refer to our response to point #1 and #3.

Are the data and the methods presented in such a way that they can be reproduced?

Yes.

Are the experiments adequately replicated and statistical analysis adequate?

Yes.

Minor comments:

Specific experimental issues that are easily addressable.

No issues.

Are prior studies referenced appropriately?

Yes.

8. Are the text and figures clear and accurate?

Figure 1D illustrates how a PRC should be obtained, but doesn't show the experimental protocol applied in the paper.

Figure 1D is a general introduction on the phase description of oscillators and phase response. It demonstrates how a perturbation can change the phase and is not supposed to represent the experimental protocol. We describe how data are analyzed and how phases are extracted in Supplementary Note 1.

9. In Figure 5B, 10 μm DAPT, the traces are already synchronized before the pulse train starts, which makes the subsequent behavior difficult to interpret.

It appears here that by chance, the samples were already almost synchronized. We notice however that the establishment of a stable rhythm with the pulses (which here is not a multiple of the natural period) supports entrainment, and is already evident when looking at the timeseries with respect to the perturbation. The temporal evolution of the instantaneous period further confirms this, showing a change in period close to ½ zeitgeber period (which is very different from the natural period of ~140 mins). This also relates to point #35, in reply to both comments we have further expanded this figure to better show the 2:1 entrainment, adding statistics on the measured period and period evolution for a zeitgeber period of 300 mins.

10. Do you have suggestions that would help the authors improve the presentation of their data and Conclusions?

The text includes several paragraphs reviewing broad principles of coarse graining and making general conclusions. This is confusing, because, as mentioned above, there is no new general advance in this paper. The interesting contributions here are specific to the applications to the segmentation clock, and the text should be focused on this aspect.

As commented above for #3 , we respectfully disagree that there is no “new general advance” in this paper. It is far from obvious that a complex ensemble of coupled oscillators implicated in embryonic development would be amenable to such coarse-graining theory. Of note, we still do not have a full understanding of neither the core oscillators in individual cells, nor what slows these down and eventually stops the oscillations, and multiple recent works suggest that both phenomena are under transient nonlinear control (e.g. our own work in Lauschke 2013). It is remarkable that despite this lack of detailed mechanistic insight, general entrainment theory can be applied to the segmentation process at the tissue level. We further show that classical entrainment theory alone is not sufficient to account for the experimental findings. Specifically, we need to account for a period change that we interpret as an internal feedback, an insight that would be impossible without our coarse-graining approach. While the results might of course be specific to the segmentation process, we think our approach motivated by coarse-graining theory and leading to new insights into the process is of general interest. We tried to make these points explicit in our conclusion.

Reviewer #1 (Significance (Required)):

Describe the nature and significance of the advance (e.g. conceptual, technical, clinical) for the field.

Description of the complex mouse segmentation clock in terms of a simple model and its PRC is an interesting, original and non-trivial result. The proposal that the segmentation clock is close to a SNIC bifurcation provides a consistent dynamical explanation of slowing behavior that has been recognized for some time, but not fully understood. This proposal also raises a hypothesis about the behavior of the underlying molecular regulatory networks, which may be tested in the future. The increase or decrease of the intrinsic period due to the zeitgeber period is not expected from theory, pointing to structures in internal biochemical feedback loops, an idea which again may be tested in the future. Also surprising from a theoretical perspective, the spatial gradient of period in the system persisted after entrainment. Although the categorization of the generic behavior is interesting, by its nature there is little from this that might give a typical developmental biologist any conclusions about pathways or molecules. The successes and limits of the theoretical description do nevertheless focus future attention on interesting behaviors.

11. Place the work in the context of the existing literature (provide references, where appropriate).

Such an analysis of the segmentation clock is based strongly on the experimental system and results in Sonnen et al., 2018, and goes well beyond it in terms of the dynamical analysis. It provisionally categorizes the mouse segmentation clock as a Class I excitable system, allowing its dynamics at a coarse grained level to be compared to other oscillatory systems. In this aspect of simplification, it is similar to approach of Riedel-Kruse et al., 2007 who used a mean-field model of oscillator coupling to explain the synchrony dynamics observed in the zebrafish segmentation clock in response to blockade of coupling pathways, thereby allowing a high-level comparison to other synchronizing systems.

It is interesting the reviewer sees similarities with the work of Riedel-Kruse et al., which uses a mean-field variable Z that corresponds to a classical approach, as described in Pikovsky’s textbook, to quantify synchronization of oscillators. In our view, while of course we work in the same context of coupled oscillators in the PSM, our approach based on perturbing and monitoring the system’s PRC in real-time provides a novel strategy to gain insight. This is evidenced by the fact that our quantifications of synchronization and insight into the PRC is the basis to exert precise control of the pace and rhythm of segmentation.

State what audience might be interested in and influenced by the reported findings.

Developmental biologists, biophysicists.

Define your field of expertise with a few keywords to help the authors contextualize your point of view. Indicate if there are any parts of the paper that you do not have sufficient expertise to evaluate.

Developmental biology, somitogenesis, dynamical systems theory, biophysics, cell signaling.

Reviewer #2 (Evidence, reproducibility and clarity (Required)):

Summary:

This is a beautifully elegant study that tests how previously published theoretical predictions about entraining nonlinear oscillators applies to a biological oscillator, the segmentation clock. The authors use a combination of state of the art experimental techniques, signal processing and analytical theory to reach a series of interesting and novel conclusions.

They show that the segmentation clock period can be entrained through Notch inhibitor (DAPT) pulses acting as an external clock (referred to as zeitgeber) using a previously developed and sophisticated microfluidic perfusion system. Pulsing DAPT every 120 to 180min can change the internal clock period while entrainment beyond this range leads to higher order coupling to the zeitgeber period, i.e. entrainment of every other pulse. They then perform entrainment experiments where the concentration of DAPT is changed to elicit a change in the strength of interaction between the internal clock and the external stimulus (referred to as zeitgeber strength); interestingly at low strength response to entrainment is more variable leading to entrainment occurring in some samples while others remain unaffected (Figure 4A); overall, higher concentration leads to faster entrainment (Figure 4C). The experimental data is then analysed using stroboscopic maps to reveal that a stable entrainment phase shift is achieved between the internal clock and the external zeitgeber. Phase response curve (PRC) analysis indicates that the system response is not sinusoidal but predominantly characterised by negative PRC, a behaviour consistent with saddle-node on invariant cycle (SNIC); it also reveals that the intrinsic period changes in a non-linear way and that this effect is reversible when external stimulation stops. Finally, a theoretical model is proposed to represent the segmentation clock as a dynamical system; this is based upon Radial Isochron Cycle with Acceleration (ERICA), an extension motivated by the PRC analysis results which are incompatible with a Radial Isochron Cycle (RIC); this model has predictive capability and could be used to design new control strategies for entrainment of the segmentation clock.

This study makes a series of key conclusions which are of particular importance in understanding the dynamic response of a biological oscillators. Firstly, given it's the characteristics of the dynamic response to entrainment, the segmentation clock is likely close to a SNIC bifurcation and this can explain the tendency for relaxation of the period over time. Secondly, the clock period was changed in a non-linear way in the direction of the zeitgeber period, a finding which is interpreted to indicate the presence of feedback of the segmentation clock onto itself, potentially via Wnt. This makes an excellent prediction that if tested experimentally would greatly improve the impact of the study. It is also noted that the entrainment of the segmentation clock does not abolish spatial periodicity and phase wave emergence suggesting that single cell oscillators can adjust to periodic perturbation while maintaining emergent properties. This is also a significant result that would need to be followed up with experiments and computation however would be best suited to a separate study.

Major comments:

12. The coarse graining is a major point that would need to be clarified since the rest of the analysis and theoretical modelling in the paper flow from this. Firstly, the interpretation of the schematic in Figure 1A on experimental data collection is not immediately obvious to the reader, lacks a clear flow between the different panels or steps (which could be numbered for example) and does not have a legend to indicate the different colour mapping.

We are grateful to the reviewer for this comment. We have implemented in Figure 1A all the changes suggested by the reviewer: we numbered the different steps and have added a colour mapping. In addition we have rephrased the caption of Figure 1A to better connect the experimental steps.

13. Secondly, Figure 2A which explicitly addresses coarse graining is not clear enough. Is the message here that by excluding the inner parts of the sample with a radial ROI, a similar dynamic response is observed over time?

Yes, indeed this is the point and we have adjusted the figure and text to explain this better. Our goal is to focus on the quantification of segmentation pace and rhythm. This is best captured by reporters such as LuVeLu, which has maximum intensity in regions where segment forms, and which dynamics is known to be strongly correlated to segmentation (Aulehla et al., 2007; Lauschke and Tsiairis et al., 20132). The global ROI is thus expected to precisely capture these segmentation and clock dynamics and we have now included more validation data and have also edited the text to make this very important point clearer:

“To perform a systematic analysis of entrainment dynamics, we first introduced a single oscillator description of the segmentation clock. We used the segmentation clock reporter LuVeLu, which shows highest signal levels in regions where segments form Aulehla_A_2007}. Hence, we reasoned that a global ROI quantification, averaging LuVeLu intensities over the entire sample, should faithfully report on the segmentation rate and rhythm, essentially quantifying 'wave arrival' and segment formation in the periphery of the sample.”

Figure 2A indeed shows that the dynamics (from the timeseries) is very similar when considering the entire field of view (global ROI) or when considering only the periphery of the 2D-assay (excluding central regions). We modified Figure 2A to clarify this point by indicating each measurement as either global ROI or global ROI minus the diameter of the excluded circular region (e.g. global ROI – 50px). We also emphasized in the caption that timeseries are obtained using global ROI, unless otherwise specified. We included a link (https://youtu.be/fRHsHYU_H2Q) in the caption to a movie of 2D-assay subjected to periodic pulses of DAPT (or DMSO) and corresponding timeseries from global ROI.

Since the inner part of the sample corresponds to the posterior side how do we interpret similarities and differences between signals with different ROIs?

As stated above, the global ROI measurements essentially capture the signal at the periphery where segments form and faithfully mirrors segmentation rate and rhythm. We have now included a comparison to the center ROI, also in response to reviewer’s comments, see our response #34.

The result shows that the period and PRC in the center matches the one found in the periphery, i.e. global ROI. We have shown previously that center and periphery differ in their oscillation phase by 2pi, i.e., one full cycle (Lauschke et al., 2013). We interpret these findings as confirmation of our analysis strategy, i.e. the global ROI allows a very reproducible, unbiased quantification that reports on segmentation clock and period.

14. A quantitative analysis of essential coarse-grained properties such as period and amplitude should be performed for different ROIs and across multiple samples. As this effectively masks any spatial differences, limitations of this approach should be clearly stated in the Discussion. For example in lines 466-470 where it is difficult to interpret the slowing down tendency and relate back to single cell level.

As outlined in our response to comment #13 and also #34, we chose an analysis that allows to determine the segmentation pace and rhythm, i.e. segment formation, which is well captured by LuVeLu signal and a global ROI analysis. We agree that a spatially resolved analysis of dynamic behaviour is important (and indeed a gradient of amplitude might be relevant in such context), but we think this is beyond the scope of the current study focused on the system level segmentation clock behaviour. We have revised the discussion as suggested by the reviewer to make this point approach and the need for future studies clearer.

15. The functional characterisation of the sample using LFNG, AXIN2 and MESP2 is unclear. The images included in Figure 2D representing expression observed when tissue explants are grown within the microfluidic chip are difficult to interpret and would require a more detailed description of anterior-posterior, pillars etc; it is also difficult to view the bright-field since it is presented as a merged image.

It is particularly difficult to see the somite boundaries for the same reason. In lines 113-117 the authors state that the global oscillation period matches the periodic boundary formation. How do we reach this conclusion from these images? What is the variability between samples?

If these two issues would be addressed it would increase confidence in the coarse graining argument and thus would strengthen the importance of the findings in the study.

We thank the reviewer for this feedback, and we have added more quantifications to address this point directly in the modified Figure 2. Importantly, we added the quantification of the rate of segmentation in multiple samples based on segment boundary formation (new Figure 2D) and compared this to the global ROI quantifications using the reporter lines LuVeLu. This data provides clear evidence that the quantification of global ROI reporter intensities closely matches the rate of morphological segment boundary formation. In addition, we show that segment formation and also Wnt-signaling oscillations (Axin2-Achilles) and the segmentation marker Mesp2 (Mesp2-GFP) are all entrained to the zeitgeber period. We have also revised the text to clarify this important validation of our quantitative approach.

In addition, we provide, in the revised Figure Suppl. 2, details of entrained samples, focusing on the segmenting regions. The brightfield and reporter channels were separated, emphasizing the segment boundaries and the expression pattern of the reporters. For ease of visualization, these samples were also re-oriented so that the tissue periphery (corresponding to anterior PSM) is at the top while the tissue center (corresponding to the posterior PSM) is at the bottom. This now additionally better shows the localization of the different reporters with respect to the segment boundary. We also included supplementary movies showing timelapse of samples expressing either Axin2-GSAGS-Achilles or Mesp2-GFP that were subjected to periodic DAPT pulses, with their respective controls.

Several minor points could be addressed to improve the manuscript and are listed below:

16. Figure 1 A the colormap and axes for the oscillatory traces should be defined

We thank the reviewer, and we have modified the figure accordingly (related to point 12). A colormap and axes for the illustrated timeseries are now included.

17. Strength of zeitgeber is not defined and there is no analytical expression provided; how does it relate to DAPT concentration? Is the fact that low DAPT concentration corresponds to weak strength expected or is it a result?

Zeitgeber strength generally refers to the magnitude of the perturbation periodically applied to an oscillator. With DAPT pulses, our expectation was that both the duration of the pulse and the drug concentration could influence the strength. Practically, the pulse duration was kept constant for all experiments and the concentration was varied. We thus expected that DAPT concentration would indeed be correlated to zeitgeber strength. We have discussed multiple evidence supporting this assumption in the main text, and this is indeed a result. In particular, as explained in the section “The pace of segmentation clock can be locked to a wide range of entrainment periods”, higher DAPT concentration gives rise to faster and better entrainment, as expected from classical theory. In the context of Arnold tongue, weaker zeitgeber strength corresponds to narrower entrainment region, which is experimentally observed (Figure 8F, showing regions where the clock is entrained).

From a modelling standpoint, Zeitgeber strength corresponds to parameter A which is the amplitude of the perturbation. Possible zeitgeber strength was inferred from the model by matching the experimental entrainment phase ϕent (Tzeit) with that obtained from the model isophases. As explained in Supplementary Note 2, we tested four concentrations of DAPT (0.5, 1, 2, and 3 uM) respectively corresponding to A values of 0.13, 0.31,0.43, 0.55. As we can see, those A values are not linear in DAPT concentrations, which is expected since multiple effects (such as saturation) can occur.

18. In some figures it looks like the amplitude of oscillations may change with DAPT concentration and hence zeitgeber strength? Is this expected?

We have not systematically analyzed the amplitude effect and have, intentionally, focused on the period and phase readout as most robust and faithful parameters to be quantified. Regarding the amplitude of LuVeLu reporter, we are cautious given that it is influenced, potentially, by the (artificial) degradation system that we included in LuVeLu, i.e. a PEST domain. This effect concerns the amplitude, but not the phase and period, explaining our strategy.

That said, we agree with the referee that DAPT concentrations might change the amplitude of oscillations. Such change could even play a role in the change of intrinsic period (in fact a similar mechanism drives overdrive suppression for cardiac oscillators, Kunysz et al., 1995). But since the change of period can be more easily measured and inferred, we prefer to directly model it instead of introducing a new hypothesis on amplitude/period coupling, at least for this first study of entrainment.

19. Figure 2A including the black area creates confusion and it is unclear which ROI is used in the rest of the study; consider moving this to a supplementary figure perhaps

We thank the reviewer for this feedback (related to point #13), and we have modified the figure accordingly. As we responded to point 13: We modified Figure 2A, by indicating each measurement as either global ROI or global ROI minus the diameter of the excluded circular region (e.g. global ROI – 50px). We also emphasized in the caption that timeseries are obtained using global ROI, unless otherwise specified.

20. What type of detrending is used in Figure 2 and throughout (include info in the figure legend)?

We used sinc-filter detrending, described and validated in detail previously (Mönke et al., 2020), as specified in Supplementary Note 1: Materials and methods > H. Data analysis > Monitoring period-locking and phase-locking: In this workflow, timeseries was first detrended using a sinc filter and then subjected to continuous wavelet transform. We thank the reviewer for pointing out that this detail is lacking in the figure captions, and we have modified the captions accordingly.

21. Figure 2D merged images are difficult to read/interpret (see major comments)

We thank the reviewer for this comment, and we have modified the figure accordingly (please see response to related point #15).

22. Kuramoto order parameter is used to quantify the level of synchrony across the different samples however it is not defined in the text. Is it also possible to assess variability in each sample? For example how quickly does entrained occur in each sample? How faithfully the peaks of expression beyond 80min (to exclude initial unsynchronised state) match with zeitgeber time? This would help make the point that weak strength leads to a more variable response which is an interesting finding.

We have now added a mathematical definition of the Kuramoto parameter in Supplementary Note 1.

A high order parameter corresponds to coherence between samples, as also elaborated in respective figure captions (e.g. in the caption for polar plots in Figure 4D).

In terms of variability in response to entrainment, we thank the reviewer for the comments, which has prompted us to perform an additional analysis, now included as Figure S13 in the Supplement.

Briefly, Author response images 2–6 and Figure 4—figure supplement 2 show how different samples get synchronized with the zeitgeber. To do this, we first represent the zeitgeber signal as a continuous uniformly increasing phase (“zeitgeber time”) with period Tzeit:ϕ˙zeit=2π/Tzeit. The initial condition for ϕzeit is chosen so that the zeitgeber phase at the moment of last pulse is matching the experimental entrainment phase for each Tzeit. We plot cos ϕ(t) for each sample (dotted lines) and the zeitgeber phase (magenta line). To quantify how well each sample is following the zeitgeber time, we compute the Kuramoto parameter: Rj(t) = |12(eiϕj(t)+eiϕzeit(t)) |. By the end of experiment most samples reach Rj 1, indicating entrainment. Most samples need 3 zeitgeber cycles to become entrained. For Tzeit=120 min the entrainment takes much longer (edge of the Arnold tongue). For Tzeit=140 min there is much variability, which can be explained by the horizontal region in the PRC around the entrainment phase. As suggested by the referee, synchronization is faster for higher DAPT concentration. So those dynamics are indeed consistent with the expectation from classical PRC theory.

Author response image 2
Tzeit = 120.
Author response image 3
Tzeit = 140.
Author response image 4
Tzeit = 150.
Author response image 5
Tzeit = 160.
Author response image 6
Tzeit = 180.

23. Do samples change period to Tzeit in similar ways – i.e. patterns over time. It looks like the kuramoto order parameter and period drop initially – why?

We do not have a direct answer as to why the Kuramoto first order parameter and the period drop for the condition the reviewer specified. It has to be noted though that because of how wavelet analysis is done (cross-correlation of the timeseries with wavelets), the period and phase determination at the boundaries of the time series are less reliable (edge effects, see Mönke et al., 2020). Because of this, we should take caution when considering data to and from the first and last pulses, respectively. This was explicitly stated in the generation of stroboscopic maps: “As wavelets only partially overlap the signal at the edges of the timeseries, resulting in deviations from true phase values (Mönke et al., 2020), the first and last pulse pairs were not considered in the generation of stroboscopic maps.”

24. In Figure 4C why is the Kuramoto order parameter already higher in the 2uM DAPT conditions at the start of the experiment?

Samples can, by chance, start synchronously and this results in a high Kuramoto first order parameter. Because of this likelihood, it is thus important to interpret the entrainment behaviour of multiple samples using various readouts, in addition to a high Kuramoto first order parameter. We investigated entrainment of the samples based on several measures: multiple samples remaining (or becoming more) synchronous (because each sample actively synchronizes with the zeitgeber), period-locking (where the pace of the samples match the pace of the zeitgeber, which can be distinct from natural pace), and phase-locking (where there is an establishment of a stable phase relationship between the samples and the zeitgeber).

25. Figure 3C and Figure S2 require statistical testing between CTRL and DAPT in each condition.

p-values were calculated for the specified conditions and were added in the caption of the figures. These values are enumerated here:

– Figure 3C

170-min 2uM DAPT (vs DMSO control): p < 0.001

– Figure S2

120-min 2uM DAPT (vs DMSO control): p = 0.064

130-min 2uM DAPT (vs DMSO control): p = 0.003

140-min 2uM DAPT (vs DMSO control): p = 0.272

150-min 2uM DAPT (vs DMSO control): p = 0.001

160-min 2uM DAPT (vs DMSO control): p < 0.001

170-min 2uM DAPT (vs DMSO control): p < 0.001

180-min 2uM DAPT (vs DMSO control): p = 0.001

To calculate p-values, two-tailed test for absolute difference between medians was done via a randomization method (Goedhart, 2019). This confirms that the period of samples subjected to pulses of DAPT is not equal to the controls, except for the 140-min condition (where the zeitgeber period is equal to the natural period, i.e. 140 mins).

26. Figure 3A gray shaded area not clearly visible on the graph.

We have decided to remove the interquartile range (IQR) in the specified figure as it does not serve a crucial purpose in this case. By removing it in Figure 3A, the timeseries of individual samples are now clearer.

27. Figure 6C colour maping of time progression is not clearly visible on the graph; the interpretation of this observation is unclear in the text and the figure.

We agree that the low quality of the image is unfortunate, and it seems that our file was greatly compressed upon submission. We have checked the proper quality of figures in the resubmitted version of the manuscript.

Regarding the interpretation of Figure 6C, we conclude that in our experiments the entrainment phase is an attractor or stable fixed point, in line with theory (Granada and Herzel, 2009; Granada et al., 2009),. We had elaborated this in the text (lines 248-252 of the submitted version of the manuscript): “at the same zeitgeber strength and zeitgeber period, faster (or slower) convergence towards this fixed point (i.e. entrainment) was achieved when the initial phase of the endogenous oscillation (φinit) was closer or farther to φent.”

28. Figure 7A circular spread not clearly visible on the graph.

Similar to point #27, we have provided a high resolution graph for the re-submission and hopefully resolved this issue.

29. Figure S7A difficult to see the difference between colours

See point #28.

30. Is it possible to compare the PRC and the plots of period over time during entrainment? The PRC is mainly negative (Figure 8A1,A2), in my understanding this means a delay, however the periods seem to decrease over time before entraining to the Tzeit (Figure 3B). Is this reflective of a decrease in Kuramoto parameter and potential de-synchronisation of single cells before re-synchronisation at Tzeit?

To address this question, we now plot the Phase response with colors indicating pulse number in new Supplementary Figure S13. While capturing the entire PRC as a function of time would require many more experiments (in particular to sample the phases far from entrainment phase), we still clearly see that the PRCs appear to translate vertically as the oscillator is being entrained, i.e. the latter time points are shifted up (down) for T_zeit = 120 (170) min, respectively.

31. Figure 8A What is the importance/meaning of the PRC being similar shape between different entrainment periods? Does this reflect that the underlying gene network is the same?

If one single gene network is responsible for oscillations, we expect from dynamical systems theory that the PRC are not only of similar shape but actually the same, independent of the entrainment period. What is surprising is that the PRC for different entrainment periods do not overlap, and the simplest explanation for this is that the intrinsic period changes with entrainment, all things being kept equal (including the underlying gene networks). This relates to the previous point since we indeed observe that the PRC “translates” vertically with the pulse number for longer periods. The change of period might be due to a long-term regulation as detailed in the discussion.

32. The spatial period gradient and wave propagation under DAPT (Figure S8) should be included in the results and not just the discussion.

We fully agree with the reviewer that both the establishment and the maintenance of a spatial phase gradient is of great interest. However, many more experiments would be required to fully quantify and understand the processes at play here, which we believe to be out of the scope of the current manuscript. To keep the focus of the paper on the global segmentation clock itself, we prefer to keep this figure in Supplement.

Reviewer #2 (Significance (Required)):

We currently do not have a detailed understanding of how biological oscillators integrate local signals from their neighbours as well global external signals to give rise to complex patterning that is important for embryonic development. Main bottlenecks that hinder our understanding are lack of real-time endogenous dynamic response together with known global inputs as well as comprehensive models that can explain emergent behaviour in a variety of tissues.

This study goes a long way in addressing these bottlenecks in the embryonic tissue responsible for somite formation, a dynamical and oscillatory system also known as the segmentation clock. Firstly, they rely on a state-of-the-art previously developed system to entrain endogenous response in live tissue explants using precise microfluidic control. They test the complete range of exogenous perturbation periods and use an existing live reporter (LuVeLu) to monitor endogenous response. They also identify higher order coupling relationships whereby every other LuVeLu peak is entrained through external stimulation.

As the stimulation system does not control but rather perturb the endogenous response, the observations from LuVeLu provide a unique opportunity in understanding input-output relationships and thus describing the dynamic response of the segmentation clock. Authors propose to study dynamic behaviour of the clock using coarse-graining and focus on describing the overall response over time while amalgamating spatial information. Appropriate coarse-graining is an important strategy in addressing complex problems and is widely used. They use sophisticated methodology such as phase response curves and Arnold tongue mapping to make several important observations. For example the nonlinear shortening and elongation of the period in response to stimulation is particularly interesting since this may indicates a feedback of the clock onto itself potentially via Wnt. Another key observation is that the spatial periodicity and phase wave activity persists in the perturbed conditions suggesting that individual single cell oscillators can adjust their behaviour to external input while retaining coordination with their neighbours. Finally, the authors go on to construct a general dynamical model of the segmentation clock and use this to conclude that the intrinsic period of the oscillator is altered and that the oscillator can be considered excitable.

This work sheds light onto mechanisms of coordination of Notch activity in assemblies of cells observed in living tissue, an area of research that is important not only for somitogenesis but also for understanding gene expression patterning in many other tissues where Notch plays a critical role, for example in the development of the neural system and organs. As a study of a real-world nonlinear oscillator this work is directly of interest to theoreticians and synthetic biology experts interested in understanding complex patterning and emergence.

Reviewer #3 (Evidence, reproducibility and clarity (Required)):

In this manuscript, authors studied the system-level responses of the somite segmentation clock by the coarse-grained theoretical-experimental approach, applying the theory of entrainment to understanding the phase responses of mouse pre-somitic mesoderm (PSM) tissues in the presence of periodic perturbation of Notch inhibitor DAPT generated by micro-fluidics technique. It was demonstrated that the segmentation clock is responsive to diverse range of the perturbation-periods from 120 to 180 min, can be period- and phase-locked, and the efficiency is dependent of the DAPT concentration (input-strength). The authors also observed two cycles of the segmentation-clock ticking in single cycles of 300 or 350 min period-perturbation, suggesting that higher order (2:1 mode) entrainment. They also applied stroboscopic maps to analysis and found that entrainment-phases are dependent of period of DAPT pulses, which is recapitulating theoretical predictions. The estimation of the phase response curve (PRC) of the segmentation clock revealed that the inferred PRC is an asymmetrical and mainly negative function, which represents characteristic features in oscillators that emerge after saddle-node on invariant cycle (SNIC) bifurcation. These results also indicated that the segmentation clock changed the intrinsic period during entrainment.

Major comments:

33. I have major concerns about the relevance of the global time-series analysis proposed in Figure 2 and conclusion about the changes of the intrinsic period during entrainment. The validity of the global time-series analysis should be carefully analyzed, because it could bring artifacts in estimated values of the intrinsic period. The authors concluded (page 3, line 172) that the period calculated by the global analysis represents similar values with the rate of segment formation, but there is no data about the quantification of the periods of segmentation, such as the frequency of Mesp2 reporter expression.

We thank the reviewer for this feedback. We have now added the quantification of the period of segment formation (new Figure 2E) and show its strong correspondence to the dynamics of reporters used (Lfng, Axin2, and Mesp2). Please see also our response to point #15 with additional comments regarding the validation of the global time-series analysis.

34. Another related issue is the presence of spatial period gradient as mentioned (page 13, line 524). One possible approach to circumvent this issue would be "local" time-series analysis; for instance, just focusing on the "putative posterior" regions that are close to source-positions of waves. Authors can re-compute and estimate PRCs by using such a method.

We thank the reviewer for this suggestion and have accordingly now included the analysis of a localized ROI at the center (center ROI) of the 2D-assays (new Figures S5-S6). We also computed the PRC from center ROIs as shown in Author response image 7. We note strong correspondence between the global ROI and the center ROI.

Author response image 7
Comparison of the PRC obtained with global ROI and center ROI, T_zeit = 130 min (left) and 170 min (right).

35. I have another major concern about the evidence of higher order entrainment shown in Figure 5. If the 1:2 entrainment is successful, we can expect that the values of observed period is close to the half of the period of pulses; However, the period shown in Figure 5B looks like 185 min longer than the half of 350 min. Is this gap due to the temporal accuracy of time-lapse movies?

We do not think the discrepancy comes from a problem of temporal accuracy as the temporal accuracy is the same for all movies and there is no reason why there would be a specific issue for this set of experiments. In addition, we have re-analyzed the data to calculate the period from the stroboscopic maps. Mathematically speaking, we take the stroboscopic map as ϕnew = (ϕold +TzeitTobserved2π) mod 2π , and use this to estimate the period of oscillation in entrained samples Tobserved , in particular inverting the formula for 1:2 entrainment we have :

TzeitTobserved=(ϕnew -ϕold )/2π + 2

The advantage of this method is that it gives a more ``instantaneous” estimation of the period.

The results are as follows:

350 10uM: 187 +- 8 min (average across entrained samples from the last zeitgeber period)

350 5uM: 193 +- 13 min (average across entrained samples from the last zeitgeber period)

300 2uM: 148 +- 8 min (averaged across entrained samples and from two last periods)

This additional analysis is in agreement with the wavelet analysis.

The reviewer is right that for 350 minutes, entrained samples show an observed period that is higher than expected, also based on this new additional analysis. The reason for this is not known. One explanation is the relatively short observation time, especially considering for pulses separated by as much as 350-minutes, i.e. only 3 pulses are applied. [We notice that for 300 minutes pulses, the period converges to 150 mins between the 3rd and the 4th pulse]. We have adjusted the text in the Results section to reflect that for 350min entrained samples, the observed period ‘approaches’ the predicted value, while for 300min entrained samples, the observed period is very close to it, i.e. 147mins In addition, we comment that the phase distribution narrows with time, another indication supporting higher order entrainment.

36. Also, authors showed the period evolution towards 1:2 locking with just one condition (350 min). Authors can show the data for multiple conditions as in Figure 3D, at least for 300 min and 325 min pulses and add the data about final entrained period with statistic analysis that supports the difference between the entrained period and the natural period (140 min).

We thank the reviewer for this feedback and have modified the figure accordingly. In particular, in Figure 5A, we have added the period evolution plot for samples subjected to 300-min periodic pulses of 2uM DAPT (or DMSO for control). Additionally, we have added Figure 5D, which plots the average period in the 300-min and 350-min conditions. We summarize the median average period here with computed p-values:

– 300-min pulses of 2uM DAPT (or DMSO for control): p-value = 0.191

CTRL: 130.39 mins

DAPT: 146.45 mins

– 350-min pulses of 5uM DAPT (or DMSO for control): p-value = 0.049

CTRL: 127 mins

DAPT: 174.86 mins

– 350-min pulses of 10uM DAPT (or DMSO for control): p-value = 0.016

CTRL: 142.82 mins

DAPT: 185.12 mins

Minor comments:

37. The authors can draw vertical lines indicating the T_zeit in Figure 3B, Figure 4B and Figure 5B in order to help comparisons between T_zeit and patterns of period (solid lines).

We thank the reviewer for this comment. We have accordingly added a horizontal line indicating Tzeit in Figures 3B, 4B, S4A, and S5A (figure panel numbers based on the submitted version of the manuscript). We similarly added a horizontal line indicating 0.5Tzeit in the period evolution plots of 300-min and 350-min conditions in Figures 5A and 5B, respectively.

38. In Figure 5A, the authors can show period evolution in the case of 300 min DAPT-pulses as shown in Figure 5B.

We thank the reviewer for this feedback (related to point #36), and we have modified the figure accordingly.

39. In Figure 6B DAPT panel, the authors can draw the points of phi_ent as shown in Figure 7A.

We thank the reviewer for this comment, and we have modified the figure accordingly.

40. In Figure 8F, authors can put the information about DAPT concentration at the right y-axis.

This is a similar comment as point #17, see above. In brief, we do not know the precise relation between the strength of the perturbation in our model and DAPT concentration, zeitgeber strength was inferred from the model by matching the experimental entrainment phase ϕent (Tzeit) with that obtained from the model isophases.

41. In Figure 8G, the PRC in the panel "170 mins" does not have any fixed point (cross sections with horizontal lines of "0" phase response). If entrainment is successful, there should be stable and unstable fixed points, but those are absent, although 170 min pulses succeeded in the entrainment as shown in Figure 3D. Authors can explain where the fixed points are.

The fixed points are indeed defined by the intersection with a horizontal line, but not with the ‘0’ line. They are found where the phase response compensates for the detuning/period mismatch, not at ‘0’ phase response. More precisely, starting with the stroboscopic map:ϕnew\ = (ϕold\ +PRC(ϕold)+ TzeitTosc2π) mod 2π , a fixed point requires that ϕnew = ϕold  ϕent.

For 1:1 entrainment we then have:

PRC(ϕent)+ TzeitTosc2π = 2π,

PRC(ϕent) =2πTosc(ToscTzeit) - period mismatch with the zeitgeber, measured in units of phase.

Note however on Figure 8G that we further observe a vertical shift of the PRC, which prompted us to propose a change of the intrinsic period with Tzeit (as explained in the text when we introduce Figures8A1-2).

Another way to visualize fixed points is offered in Figure 16 D-E, where we plot the inferred corrected PTC and the stroboscopic maps: there, fixed points correspond to intersections with the diagonal.

Reviewer #3 (Significance (Required)):

Although the phase-analysis has been widely applied to various biological systems, such as circadian clocks, cardiac tissues and neurons, this paper represents the first detailed experimental analysis of the segmentation clock based on the theory of phase dynamics. The major results are inline with theoretical predictions, whereas the suggestion about the SNIC bifurcation is attractive not only to the theoretical researchers but also to the experimental biologists; it has been believed that the segmentation clock consists of negative-feedback oscillator that emerge by Hopf bifurcation, whereas this paper proposes another possibility of the molecular network structure for the clockwork. This issue is related to recently proposed hypothesis about the excitable system in the segmentation clock based on the Yap signaling (Hubaud et al. Cell 171, 668 (2017)). However, unfortunately, discussion about detailed molecular networks are not abundant.

42. Thus, maybe the main readers are computational biologists and systems biologists.

We thank the reviewer for his/her significance comment. We have added comments on the bifurcation structure of the segmentation clock and on excitable systems in the discussion. While our focus is on coarse-graining so that we do not and cannot infer precise molecular details, we can still infer some properties of the underlying networks. In particular we now cite several papers explaining how systems with tunable periods/excitable are indicative of the interplay between positive and negative feedbacks. We think those considerations are of interest to a broad range of biologists interested in connecting experiments to theory.

https://doi.org/10.7554/eLife.79575.sa2

Article and author information

Author details

  1. Paul Gerald Layague Sanchez

    European Molecular Biology Laboratory (EMBL), Developmental Biology Unit, Heidelberg, Germany
    Present address
    Department of Genetics and Evolution, University of Geneva, Geneva, Switzerland
    Contribution
    Conceptualization, Formal analysis, Investigation, Writing – original draft
    Contributed equally with
    Victoria Mochulska
    Competing interests
    No competing interests declared
    ORCID icon "This ORCID iD identifies the author of this article:" 0000-0001-6213-8927
  2. Victoria Mochulska

    McGill University, Montreal, Canada
    Contribution
    Formal analysis, Investigation, Methodology, Writing – review and editing
    Contributed equally with
    Paul Gerald Layague Sanchez
    Competing interests
    No competing interests declared
    ORCID icon "This ORCID iD identifies the author of this article:" 0000-0002-1765-6132
  3. Christian Mauffette Denis

    McGill University, Montreal, Canada
    Contribution
    Formal analysis, Investigation
    Competing interests
    No competing interests declared
    ORCID icon "This ORCID iD identifies the author of this article:" 0000-0003-1763-8249
  4. Gregor Mönke

    European Molecular Biology Laboratory (EMBL), Developmental Biology Unit, Heidelberg, Germany
    Present address
    Ernst Strüngmann Institute for Neuroscience, Frankfurt, Germany
    Contribution
    Conceptualization, Methodology
    Competing interests
    No competing interests declared
    ORCID icon "This ORCID iD identifies the author of this article:" 0000-0002-3521-715X
  5. Takehito Tomita

    European Molecular Biology Laboratory (EMBL), Developmental Biology Unit, Heidelberg, Germany
    Contribution
    Formal analysis, Methodology
    Competing interests
    No competing interests declared
  6. Nobuko Tsuchida-Straeten

    European Molecular Biology Laboratory (EMBL), Developmental Biology Unit, Heidelberg, Germany
    Present address
    Universitätsklinikum Heidelberg, Heidelberg, Germany
    Contribution
    Investigation
    Competing interests
    No competing interests declared
  7. Yvonne Petersen

    European Molecular Biology Laboratory (EMBL), Transgenic Service, Heidelberg, Germany
    Contribution
    Methodology
    Competing interests
    No competing interests declared
  8. Katharina Sonnen

    European Molecular Biology Laboratory (EMBL), Developmental Biology Unit, Heidelberg, Germany
    Present address
    Hubrecht Institute, Utrecht, Netherlands
    Contribution
    Conceptualization, Investigation, Methodology
    Competing interests
    No competing interests declared
    ORCID icon "This ORCID iD identifies the author of this article:" 0000-0002-2902-9419
  9. Paul François

    McGill University, Montreal, Canada
    Contribution
    Conceptualization, Formal analysis, Supervision, Writing – original draft
    For correspondence
    paul.francois2@mcgill.ca
    Competing interests
    No competing interests declared
    ORCID icon "This ORCID iD identifies the author of this article:" 0000-0002-2223-839X
  10. Alexander Aulehla

    European Molecular Biology Laboratory (EMBL), Developmental Biology Unit, Heidelberg, Germany
    Contribution
    Conceptualization, Resources, Supervision, Funding acquisition, Investigation, Writing – original draft
    For correspondence
    aulehla@embl.de
    Competing interests
    No competing interests declared
    ORCID icon "This ORCID iD identifies the author of this article:" 0000-0003-3487-9239

Funding

European Molecular Biology Laboratory

  • Alexander Aulehla

H2020 European Research Council (866537)

  • Alexander Aulehla

Deutsche Forschungsgemeinschaft (331351713)

  • Alexander Aulehla

Natural Sciences and Engineering Research Council of Canada (NFRFE-2020-00703)

  • Paul François

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

Acknowledgements

This work is presented, in part, in the doctoral thesis of Paul Gerald Layague Sanchez (https://doi.org/10.11588/heidok.00029209) submitted to the University of Heidelberg, Germany. We thank Mogens Høgh Jensen, Jonas Juul, Mathias Heltberg, Sandeep Krishna, Adrián Granada, and Istvan Kiss for in-depth discussions on Arnold tongues and PRC; Justin Crocker, Anne Ephrussi, Lars Hufnagel, Jelle Scholtalbers, Ulrich Schwarz, Stefano Vianello, Nachiket Satish Shembekar, and Ramesh Utharala, and all the members of the Aulehla group for the discussions and feedback. We are very grateful to Stefano Vianello for also designing some of the figures. We thank Antonio Politi for providing the Pipeline Constructor macro for automated live microscopy. We also acknowledge the assistance from Maximilian Beckers, Aliaksandr Halavatyi, and Oleksandr Maistrenko in writing scripts for data analysis. We are thankful to Ricardo Henriques for kindly sharing the template that was used to format this preprint. Following EMBL core facilities are acknowledged for their support of this work: the EMBL Advanced Light Microscopy Facility (ALMF), the EMBL Laboratory Animal Resources (LAR), the EMBL Mechanical Design Office, and the EMBL Mechanical Workshop. The work was supported by the European Molecular Biology Laboratory (EMBL), by the German Research Foundation/DFG (CRC/SFB 1324, project number 331351713) to AA, the National Research Council of Canada, Discovery Grant program to Paul François, by the New Frontiers in Research Fund, Exploration program to PF and AA and the European Research Council (ERC) via an ERC Starting Grant (agreement no. 633943) and Consolidator Grant (agreement no. 866537) to Alexander Aulehla.

Ethics

All animal experiments were conducted under veterinarian supervision and after project approval by European Molecular Biology Laboratory, following the guidelines of the European Commission, Directive 2010/63/EU and AVMA Guidelines 2007.

Senior Editor

  1. Aleksandra M Walczak, CNRS LPENS, France

Reviewing Editor

  1. Sandeep Krishna, National Centre for Biological Sciences­‐Tata Institute of Fundamental Research, India

Publication history

  1. Preprint posted: October 21, 2021 (view preprint)
  2. Received: April 19, 2022
  3. Accepted: August 23, 2022
  4. Version of Record published: October 12, 2022 (version 1)

Copyright

© 2022, Sanchez, Mochulska et al.

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

Metrics

  • 467
    Page views
  • 70
    Downloads
  • 0
    Citations

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

Download links

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

Downloads (link to download the article as PDF)

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

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

  1. Paul Gerald Layague Sanchez
  2. Victoria Mochulska
  3. Christian Mauffette Denis
  4. Gregor Mönke
  5. Takehito Tomita
  6. Nobuko Tsuchida-Straeten
  7. Yvonne Petersen
  8. Katharina Sonnen
  9. Paul François
  10. Alexander Aulehla
(2022)
Arnold tongue entrainment reveals dynamical principles of the embryonic segmentation clock
eLife 11:e79575.
https://doi.org/10.7554/eLife.79575

Further reading

    1. Physics of Living Systems
    Rakesh Das, Takahiro Sakaue ... Tetsuya Hiraiwa
    Research Article

    Spatial organization of chromatin plays a critical role in genome regulation. Previously, various types of affnity mediators and enzymes have been attributed to regulate spatial organization of chromatin from a thermodynamics perspective. However, at the mechanistic level, enzymes act in their unique ways and perturb the chromatin. Here, we construct a polymer physics model following the mechanistic scheme of Topoisomerase-II, an enzyme resolving topological constraints of chromatin, and investigate how it affects interphase chromatin organization. Our computer simulations demonstrate Topoisomerase-II's ability to phase separate chromatin into eu- and heterochromatic regions with a characteristic wall-like organization of the euchromatic regions. We realized that the ability of the euchromatic regions to cross each other due to enzymatic activity of Topoisomerase-II induces this phase separation. This realization is based on the physical fact that partial absence of self-avoiding interaction can induce phase separation of a system into its self-avoiding and non-self-avoiding parts, which we reveal using a mean-field argument. Furthermore, motivated from recent experimental observations, we extend our model to a bidisperse setting and show that the characteristic features of the enzymatic activity driven phase separation survive there. The existence of these robust characteristic features, even under the non-localized action of the enzyme, highlights the critical role of enzymatic activity in chromatin organization.

    1. Cell Biology
    2. Physics of Living Systems
    Nicola Bellotto, Jaime Agudo-Canalejo ... Victor Sourjik
    Research Article

    Inside prokaryotic cells, passive translational diffusion typically limits the rates with which cytoplasmic proteins can reach their locations. Diffusion is thus fundamental to most cellular processes, but the understanding of protein mobility in the highly crowded and non-homogeneous environment of a bacterial cell is still limited. Here we investigated the mobility of a large set of proteins in the cytoplasm of Escherichia coli, by employing fluorescence correlation spectroscopy (FCS) combined with simulations and theoretical modeling. We conclude that cytoplasmic protein mobility could be well described by Brownian diffusion in the confined geometry of the bacterial cell and at the high viscosity imposed by macromolecular crowding. We observed similar size dependence of protein diffusion for the majority of tested proteins, whether native or foreign to E. coli. For the faster-diffusing proteins, this size dependence is well consistent with the Stokes-Einstein relation once taking into account the specific dumbbell shape of protein fusions. Pronounced subdiffusion and hindered mobility are only observed for proteins with extensive interactions within the cytoplasm. Finally, while protein diffusion becomes markedly faster in actively growing cells, at high temperature, or upon treatment with rifampicin, and slower at high osmolarity, all of these perturbations affect proteins of different sizes in the same proportions, which could thus be described as changes of a well-defined cytoplasmic viscosity.