Research ArticlePHYSICS

First-order synchronization transition in a large population of strongly coupled relaxation oscillators

See allHide authors and affiliations

Science Advances  23 Sep 2020:
Vol. 6, no. 39, eabb2637
DOI: 10.1126/sciadv.abb2637


Onset and loss of synchronization in coupled oscillators are of fundamental importance in understanding emergent behavior in natural and man-made systems, which range from neural networks to power grids. We report on experiments with hundreds of strongly coupled photochemical relaxation oscillators that exhibit a discontinuous synchronization transition with hysteresis, as opposed to the paradigmatic continuous transition expected from the widely used weak coupling theory. The resulting first-order transition is robust with respect to changes in network connectivity and natural frequency distribution. This allows us to identify the relaxation character of the oscillators as the essential parameter that determines the nature of the synchronization transition. We further support this hypothesis by revealing the mechanism of the transition, which cannot be accounted for by standard phase reduction techniques.


Since C. Huygens’s discovery of synchronization in coupled pendulum clocks in 1665 (1), emergent synchronization of oscillating units has been observed in a myriad of natural systems, including firing neurons (2, 3), contracting cardiomyocytes (4), quorum-sensing bacteria (5), beating cilia (6), rainfall extremes (7), and neutrino oscillations (8). In addition, synchronization underpins the dynamics of a variety of technological systems such as power grids (9), bridge instabilities (10), traffic patterns (11), and lasers (12, 13). The process of synchronization can be interpreted from a statistical mechanics perspective as a nonequilibrium phase transition, where a synchronized state emerges from an incoherent one as the coupling between the oscillators is increased. The seminal works of Winfree and Kuramoto (14, 15) and subsequent experiments (16) have shown that in populations of weakly coupled phase oscillators with a unimodal natural frequency distribution, such a transition proceeds continuously and reversibly in a second-order fashion. Under special conditions, the synchronization transition can also be of first-order type (1721). In this case, synchronization abruptly ensues at a critical coupling strength K but disappears below a coupling strength K < K, resulting in hysteresis (22, 23). Such a discontinuous transition has been hypothesized to play a role in the onset of anesthesia-induced unconsciousness (24), epileptic seizures (25), acoustical signal transduction in the cochlea (26), hypersensitivity in chronic pain (27), and memory processes (28).

The simplifying assumptions underlying the Kuramoto phase oscillator model, however, render it unsuitable for understanding systems where strongly coupled relaxation oscillators prevail. To experimentally investigate the onset of synchronization in large ensembles with N = 200 and N = 1000 relaxation oscillators, we use the well-studied Belousov-Zhabotinsky (BZ) chemical reaction (29). Despite intrinsic differences in the underlying microscopic mechanisms with biological neurons, this chemical reaction shows qualitatively identical emergent behavior: spiked slow-fast oscillations in the time traces of the concentrations and phase-dependent excitable response to external perturbations (3034).


Chemical micro-oscillators

We synthesized individual oscillatory units from ion-exchange resin beads saturated with a photosensitive BZ reaction catalyst, ruthenium(II)-tris(2,2′-bipyridine-dimethylene)-chloride (30, 3537). The particles were then immobilized under a hydrogel layer on an acrylic plate and immersed in a catalyst-free reaction solution. This resulted in a reservoir of uncoupled chemical micro-oscillators. The oscillators can be coupled together photochemically with the experimental setup presented in Fig. 1A. During an oscillation cycle, the catalyst varies periodically between its fluorescing and nonfluorescing oxidation states. The phase of each oscillator can thus be monitored by optically measuring the oxidized catalyst concentration via its fluorescence intensity fi with a complementary metal-oxide semiconductor (CMOS) camera. On the basis of these intensities, an individual photochemical feedbackIi(t)=I0+KΣj=1NAij[fj(t)fi(t)](1)is calculated and projected on each photosensitive micro-oscillator with a spatial light modulator. The natural frequencies ωi of all oscillators are measured at the start of each experiment under a uniform background light of intensity I0. Subsequently, a suitable subset of the oscillators is selected to obtain a desired frequency distribution. Different network connectivities can be implemented by choosing an appropriate adjacency matrix Aij. Figure 1B depicts a typical camera image of the fluorescing bead reservoir during an experiment together with a chosen network graph superimposed to highlight the connected oscillators. During each experiment, the coupling strength K is cycled from low to high values and back. We monitor the synchronization level using the Kuramoto order parameter (15)R=1N|Σj=1Neiϕj(t)|t(2)where ⟨…⟩t denotes the time average and ϕj(t) represents the phase of the j-th oscillator, which is calculated by linear interpolation between consecutive firing events. The order parameter R ranges from 0, when the phases are incoherent, to 1, where all phases align perfectly. The onset of synchronization can be captured from the dependence of the order parameter on the coupling strength.

