Rhythmic synchronization and hybrid collective states of globally coupled oscillators

Macroscopic rhythms are often signatures of healthy functioning in living organisms, but they are still poorly understood on their microscopic bases. Globally interacting oscillators with heterogeneous couplings are here considered. Thorough theoretical and numerical analyses indicate the presence of multiple phase transitions between different collective states, with regions of bi-stability. Novel coherent phases are unveiled, and evidence is given of the spontaneous emergence of macroscopic rhythms where oscillators’ phases are always found to be self-organized as in Bellerophon states, i.e. in multiple clusters with quantized values of their average frequencies. Due to their rather unconditional appearance, the circumstance is paved that the Bellerophon states grasp the microscopic essentials behind collective rhythms in more general systems of interacting oscillators.


Results
We start by considering a frequency-weighted Kuramoto model of N globally coupled oscillators: where θ i , ω i , and κ i are respectively the instantaneous phase of the i th oscillator, its natural frequency, and the strength of its coupling to the other oscillators. Here ω i and κ i are randomly chosen variables from a joint probability density of the form G(ω, κ) = g(ω)Γ(κ|ω), where g(ω) is the natural frequency distribution -hereafter assumed symmetric and unimodal -and Γ(κ|ω) is an additional density describing the type of the intrinsic frequency-coupling correlations 36 .
Suppose now that the couplings κ i take binary (for simplicity, integer) values of opposite signs, i.e. κ 1 < 0 and κ 2 > 0. Oscillators are then grouped into two populations 37,38 : contrarians (opposing the system's beat) and conformists (attempting to follow the global rhythm). In this scenario, we further specialize on the following three cases: 1 , , 1 2 κ ω ω ω δ ω ω δ , where the subscripts relate to specific cases), Θ ⋅ ( ) is a Heaviside step function, δ ⋅ ( ) is the Kronecker delta. The three densities (2)-(4) reflect three different realistic strategies where contrarians are gradually flipped into conformists (Initially all oscillators are contrarians): in Case I randomly chosen contrarians are progressively flipped into conformists, in Case II contrarians are ranked according to the magnitude of their natural frequencies |ω i | and then flipped from the largest |ω i | to a given threshold ω 0 , and Case III is the opposite of Case II. As the control parameter p is adiabatically tuned, the system (1) typically undergoes a series of PT's to different coherent states.
The rigorous treatment of model (1) for all the Cases I-III consists of i) performing a linear stability analysis for detecting the thresholds (the forward critical fraction σ p c , σ = 1, 2, 3) at which the incoherent state loses stability in the three Cases I-III, and ii) adopting the Kuramotos self-consistent method 9 for identifying all possible coherent stationary states of the three models, as well as the backward critical thresholds σ p b at which such coherent states lose their stability. The predicted behaviors are then compared with the numerical solutions obtained integrating directly Eq. (1) (Unless otherwise specified, the following stipulations are chosen: (i) the strength of couplings for conformist oscillators is kept fixed to κ 2 = 5; (ii) a Lorentzian frequency distribution g(ω) = γ/π/ [(ω − ω 0 ) 2 + γ 2 ] with γ = 0.05 is adopted; (iii) numerical integrations are performed with a fourth-order Runge-Kutta method with integration time step Δt = 0.01; (iv) the initial conditions for the phase variables are randomly taken; (v) the typical number of oscillators in the ensemble is N = 5 × 10 4 ; (vi) a sufficiently long time interval (much larger than the oscillation period T f = 2π/Ω f ) is used for the average of the order parameter). In the following, we report the details and results of the linear stability analysis and mean-field theory.
Linear stability analysis. In the mean-field form, Eq. (1) can be rewritten as: where r(t) and Ψ(t) are amplitude and phase of the Kuramoto order parameter Here, Z(t) can be interpreted as the collective rhythm produced by the whole population of oscillators. The amplitude 0 ≤ r ≤ 1 measures the phase coherence in the system and Ψ gives the average phase.
In the thermodynamic limit (N → ∞), a density function ρ(θ, t|ω, κ) can be defined, where ρ dθ denotes the fraction of oscillators with natural frequency ω and coupling strength κ, whose phases have values between θ and 1 1 1 n n The right-hand side of Eq. (9) defines a linear operator A as From Eq. (10), it is obvious that the higher Fourier harmonics are neutrally stable to the perturbation, and the stability of the incoherent state depends on the spectrum of Eq. (9). The spectrum of A has both continuous and discrete sets. Following refs 10,11 , the continuous part of the spectrum is the set {−iω : ω ∈ Support(g)}, which is the whole imaginary axis for a Lorentzian frequency distribution (FD). Therefore, the incoherent state is either unstable or neutrally stable. As for the discrete part of the spectrum, one has to seek solutions of the form c 1 (ω, κ, t) = b(ω, κ)e λt , so that the characteristic equation ∂ t c 1 = A(ω, κ)c 1 given in Eq. (11) takes now the form We can then evaluate the latter integral equation self-consistently by defining the auxiliary function  so that Eq. (12) can be solved for b(ω, κ), yielding }. Inserting now the expression for b into Eq. (13), we are led to the characteristic equation where λ is the complex eigenvalue of A except for the points in the set {−iω}. Notice that Eq. (14) relates implicitly p (or ω 0 ), which serves as the control parameter, with the eigenvalue λ. Since the real part of λ determines the stability of the incoherent state, we write Eq. (14) into two equations by letting λ = x + iy, i.e., Based on Eq. (15), one can determine the stability of the incoherent state, and obtain the critical proportion of conformists (p c ) for the forward phase transition.  In contrast with the classical Kuramoto model 10,11 , x 1 has not necessarily to be positive, as κ 1 < 0. As p 1 increases, if x 1 changes from negative to positive, the incoherent state will lose its stability. Imposing the critical condition x 1 → 0, one obtains the critical proportion of conformists for the forward PT as where y j are determined by Eq. (15) in the limit x 1 → 0, and y j cannot be zero. Eq. (15) may have more than one root as x 1 → 0. sup j means that we choose the jth root y j which makes the product |y j |g(y j ) maximal, so that p c 1 corresponds to the foremost critical point for the onset of synchronization. Specifically, for a Lorentzian FD (Unless otherwise specified, the following stipulations are chosen: (i) the strength of couplings for conformist oscillators is kept fixed to κ 2 = 5; (ii) a Lorentzian frequency distribution g(ω) = γ/π/[(ω − ω 0 ) 2 + γ 2 ] with γ = 0.05 is adopted; (iii) numerical integrations are performed with a fourth-order Runge-Kutta method with integration time step Δt = 0.01; (iv) the initial conditions for the phase variables are randomly taken; (v) the typical number of oscillators in the ensemble is N = 5 × 10 4 ; (vi) a sufficiently long time interval (much larger than the oscillation period T f = 2π/Ω f ) is used for the average of the order parameter), one gets y j = {0, ±γ}.
1 , which implies that the critical point for the forward PT is uniquely determined by the coupling strengths. In panel (a) of Fig. 1 Here, the symbol P.V. means the Cauchy principal-value integration within the real line. Considering g(ω) = g(−ω) and Γ 2 (κ|ω) = Γ 2 (κ| − ω), a pair of y c 2 with opposite sign might emerge together. For a Lorentzian FD, one eventually obtains 2 . First, one sees that y c 2 does not exist if κ 2 < 4, suggesting that x 2 → 0 is self-contradictory. Therefore, the real part of λ 2 is either positive or negative when κ 2 < 4. Obviously, x 2 > 0 is physically unreasonable. Therefore, x 2 is supposed to be negative in this case, which implies that the coherent state will not emerge no matter how large the proportion of conformist is. Then, when κ 2 > 4, one can obtain the critical proportion of conformists for forward PT, which is reported in panel (b) of Fig. 1 3 . This can be proved briefly as follows. Because of z 2 z 3 = 1, one has , one straightforwardly obtains Mean-field theory. In order to unveil all possible coherent states in the system (as the proportion of conformists p increases), one has to rely on self-consistence analysis. For stationary coherent state, the mean-field phase Ψ rotates uniformly with frequency Ω, i.e., Ψ(t) = Ωt + Ψ(0). Without loss of generality, one can set Ψ = 0 after an appropriate time shift. In the rotating frame with frequency Ω, the phase difference φ i = θ i − Ψ is introduced, and the model can be written as When replacing summation by integration, the contribution of the phase-locked oscillators to r reads where conformists take the positive sign and contrarians take the negative sign in the first integral. This is because, in a stationary coherent state, conformists attempt to follow the global rhythm of the system (and therefore cos φ i > 0 is always the case), whereas contrarians try to oppose the system's mean-field featuring cos φ i < 0. In contrast to phase-locked oscillators, the drifting oscillators cannot be entrained by the mean-field. Self-consistently, they form a stationary distribution on the circle 11 , i.e., ∂ρ/∂t = 0, and the constant value of the order parameter must be consistent with that implied by Eq. (22). Then, the distribution of the drifting oscillators in the rotating frame is given by s in 2 2 SCIentIfIC RepoRTs | (2018) 8:12950 | DOI:10.1038/s41598-018-31278-9 It is then easy to obtain that, for drifting oscillators, The drifting oscillators have no contributions to the real part of r. However, their contributions to the imaginary part of r should not be neglected. The closed form of self-consistence equations for the real and imaginary parts of r are 2 Taken together, Eqs (24) and (25) provide a closed equation for the dependence of the magnitude r and the frequency of the mean-field Ω on p. We notice that Ω = 0 is always a trivial solution of Eq. (25), corresponding to the π-state reported in refs 37,38 . In this state, the conformist and contrarian oscillators converge to a partially synchronized state where they both satisfy a stationary distribution of phases, and the phase difference between these two clusters is always δ = π. Since Ω = 0 may not be the only solution, there could be more than one value for Ω that satisfies the phase balance equation. Ω ≠ 0 corresponds to the travelling wave (TW) state, in which the two clusters always maintain a constant separation δ π ≠ , and rotate with the same frequency along the unit circle, i.e., they are relatively static with each other.

