## INTRODUCTION

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 (*17*–*21*). 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 (*30*–*34*).

## RESULTS AND DISCUSSION

### 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*, *35*–*37*). 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 *f _{i}* with a complementary metal-oxide semiconductor (CMOS) camera. On the basis of these intensities, an individual photochemical feedback

(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

*I*

_{0}. 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

*A*. 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

_{ij}*K*is cycled from low to high values and back. We monitor the synchronization level using the Kuramoto order parameter (

*15*)

(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.

### 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.

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).

### 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 is

(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).

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 *n*_{1} other oscillators, whose phases are in the excitable interval. As a result, the mean field amplitude

contains a spike amplified by a factor *n*_{1}, which will, in turn, induce even more (*n*_{2} > *n*_{1}) 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.

The situation at low

$\overline{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 *K* ≪ *K*_{↑}; 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.

## MATERIALS AND METHODS

### 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)_{3}Cl_{2}) 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 ([H_{2}SO_{4}] = 0.77 M, [NaBrO_{3}] = 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,

${\overline{\mathrm{\omega}}}_{i}(t)={\stackrel{\xb7}{\mathrm{\varphi}}}_{i}(t)$. As we are interested in the time evolution of

${\overline{\mathrm{\omega}}}_{i}$ on time scales comparable with τ_{coup}, we convolve

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, *u _{i}* and

*v*, which are proportional to the concentrations of the HBrO

_{i}_{2}and

reaction intermediates, respectively. Their time evolution is given by

$${\mathrm{\u03f5}}_{1}{\stackrel{\xb7}{u}}_{i}={I}_{i}+(\frac{\mathrm{\alpha}{\mathit{qv}}_{i}}{{\mathrm{\u03f5}}_{3}+1-{v}_{i}}+\mathrm{\beta})\frac{\mathrm{\mu}-{u}_{i}}{\mathrm{\mu}+{u}_{i}}+\mathrm{\gamma}{\mathrm{\u03f5}}_{2}{w}_{i,\mathit{ss}}^{2}+(1-{v}_{i}){w}_{i,\mathit{ss}}-{u}_{i}^{2}-{u}_{i}$$(4)

$${\stackrel{\xb7}{v}}_{i}=2{I}_{i}+(1-v){w}_{i,\mathit{ss}}-\frac{\mathrm{\alpha}{\mathit{qv}}_{i}}{{\mathrm{\u03f5}}_{3}+1-{v}_{i}}$$(5)where

$${w}_{i,\mathit{ss}}=\frac{1}{4\mathrm{\gamma}{\mathrm{\u03f5}}_{2}}(\sqrt{16\mathrm{\gamma}{\mathrm{\u03f5}}_{2}{u}_{i}+{v}_{i}^{2}-2{v}_{i}+1}+{v}_{i}-1)$$(6)represents the steady-state concentration of

${\text{HBrO}}_{2}^{+}$and

$${I}_{i}={I}_{0}+K{\displaystyle \underset{j=1}{\overset{N}{\Sigma}}}{L}_{\mathit{ij}}[{v}_{j}+{\mathrm{\xi}}_{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, *I*_{0} is the background light intensity, *K* is the coupling strength, *L _{ij}* 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” (*u _{i}*) and an “inhibitor” (

*v*) whose time evolution is given by

_{i}(8)

$${\stackrel{\xb7}{v}}_{i}={u}_{i}+a$$(9)where *a* and ϵ represent dynamical parameters and the feedback

(10)acts additively on each node. The symbols *K*, *L _{ij}*, and ξ

*(*

_{i}*t*) have the same meanings as in Eq. 7. All the parameter values are given in tables S1 and S2.

**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.