Fig. 1 Coupled photochemical oscillators.

(A) Experimental setup. A thermostatted open reactor, hosting 2600 chemical oscillators, is spectrophotometrically monitored in fluorescence light (λ > 550 nm) with a CMOS camera. Recorded light intensities fi determine the photochemical feedback Ii, which is applied with a spatial light modulator. (B) Camera image of the fluorescing oscillator reservoir. The connectivity between oscillators is overlaid in blue. (C) Hysteresis loop of the Kuramoto order parameter R in the case of N = 1000 all-to-all coupled oscillators with unimodally distributed natural frequencies (SD, σω = 0.012 rad/s).

Role of natural frequency distribution

An experimentally recorded order parameter curve is shown in Fig. 1C for the case of N = 1000 globally coupled oscillators with a unimodal distribution of natural frequencies. Upon increasing K, there is an abrupt transition to a highly synchronized state at a critical value of the coupling strength, K = 0.4. Once this phase is formed, it remains stable until the coupling strength K is decreased below K to a smaller value, K = 0.1. This hysteresis cycle is characteristic of a first-order phase transition. The evolution of the instantaneous frequencies for each oscillator (Fig. 2) confirms the hysteretic nature of the transition: For a time-reversal symmetric protocol of the coupling strength K (Fig. 2A), the oscillators’ frequencies evolve asymmetrically with respect to time reversal (Fig. 2B). Moreover, the fluorescence time plots (Fig. 2, C to E) indicate that the emergence of an in-phase synchronized state is preceded by the formation of antiphase clusters. In addition, the frequency of the oscillators in the in-phase synchronized state is approximately given by the frequency of the fastest unperturbed oscillators. This contrasts with the Kuramoto model, which predicts that the oscillators phase-lock at a frequency equal to the average natural frequency of the entire population.

Fig. 2 Experimental observation of hysteresis in the case of N = 1000 all-to-all coupled oscillators with normally distributed natural frequencies.

(A) Time protocol for the coupling strength K. (B) Time evolution of the instantaneous frequencies (ω¯i). The color of each line corresponds to the natural frequency of the nodes (ωi), respectively. The oscillators synchronize in-phase at K = 0.4 but transition back to incoherence at K = 0.1 < K (see also Fig. 1C). (C to E) Fluorescence value (fi) plots for clustering, synchronized, and incoherent states as observed during the cycle. The oscillators are indexed (i) in order of increasing natural frequencies. Synchronization from incoherent initial conditions (E) proceeds with the formation of antiphase clusters (C), which delay the onset of synchronization (D) leading to hysteretic behavior. The antiphase clustering state is characterized by a higher degree of frequency coherence for low-frequency oscillators.

