Arnold tongue entrainment reveals dynamical principles of the embryonic segmentation clock
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 coarsegrained approach, in which theory guides the experimental identification of the segmentation clock entrainment responses. We demonstrate period and phaselocking of the segmentation clock across a wide range of entrainment parameters, including higherorder 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 nonlinear oscillator close to an infinite period bifurcation and suggests the presence of longterm feedbacks. Combined, this coarsegrained theoreticalexperimental 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 coarsegrained 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 coarsegrained 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.sa0Introduction
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 coarsegrain degrees of freedom to obtain scalefree 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; ProulxGiraldeau 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, macrolevel 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 coarsegrained 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; DiazCuadros et al., 2020; Matsuda et al., 2020).
Strikingly, as individual oscillating cells are coupled to their neighbors via NotchDelta signaling, the oscillations occur synchronized and wavelike activity patterns appear to periodically sweep along the embryonic anteriorposterior axis (Maroto et al., 2005; Masamizu et al., 2006; Aulehla et al., 2008; Herrgen et al., 2010; Webb et al., 2016; YoshiokaKobayashi 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; ElSherif et al., 2012; Brena and Akam, 2013; Schönauer et al., 2016).
In this work, we coarsegrain these underlying complexities and take a dynamical systems, macroperspective 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.
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 ${\varphi}_{ent}$ at the time of the zeitgeber is constant (technically, a fixed point of the Poincaré return map [Cross and Siggia, 2005]). This defines periodlocking (also termed modelocking) (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, ${T}_{zeit}{T}_{osc}$) between the zeitgeber and the freerunning 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 ($\u03f5$) versus the corresponding detuning (or simply ${T}_{zeit}$, if ${T}_{osc}$ is constant): these maps are more commonly known as Arnold tongues (Figure 1E). Arnold tongues predict the period and phaselocking 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$ periodlocking. In that case, the instantaneous period of the oscillator matches $m/n$ the period of the zeitgeber (${T}_{zeit}$). 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 shockdriven expression of the ligand Delta was used to modulate the segmentation of PSM in zebrafish (SozaRied 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 microfluidicsbased 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 quasitwodimensional in vitro segmentation assay (hereafter referred to as a 2Dassay), which recapitulates segmentation clock dynamics and PSM patterning (Lauschke et al., 2013), the microfluidicsentrainment 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 coarsegrained, 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.
Hence, we conclude that LuVeLu global ROI timeseries analysis provides a valid coarsegrained 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, coarsegrained readout for segmentation clock pace and rhythm, we next analyzed whether pace and rhythm can be experimentally tuned using microfluidicsbased 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 freerunning, natural period or ${T}_{osc}$. 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 ${T}_{osc}$ (Figure 3A–C, Figure 3—figure supplement 1), the segmentation clock rhythm in DAPTentrained samples closely adjusted to ${T}_{zeit}$ 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.
Notably, periodlocking 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.
We also tested the effect of varying DAPT concentration, as means of changing zeitgeber strength ($\u03f5$ in Figure 1), on entrainment dynamics. Synchronization theory predicts that zeitgeber strength correlates with the time it takes to reach periodlocking (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 periodlocking 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 inphase 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 RiedelKruse et al., 2007).
Higherorder entrainment of the segmentation clock
In theory, a nonlinear oscillator should be amenable not only to 1:1 entrainment, but also to higherorder $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 higherorder 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 higherorder 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 ${T}_{zeit}=300$ min and approaching 175 min for ${T}_{zeit}=350$ min (Figure 5C). We also found evidence for phaselocking during 2:1 entrainment, that is, a narrowing of the phase distribution among the individual samples (Figure 5D).
The segmentation clock establishes a stable entrainment phase relative to the zeitgeber
We next analyzed the phaselocking behavior between the segmentation clock and the zeitgeber. According to dynamical systems theory, phaselocking, that is, the entrainment phase (${\varphi}_{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 ${\varphi}_{ent}$ as a stable fixed point that lies on the diagonal, where there is phaselocking. 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 ${\varphi}_{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 (${\varphi}_{init}$) was closer (or farther) to ${\varphi}_{ent}$.
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 $\pi $) within the range of permissible zeitgeber entrainment periods, ${T}_{zeit}$ (Aschoff and Pohl, 1978; Granada et al., 2013).
To test these theoretical predictions, we quantified ${\varphi}_{ent}$ across a wide range of detuning, that is, from 120 to 180 min zeitgeber period. We found that ${\varphi}_{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, ${\varphi}_{ent}$ systematically shifted from $\sim \pi /2$ to $\sim 3\pi /2$ (Figure 7A–B), spanning a range of almost $\pi $ (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.
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 2Dassays (Tsiairis and Aulehla, 2016; SozaRied et al., 2014). We were then curious whether or not the dynamics from our coarsegrained analysis match that for oscillations in the center of the 2Dassays, 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 2Dassays (‘center ROI’, radius: 25 pixels, 1 pixel = 1.38 µm). Upon comparison, we noted that the oscillations at the center of 2Dassays adjusted its rhythm to the zeitgeber with strong correspondence to the coarsegrained 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 phaselocking 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 2Dassays (i.e. corresponding to the posterior PSM), the entrainment phase at ${T}_{zeit}=130$ min is different to the entrainment phase at ${T}_{zeit}=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 $\u03f5$ (Pikovsky et al., 2001) and write:
where ${\varphi}_{\text{pulse}}$ is the phase of the segmentation clock on the cycle at the moment of the pulse, and ${\varphi}_{\text{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 $\u03f5$ 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 $\u03f5$ 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 ${T}_{zeit}$. The sequence of phases at each pulse is then given by the stroboscopic map, introduced in Figure 6A:
and 1:1 entrainment occurs when this stroboscopic map converges toward a single fixed point ${\varphi}_{ent}$ for a given ${T}_{zeit}$. When ${T}_{zeit}$ is varied, different phases of entrainment ${\varphi}_{ent}({T}_{zeit})$ are observed, here plotted in Figure 7B.
To derive the PRC directly from the experimental stroboscopic maps, we invert Equation 2 into:
By estimating ${\varphi}_{n},{\varphi}_{n+1}$ for a given ${T}_{zeit}$, 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 ${\varphi}_{n}$ , even far from the entrainment phase ${\varphi}_{ent}$. Notice, however, that phases far from ${\varphi}_{ent}$ will be only sampled over the first few pulses (so with possibly much noise) while conversely, ${\varphi}_{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 ${T}_{zeit}$ is changed.
Such vertical shifts in PRCs as a function of entrainment period have been previously observed in socalled ‘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 ${T}_{osc}$). 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 ${\varphi}_{ent}({T}_{zeit})$ 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 ${T}_{zeit}=140$ min to be the ‘reference’ curve, since this ${T}_{zeit}$ 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, periodindependent PRC in Figure 8A2 (see Appendix 1, section Theoretical methods). From there, we build upon the simplest nonlinear phaseoscillator, that is, the classical Radial Isochron Cycle (RIC). To modulate its sineform 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 $\lambda $). Both high values of $\lambda $ 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 $\u03f5=0.4$, high value of $\lambda =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, coarsegrained 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 coarsegraining, entrainment approach to gain new insights into the properties of the segmentation clock from dynamical systems’ perspective. We modelock the segmentation clock to various entrainment periods and use the information about the dynamic phaselocking behavior to derive the somite segmentation clock PRC from the experimental data.
A coarsegraining 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 systemslevel single oscillator phase description, we observe a consistent entrainment behavior, that is, periodlocking with convergence toward a welldefined entrainment phase, as predicted from the oscillator phase response theory. It is important to point out that theoretical predictions are nontrivial 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 coarsegraining 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 constantsign 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 saddlenode 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 relaxationtype 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 $\mathrm{\infty}$ to $+\mathrm{\infty}$ in a finite time before resetting to $\mathrm{\infty}$ (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, constantsign 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 relaxationtype 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 longterm 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 systemlevel change in segmentation clock rate is not understood (see also below for macro versus microscale 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 (JutrasDubé et al., 2020), reflecting network modularity. Such modularity and robustness with changes in regulatory logic is a hallmark of developmental plasticity (WestEberhard, 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 higherorder 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 phaselocked, 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 macroscale behavior relates to the underlying cellular scale oscillations.
Conclusions
Our work demonstrates how, despite all the molecular and functional complexities, coarsegraining 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 coarsegrained level.
We also aim to illustrate the potential of an integrated, theoreticalexperimental 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 experimentalobservational 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 macroscale, 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’ theoreticalexperimental 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
Materials and methods
Mouse lines
Request a detailed protocolFor 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 Mesp2GFP mouse line generated by Morimoto et al., 2006. Axin2linkerAchilles, a mouse line expressing Axin2 fused with a GSAGS linker to Achilles, a fastmaturing variant of YFP (YoshiokaKobayashi et al., 2020), was generated inhouse. To generate the knockin 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: linkerAchillesloxPPGK NeoloxP. The selection cassette was flanked by loxP sites for eventual Cremediated excision. Axin2linkerAchilles knockin 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 (RIKENCBS) 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 protocolOn 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% CO_{2} until use.
Mouse dissection and embryo recovery
Request a detailed protocolFemale 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 reporterofinterest (e.g. LuVeLu) using a stereo fluorescence microscope. The tailbud was cut from the rest of the tail and was immediately transferred to preequilibriated 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 microfluidicsbased entrainment of 2Dassays
We used in this study a microfluidicsbased 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 microfluidicsbased 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 twentyfour 1 cm plugs made from cut PDMSfilled PTFE tubing for the sample inlets and unused drug/medium inlets. These were all sterilized under UV for at least 20 min.
Fibronectincoating of PDMS device and overnight soaking in buffer
Request a detailed protocolWhile waiting, 5 mL of PenStrep (Gibco, 15140122) 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. PDMSfilled PTFE plugs were also immersed in PBS + PenStrep. To prepare the fibronectin solution, 280 µL of fibronectin (SigmaAldrich, 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 LuerLok REF 309658). A UVirradiated needle, which was earlier inserted into a 1meter 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 protocolFor a microfluidicsbased experiment with periodic pulses of drug, four 10 mL syringes (diameter: 14.5 mm, BD LuerLok 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 (SigmaAldrich, D59425MG). To prepare 10 mM stock of DAPT, 5 mg DAPT (MW = 432.46 g/mol) was dissolved in 1156.2 µL DMSO (SigmaAldrich, D8418). This was aliquoted and stored at –20°C until use.
Degassing drug/medium and PDMS device
Request a detailed protocolAfter 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 protocolBefore mouse dissection and recovery of mouse embryos, syringes containing the degassed drug/medium were each connected to a UVirradiated 3 m PTFE tubing via attached syringe needle, and were carefully mounted on programmable pumps (World Precision Instruments, AL400) 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% CO_{2}. After recovery of embryonic tissues, using a pipette (i.e. P200 for 2Dassays 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 PDMSfilled 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 lintfree wipes (KimberlyClark). The PDMS device, with more than 1 m of attached tubings, was carefully placed inside a preequilibriated 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 Ushaped 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.
Setting up automated pumping
Request a detailed protocolEntrainment 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, AL400). Diameter of syringe (14.5 mm for 10 mL syringe, BD LuerLok 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 ${T}_{zeit}$ = 170 min.
Confocal microscopy
Request a detailed protocolFor most of the experiments here, samples were imaged in an LSM 780 laserscanning microscope (Carl Zeiss Microscopy) fitted with an incubation chamber (EMBL Mechanical Workshop). A PlanApochromat ×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 zstacks (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 zstack, 32 pixels × 32 pixels (pixel size = 22.14 µm). Imaging and automated pumping of drug/medium through the microfluidics device were started simultaneously.
Data analysis
Details on the coarsegraining strategy and analyses of entrainment dynamics are specified here.
Extracting timeseries from global intensity of 2Dassays for subsequent analyses
Request a detailed protocolTo 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). Zstacks were first projected based on maximum intensity (in Fiji: Image > Stacks > Z Project). Subsequently, the timeseries was obtained by plotting the zaxis profile (in Fiji: Image > Stacks > Plot Zaxis 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 2Dassays, 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 reoriented (and registered, if necessary) 2Dassay timelapse file (in Fiji: Edit > Selection > Specify). Timeseries was extracted after specifying the center ROI (in Fiji: Image > Stacks > Plot Zaxis 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
(see also RiedelKruse et al., 2007), the first Kuramoto order parameter is defined by
Monitoring periodlocking and phaselocking
Request a detailed protocolEntrainment was evaluated based on periodlocking and phaselocking 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 Pythonbased 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). Timeresolved frequency analysis was done by crosscorrelating 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 ($\mathrm{\Delta}\varphi $) was defined as
where $t$ is the time of perturbation and ${T}_{zeit}$ is the zeitgeber period (i.e. one cycle after time = $t$). $\varphi (t)$ and $\varphi (t+{T}_{zeit})$ 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 (${\varphi}_{ent}$), was determined considering only phaselocked 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 $\pi /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:
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), scikitimage (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 pvalues, twotailed test for absolute difference between medians was done via a randomization method using PlotsOfDifferences (Goedhart, 2019). While not recommended, pvalues 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 filterbased 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 2Dassays
Reoriented 2Dassays (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 2Dassays. 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, reoriented (and registered, if necessary) timelapse movies of 2Dassays 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 ${T}_{zeit}=(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 ${T}_{zeit}=(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 3DE.
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 $\u03f5$ (Appendix 1—figure 3B).
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 $\lambda $, so that the corresponding equation of the ellipse is
combined with radial equidistant isochrons, $\dot{\theta}=1$, where $\theta $ 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:
where ${\dot{r}}_{cycle}=\lambda \frac{rx}{{(r+\lambda y)}^{2}}$ (again assuming $\dot{\theta}=1$). This leads to the following equations in cartesian coordinates:
$\dot{y}=x+\left({\displaystyle \frac{r}{r+\lambda y}}r{\displaystyle \frac{\lambda x}{{(r+\lambda y)}^{2}}}\right).y=d{y}_{\lambda}(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 $\u03f5$ toward the right can be computed in a straightforward way and is:
As said above, the effect of introducing $\lambda $ is to smoothly flatten the PRC over half of the cycle, as illustrated in Appendix 1—figure 4A.
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
This will keep isochrons radial but will change their spacing. With the elliptic limit cycle, this leads to the differential equation in cartesian coordinates
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}_{\ast}>1$ for one sector and $s=1$ otherwise. The sped up sector is centered at angle $\alpha $ and has width $\beta $, so that the modified period of the cycle is ${T}_{{s}_{*}}=\beta /{s}_{*}+(2\pi \beta )$. It is also convenient to define the rescaled angular velocity ${\omega}_{{s}_{*}}=2\pi /{T}_{{s}_{*}}=\frac{2\pi {s}_{*}}{\beta +(2\pi \beta ){s}_{*}}$. From there, the phase of the cycle as a function of the angle $\theta $ in the plane is the simple linear transformation (defining ${\theta}_{0}=\alpha \beta /2$):
for $0<\theta <{\theta}_{0}$,
for ${\theta}_{0}<\theta <{\theta}_{0}+\beta$, and
for ${\theta}_{0}+\beta <\theta <2\pi$. Those functions ensure that the rate of phase evolution in sector ${\theta}_{0}<\theta <{\theta}_{0}+\beta$ is $1/{s}_{*}$ times the rate in the other sectors (compare $\theta $ coefficients in Equations 1416), that angle 0 in $\theta $ is phase ${\varphi}_{s}=0$, and that phase is continuous so that ${\varphi}_{s}({\theta}_{0})={\omega}_{{s}_{*}}{\theta}_{0}$ and ${\varphi}_{s}({\theta}_{0}+\beta )=({\theta}_{0}+\beta /{s}_{*}){\omega}_{{s}_{*}}$. Notice also that for $\theta =2\pi $ we get from Equation 16 ${\varphi}_{s}(2\pi )=\left(\frac{\beta}{{s}_{*}}+2\pi \beta \right)\frac{2\pi {s}_{*}}{\beta +(2\pi \beta ){s}_{*}}=2\pi $ 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 $\lambda ,{s}_{*}$, $\alpha $, and $\beta $ 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 ${\theta}_{i}$ linearly spaced on the interval $[0,2\pi ]$. We then compute the corresponding position on the ERIC limit cycle, and compute numerically the angle ${\theta}_{i}^{\to \u03f5}$ of the point at distance $\u03f5$ on the right. By definition of ERICA, the PRC for each index $i$ is then
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 $\u03f5$, $\lambda $, ${s}_{*}$, sped up sector location $\alpha $ and width $\beta $, and the location of the zerophase reference point on the cycle ${\varphi}_{0}$. We minimize ${\chi}^{2}$:
where $PR{C}_{i}^{corr}$ are data points of the corrected experimental PRC and $\sigma $ 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: $\u03f5=0.43$, $\lambda =0.53$, ${s}_{*}=5.64$, $\alpha =1.51\pi $, $\beta =1.15\pi $, ${\varphi}_{0}=1.17\pi $.
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 ${T}_{osc}$ changes linearly with ${T}_{zeit}$.
To see this, we start with the classical condition for stability of the entrainment phase (see, e.g., Izhikevich, 2006),
or equivalently:
Now we also have by definition of the entrainment phase ${\varphi}_{ent}$:
where the factor $\frac{2\pi}{{T}_{osc}}$ is the unit conversion from time to phase. This equation defines in an implicit way a function ${\varphi}_{ent}({T}_{zeit})$, 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:
This allows to implicitly define the curve ${\varphi}_{ent}({T}_{zeit})$ with the help of the $PRC$:
There is a geometric interpretation here: the curve ${\varphi}_{ent}({T}_{zeit})$ is in fact a (rescaled) $+\pi /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 $\varphi $ and $PRC(\varphi ,PRC(\varphi ))\to (PRC(\varphi ),\varphi )$). 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\pi /4=\pi /2$). For instance, if the PRC is sinusoidal, ${\varphi}_{ent}({T}_{zeit})$ 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 $\pi /2$ is compared with the curve ${\varphi}_{ent}({T}_{zeit})$ 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 ${\varphi}_{ent}({T}_{zeit})$ , for example, combining Equations 2023 and counting time in units of the intrinsic period ${T}_{osc}$, $2\pi /{T}_{osc}=1,$ we get
This defines an absolute, minimum slope of the curve ${\varphi}_{ent}({T}_{zeit})$. 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 ${T}_{osc}$ in fact depends on the zeitgeber. We then get a generalized version of Equation 22 with changing intrinsic period
giving the new implicit definition
Now we see that when $\frac{d{T}_{osc}}{d{T}_{zeit}}\frac{{T}_{zeit}}{{T}_{osc}}\sim 1$, the righthand side of Equation 26 can be very small, so that one can get $\frac{d{\varphi}_{ent}}{d{T}_{zeit}}\sim 0$, 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 $\frac{d{T}_{osc}}{d{T}_{zeit}}\frac{{T}_{zeit}}{{T}_{osc}}=\frac{d\mathrm{ln}{T}_{osc}}{d\mathrm{ln}{T}_{zeit}}\sim 1$ could be indicative of a mechanism where the oscillator adapts its intrinsic period (here to the zeitgeber period).
Practically, we obtain the ${T}_{osc}({T}_{zeit})$ 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 3DE). The values of ${T}_{osc}$ 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 $\sim 120$ min and above $\sim 200$ min which indicates that $\frac{d{T}_{osc}}{d{T}_{zeit}}$ 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 cutoff for ${T}_{osc}$ 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 (MorenoRisueno et al., 2010). With such cutoff, 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 $\u03f5=0.43$. Assuming that the DAPT concentration is reflected in the strength of perturbation, we can find crosssections of the isophases at different values of $\u03f5$ to have an excellent agreement with experiments, as illustrated in Figure 8—figure supplement 3G. We take $\u03f5=0.55$ for 3.0 µM, $\u03f5=0.31$ for 1.0 µM, and $\u03f5=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 ${\varphi}_{ent}({T}_{zeit})$ 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/EMBLfiles).
References

Collective modes of coupled phase oscillators with delayed couplingPhysical Review Letters 108:204101.https://doi.org/10.1103/PhysRevLett.108.204101

Phase relations between a circadian rhythm and its zeitgeber within the range of entrainmentDie Naturwissenschaften 65:80–84.https://doi.org/10.1007/BF00440545

BookSynchronization: From Simple to ComplexSpringer Berlin Heidelberg. GoogleBooksID y6_UDAEACAAJ.

Tuning the phase of circadian entrainmentJournal of the Royal Society, Interface 12:20150282.https://doi.org/10.1098/rsif.2015.0282

BookPeriodic perturbation of the BZreaction in a CSTR: chemical resonance, entrainment and quasiperiodic behaviorIn: Rensing L, Jaeger NI, editors. Temporal Order Springer Series in Synergetics. Springer. pp. 116–121.https://doi.org/10.1007/9783642703324

BookScaling and Renormalization in Statistical PhysicsCambridge: Cambridge University Press.https://doi.org/10.1017/CBO9781316036440

Mode locking the cell cyclePhysical Review. E, Statistical, Nonlinear, and Soft Matter Physics 72:021910.https://doi.org/10.1103/PhysRevE.72.021910

Segmental patterning of the vertebrate embryonic axisNature Reviews. Genetics 9:370–382.https://doi.org/10.1038/nrg2320

Methods in Enzymology1–27, Chapter 1 Phase Response Curves: Elucidating the Dynamics of Coupled Oscillators, Methods in Enzymology, Academic Press, 10.1016/S00766879(08)038019.

Intercellular coupling regulates the period of the segmentation clockCurrent Biology 20:1244–1253.https://doi.org/10.1016/j.cub.2010.06.034

Numerical parameter space compression and its application to biophysical modelsBiophysical Journal 118:1455–1465.https://doi.org/10.1016/j.bpj.2020.01.023

Signalling dynamics in vertebrate segmentationNature Reviews. Molecular Cell Biology 15:709–721.https://doi.org/10.1038/nrm3891

Matplotlib: a 2d graphics environmentComputing in Science & Engineering 9:90–95.https://doi.org/10.1109/MCSE.2007.55

Perturbation theory of selfoscillating system with a periodic perturbationProgress of Theoretical Physics 61:45–53.https://doi.org/10.1143/PTP.61.45

BookDynamical Systems in Neuroscience: The Geometry of Excitability and BurstingCambridge, MA, USA: MIT Press.https://doi.org/10.7551/mitpress/2526.001.0001

Stroboscopic phase portrait of a forced nonlinear oscillatorProgress of Theoretical Physics 61:54–73.https://doi.org/10.1143/PTP.61.54

Adjoint method provides phase response functions for delayinduced oscillationsPhysical Review Letters 109:044101.https://doi.org/10.1103/PhysRevLett.109.044101

Evolutionary plasticity of segmentation clock networksDevelopment 138:2783–2792.https://doi.org/10.1242/dev.063834

Overdrive suppression of spontaneously beating chick heart cell aggregates: experiment and theoryThe American Journal of Physiology 269:H1153–H1164.https://doi.org/10.1152/ajpheart.1995.269.3.H1153

BookChemical Oscillations, Waves, and TurbulenceCourier Corporation. GoogleBooksID.

Synchronised cycling gene oscillations in presomitic mesoderm cells require cellcell contactThe International Journal of Developmental Biology 49:309–315.https://doi.org/10.1387/ijdb.041958mm

ConferenceData Structures for Statistical Computing in PythonPython in Science Conference.https://doi.org/10.25080/Majora92bf192200a

Theory in biology: figure 1 or figure 7?Trends in Cell Biology 25:723–729.https://doi.org/10.1016/j.tcb.2015.10.007

BookV. An oscillator model for biological clocksIn: Rudnick D, editors. Rhythmic and Synthetic Processes in Growth. Princeton University Press. pp. 75–110.https://doi.org/10.1515/9781400876167006

Untangling the hairball: fitnessbased asymptotic reduction of biological networksBiophysical Journal 113:1893–1906.https://doi.org/10.1016/j.bpj.2017.08.036

A circadian surface of entrainment: varying T, τ, and photoperiod in Neurospora crassaJournal of Biological Rhythms 25:318–328.https://doi.org/10.1177/0748730410379081

Fiji: an opensource platform for biologicalimage analysisNature Methods 9:676–682.https://doi.org/10.1038/nmeth.2019

ConferenceThe NumPy Array: A Structure for Efficient Numerical ComputationComputing in Science & Engineering. pp. 22–30.

A microfluidics approach for the functional investigation of signaling oscillations governing somitogenesisJournal of Visualized Experiments 1:e2318.https://doi.org/10.3791/62318

SoftwareSeabornV0.5.

BookDevelopmental Plasticity and EvolutionNY: Oxford University Press.https://doi.org/10.1093/oso/9780195122343.001.0001

Biological rhythms and the behavior of populations of coupled oscillatorsJournal of Theoretical Biology 16:15–42.https://doi.org/10.1016/00225193(67)900513
Decision letter

Sandeep KrishnaReviewing Editor; National Centre for Biological Sciences‐Tata Institute of Fundamental Research, India

Aleksandra M WalczakSenior 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.sa1Author 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(\theta )=\text{}{\theta}_{\text{new}}\theta$
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 quasiperiodic regime since all the points in the circle will be sampled. Since no quasiperiodic response was demonstrated, the claim of entrainment is not convincing.
While, in principle, PRC can be indeed obtained from responses in the “quasiperiodic” 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.
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 S6S7 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. perioddoubling, higher harmonics (Figure 5), quasiperiodicity, 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: 1D oscillators with resetting (such as the quadratic fireandintegrate 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 theorydriven 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 coarsegraining 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 counterintuitive 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 redone, 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 coarsegraining 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 coarsegraining approach. While the results might of course be specific to the segmentation process, we think our approach motivated by coarsegraining 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 nontrivial 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 RiedelKruse et al., 2007 who used a meanfield 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 highlevel comparison to other synchronizing systems.
It is interesting the reviewer sees similarities with the work of RiedelKruse et al., which uses a meanfield 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 realtime 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 saddlenode on invariant cycle (SNIC); it also reveals that the intrinsic period changes in a nonlinear 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 nonlinear 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 2Dassay (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 2Dassay 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 coarsegrained 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 466470 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 anteriorposterior, pillars etc; it is also difficult to view the brightfield since it is presented as a merged image.
It is particularly difficult to see the somite boundaries for the same reason. In lines 113117 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 Wntsignaling oscillations (Axin2Achilles) and the segmentation marker Mesp2 (Mesp2GFP) 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 reoriented 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 Axin2GSAGSAchilles or Mesp2GFP 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 ${\varphi}_{\mathrm{e}\mathrm{n}\mathrm{t}\text{}}({T}_{\mathrm{z}\mathrm{e}\mathrm{i}\mathrm{t}})$ 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 sincfilter 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 periodlocking and phaselocking: 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 ${T}_{\mathrm{z}\mathrm{e}\mathrm{i}\mathrm{t}}:$${\dot{\varphi}}_{\mathrm{z}\mathrm{e}\mathrm{i}\mathrm{t}}=2\pi /{T}_{\mathrm{z}\mathrm{e}\mathrm{i}\mathrm{t}}.$ The initial condition for $\varphi}_{\mathrm{z}\mathrm{e}\mathrm{i}\mathrm{t}$ is chosen so that the zeitgeber phase at the moment of last pulse is matching the experimental entrainment phase for each $T}_{\mathrm{z}\mathrm{e}\mathrm{i}\mathrm{t}$. We plot $cos\text{}\varphi (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: ${R}_{j}(t)\text{}=\text{}\frac{1}{2}({e}^{i{\varphi}_{j}(t)}+{e}^{i{\varphi}_{\mathrm{z}\mathrm{e}\mathrm{i}\mathrm{t}}(t)})\text{}$. By the end of experiment most samples reach ${R}_{j}\text{}\approx 1$, indicating entrainment. Most samples need $\sim 3$ zeitgeber cycles to become entrained. For ${T}_{\mathrm{z}\mathrm{e}\mathrm{i}\mathrm{t}}=120$ min the entrainment takes much longer (edge of the Arnold tongue). For ${T}_{\text{zeit}}=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.
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 (crosscorrelation 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), periodlocking (where the pace of the samples match the pace of the zeitgeber, which can be distinct from natural pace), and phaselocking (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.
pvalues were calculated for the specified conditions and were added in the caption of the figures. These values are enumerated here:
– Figure 3C
170min 2uM DAPT (vs DMSO control): p < 0.001
– Figure S2
120min 2uM DAPT (vs DMSO control): p = 0.064
130min 2uM DAPT (vs DMSO control): p = 0.003
140min 2uM DAPT (vs DMSO control): p = 0.272
150min 2uM DAPT (vs DMSO control): p = 0.001
160min 2uM DAPT (vs DMSO control): p < 0.001
170min 2uM DAPT (vs DMSO control): p < 0.001
180min 2uM DAPT (vs DMSO control): p = 0.001
To calculate pvalues, twotailed 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 140min 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 248252 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 resubmission 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 desynchronisation of single cells before resynchronisation 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 longterm 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 realtime 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 stateoftheart 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 inputoutput relationships and thus describing the dynamic response of the segmentation clock. Authors propose to study dynamic behaviour of the clock using coarsegraining and focus on describing the overall response over time while amalgamating spatial information. Appropriate coarsegraining 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 realworld 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 systemlevel responses of the somite segmentation clock by the coarsegrained theoreticalexperimental approach, applying the theory of entrainment to understanding the phase responses of mouse presomitic mesoderm (PSM) tissues in the presence of periodic perturbation of Notch inhibitor DAPT generated by microfluidics technique. It was demonstrated that the segmentation clock is responsive to diverse range of the perturbationperiods from 120 to 180 min, can be period and phaselocked, and the efficiency is dependent of the DAPT concentration (inputstrength). The authors also observed two cycles of the segmentationclock ticking in single cycles of 300 or 350 min periodperturbation, suggesting that higher order (2:1 mode) entrainment. They also applied stroboscopic maps to analysis and found that entrainmentphases 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 saddlenode 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 timeseries analysis proposed in Figure 2 and conclusion about the changes of the intrinsic period during entrainment. The validity of the global timeseries 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 timeseries 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" timeseries analysis; for instance, just focusing on the "putative posterior" regions that are close to sourcepositions of waves. Authors can recompute 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 2Dassays (new Figures S5S6). 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.
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 timelapse 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 reanalyzed the data to calculate the period from the stroboscopic maps. Mathematically speaking, we take the stroboscopic map as ${\varphi}_{\mathrm{n}\mathrm{e}\mathrm{w}\text{}}=\text{}({\varphi}_{\mathrm{o}\mathrm{l}\mathrm{d}\text{}}+\frac{{T}_{\mathrm{z}\mathrm{e}\mathrm{i}\mathrm{t}}}{{T}_{\mathrm{o}\mathrm{b}\mathrm{s}\mathrm{e}\mathrm{r}\mathrm{v}\mathrm{e}\mathrm{d}}}2\pi )\text{}mod\text{}2\pi \text{},$ and use this to estimate the period of oscillation in entrained samples $T}_{\mathrm{o}\mathrm{b}\mathrm{s}\mathrm{e}\mathrm{r}\mathrm{v}\mathrm{e}\mathrm{d}$ , in particular inverting the formula for 1:2 entrainment we have :
$\frac{{T}_{\mathrm{z}\mathrm{e}\mathrm{i}\mathrm{t}}}{{T}_{\mathrm{o}\mathrm{b}\mathrm{s}\mathrm{e}\mathrm{r}\mathrm{v}\mathrm{e}\mathrm{d}}}$=$(\varphi}_{\mathrm{n}\mathrm{e}\mathrm{w}\text{}$${\varphi}_{\mathrm{o}\mathrm{l}\mathrm{d}\text{}})/2\pi$ + 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 350minutes, 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 300min periodic pulses of 2uM DAPT (or DMSO for control). Additionally, we have added Figure 5D, which plots the average period in the 300min and 350min conditions. We summarize the median average period here with computed pvalues:
– 300min pulses of 2uM DAPT (or DMSO for control): pvalue = 0.191
CTRL: 130.39 mins
DAPT: 146.45 mins
– 350min pulses of 5uM DAPT (or DMSO for control): pvalue = 0.049
CTRL: 127 mins
DAPT: 174.86 mins
– 350min pulses of 10uM DAPT (or DMSO for control): pvalue = 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 300min and 350min conditions in Figures 5A and 5B, respectively.
38. In Figure 5A, the authors can show period evolution in the case of 300 min DAPTpulses 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 yaxis.
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 ${\varphi}_{\mathrm{e}\mathrm{n}\mathrm{t}\text{}}({T}_{\mathrm{z}\mathrm{e}\mathrm{i}\mathrm{t}})$ 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:${\varphi}_{\text{new}}=\text{}({\varphi}_{\text{old}}+PRC({\varphi}_{\text{old}})+\text{}\frac{{T}_{\text{zeit}}}{{T}_{\text{osc}}}2\pi )\text{}mod\text{}2\pi \text{},$ a fixed point requires that $\varphi}_{\mathrm{n}\mathrm{e}\mathrm{w}\text{}}=\text{}{\varphi}_{\mathrm{o}\mathrm{l}\mathrm{d}\text{}}\equiv \text{}{\varphi}_{\mathrm{e}\mathrm{n}\mathrm{t}$.
For 1:1 entrainment we then have:
$PRC({\varphi}_{\mathrm{e}\mathrm{n}\mathrm{t}})+\text{}\frac{{T}_{\mathrm{z}\mathrm{e}\mathrm{i}\mathrm{t}}}{{T}_{\mathrm{o}\mathrm{s}\mathrm{c}}}2\pi \text{}=\text{}2\pi$,
$PRC({\varphi}_{\mathrm{e}\mathrm{n}\mathrm{t}})\text{}=\frac{2\pi}{{T}_{\mathrm{o}\mathrm{s}\mathrm{c}}}({T}_{\mathrm{o}\mathrm{s}\mathrm{c}}{T}_{\mathrm{z}\mathrm{e}\mathrm{i}\mathrm{t}})$  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 $T}_{\mathrm{z}\mathrm{e}\mathrm{i}\mathrm{t}$ (as explained in the text when we introduce Figures8A12).
Another way to visualize fixed points is offered in Figure 16 DE, 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 phaseanalysis 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 negativefeedback 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 coarsegraining 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.sa2Article and author information
Author details
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 (NFRFE202000703)
 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 indepth 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
 Aleksandra M Walczak, CNRS LPENS, France
Reviewing Editor
 Sandeep Krishna, National Centre for Biological Sciences‐Tata Institute of Fundamental Research, India
Version history
 Preprint posted: October 21, 2021 (view preprint)
 Received: April 19, 2022
 Accepted: August 23, 2022
 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

 1,379
 Page views

 174
 Downloads

 2
 Citations
Article citation count generated by polling the highest count across the following sources: Crossref, PubMed Central, Scopus.
Download links
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)
Further reading

 Microbiology and Infectious Disease
 Physics of Living Systems
Microsporidia are eukaryotic, obligate intracellular parasites that infect a wide range of hosts, leading to health and economic burdens worldwide. Microsporidia use an unusual invasion organelle called the polar tube (PT), which is ejected from a dormant spore at ultrafast speeds, to infect host cells. The mechanics of PT ejection are impressive. Anncaliia algerae microsporidia spores (3–4 μm in size) shoot out a 100nmwide PT at a speed of 300 μm/s, creating a shear rate of 3000 s^{1}. The infectious cargo, which contains two nuclei, is shot through this narrow tube for a distance of ∼60–140 μm (Jaroenlak et al, 2020) and into the host cell. Considering the large hydraulic resistance in an extremely thin tube and the lowReynoldsnumber nature of the process, it is not known how microsporidia can achieve this ultrafast event. In this study, we use Serial BlockFace Scanning Electron Microscopy to capture 3dimensional snapshots of A. algerae spores in different states of the PT ejection process. Grounded in these data, we propose a theoretical framework starting with a systematic exploration of possible topological connectivity amongst organelles, and assess the energy requirements of the resulting models. We perform PT firing experiments in media of varying viscosity, and use the results to rank our proposed hypotheses based on their predicted energy requirement. We also present a possible mechanism for cargo translocation, and quantitatively compare our predictions to experimental observations. Our study provides a comprehensive biophysical analysis of the energy dissipation of microsporidian infection process and demonstrates the extreme limits of cellular hydraulics.

 Computational and Systems Biology
 Physics of Living Systems
The adaptive dynamics of evolving microbial populations takes place on a complex fitness landscape generated by epistatic interactions. The population generically consists of multiple competing strains, a phenomenon known as clonal interference. Microscopic epistasis and clonal interference are central aspects of evolution in microbes, but their combined effects on the functional form of the population’s mean fitness are poorly understood. Here, we develop a computational method that resolves the full microscopic complexity of a simulated evolving population subject to a standard serial dilution protocol. Through extensive numerical experimentation, we find that stronger microscopic epistasis gives rise to fitness trajectories with slower growth independent of the number of competing strains, which we quantify with powerlaw fits and understand mechanistically via a random walk model that neglects dynamical correlations between genes. We show that increasing the level of clonal interference leads to fitness trajectories with faster growth (in functional form) without microscopic epistasis, but leaves the rate of growth invariant when epistasis is sufficiently strong, indicating that the role of clonal interference depends intimately on the underlying fitness landscape. The simulation package for this work may be found at https://github.com/nmboffi/spin_glass_evodyn.