Case I.
Substituting Γ 1 (κ) into Eqs (24) and (25) yields   From Eq. (29), one can extract the critical proportion of conformists for backward PT as well as the π-state 37,38 theoretically. Interestingly, this solution is independent of the specific form of g(ω) as long as g(ω) is symmetric and centered at 0, which is like the case of two-cluster synchronous state in ref. 17 . The critical proportion of conformists for backward PT (p b 1 ) is the minimum p 1 which satisfies Eq. (29), thus p b 1 can be determined by setting dp 1 /dr = 0 in Eq. (29). When |κ 1 When |κ 1 | > κ 2 , the equation for p b 1 turns out to be tedious so we do not show the exact results here. In Fig. 2(a-c), we have plotted the phase diagram for typical parameters in Case I, as well as the solutions of Eq. (29), which perfectly coincide with the numerical results. In Case I, the typical coherent state is the πstate. In this state, there are four coherent clusters in the system, including a pair of conformist clusters and a pair of contrarian clusters. They are all static and their phases are locked with zero average speed, thus the order parameter is a fixed point on the complex plane. Moreover, the average phase of all conformists and that of all contrarians maintain a constant difference of π. During the backward transition, as p 1 decreases, contrarian clusters first begin to desynchronize, while conformist clusters still keep synchronized. Only when p 1 becomes small enough, the conformist clusters begin to desynchronize and the system finally returns back to the incoherent state. For the case of α 1,2 < 1, Ω 1 is not supposed to be 0, the solution of Eqs (26) and (27) can be solved numerically, corresponding to the TW state. Particularly, in the limit case r → 0 + , one can obtain p c 1 again, which is exactly the same as Eq. (17): (1) as the proportion of conformists p increases. Lorentzian FD with γ = 0.05. From top to bottom, the three rows correspond to Cases I-III, respectively, while from left to right, the three columns correspond to Q < 1, Q = 1, and Q > 1, respectively (Q = |κ 1 |/κ 2 ). Both the forward (red circles) and the backward (blue squares) transitions are studied in an adiabatic way, and the black (dashed pinkish red) curves correspond to theoretical predictions of the stable (unstable) stationary coherent states, including the π-state and the TW states. Other typical non-stationary coherent states, such as the strange π-state, the B-state, and the hybrid-B state, can also be observed in broad parameter ranges.
SCIentIfIC RepoRTs | (2018) 8:12950 | DOI:10.1038/s41598-018-31278-9 where Ω c 1 is the critical mean-field frequency. By Taylor expansion of Eq. (27), we find that Ω c 1 satisfies the following balance equation, when the incoherent state loses its stability, i.e., a stationary TW state will emerge when the system get synchronized. Although a linear stability analysis to this TW state is hard to be performed, its stability can still be studied through numerical simulations. It is found that the TW state predicted by the mean-field theory turns out to be unstable. Taking a Lorentzian FD, Eqs (26) and (27) can be simplified to  Numerical simulations suggest that this bifurcation is supercritical and unstable. Therefore, above p c 1 the incoherent state (α 1,2 < 0) loses its stability and the system jumps to π-state (α 1,2 > 0) predicted by Eq. (29) because the TW state is unstable, i.e., a first-order synchronization transition takes place. In our numerical studies, it is found that in the backward PT, however, the π-state does not always return to incoherence directly as shown in Fig. 2(a-c), and rather bifurcates continuously (for small enough κ 1 ) towards a novel stationary state, here called the strange π-state. For example, in Fig. 3(a), we plot the phase diagram for κ 1 = −2 and κ 2 = 5. As p 1 decreases, when the π-state loses stability, instead of returning back to incoherence state directly, the system bifurcates to the strange π-state. In Fig. 3(b), we show the microscopic characterization of this state. It is found that in this state, the contrarian clusters have desynchronized, while the conformist clusters still maintain complete coherence. The contrarians maximize their distance from the conformists' clusters, resulting in an averaged phase difference of π. Physically, this can be understood as follows. When |κ 1 | is significantly smaller than κ 2 , contrarians are less affected by the mean field, and thus are easier to desynchronize as p decreases. Furthermore, as Ω 1 = 0, we can predict this kind of state through mean-field method. Now, one should remove the contribution of contrarian clusters in Eq. (26), and get κ As shown in Fig. 3(a), the theoretical predictions agree perfectly with the numerical results. With Lorentzian FD, the forward threshold p c 1 (for the transition from incoherence) always exceeds the backward one p b 1 where π or strange-π states lose their stability. With uniform distributions, instead, one can trigger (by appropriately choosing κ 1 ) the forward transition so that < p p c b 1 1 . For example, we focus on the case of uniform distribution g(ω) = 1/2, ω ∈ (−1, 1), and |κ 1 | = κ 2 = 5. From Eqs (30) and (31), the exact expression for critical proportion of conformists is π π = + < = .  Fig. 4(a)]. In this regime, oscillators split into four coherent clusters, two for each population. Like π-states, different clusters maintains a phase difference of π between each other, but oscillators within each clusters are neither phase-nor frequency-locked [ Fig. 4(b)]. In fact, they evolve with different periodic patterns [ Fig. 4(c)], and correlate with each other so that their average frequencies lock to binary constants ±Ω f , forming a staircase structure [ Fig. 4(b3)]. Henceforth, both populations self-organize their phases as B's 23,24 , resulting in a macroscopic rhythmic behaviour [ Fig. 4(d)]. Numerically, we find that as p 1 increases in the regime from p c 2. Case II.
In this case, we substitute Γ 2 (κ|ω) into Eqs (24) and (25) and obtain  Using the same methods developed in Case I, one can identify the entire π-state. Interestingly, the result is exactly the same as Eq. (29). For the case of |κ 1 r| < 1 and κ 2 r < 1, Ω 2 is not supposed to be 0, the solution of Eqs (34) and (35) can be solved numerically, corresponding to the TW state. Particularly, in the limit case r → 0 + , one can obtain p c 1 2 1 As for the TW state, Eqs (34) and (35) can be simplified to  Numerical results highlight that the bifurcation is supercritical and the TW state is unstable, like the TW state in Case I. In Fig. 2(d-f), we have plotted the phase diagram in Case II, as well as the solutions of Eq. (36). However, the predictions of π-state made by the mean-field theory do not coincide with the numerical results, suggesting that the stationary π-state is unstable, which has been further verified by numerical method. In Case II, forward thresholds p c 2 smaller than p b 2 commonly appears, leading generically to the emergence of non-stationary coherent states. Typically, the system undergoes two PT's: one at p c 2 , where the incoherent state loses stability and a B-state discontinuously emerges with hysteresis, and another at > ′ p p c c 2 2 , where the Bellerophon state becomes unstable and bifurcates to a novel rhythmic state with a PT whose type changes from discontinuous to continuous as the coupling ratio Q ≡ |κ 1 |/κ 2 crosses one [ Fig. 2(d-f)]. Physically, Q denotes the ratio of entrainment received by the contrarians and conformists from the mean-field. In Fig. 5, one such novel rhythmic state is illustrated. In this state, there are six (two pairs of contrarians' and one pair of conformists') coherent clusters [ Fig. 5(a) Being a superposition of odd and even B's, we refer to such a macroscopic rhythmic state as hybrid-Bellerophon state. To conclude, there are two synchronization PT's in Case II. Further numerical simulations show that when κ 2 is a constant, with the increasing of |κ 1 |, the width of hysteresis of the first PT almost remains the same, but the width of hysteresis of the second PT becomes smaller and smaller. When Q = |κ 1 |/κ 2 < 1, the two staircases of Bellerophon state are one in the forward PT and the other in the backward PT. When Q = |κ 1 |/κ 2 = 1, the two staircases start joining with each other. Until Q = |κ 1 |/κ 2 > 1, the two staircases coincides with each other and the hysteresis disappears eventually. Although the second PT becomes a  Fig. 2(d): distribution of θ i (a1), θ  i (a2), and θ 〈 〉  i (a3) vs. ω i . One can identify six (two conformist and four contrarian) coherent clusters of oscillators, whose phases are self-organized as in pure B's and in oscillating π-states. (b) The order parameters for all oscillators with positive (blue) and negative (red) frequencies, and their frequency average (pinkish red). The insets plot the typical rhythmic behaviors of the order parameter restricted to the clusters of oscillators in B (b1) and in oscillating-π (b2) states. Panels (b3 and b4) report the rhythmic evolution in time of the global order parameters r(t) and Ψ(t), respectively. state to π-state. When Q > 1, the results turn out to be the same as Q = 1, except that the TW state can be stable and it vanishes at > ′ p p c b 3 3 with a discontinuous transition to the π-state as p 3 increases. Further numerical simulations show that if |κ 1 | is large enough, the Bellerophon state will emerge from subcritical bifurcation to supercritical bifurcation, which causes the first PT change from a discontinuous one to a continuous one. In conclusion, with the increasing of |κ 1 |, the first PT changes from first-order to second-order, meanwhile the second PT emerges and changes from the second-order to the first-order.

