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 coarse-grained approach, in which theory guides the experimental identification of the segmentation clock entrainment responses. We demonstrate period- and phase-locking of the segmentation clock across a wide range of entrainment parameters, including higher-order coupling. These quantifications allow to derive the phase response curve (PRC) and Arnold tongues of the segmentation clock, revealing its essential dynamical properties. Our results indicate that the somite segmentation clock has characteristics reminiscent of a highly non-linear oscillator close to an infinite period bifurcation and suggests the presence of long-term feedbacks. Combined, this coarse-grained theoretical-experimental approach reveals how we can derive simple, essential features of a highly complex dynamical system, providing precise experimental control over the pace and rhythm of the somite segmentation clock.
Editor's evaluation
This manuscript undertakes a fundamental study of the entrainment of the mouse segmentation clock to an external oscillating signal. This provides insight into both the nonlinear properties of the clock as well as how it might integrate global signals. By combining experimental measurement of the finite pulse phase response curve with theoretical analysis using tools from nonlinear dynamics, this work produces a compelling coarse-grained model of the segmentation clock, albeit in an in vitro setting. In future, it will be of great interest to see how the predictions of this coarse-grained view, in particular the consequences of operating near a bifurcation, hold up in an in vivo system. This work is likely to be of interest to developmental biologists, synthetic biology experts, theoreticians and others interested in the emergence of complex patterning in living systems.
https://doi.org/10.7554/eLife.79575.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 coarse-grain degrees of freedom to obtain scale-free theories, allowing to define universality classes independent of the precise details of interactions (Cardy, 1996; Bialek, 2015). Another recent example is the parameter space compression theory, showing how complex systems (in biology or physics) can be typically reduced to simpler descriptions with few parameters (Machta et al., 2013; Proulx-Giraldeau et al., 2017; Hsu et al., 2020).
Going one step further, this suggests that one might be able to study – and control – complex systems provided we identify the essential, macro-level behavior. This is possible because only a limited number of universal descriptions exist, with defining behaviors and properties, that do not depend on the detailed implementation. A central challenge that remains is to implement these theoretical ideas to the experimental study of biological complexity.
Here, we develop a coarse-grained approach combining theory and experiments to study a cellular oscillator ensemble that constitutes the embryonic somite segmentation clock. Functionally, this clock controls the periodic formation of somites, the precursors of the vertebral column, and other tissues (Palmeirim et al., 1997; Krol et al., 2011). Molecularly, the segmentation clock comprises the oscillatory activity of several major signaling pathways, such as the Notch, Wnt, and Fgf signaling pathways, which show oscillatory dynamics with a period matching somite formation, that is, ∼2 hr in mouse embryos (Dequéant and Pourquié, 2008; Hubaud and Pourquié, 2014; Aulehla et al., 2003; Dequéant et al., 2006; Aulehla et al., 2008). More recently, segmentation clock oscillations with a period of ∼5 hr have been identified in human induced pluripotent stem cells differentiated into paraxial mesoderm, identifying a set of ∼200 oscillating genes, including targets of Notch and Wnt signaling (Chu et al., 2019; Diaz-Cuadros et al., 2020; Matsuda et al., 2020).
Strikingly, as individual oscillating cells are coupled to their neighbors via Notch-Delta signaling, the oscillations occur synchronized and wave-like activity patterns appear to periodically sweep along the embryonic anterior-posterior axis (Maroto et al., 2005; Masamizu et al., 2006; Aulehla et al., 2008; Herrgen et al., 2010; Webb et al., 2016; Yoshioka-Kobayashi et al., 2020; Rohde et al., 2021).
Adding to the complexity, these periodic spatiotemporal wave patterns are linked to an underlying spatial period gradient along the embryonic axis, that is, signaling dynamics in cells close to the posterior of the embryo oscillate faster compared to those in cells located more anteriorly. Such a period gradient linked to the segmentation clock has been identified in several species (Gomez et al., 2008; Morelli et al., 2009; Oates et al., 2012; Lauschke et al., 2013; Soroldoni et al., 2014; Tsiairis and Aulehla, 2016; Falk et al., 2022) and also in in vitro assays culturing intact or even dissociated presomitic mesoderm (PSM) (Lauschke et al., 2013; Tsiairis and Aulehla, 2016).
Of note, an analogous oscillatory system was also described during segmentation in arthropods (Clark et al., 2019) and while distinct at molecular level, it also exhibits spatiotemporal wave patterns traversing the embryo axis, again with indication of a period gradient (Sarrazin et al., 2012; El-Sherif et al., 2012; Brena and Akam, 2013; Schönauer et al., 2016).
In this work, we coarse-grain these underlying complexities and take a dynamical systems, macro-perspective on the segmentation clock, studying it as a single phase oscillator (Figure 1A). We build on the theory of synchronization and entrainment (see below) to first perform a systematic experimental characterization of its response to perturbation. We compare the outcome to qualitative and quantitative theoretical predictions. In turn, these experimental quantifications allow to derive a phase response curve (PRC) that uniquely characterizes the dynamical properties of the segmentation clock. This new insight provides the means to understand – and control – the timing of a complex embryological patterning process.
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 at the time of the zeitgeber is constant (technically, a fixed point of the Poincaré return map [Cross and Siggia, 2005]). This defines period-locking (also termed mode-locking) (Pikovsky et al., 2001). Entrainment is not always manifested and conditions for its existence can be derived. Quantitatively, when entrainment occurs, the zeitgeber induces a periodic phase perturbation (or response) of the entrained oscillator, which exactly compensates the detuning (or period mismatch, ) between the zeitgeber and the free-running oscillator. For this reason, when the detuning is very small, a weak external perturbation is enough to entrain an oscillator. Conversely, if the detuning is big, a strong signal and associated response is required for entrainment. One can then plot the minimal strength of the zeitgeber () versus the corresponding detuning (or simply , if is constant): these maps are more commonly known as Arnold tongues (Figure 1E). Arnold tongues predict the period- and phase-locking behavior in oscillatory systems as different as electrical circuits (Balanov, 2008), oscillatory chemical reactions (Ito, 1979; Kai and Tomita, 1979; Buchholz et al., 1985; Makki et al., 2014), or living systems such as circadian rhythms (Granada et al., 2013).
Lastly, more complex patterns of entrainment can be observed: for instance, stable phase relationships can be established where the entrained oscillator goes through cycles for every cycles of the external signal, defining period-locking. In that case, the instantaneous period of the oscillator matches the period of the zeitgeber (). Corresponding Arnold tongues can be obtained, leading to a rich structure for entrainment in parameter space (Figure 1F–H).
Previously, periodic activation of Notch signaling via heat shock-driven expression of the ligand Delta was used to modulate the segmentation of PSM in zebrafish (Soza-Ried et al., 2014). In the cited study, the readout to assess the effect of the periodic perturbation relied on somite size and morphology. While this gave important insight into how morphological segmentation of the PSM and somite length are affected by a zeitgeber, experimental investigation on the underlying signaling dynamics upon periodic perturbation remains limited.
In this work, to experimentally apply the theory of synchronization to the segmentation clock, we make use of a microfluidics-based entrainment setup, which we had established previously in the lab (Figure 1—figure supplement 1; Sonnen et al., 2018).
We showed before that using a quasi-two-dimensional in vitro segmentation assay (hereafter referred to as a 2D-assay), which recapitulates segmentation clock dynamics and PSM patterning (Lauschke et al., 2013), the microfluidics-entrainment approach allowed us to take control of Notch and Wnt signaling oscillations, providing direct functional evidence that the oscillation phase shift between Wnt and Notch signaling is critical for PSM patterning (Sonnen et al., 2018).
Results
A coarse-grained, single oscillator description of the segmentation clock
To perform a systematic analysis of entrainment dynamics, we first introduced a single oscillator description of the segmentation clock. We used the segmentation clock reporter LuVeLu, which shows highest signal levels in regions where segments form (Aulehla et al., 2008). Hence, we reasoned that a global region of interest (ROI) quantification, averaging LuVeLu intensities over the entire sample, should faithfully report on the segmentation rate and rhythm, essentially quantifying ‘wave arrival’ and segment formation in the periphery of the sample. Indeed, our validation confirmed that the oscillation period based on the global ROI analysis using the LuVeLu reporter closely matched the rate of morphological segment boundary formation (Figure 2D, Figure 2—figure supplement 1). In addition, global ROI quantifications using additional reporters to measure Wnt signaling oscillations (i.e. Axin2, Figure 2C–D, Figure 2—figure supplement 1) and the segmentation marker Mesp2 (Figure 2D, Figure 2—figure supplement 1) also showed close correspondence with LuVeLu global ROI measurements of segmentation clock period.
Hence, we conclude that LuVeLu global ROI timeseries analysis provides a valid coarse-grained quantification of the pace (period ∼ rate of segment formation) and rhythm (phase) of the segmentation clock.
The pace of the segmentation clock can be locked to a wide range of entrainment periods
Having established a quantitative, coarse-grained readout for segmentation clock pace and rhythm, we next analyzed whether pace and rhythm can be experimentally tuned using microfluidics-based entrainment (Figure 1A, Figure 1—figure supplement 1).
First, we tested whether the segmentation clock can be entrained to periods different and far from the endogenous period of ∼140 min, which we refer to as its free-running, natural period or . 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 (Figure 3A–C, Figure 3—figure supplement 1), the segmentation clock rhythm in DAPT-entrained samples closely adjusted to 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, period-locking was less precise (i.e. higher standard deviation as shown in Table 1) with 120 and 180 min zeitgeber periods, a possible indication that the limit of entrainment range is approached at these conditions.
We also tested the effect of varying DAPT concentration, as means of changing zeitgeber strength ( in Figure 1), on entrainment dynamics. Synchronization theory predicts that zeitgeber strength correlates with the time it takes to reach period-locking (Pikovsky et al., 2001; Granada and Herzel, 2009). To test this prediction, we entrained samples with periodic DAPT pulses at fixed intervals of 170 min and varied drug concentration. We indeed found that the time needed to show period-locking was shortened in samples using higher DAPT concentrations (Figure 4A–B, Figure 4—figure supplement 1A). Additionally, as expected, higher drug concentrations also resulted in more robust entrainment, indicated by the quantifications of the first Kuramoto order parameter , a measure for in-phase synchrony between samples (Figure 4C–D, Figure 4—figure supplement 1B, Figure 4—figure supplement 2) (see definition in Appendix 1, notice that is the modulus of the mean field variable used in Riedel-Kruse et al., 2007).
Higher-order entrainment of the segmentation clock
In theory, a non-linear oscillator should be amenable not only to 1:1 entrainment, but also to higher-order entrainment, in which cycles of the endogenous oscillation lock to cycles of the zeitgeber (Pikovsky et al., 2001; Mitarai et al., 2013). In practice, demonstration of higher-order entrainment is challenging due to the narrow permissive parameter region, that is, Arnold tongues are progressively narrower away from the 1:1 regime. Strikingly, we found experimental evidence for higher-order 2:1 entrainment (Figure 5). Samples entrained with either 300 or 350 min pulses showed evidence of 2:1 entrainment (Figure 5A–C), that is, the segmentation clock oscillated twice per each zeitgeber pulse and hence the segmentation clock rhythm adjusted to a period very close to 150 min for min and approaching 175 min for min (Figure 5C). We also found evidence for phase-locking during 2:1 entrainment, that is, a narrowing of the phase distribution among the individual samples (Figure 5D).
The segmentation clock establishes a stable entrainment phase relative to the zeitgeber
We next analyzed the phase-locking behavior between the segmentation clock and the zeitgeber. According to dynamical systems theory, phase-locking, that is, the entrainment phase (), 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 as a stable fixed point that lies on the diagonal, where there is phase-locking. Plotting the stroboscopic maps shows that in samples entrained to periodic pulses of DAPT the segmentation clock phase dynamics converge to a region on the diagonal, which marks (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 () was closer (or farther) to .
The phase of entrainment varies according to zeitgeber period
One fundamental property of entrained oscillatory systems is that its phase of entrainment varies as a function of the detuning, that is, the period mismatch between the endogenous oscillator and the entrainment period (Pittendrigh and Bruce, 2015; Aschoff and Pohl, 1978; Rémi et al., 2010; Bordyugov et al., 2015; Granada et al., 2013). In addition, theoretical studies have supported a ‘180 degree’ rule (Granada et al., 2013), stating that the phase of entrainment varies by half a cycle (180° or ) within the range of permissible zeitgeber entrainment periods, (Aschoff and Pohl, 1978; Granada et al., 2013).
To test these theoretical predictions, we quantified across a wide range of detuning, that is, from 120 to 180 min zeitgeber period. We found that , 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, systematically shifted from to (Figure 7A–B), spanning a range of almost (Figure 7B–C). Hence, as predicted by theory, our results show that the segmentation clock phase of entrainment varies as a function of the detuning relative to zeitgeber pulses.
It was previously documented that the period of tissue segmentation matches the period of oscillations localized in the posterior PSM, that is, center of of the 2D-assays (Tsiairis and Aulehla, 2016; Soza-Ried et al., 2014). We were then curious whether or not the dynamics from our coarse-grained analysis match that for oscillations in the center of the 2D-assays, and if the behavior of oscillations in this localized region also exhibits dependence to zeitgeber period. To study this, we compared the timeseries acquired from (a) global ROI spanning the entire field of view and from (b) a circular ROI restricted at the center of the 2D-assays (‘center ROI’, radius: 25 pixels, 1 pixel = 1.38 µm). Upon comparison, we noted that the oscillations at the center of 2D-assays adjusted its rhythm to the zeitgeber with strong correspondence to the coarse-grained oscillator (Figure 7—figure supplement 2), regardless of whether the segmentation clock was sped up (Figure 7—figure supplement 2A–B) or slowed down (Figure 7—figure supplement 2C–D). Moreover, the phase-locking between the zeitgeber and either the segmentation clock or the oscillations at the center of the samples are similar (Figure 7—figure supplement 2B1–B2, D1–D2, Figure 7—figure supplement 3). Concomitantly, for oscillations in the center of 2D-assays (i.e. corresponding to the posterior PSM), the entrainment phase at min is different to the entrainment phase at min (Figure 7—figure supplement 3B), here too exemplifying the effect of zeitgeber period on the phase of entrainment.
PRCs derivation and period change
Building on our finding that the phase of entrainment varies as a function of detuning (Figure 7B), our goal was to extract the quantitative information embedded in this dynamic behavior. This is possible since the phase of entrainment dynamics reflect, at a quantitative level, the fundamental properties of a dynamical system that responds to external perturbations. This behavior, in turn, can usually be captured with a single function, the PRC (Izhikevich, 2006; Granada et al., 2013). The PRC describes the change of phase induced by a perturbation, and a priori depends on both the nature of the perturbation received and the phase of the cycle. Because of the direct dependence of the phase of entrainment dynamics on the PRC, our goal was to use the experimental data to gain insight into the segmentation clock PRC.
We follow Izhikevich, 2006, to define a generalized PRC not tied to the shape or duration of the perturbation. Without loss of generality, we then model perturbations in form of pulses of amplitude (Pikovsky et al., 2001) and write:
where is the phase of the segmentation clock on the cycle at the moment of the pulse, and the phase after the pulse (which can be defined even if the system transiently moves outside of the limit cycle, via isochrons, see Winfree, 1967). Amplitude is a dimensionless parameter, quantifying zeitgeber strength in the mathematical formalism we develop below (see Appendix 1, section Arnold tongue computations to see how we relate to DAPT concentrations). If, following a perturbation, the oscillator relaxes quickly toward the limit cycle, one can use the PRC to compute response to periodic pulses with period . 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 for a given . When is varied, different phases of entrainment are observed, here plotted in Figure 7B.
To derive the PRC directly from the experimental stroboscopic maps, we invert Equation 2 into:
By estimating for a given , 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 , even far from the entrainment phase . Notice, however, that phases far from will be only sampled over the first few pulses (so with possibly much noise) while conversely, 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 is changed.
Such vertical shifts in PRCs as a function of entrainment period have been previously observed in so-called ‘overdrive suppression’ for cardiac cell oscillators (Kunysz et al., 1995). Shifts occur when the entrainment signal impacts the biochemical control of the system, leading to a change of intrinsic period (here ). 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 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 min to be the ‘reference’ curve, since this coincides with the natural period of the process. We can then estimate for each entrainment period the period shift most consistent with the data, to obtain a single, period-independent PRC in Figure 8A2 (see Appendix 1, section Theoretical methods). From there, we build upon the simplest non-linear phase-oscillator, that is, the classical Radial Isochron Cycle (RIC). To modulate its sine-form PRC, which is incompatible with our data, we perturb it into an Elliptic Radial Isochron Cycle with Acceleration or ERICA.
ERICA is designed to have radial isochrons, meaning that the phase (and the PRC) can be analytically computed in the entire plane. ERICA then allows for the speeding up of the cycle for some angle range (parameter ) or for changes of the limit cycle shape into an ellipse of increasing eccentricity (parameter ). Both high values of and 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 , high value of (indicating a strong elliptical shape), and high value of over half a cycle (indicative of a relaxation type dynamics in with a slow decay phase followed by a fast reset). Figure 8C compares the PRC of this optimized model with the multiple data points, showing excellent agreement. In addition, we combined the ERICA framework with a simple fit for intrinsic period change (Figure 8D) to account for the experimental phase transition curves, as illustrated in Figure 8G (see also Figure 8—figure supplement 2).
With the ERICA model at hand, we derived numerically the Arnold tongues of the system and the phase/detuning curves for all entrainment parameters (Figure 8E–F). We notice that the Arnold tongues are heavily skewed toward the right, meaning that the system can be more easily slowed down than sped up, consistent with the negative PRC shape. Remarkably, while we build the model using only one entrainment drug concentration, we can also largely explain data obtained at other concentrations (Figure 8F). In particular, the Arnold tongues/our model predict a specific change of the entrainment phase as the entrainment strength is varied (Figure 8—figure supplement 3G). The comparison to experimental data (Figure 8—figure supplement 4) shows that this prediction is, qualitatively, verified. More generally, having the PRC and a minimal, coarse-grained model that captures the essential dynamical features of the segmentation clock during entrainment, including the change in intrinsic period, enables predictable control over the pace and the rhythm of the segmentation clock using the entrainment strategy.
Discussion
In this work, we used a coarse-graining, entrainment approach to gain new insights into the properties of the segmentation clock from dynamical systems’ perspective. We mode-lock the segmentation clock to various entrainment periods and use the information about the dynamic phase-locking behavior to derive the somite segmentation clock PRC from the experimental data.
A coarse-graining approach captures essential dynamical features using a simple ‘one oscillator’ phase description
Given the complexity underlying the somite segmentation clock, comprised of several, interconnected, signaling pathways with countless molecular interactions, it was a priori unclear whether a simple ‘one oscillator’ phase description and a perturbation approach would capture the essential dynamical characteristics at the systems level. Another potential difficulty arises from the fact that cellular oscillators define a phase gradient, controlling the segmentation pattern (Lauschke et al., 2013). For all those reasons, one could imagine that no single phase could describe the system’s behavior, a general concern for systems of interacting oscillators already mentioned by Winfree, 1967.
Remarkably, one key finding of this work is that indeed, using a systems-level single oscillator phase description, we observe a consistent entrainment behavior, that is, period-locking with convergence toward a well-defined entrainment phase, as predicted from the oscillator phase response theory. It is important to point out that theoretical predictions are non-trivial and quantitative, such as the higher order entrainment (Figure 5) and the dependence of the entrainment phase on detuning (Figure 7A–B). Given these experimental findings we conclude that the coarse-graining approach and the description of the entire system using a single oscillator phase is justified and enables to extract the essential dynamical properties, which are captured by a mathematical model including the system’s PRC.
Asymmetrical segmentation clock PRC compatible with SNIC bifurcation
Insight into the system’s PRC allows to infer about the nature and characteristics of the oscillatory system, independent of its specific molecular realization (Izhikevich, 2006). For instance, excitable systems can be functionally categorized with their PRC, that is, in neuroscience, Class I excitable systems have PRCs with constant sign (Izhikevich, 2006). The canonical example of (Class I) oscillators with constant-sign PRCs is the quadratic fire and integrate model (Izhikevich, 2006), which can be used to model neuronal spiking.
In contrast, Class II excitable systems exhibit sinusoidal PRCs, examples include many biological systems such as the circadian clock (Granada et al., 2009).
This difference in PRC reflects fundamental properties of the underlying dynamics, including the associated bifurcations: Class I systems are associated with infinite period bifurcations, a typical example being the saddle-node on invariant cycle (SNIC) bifurcation. For such bifurcations, the period can be tuned and can become arbitrarily long depending on the proximity to the bifurcation, and oscillations present relaxation-type dynamics (i.e. combination of slow dynamics with fast resetting). The simplest example, the quadratic fire and integrate model, is a one variable oscillator, where goes from to in a finite time before resetting to (Izhikevich, 2006). The PRC for this oscillator is ‘naturally’ of constant sign because there is only one variable, which can only advance (or delay) in response to a constant signal. On the other hand, Class II systems are associated with Hopf bifurcations (Izhikevich, 2006), have a fixed period, and more sinusoidal dynamics close to the bifurcation.
Remarkably, the segmentation clock PRC that we computed here shares features of systems close to an infinite period bifurcation, that is, a highly asymmetrical, constant-sign PRC with an extended flat region. Consistent with this, the analytical ERICA model that we developed to capture the experimental PRCs in a wide range of conditions exhibits relaxation-type dynamics (Figure 8B, slow phase in white and fast resetting phase in gray). It is nevertheless important to point out that a given PRC can result from different oscillators, that is, a constant sign does not prove the existence of a SNIC oscillator. For instance, delayed oscillators close to a Hopf bifurcation can show PRCs with an (almost) constant sign (Kotani et al., 2012). Infinite period can also be obtained with saddle homoclinic bifurcations (Izhikevich, 2006). However, such models are mathematically more complex, for instance delays include long-term memory in the system, and homoclinic bifurcations require the existence of other attractors outside of the cycle.
In agreement with SNIC oscillators, the segmentation clock does show a tunable period and a slowing down behavior at multiple levels: it is known that the overall rate of segmentation slows down over developmental time, a feature described in several species. In addition, the oscillations also slow down along the embryonic axis, that is, a period gradient is present within the PSM, reflecting tunability of oscillations at cellular level, as PSM cells progress toward differentiation and eventually stop oscillations. How the tunability at cellular level is linked to the system-level change in segmentation clock rate is not understood (see also below for macro- versus micro-scale comparison). However, our findings of mostly negative PRCs might reflect that the segmentation clock has a natural bias of slowing down the oscillations.
From a network standpoint, biological networks with tunable periods have been associated with specific structural features, combining negative and positive feedback loops (Tsai et al., 2008). Positive feedback in particular puts oscillators close to a multistable regime, leading to relaxation oscillations (Ferrell et al., 2011). SNICs appear naturally when the regulatory logic of a system moves from a negative feedback oscillator to multistability (Jutras-Dubé et al., 2020), reflecting network modularity. Such modularity and robustness with changes in regulatory logic is a hallmark of developmental plasticity (West-Eberhard, 2003) and hence one would anticipate to find SNIC oscillators to be abundant in developing systems. First examples are indeed emerging, such as during Caenorhabditis elegans development (Meeuse et al., 2020).
Findings not predicted by the theory of PRC
We also made several unexpected findings, not predicted by the theory of PRC.
First, we found evidence that during entrainment, not only does the segmentation clock adjust its observed period to the zeitgeber pulses, but also changes its intrinsic period in the direction of the zeitgeber rhythm. Hence, during entrainment that slows down the clock (i.e. 170 min), we find evidence that the intrinsic period lengthens (not just the observed rhythm), while during entrainment that speeds up the clock (i.e. 120 min), we find evidence that the intrinsic period decreased. Again, such a change of intrinsic period is not predicted in entrained systems. The simplest explanation could be that the drug pulses change the period of the oscillator by changing some biochemical parameter in the system, similar to overdrive suppression in cardiac cells (Kunysz et al., 1995). However, as stated above, we do not find evidence for a consistent slowing down or speeding up effect: the intrinsic period is decreased for short entrainment cycles and increased for longer entrainment cycles. This rather suggests the existence of feedbacks of the clock on itself, leading to higher-order adaptations beyond the rapid Notch phase response. One can only speculate on the mechanisms underlying such adaptation, but this is compatible with the idea that two interacting oscillators control the intrinsic period. Here, it is possible that entraining Notch with a zeitgeber modifies the Wnt oscillation period, which in turn feeds back on Notch on a longer time scale. So the internal period change that we see might in fact come from the induced change of Wnt period. We have shown previously that Wnt and Notch oscillators are coupled but are not phase-locked, with functional impact on tissue patterning (Sonnen et al., 2018). Alternatively, a similar role could also be played by the long inter- or intracellular delays in the system, postulated in multiple theoretical works (Lewis, 2003; Morelli et al., 2009; Ares et al., 2012). Such delays could effectively couple multiple cycles, changing clock parameters (such as the period) beyond the instantaneous phase response. More experimental and theoretical work is needed to explore these ideas.
Second, a striking outcome we obtained was that even after entrainment, a spatial period gradient and phase waves emerged over time (Appendix 1—figure 1). Put differently, while the overall system is entrained, as evidenced by a control of the timing of morphological segmentation and of segmentation clock rhythm, the underlying cellular oscillators show a divergent, yet, coherent response, that is, a spatial period gradient emerges within the tissue. A recent entrainment model of coupled PSM oscillators (Juul et al., 2018) does not account for the emergence of a period gradient during entrainment. This unexpected behavior hence awaits further experimental and theoretical analysis. In particular, the role of intra- and intercellular coupling needs to be addressed to reveal how the macro-scale behavior relates to the underlying cellular scale oscillations.
Conclusions
Our work demonstrates how, despite all the molecular and functional complexities, coarse-graining and theory can be used to effectively take control of complex biological processes. A molecular mechanism is not needed to exert control, as long as we have a mathematical one, one that captures the essential features of a system at a meaningful coarse-grained level.
We also aim to illustrate the potential of an integrated, theoretical-experimental approach to complex biological systems: while from theoretical viewpoint, the fact that entrainment phase varies as detuning is altered is a mathematical necessity and hence ‘obvious’, this outcome is not at all intuitive, a priori, from an experimental-observational viewpoint. Theory guides experimentation, that is, Figure 1 Theory (Phillips, 2015), leading to new insight and understanding of a complex biological phenomenon.
Future studies are needed to gain further insight into the response to other perturbation regimes. In addition, it will be essential to perform the analysis both at the macro-scale, as done in this study, but also extending it to the cellular scale, that is, taking spatial differences, such as the period gradient, fully into account. We consider that such a dynamical systems’ theoretical-experimental approach is a promising and, as we show here, feasible way forward with the goal in categorizing the segmentation clock in its universality class.
Materials and methods
Please refer to the Appendix 1 for detailed materials and methods.
Appendix 1
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 Mesp2-GFP mouse line generated by Morimoto et al., 2006. Axin2-linker-Achilles, a mouse line expressing Axin2 fused with a GSAGS linker to Achilles, a fast-maturing variant of YFP (Yoshioka-Kobayashi et al., 2020), was generated in-house. To generate the knock-in alleles, we targeted the stop codon of endogenous Axin2 locus with vector containing the reporter sequence coding for Achilles and a selection cassette. This targeting vector was constructed as follows: linker-Achilles-loxP-PGK Neo-loxP. The selection cassette was flanked by loxP- sites for eventual Cre-mediated excision. Axin2-linker-Achilles knock-in reporter line was generated by standard gene targeting techniques using R1 embryonic stem cells. Briefly, chimeric mice were obtained by C57BL/6 blastocyst injection and then outbred to establish the line through germline transmission. The Achilles/pRSETB plasmid was a gift from the lab of Atsushi Miyawaki at RIKEN Center for Brain Science (RIKEN-CBS) in Japan. Mice were kept in an outbred background and were housed in the EMBL Laboratory Animal Resources (LAR). All animal experiments were conducted under veterinarian supervision and after project approval by European Molecular Biology Laboratory, following the guidelines of the European Commission, Directive 2010/63/EU and AVMA Guidelines 2007.
Media preparation
Request a detailed 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% CO2 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 reporter-of-interest (e.g. LuVeLu) using a stereo fluorescence microscope. The tailbud was cut from the rest of the tail and was immediately transferred to pre-equilibriated culture medium that was supplemented with HEPES (170 µL of 1 M HEPES in 10 mL culture medium). These isolated embryonic tissues were then directly loaded in the microfluidics device.
Preparations for microfluidics-based entrainment of 2D-assays
We used in this study a microfluidics-based experimental entrainment platform with a general protocol elaborated in a recent publication (van Oostrom et al., 2021). The PDMS microfluidics device was made using standard soft lithography (El Debs et al., 2012) and as previously described (Sonnen et al., 2018). The ratio of Sylgard 184 silicone elastomer base (Dow) to curing agent (Dow) was 9:1 (w/w). PDMS chip was attached to cover glass (70 mm × 70 mm, 1.5H, Marienfeld 0107999 098) via plasma bonding. UV irradiation of PTFE tubing and PDMS device PTFE tubing (inner diameter: 0.6 mm, APT AWG24T) and syringe needles (22 G 1 1/4 – Nr. 12), as summarized in Appendix 1—table 3, were prepared a day prior to actual microfluidics-based entrainment experiment. In addition to one PDMS microfluidics device (Figure 1—figure supplement 1), each experiment required four 3 m PTFE tubing (each with syringe needle inserted in one end) for the drug/medium inlets, two 1 m PTFE tubing (each with syringe needle inserted in one end) for the outlets, and twenty-four 1 cm plugs made from cut PDMS-filled PTFE tubing for the sample inlets and unused drug/medium inlets. These were all sterilized under UV for at least 20 min.
Fibronectin-coating of PDMS device and overnight soaking in buffer
Request a detailed protocolWhile waiting, 5 mL of PenStrep (Gibco, 15140-122) was added to 495 mL of ×1 PBS (PBS + PenStrep); 5.6 mL of the buffer was set aside to prepare fibronectin solution, while the rest was poured in a glass dish. After UV irradiation, the sterilized PDMS device was immersed in the PBS + PenStrep and bubbles were removed by flushing the channels with buffer. PDMS-filled PTFE plugs were also immersed in PBS + PenStrep. To prepare the fibronectin solution, 280 µL of fibronectin (Sigma-Aldrich, F1141) was added to the set aside 5.6 mL PBS + PenStrep. At least 2.5 mL of fibronectin solution was loaded into a 3 mL syringe (diameter: 8.66 mm, BD Luer-Lok REF 309658). A UV-irradiated needle, which was earlier inserted into a 1-meter PTFE tubing, was attached to the filled syringe. The tubing was then inserted into an outlet in the PDMS device, carefully avoiding introduction of bubbles. Fibronectin was flowed (flow rate: 50 µL/hr) into the PDMS device at room temperature overnight.
Preparation of syringes containing drug/medium
Request a detailed protocolFor a microfluidics-based experiment with periodic pulses of drug, four 10 mL syringes (diameter: 14.5 mm, BD Luer-Lok REF 300912) were filled with either the drug, DMSO control, or culture medium (see Appendix 1—table 2). Components of solution in each of these syringes are specified in Appendix 1—table 4. Drug used in this study was DAPT (Sigma-Aldrich, D5942-5MG). To prepare 10 mM stock of DAPT, 5 mg DAPT (MW = 432.46 g/mol) was dissolved in 1156.2 µL DMSO (Sigma-Aldrich, D8418). This was aliquoted and stored at –20°C until use.
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 UV-irradiated 3 m PTFE tubing via attached syringe needle, and were carefully mounted on programmable pumps (World Precision Instruments, AL-400) next to the microscope. Gas in the tubing was displaced with drug/medium by careful pushing of the syringes’ plunger. Pumps were turned on and flow rate was set to 900 μL/hr. The microscope was then equilibriated to 37°C and 5% CO2. After recovery of embryonic tissues, using a pipette (i.e. P200 for 2D-assays and P1000 for intact PSM), each sample was carefully loaded into the microfluidics device, which was already coated with fibronectin, immersed in a buffer of PBS and PenStrep, and degassed. Each sample inlet was plugged with a PDMS-filled PTFE tubing immediately after sample loading. Unused inlets, if any, were also plugged. Flow rate of drug/medium was set to 20 μL/hr and the tubings were carefully inserted into the drug/medium inlets in the PDMS device while it is immersed in buffer. The tubings connected to the syringes with medium were inserted first, and the tubing connected to the syringe with the drug was inserted last. A 15 min timer was started after insertion of the drug tubing. PDMS device was then removed from the buffer and excess liquid on the cover glass was removed carefully with lint-free wipes (Kimberly-Clark). The PDMS device, with more than 1 m of attached tubings, was carefully placed inside a pre-equilibriated microfluidics holder (EMBL Mechanical Workshop, Appendix 1—figure 2) customized to fit the 70 mm × 70 mm cover glass (Marienfeld 0107999 098) and some PTFE tubing. Putting some of the tubing inside the microfluidics holder was necessary to equilibriate the drug/medium to desired environmental conditions before they are perfused into the microfluidics device. The cover glass was secured in place with grease and a U-shaped metal clamp. The microfluidics holder was then carefully mounted on the stage of the microscope, and the end of the outlet tubings were placed in a beaker. Fifteen minutes after insertion of the drug tubing in the PDMS device, each pump was tilted on its side and was equilibriated for another 15 min. Afterward, the pump mounting the syringes with the drug and DMSO control was turned off, and the pump mounting the syringes with the medium was set to flow rate of 60 µL/hr. Samples were equilibrated at these conditions for at least 30 min before start of imaging/entrainment.
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, AL-400). Diameter of syringe (14.5 mm for 10 mL syringe, BD Luer-Lok REF 300912) was accordingly set to match a defined flow rate and a specified volume of solution to be perfused with the duration of perfusion. Standard pumping programs of medium and drug/control are summarized in Appendix 1—table 5, respectively, considering an entrainment experiment with = 170 min.
Confocal microscopy
Request a detailed protocolFor most of the experiments here, samples were imaged in an LSM 780 laser-scanning microscope (Carl Zeiss Microscopy) fitted with an incubation chamber (EMBL Mechanical Workshop). A Plan-Apochromat ×20 air objective with a numerical aperture (NA) of 0.8 (Carl Zeiss Microscopy) was used for imaging, and the zoom was set to 0.6. Three z-stacks (spacing: 8 µm) were scanned for each sample every 10 min to acquire timelapse movies. Imaging of multiple samples in multiple locations was done with a motorized stage, controlled using Zen Black software (Carl Zeiss Microscopy), and automated using a VBA macro developed by Antonio Politi (Politi et al., 2018), which is available at https://github.com/manerotoni/mypic. The dimension of the images was 512 pixels × 512 pixels, with a pixel size of 1.38 µm and bit depth of 16 bit. Detection of drug pulses, using Cascade Blue (excited with 405 nm) added to the solution, was also done every 10 min with lower image resolution: 1 z-stack, 32 pixels × 32 pixels (pixel size = 22.14 µm). Imaging and automated pumping of drug/medium through the microfluidics device were started simultaneously.
Data analysis
Details on the coarse-graining strategy and analyses of entrainment dynamics are specified here.
Extracting timeseries from global intensity of 2D-assays for subsequent analyses
Request a detailed 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). Z-stacks were first projected based on maximum intensity (in Fiji: Image > Stacks > Z Project). Subsequently, the timeseries was obtained by plotting the z-axis profile (in Fiji: Image > Stacks > Plot Z-axis Profile). Timeseries of replicate samples were compiled in a .txt file for subsequent analyses. To extract timeseries corresponding to signaling oscillations localized at the center of 2D-assays, samples that were centered in the field of view were considered. A 50 pixel × 50 pixel oval ROI was specified at the center (x=256 and y=256) of a re-oriented (and registered, if necessary) 2D-assay timelapse file (in Fiji: Edit > Selection > Specify). Timeseries was extracted after specifying the center ROI (in Fiji: Image > Stacks > Plot Z-axis Profile). Timeseries of replicate samples were compiled in a .txt file for subsequent analyses. Definition of Kuramoto order parameter We follow the definition from Pikovsky et al., 2001: introducing the mean field variable
(see also Riedel-Kruse et al., 2007), the first Kuramoto order parameter is defined by
Monitoring period-locking and phase-locking
Request a detailed protocolEntrainment was evaluated based on period-locking and phase-locking of the signaling oscillations to the periodic drug pulses. Oscillatory components were extracted from timeseries using a wavelet analysis workflow that was developed by Gregor Mönke, which was recently implemented as a Python-based standalone software (Mönke et al., 2020c, available at https://github.com/tensionhead/pyBOAT; Mönke, 2020a). In this workflow, timeseries was first detrended using a sinc filter and then subjected to continuous wavelet transform. The Morlet wavelet is used as the mother wavelet (Mönke et al., 2020c). Time-resolved frequency analysis was done by cross-correlating the signal to wavelet functions of known frequencies, generating a power spectrum. A high power score was assigned to wavelets that correlated well with the signal relative to white noise. Instantaneous period and phase were extracted upon evaluation of the power spectrum along a ridge tracing wavelet with maximum power for every timepoint. Phase dynamics of signaling oscillations, upon subjecting them to periodic perturbation, were analyzed using stroboscopic maps (Balanov, 2008; Cross and Siggia, 2005; Isomura et al., 2017). Briefly, the phase difference () was defined as
where is the time of perturbation and is the zeitgeber period (i.e. one cycle after time = ). and 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 (), was determined considering only phase-locked samples (where the difference between a sample’s phase at the time of the final drug pulse considered and its phase one drug pulse before is less than ). 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 is the first Kuramoto order parameter. As wavelets only partially overlap the signal at the edges of the timeseries, resulting in deviations from true phase values (Mönke et al., 2020c), the first and last pulse pairs were not considered in the generation of stroboscopic maps. Polar plots were also generated summarizing the instantaneous phase of replicate samples and their first Kuramoto order parameter as shown in Figure 4D, Figure 5D, and Figure 4—figure supplement 1B. The Python code is available as a Jupyter notebook (.ipynb) at https://github.com/PGLSanchez/EMBL_OscillationsAnalysis/tree/master/EntrainmentAnalysis (Sanchez, 2022; copy archived at swh:1:rev:bd4a258570e7b7b5f0c179bc976f0d45988efa1a). This code uses Matplotlib (Hunter, 2007), NumPy (van der Walt et al., 2011; Harris et al., 2020), pandas (McKinney, 2010), scikit-image (van der Walt et al., 2014), SciPy (Virtanen et al., 2020), and seaborn (Waskom et al., 2014).
Detrended timeseries of replicate samples were in some part of the study represented as a heatmap using PlotTwist (Goedhart, 2020), as shown in Figure 4A. Average periods (from 650 to 850 min after start of experiment) were meanwhile plotted using PlotsOfData (Postma and Goedhart, 2019), as shown in Figure 3C–D and Figure 3—figure supplement 1. These apps are available at https://huygens.science.uva.nl/PlotTwist/ and at https://huygens.science.uva.nl/PlotsOfData/, respectively. To calculate p-values, two-tailed test for absolute difference between medians was done via a randomization method using PlotsOfDifferences (Goedhart, 2019). While not recommended, p-values were also calculated similarly for conditions with n<10. PlotsOfDifferences is available at https://huygens.science.uva.nl/PlotsOfDifferences/.
Generating wavelet movies
Period and phase wavelet movies were generated using the Wavelet Processing and Export workflow developed by Gregor Mönke, which runs on the EMBL cluster and is implemented in Galaxy with technical assistance from Jelle Scholtalbers (EMBL Genome Biology Computational Support). This workflow extracts timeseries of every pixel in a timelapse movie and subjects them to sinc filter-based detrending and subsequent continuous wavelet transform (Mönke et al., 2020c). This results in extraction of instantaneous period and phase of each pixel, recovering period and phase wavelet movies corresponding to the input timelapse movie (for scheme, please refer to Appendix 1—figure 1A). The workflow is available at https://github.com/tensionhead/SpyBOAT, (Mönke, 2020b).and can be used via a public Galaxy server at https://usegalaxy.eu/root?tool_id=/spyboat. Settings used to generate wavelet movies in this study were: sigma of 8.0, sample interval of 10 min, period range from 100 to 250 min, and number of periods analyzed of 151.
Examining period gradient in 2D-assays
Re-oriented 2D-assays (dorsal side up) and their corresponding period and phase wavelet movies were used for the analyses. Temporal evolution of the period gradient during the course of entrainment experiments was evaluated from the period wavelet movies of the 2D-assays. As the period wavelet movies also contained wavelet transformations for pixels in the background, a binary mask was first created to differentiate pixels corresponding to signal. To create the mask, re-oriented (and registered, if necessary) timelapse movies of 2D-assays were blurred using a Gaussian blur (in Fiji: Process > Filters > Gaussian Blur [sigma radius: 8, scaled units in µm]). Then, signal was specified by thresholding (in Fiji: Image > Adjust > Threshold [Default method, dark bakground]). After thresholding, pixels corresponding to the signal were assigned a value of 255, while those corresponding to the background were assigned a value of 0. If opposite, the values were inverted (in Fiji: Edit > Invert). Then, the binary mask (signal = 1 and background = 0) was created by dividing all values by 255 (in Fiji: Process > Math > Divide [value: 255]), and was used to mask the period wavelet movie.
Theoretical methods
PRC correction and computation
PRCs at different entrainment periods appear shifted with respect to one another. To correct for this, we fitted the PRCs at entrainment periods 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 for min, respectively. The data points shown in Figure 8C are used for MC optimization (see next section). Our choice was later validated by recomputation of PRCs and comparison with data in Figure 8G. We also show the corresponding PTC and stroboscopic maps in Figure 8—figure supplement 3D-E.
Model
Our model (Appendix 1—figure 3A) is based on simple modifications of the classical Poincaré oscillator, otherwise called RIC model. The main motivation is to find the simplest possible model able to reproduce the general shape of the experimental PRC, with both a flat and a negative region, while being analytical with simple isochrons. The full description and analytical calculations for this model will be described elsewhere, here we outline the basic equations and some general properties of the model. For PRC computations, we will assume that a perturbation is a horizontal shift of magnitude (Appendix 1—figure 3B).
To flatten the PRC for half the cycle, we first modify the limit cycle of the standard RIC model into an ellipse, keeping the origin as one of the foci. This is done through the introduction of a parameter , so that the corresponding equation of the ellipse is
combined with radial equidistant isochrons, , where is the polar angle. This is however not enough to define the entire flow in the plane since one needs to specify how 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 (again assuming ). This leads to the following equations in cartesian coordinates:
(9)
We name this model Elliptic Radial Isochron Cycle, or ERIC. An important feature of ERIC is that the angle in the plane is the phase of the oscillator, in particular since the phase is defined by the planar angle, the PRC following a horizontal perturbation of size toward the right can be computed in a straightforward way and is:
As said above, the effect of introducing is to smoothly flatten the PRC over half of the cycle, as illustrated in Appendix 1—figure 4A.
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’ 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 functions linear by piece, that is, for one sector and otherwise. The sped up sector is centered at angle and has width , so that the modified period of the cycle is . It is also convenient to define the rescaled angular velocity . From there, the phase of the cycle as a function of the angle in the plane is the simple linear transformation (defining ):
for ,
for , and
for . Those functions ensure that the rate of phase evolution in sector is times the rate in the other sectors (compare coefficients in Equations 14-16), that angle 0 in is phase , and that phase is continuous so that and . Notice also that for we get from Equation 16 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 , , and independently, for various choices of sped up sectors.
For optimization purpose, it is easier to numerically compute the PRC in the following way. We consider an ensemble of angles linearly spaced on the interval . We then compute the corresponding position on the ERIC limit cycle, and compute numerically the angle of the point at distance on the right. By definition of ERICA, the PRC for each index 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 , , , sped up sector location and width , and the location of the zero-phase reference point on the cycle . We minimize :
where are data points of the corrected experimental PRC and is the noise, estimated as the root mean square error of a Fourier fit to data (see Figure 8A). We use Markov chains to generate distributions of parameters and take the average values to plot the model limit cycle and PRC. The optimized parameter values are: , , , , , .
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 ( 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 changes linearly with .
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 :
where the factor is the unit conversion from time to phase. This equation defines in an implicit way a function , 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 with the help of the :
There is a geometric interpretation here: the curve is in fact a (rescaled) rotation of the PRC! This can be seen here because Equation 23 combines a mirror image of the PRC along the axis (minus sign) with a mirror image along the first diagonal (power −1 which corresponds to the inversion function exchanging and ). 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 ). For instance, if the PRC is sinusoidal, 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 is compared with the curve 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 , for example, combining Equations 20-23 and counting time in units of the intrinsic period , we get
This defines an absolute, minimum slope of the curve . 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 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 , the right-hand side of Equation 26 can be very small, so that one can get , 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 could be indicative of a mechanism where the oscillator adapts its intrinsic period (here to the zeitgeber period).
Practically, we obtain the dependency by fitting the stroboscopic maps and entrainment phase from all experiments with 2.0 µM DAPT, using the PTC of our optimized model (Figure 8—figure supplement 3D-E). The values of 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 min and above min which indicates that is becoming smaller around those zeitgeber periods. More biologically, this indicates that the system does not adjust to any entrainment period. Again some arbitrary choices have to be made, but based on those constraints and our own interpolation, a cut-off for of ±20% of the natural period of 140 min was assumed, and is further consistent with the relative period changes observed experimentally in zebrafish segmentation clock mutants (Herrgen et al., 2010) and theoretically derived from delayed coupled models (Moreno-Risueno et al., 2010). With such cut-off, the obtained entrainment range of zeitgeber periods is from 118 to 200 min.
Computation of Arnold tongue and isophases in Figure 8F was done by iterations of the return map Equation 2 and detection of fixed points, using as the intrinsic period the interpolation and extrapolation shown in Figure 8—figure supplement 3F.
Lastly, we can use this model to compute the entrainment phase for different DAPT concentrations. For the main concentration of 2.0 µM DAPT, we have the optimized amplitude of perturbation . Assuming that the DAPT concentration is reflected in the strength of perturbation, we can find cross-sections of the isophases at different values of to have an excellent agreement with experiments, as illustrated in Figure 8—figure supplement 3G. We take for 3.0 µM, for 1.0 µM, and for 0.5 µM DAPT, as shown on the Arnold tongue in Figure 8F. We capture both the presence/absence of entrainment and the relation. In particular, we see that for lower DAPT concentrations, the range of entrainment is smaller but we still get plateaus in entrainment phase consistent with the change of period we computed.
Data availability
Source data and Python codes are made available via GitHub (https://github.com/PGLSanchez/EMBL_OscillationsAnalysis [copy archived at swh:1:rev:bd4a258570e7b7b5f0c179bc976f0d45988efa1a] and https://github.com/PGLSanchez/EMBL-files).
References
-
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. Google-Books-ID 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 BZ-reaction in a CSTR: chemical resonance, entrainment and quasi-periodic behaviorIn: Rensing L, Jaeger NI, editors. Temporal Order Springer Series in Synergetics. Springer. pp. 116–121.https://doi.org/10.1007/978-3-642-70332-4
-
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/S0076-6879(08)03801-9.
-
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 self-oscillating 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 delay-induced 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. Google-Books-ID.
-
Synchronised cycling gene oscillations in presomitic mesoderm cells require cell-cell 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/Majora-92bf1922-00a
-
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/9781400876167-006
-
Untangling the hairball: fitness-based 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 open-source platform for biological-image 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/0022-5193(67)90051-3
Article 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 (NFRFE-2020-00703)
- Paul François
The funders had no role in study design, data collection and interpretation, or the decision to submit the work for publication.
Acknowledgements
This work is presented, in part, in the doctoral thesis of Paul Gerald Layague Sanchez (https://doi.org/10.11588/heidok.00029209) submitted to the University of Heidelberg, Germany. We thank Mogens Høgh Jensen, Jonas Juul, Mathias Heltberg, Sandeep Krishna, Adrián Granada, and Istvan Kiss for in-depth discussions on Arnold tongues and PRC; Justin Crocker, Anne Ephrussi, Lars Hufnagel, Jelle Scholtalbers, Ulrich Schwarz, Stefano Vianello, Nachiket Satish Shembekar, and Ramesh Utharala, and all the members of the Aulehla group for the discussions and feedback. We are very grateful to Stefano Vianello for also designing some of the figures. We thank Antonio Politi for providing the Pipeline Constructor macro for automated live microscopy. We also acknowledge the assistance from Maximilian Beckers, Aliaksandr Halavatyi, and Oleksandr Maistrenko in writing scripts for data analysis. We are thankful to Ricardo Henriques for kindly sharing the template that was used to format this preprint. Following EMBL core facilities are acknowledged for their support of this work: the EMBL Advanced Light Microscopy Facility (ALMF), the EMBL Laboratory Animal Resources (LAR), the EMBL Mechanical Design Office, and the EMBL Mechanical Workshop. The work was supported by the European Molecular Biology Laboratory (EMBL), by the German Research Foundation/DFG (CRC/SFB 1324, project number 331351713) to AA, the National Research Council of Canada, Discovery Grant program to Paul François, by the New Frontiers in Research Fund, Exploration program to PF and AA and the European Research Council (ERC) via an ERC Starting Grant (agreement no. 633943) and Consolidator Grant (agreement no. 866537) to Alexander Aulehla.
Ethics
All animal experiments were conducted under veterinarian supervision and after project approval by European Molecular Biology Laboratory, following the guidelines of the European Commission, Directive 2010/63/EU and AVMA Guidelines 2007.
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,943
- views
-
- 238
- downloads
-
- 17
- citations
Views, downloads and citations are aggregated across all versions of this paper published by eLife.
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
-
- Physics of Living Systems
Many proteins have been recently shown to undergo a process of phase separation that leads to the formation of biomolecular condensates. Intriguingly, it has been observed that some of these proteins form dense droplets of sizeable dimensions already below the critical concentration, which is the concentration at which phase separation occurs. To understand this phenomenon, which is not readily compatible with classical nucleation theory, we investigated the properties of the droplet size distributions as a function of protein concentration. We found that these distributions can be described by a scale-invariant log-normal function with an average that increases progressively as the concentration approaches the critical concentration from below. The results of this scaling analysis suggest the existence of a universal behaviour independent of the sequences and structures of the proteins undergoing phase separation. While we refrain from proposing a theoretical model here, we suggest that any model of protein phase separation should predict the scaling exponents that we reported here from the fitting of experimental measurements of droplet size distributions. Furthermore, based on these observations, we show that it is possible to use the scale invariance to estimate the critical concentration for protein phase separation.
-
- Computational and Systems Biology
- Physics of Living Systems
Explaining biodiversity is a fundamental issue in ecology. A long-standing puzzle lies in the paradox of the plankton: many species of plankton feeding on a limited variety of resources coexist, apparently flouting the competitive exclusion principle (CEP), which holds that the number of predator (consumer) species cannot exceed that of the resources at a steady state. Here, we present a mechanistic model and demonstrate that intraspecific interference among the consumers enables a plethora of consumer species to coexist at constant population densities with only one or a handful of resource species. This facilitated biodiversity is resistant to stochasticity, either with the stochastic simulation algorithm or individual-based modeling. Our model naturally explains the classical experiments that invalidate the CEP, quantitatively illustrates the universal S-shaped pattern of the rank-abundance curves across a wide range of ecological communities, and can be broadly used to resolve the mystery of biodiversity in many natural ecosystems.