We investigated the robustness of the first-order synchronization transition with regard to the frequency distribution in an all-to-all coupled network of BZ oscillators with frequencies drawn from a bimodal distribution (18, 38). Incidentally, this reveals the hierarchy of the emergent synchronization dynamics (Fig. 3). While cycling the coupling strength K up and down (Fig. 3A), the evolution of the frequencies is asymmetric in time, which again indicates hysteretic behavior (Fig. 3B). At the beginning of the experiment, the coupling strength is small, and the oscillators are desynchronized and incoherent. Upon a slight increase of the coupling strength, the fluorescence time plot (Fig. 3C) shows the presence of approximately antiphase clusters (α) in the subpopulations associated with fast and slow intrinsic frequencies. Oscillator heterogeneity induces intercluster switching, so there is no perfect frequency synchronization (36). Before the onset of global synchronization, the low-frequency subpopulation achieves in-phase synchronization, while the high-frequency group remains approximately antiphase but displays an average increase in instantaneous frequency (β). Once the synchronized state is established at t ≈ 5300 s (K = 0.72, blue dashed line), the population oscillates with the natural frequency of the fastest oscillators (γ). This regime is characterized by almost perfect phase alignment, with the fast oscillators entraining the entire population (Fig. 3D). The destabilization of the synchronized state at t ≈ 7800 s (K = 0.56, orange dashed line) is initiated by the loss of frequency coherence in the slower subpopulation (δ). At the end of the experiment, we recover the fully incoherent state that is also observed in the beginning for very low (K < 0.2) coupling strengths (Fig. 3E).

Fig. 3 Experimental observation of hysteresis for a bimodal frequency distribution.

(A) Time protocol for coupling strength K. (B) Time evolution of the instantaneous frequencies (ω¯i) of N = 200 oscillators. The color of each line encodes the natural frequency of the nodes (ωi), respectively. The oscillators synchronize in-phase at K = 0.72 upon increasing K but transition back to incoherence at K = 0.56 < K upon decreasing K. (C to E) Individual fluorescence values for clustering (C), synchronized (D), and incoherent states (E) during the experiment. The oscillators are indexed (i) in order of increasing natural frequencies. Synchronization from incoherent initial conditions proceeds with the formation of approximately antiphase clusters, which delay the onset of synchronization leading to hysteresis.

Synchronization mechanism for relaxation oscillators

To gain more insight in the mechanism of the observed first-order synchronization transition, we performed numerical simulations with an established model of the BZ chemical kinetics (30, 39). Figure 4 shows the comparison of the hysteretic order parameter curves between experiments and simulations in the case of globally coupled oscillators with unimodal and bimodal distributions. In both cases, the ascending branch of the hysteresis loop is characterized by a persistence of low order parameter values. The detailed inspection of the collective node dynamics in the case of unimodal (Fig. 2 and fig. S1) and bimodal (Fig. 3) frequency distributions reveals that the abrupt emergence of an in-phase synchronized state is preceded by the formation of antiphase clusters. The antiphase synchronized state effectively suppresses the onset of in-phase synchronization, resulting in hysteretic behavior. The critical coupling strength for in-phase synchronization corresponds to the point where the antiphase state vanishes.

Moreover, we found that the collective transition for many oscillators can be described by a reduced model of just two coupled identical oscillators. In this case, direct simulations reveal that, up to a certain critical coupling strength, both in-phase and antiphase states coexist (fig. S2). Beyond a critical point, the antiphase state becomes unstable. An estimate for the critical coupling strength, K* = 4.3 × 10−3, agrees quantitatively in the case of a unimodal distribution of natural frequencies and qualitatively for a bimodal distribution, where the effects of frequency distribution are more pronounced. The peculiarities of the latter case (the sequence of α, β, γ, and δ states) are also accurately reproduced by the simulations in both the order parameter curves and fluorescence intensity plots (fig. S3).

Network topology