Discussion
While in the classical Kuramoto model, phase oscillators are globally and homogeneously coupled to the mean field, in some practical situations such a coupling may instead be heterogeneous. In ref. 23 , a frequency-weighted Kuramoto model was studied, and a novel quantized, time-dependent, and clustered state (the Bellerophon state) was revealed. In ref. 24 , Bellerophon states were again observed in coupled conformist and contrarian oscillators.
Here, we combined frequency-weighted coupling strengths with positive and negative interactions. As compared with refs 23,24 , such a higher order of heterogeneity gives rise to a much richer scenario, which includes different regions of bi-stability. Furthermore, we identified the multiple phase transition character of the root to synchronization, and gave evidence of the generic emergence of macroscopic rhythmic regimes where oscillators self-organize as in Bellerophon states. Together with the Bellerophon states of refs 23,24 , novel collective phases (namely, the strange πand the hybrid B-states) were revealed. Unlike precedent studies, where periodic synchronization behaviors were observed in the presence of an external periodic driving [40][41][42] , the collective rhythms here reported emerge spontaneously, as soon as the forward critical threshold precedes the backward one.
The strange πand hybrid B-states observed in the present model can be regarded as a transitional state between the incoherent state (full asynchrony) and the π-state (full synchrony). On the one hand, the control parameter is not strong enough to completely entrain the system into the π-state, on the other hand it is large enough to maintain certain correlations among the instantaneous frequencies of oscillators. As a compromise of this competition, the instantaneous frequencies of oscillators are not locked but their average frequencies are locked to certain constant.
The current model of conformist and contrarian may describe neuron systems, political systems, or economical systems. The rich synchronization phenomena and the nontrivial coherent states observed in this work could help us better understand the collective behaviors in such systems. For instance, diverse neuro-degenerative disorders, like Parkinson's disease, have been proved to be correlated to the spontaneous emergence of global neuronal oscillations, which makes the present model appealing for a better understanding of their dynamical origins. Moreover, due to their ubiquitous appearance, our study paves the possibility that the microscopic, quantized features of B's are actually the fundamental building blocks behind spontaneous emergence of collective rhythms in more general systems of interacting oscillators, opening a theoretical challenge in the study of rhythmic synchronization.