To further demonstrate the robustness of the first-order synchronization transition for relaxation oscillators, we investigated the role of network connectivity. We chose two paradigmatic random networks, the Barabási-Albert and the Erdős-Rényi graphs, where the natural frequencies depended linearly on the corresponding node degrees (23). Our experiments and simulations with chemical relaxation oscillators show that there is a discontinuous first-order transition to in-phase synchronization, with hysteresis irrespective of the chosen network connectivity models (fig. S4). This suggests that the occurrence of abrupt synchronization in the case of relaxation oscillators depends only weakly on the underlying network topology. Moreover, a close inspection of the collective node dynamics again shows the presence of approximately antiphase clusters suppressing the onset of in-phase synchronization (figs. S5 and S6).

Relaxation dynamics and time scale separation

To validate our hypothesis on the role of the relaxation character determining the type of synchronization transition, we use the FitzHugh-Nagumo (FHN) model (33), a canonical model for relaxation oscillations in the context of neuronal excitability. Varying the time scale separation, parameter ϵ allows for tuning between harmonic and slow-fast relaxation oscillations (Fig. 5). Every oscillator can be characterized by its phase response curve (PRC), which encodes the resultant phase change Δϕ due to a short perturbation applied at a phase ϕ (33). We observe that the PRC evolves from a linear to a nonlinear dependence on the perturbation amplitude A: While for negligible time scale separation, the PRC is roughly sinusoidal (Fig. 5A), for strong time scale separation, it is a discontinuous function whose jump point ϕ*(A) shifts to smaller phases with increasing perturbation strength A (Fig. 5B), thus enlarging the excitable interval. A simple approximation for the PRC isP(ϕ,A)={0ϕ<ϕ*(A)2πϕϕϕ*(A)(3)

This PRC encodes phase-dependent excitability: At early phases, an oscillator is refractory, as perturbations do not affect it. However, in the excitable window with phases above ϕ*(A), perturbations immediately trigger a new spike (31, 32). We observe identical behavior for our chemical oscillators (Fig. 5C).

Simulations of a globally coupled network of N = 200 FHN oscillators indicate the presence of a continuous transition without hysteresis for negligible time scale separation (ϵ = 0.3) and a discontinuous transition with hysteresis for a pronounced time scale separation (ϵ = 10) (Fig. 5, D and E). In agreement with the experiments, individual dynamics of oscillators in the bistable region reveal antiphase states and in-phase states during the up- and down-sweep, respectively. This behavior is further confirmed in a reduced two-oscillator model, which shows the presence of stable antiphase states only in the case of pronounced time scale separation (fig. S7).

Fig. 4 First-order synchronization transitions in all-to-all coupled oscillator populations.

Hysteresis curves for the order parameter (R) with increasing (blue) or decreasing (orange) coupling strength in chemical experiments (A and B) and numerical simulations (C and D). We consider N = 200 globally coupled oscillators with unimodal (A and C) or bimodal (B and D) natural frequency distributions. Because of finite-size effects, RO(1/N) in the unsynchronized phase (low K).

The mechanism underlying the hysteretic synchronization transition for relaxation oscillators is thus revealed to be deeply rooted in the nonlinear behavior of the PRC (Fig. 5F). In an initially incoherent population of coupled oscillators at low K, a firing event from an oscillator (the pacemaker) triggers firing events in n1 other oscillators, whose phases are in the excitable interval. As a result, the mean field amplitude Kv¯(t) contains a spike amplified by a factor n1, which will, in turn, induce even more (n2 > n1) oscillators to fire, since the excitable interval is now even larger (see Fig. 5, B and C). Thus, synchronization for relaxation oscillators can be seen to be driven by pacemakers, which are abundant due to the random initial conditions.

Fig. 5 Relaxation character determines the order of the synchronization transition.

(A and B) PRCs of the FHN model for time scale separation parameter ϵ = 0.3 (A) with perturbation strengths of 0.4, 0.6, and 0.8 (light-orange, purple, and light blue, respectively) and ϵ = 10 (B) with perturbation strengths of 0.5, 1.0, and 1.5 (light orange, purple, and light blue, respectively). The phase response curves are normalized by the respective perturbation strengths in (A). (C) Phase response curves determined from chemical experiments (dots) together with functions fitted with model (3) for perturbing light intensities of 0.01 (light orange), 0.06 (purple), and 0.25 mW cm−2 (light blue). (D and E) Order parameter curves in the case of N = 200 globally coupled FHN oscillators with unimodally distributed natural frequencies for ϵ = 0.3 (D) and ϵ = 10 (E). (F) Schematic representation of the mechanism of first-order synchronization via antiphase cluster formation: time evolution of two oscillators from different subpopulations (v1,v2) and their excitable intervals (hatched regions) together with the mean amplitude (v¯) in antiphase (★) and in-phase synchronized (◆) states.

The situation at low v¯, where the excitable window of oscillators is short, will generally lead to the emergence of synchronized subpopulations that are effectively uncoupled. The simplest example is the antiphase state that is formed by two subpopulations with opposing phases: A firing event occurring in one subpopulation may entrain oscillators with similar phases belonging to the same subpopulation; however, the event cannot illicit spiking of the oscillators belonging to the other subpopulation, since they are in the refractory interval and cannot be excited.

The width of the excitable interval grows with increasing coupling strength K until it covers most of the oscillation cycle (ϕ*(A) < π) at the transition point (K), as illustrated in the bottom panel of Fig. 5F. In this situation, oscillators residing in one subpopulation can trigger responses in oscillators belonging to the other subpopulation, thus giving rise to in-phase synchronization across the entire system.

In the K-decreasing branch, the in-phase synchronized state can persist for KK; the already synchronized oscillators require only a small excitable interval to maintain their synchronized state. Last, the loss of synchronization at K happens almost instantaneously, with the exact transition point depending on the frequency distribution.


Our theoretical and experimental analysis points to the fact that simple phase models, while analytically tractable, can fail to capture emergent phenomena in ensembles of strongly coupled relaxation oscillators correctly. In the systems studied here, the type of synchronization transition is much more sensitive to the relaxation character than the frequency distribution or network connectivity as in the case of the weakly coupled phase oscillator model. Because of the ubiquitous nature of these systems, we expect the presented mechanism to play an essential role in understanding further oscillatory systems, such as next-generation neuromorphic photonic devices (40) and optogenetically addressable neural networks (41), as well as in medical therapies of Tinnitus and Parkinson’s disease based on neural desynchronization strategies (2, 3).


Preparation of the chemical oscillators

Cation-exchange resin beads (75 to 150 μm in diameter; DOWEX WX4 100-200) were sieved to obtain a narrow size distribution (106 to 112 μm). One gram of sifted beads was placed in 5 ml of water, and under constant vortex mixing, 15 ml of ruthenium(II)-tris(2,2′-bipyridine-dimethylene)-chloride (Ru(dmpy)3Cl2) catalyst solution (1.66 mM) was added slowly over the course of 10 min. Mixing was continued for 48 hours until a homogeneous (verified by color saturation measurements) catalyst loading of resin (2.5 × 10−5 mol/g) was achieved in all beads. The catalyst-soaked beads were placed on a drilled acrylic plate with a grid of 64 × 44 cylindrical wells (diameter, 200 μm; depth, 150 μm; separation, 400 μm) and then evenly distributed with a fine brush by applying a water surfactant solution (0.5% Triton X-100). After 3 hours, the beads were sealed with liquid silica hydrogel that solidified more than 30 min (30). The chemical oscillations are started when the acrylic plate is submerged in a Belousov-Zhabotinsky reaction solution ([H2SO4] = 0.77 M, [NaBrO3] = 0.51 M, [NaBr] = 0.08 M, [malonic acid] = 0.16 M).

Data analysis

During an experimental run, the coupling strength K is cycled from low to high values and back. Each value of K is maintained for τcoup = 550 s at the end of which the Kuramoto order parameter R is determined by averaging over the past interval τav = 100 s according to Eq. 3. Note that τcoup is large enough to ensure that the oscillators reach the (de-)synchronized steady state at the end of the respective coupling stage. Because of the phase-resetting nature of the relaxation oscillators, the oscillators synchronized within one or two periods after K is reached (characteristic time scale of 1/ωi~100 s). In a similar manner, desynchronization happens abruptly with the order parameter relaxing to its steady-state value on a characteristic time scale, which is inversely proportional to the spread of natural frequency (1/σω~150 s).

The instantaneous frequency of each chemical oscillator is computed directly from its temporal phase, ω¯i(t)=ϕ·i(t). As we are interested in the time evolution of ω¯i on time scales comparable with τcoup, we convolve ω¯i(t) with a Gaussian of width σt = 15 s in Figs. 2 and 3 and figs. S1, S5, and S6.

Numerical simulations

For simulating the chemical kinetics of the BZ reaction, we use the nondimensionalized Zhabotinsky-Buchholtz-Kiyatkin-Epstein (ZBKE) model, which is further modified to account for photochemical effects due to coupling with light (30, 39). The state of the i-th node is given by two variables, ui and vi, which are proportional to the concentrations of the HBrO2 and Ru(dmbpy)33+ reaction intermediates, respectively. Their time evolution is given byϵ1u·i=Ii+(αqviϵ3+1vi+β)μuiμ+ui+γϵ2wi,ss2+(1vi)wi,ssui2ui(4)v·i=2Ii+(1v)wi,ssαqviϵ3+1vi(5)wherewi,ss=14γϵ2(16γϵ2ui+vi22vi+1+vi1)(6)represents the steady-state concentration of HBrO2+ andIi=I0+KΣj=1NLij[vj+ξj(t)](7)is the light intensity projected on the i-th node. In Eqs. 4 to 7, ϵ1, ϵ2, ϵ3, α, β, γ, μ, and q are kinetic and time scale parameters, I0 is the background light intensity, K is the coupling strength, Lij is the Laplacian matrix of the corresponding coupling network, while ξi(t) denotes white Gaussian noise.

We use the FHN model (33) to study the influence of the relaxation character of the oscillators on the order of the synchronization transition. In a similar fashion to the ZBKE model, each oscillator is characterized by two dynamical variables: an “activator” (ui) and an “inhibitor” (vi) whose time evolution is given byu·i=ϵ(uiui33vi+Ii)(8)v·i=ui+a(9)where a and ϵ represent dynamical parameters and the feedbackIi=KΣj=1NLij[uj+ξj(t)](10)acts additively on each node. The symbols K, Lij, and ξi(t) have the same meanings as in Eq. 7. All the parameter values are given in tables S1 and S2.


Supplementary material for this article is available at

This is an open-access article distributed under the terms of the Creative Commons Attribution-NonCommercial license, which permits use, distribution, and reproduction in any medium, so long as the resultant use is not for commercial advantage and provided the original work is properly cited.


Acknowledgments: We acknowledge discussions with M. Bär, S. Yanchuk, and W. J. A. Martin. We thank U. Künkel for support in the preparation of the experiments. Funding: J.F.T. and H.E. thank SFB 910 and GRK 1558. D.C. and J.F.T. thank DAAD RISE 2017, and D.C. thanks Trinity College, Cambridge for Trinity Summer Studentship Scheme 2017. Author contributions: D.C. and J.F.T. devised the study and did experimental and numerical work. D.C., J.F.T., E.A.M., and H.E. discussed the results, commented extensively on the manuscript at all stages of preparation, and jointly wrote the manuscript. Competing interests: The authors declare that they have no competing interests. Data and materials availability: All data needed to evaluate the conclusions in the paper are present in the paper and/or the Supplementary Materials. Additional data related to this paper may be requested from the authors.

Stay Connected to Science Advances

Navigate This Article