Hybrid-type synchronization transitions: where marginal coherence, scale-free avalanches, and bistability live together

The human cortex is never at rest but in a state of sparse and noisy neural activity that can be detected at broadly diverse resolution scales. It has been conjectured that such a state is best described as a critical dynamical process -- whose nature is still not fully understood -- where scale-free avalanches of activity emerge at the edge of a phase transition. In particular, some works suggest that this is most likely a synchronization transition, separating synchronous from asynchronous phases. Here, by investigating a simplified model of coupled excitable oscillators describing the cortex dynamics at a mesoscopic level, we investigate the possible nature of such a synchronization phase transition. Within our modeling approach, we conclude that -- in order to reproduce all key empirical observations, such as scale-free avalanches and bistability, on which fundamental functional advantages rely -- the transition to collective oscillatory behavior needs to be of an unconventional hybrid type, with mixed features of type-I and type-II excitability, opening the possibility for a particularly rich dynamical repertoire.


I. INTRODUCTION
Neurons in the cerebral cortex fire in a rather irregular and sparse way, even in the absence of external stimuli or tasks [1][2][3]. This activity is also manifested at large scales in the so-called resting-state networks [4,5]. Understanding the origin and functionality of such an energetically-costly fluctuating "ground state" is a fundamental question in neuroscience, essential to shed light on how the cortex processes and transmits information [6][7][8][9][10].
Two twin sides of spontaneous neuronal activity are synchronization and avalanches. Depending mostly on cortical regions and functional states, diverse synchronization levels across a continuum spectrum are observed. Both synchronous and asynchronous states are retained to be essential for diverse aspects of information processing; e.g., neuronal synchronization is at the root of collective oscillatory rhythms, a crucial aspect for information transmission between distant areas [11][12][13], while asynchronous states also play key roles for information coding [14]. Accumulating evidence -including results from the Human Brain Project [15]-suggests that the ground state of spontaneous activity of a healthy cortex lies close to the edge of a synchronization phase transition, neither exceedingly synchronous nor overly incoherent, allowing for transient and flexible levels of neural coherence, as well as a very-rich dynamical repertoire [10,15,16]. Indeed, abnormalities in the synchronization level are linked to pathologies such as e.g. Parkinsonian disease (excess) and autism (deficit) [17][18][19].
These empirical observations triggered the development of the "criticality hypothesis", conjecturing that the cortex might extract crucial functional advantages -e.g., large sensitivity and optimal computational capabilitiesby operating close to a critical point [20,23,24,[29][30][31][32]. However, there is still no consensus on what type of phase transition is required for such a critical behavior; a critical branching process, as mentioned above, is the most common and straightforward interpretation, but alternative scenarios have been put forward [27,[33][34][35][36]. Among these, a Landau-Ginzburg theory has been recently proposed to describe cortical dynamics at a mesoscopic level suggesting that the cortex might operate in a regime close to the edge of a synchronization phase transition at which, quite remarkably, scale-free avalanches emerge in concomitance with incipient oscillations [37]. This scenario that had been previously suggested in the literature [35,[37][38][39][40][41][42] and that is supported by empirical evidence [15,[43][44][45]. In spite of these advances, a minimal model -simpler than the analytically un-tractable Landau-Ginzburg theory-capturing the gist of such a synchronization phase transition and allowing for indepth theoretical understanding is still missing.
Here we pose the following questions: can scale-free avalanches possibly occur at the critical point of the canonical model for phase synchronization? If not, which is the minimal model for synchronization able to accommodate them? What type of phase transition does it exhibit? Answering these questions will pave the road for the more ambitious goal of constructing a statistical mechanics of cortical networks, shedding light on the collective states -with crucial functional roles-that they sustain, as well as advancing our general understanding of phase transitions.

II. LACK OF SCALE-FREE AVALANCHES IN THE ANNEALED KURAMOTO MODEL
We start by analyzing the canonical model for phase synchronization, customarily used in neuroscience [46][47][48][49] and other fields [50]: the Kuramoto model [50][51][52][53]. However, given that -as revealed by the Landau-Ginzburg model [37]-node heterogeneity is not an essential ingredient to generate avalanches, and that, on the other hand, noise is inherent to neural dynamics, we consider here the annealed version of the Kuramoto model with homogeneous frequencies [49,50], i.e.
where the phase ϕ j (t) describes the dynamical state of the j-th oscillatory unit, with j ∈ [1, N ], ω is a common intrinsic frequency, η j (t) a zero-mean unit-variance Gaussian white noise with amplitude σ, and J is the coupling strength with all the neighbours on, e.g., a fullyconnected network [50][51][52][53][54]. As is well known, Eq.(1) exhibits a synchronization phase transition where the synchronization (Kuramoto) order parameter, Z = e iϕ , experiences a (supercritical) Hopf bifurcation from a fixed point to a limit cycle, revealing the emergence of collective oscillations [51,53].
In what follows, we set to analyze if this model (and related models) can exhibit scale-free outbursts of activity. For this purpose, we define and measure avalanches following the standard protocol in neuroscience [20,23], which relies on the identification of individual-unit "spikes" (in the case of models consisting of oscillators "spikes" can be identified in an effective way as crossings over a given phase value) as well as on the definition of a time discretization window, needed to cluster close-in-time spikes together [20] (details of the protocol are carefully explained in Appendix A; see also [55]). Extensive computational simulations (reported in detail in Appendix B) reveal that neither at the critical point of Eq.(1) nor around it scale-free avalanches can be found; i.e. P (S) and P (T ) always show exponential decays, even if other standard quantities are known to exhibit scaling for different versions of the Kuramoto model [56][57][58].

III. HYBRID-TYPE SYNCHRONIZATION EMERGES FROM EXCITABLE OSCILLATORS
To search for a better suited minimal model for synchronization with scale-free avalanches, we scrutinize the Landau-Ginzburg theory (LG) in [37], known to exhibit these two features in concomitance. In a nutshell, the LG model consists of a set of diffusively coupled units, each of which represents a mesoscopic region of the cortex, and is described by a set of two dynamical equations for the local density of: (i) neuronal activity and (ii) available synaptic resources, respectively (see Appendix C for a more detailed presentation of the model). A salient feature is that as the control parameter (which can be taken to be the maximum possible level of synaptic resources) is increased, each isolated individual unit can experience an infinite-period bifurcation from a low-activity fixed point to a limit cycle -with zero-frequency and fixedamplitude at the bifurcation point-where both activity and synaptic resources oscillate in a "spike-like" manner, i.e., with a phase-dependent angular velocity, in contrast with the sinusoidal oscillation in the Kuramoto model of Eq.1, suggesting that a different minimal model is needed to describe individual units. Moreover, isolated units can produce spikes even when they are slightly below threshold owing to the effect of noise; in other words, they behave as type-I excitable units [59,60] (see Appendix D for a brief account of excitability types). Importantly, for sufficiently large values of the control parameter, the coupled oscillatory units become synchronized and, at the edge of the synchronization phase transition, scalefree avalanches emerge [37]. Finally, let us recall that a similar phenomenology is obtained in a variant of the LG model, including inhibitory neural populations [37], a well-recognized key player in the generation of neural rhythms [61].
Guided by these observations, we consider a set of coupled oscillatory units, each of them represented by the canonical form of type-I excitable units (see Appendix D) or "active rotors" as defined bẏ where ω and a are parameters [62,63]. For a > ω the deterministic dynamics of each isolated unit exhibits a stable fixed point at ϕ * = − arcsin(ω/a), as well as a saddle point with the opposite sign. The addition of a stochastic term ση(t) can induce fluctuations beyond the saddle, thus generating large excursions of the phase before relaxing back to its equilibrium. On the other hand, for a < ω the system oscillates with phase-dependent angular velocity and, as the "saddle-node into invariant FIG. 1. Illustrative diagrams sketching in a qualitative way the phase diagram -underlining the main bifurcations and phasesas analytically obtained. It reveals the existence of synchronous and asynchronous states, and a collectively-excitable regime, separated by diverse types of bifurcations for the collective order parameter. The central triangular-shaped region (whose size has been increased for clarity) describes a bistability regime; its vertices correspond to codimension-2 bifurcations (see main text). Finally, there is also a homoclinic line (Hom) linking the SNL and the BT points. (b) Sketch of the different phases/regimes represented in terms of their corresponding complex Kuramoto order parameter, Z; a point represents each regime in the unit circle (since |Z| ≤ 1) with filled circles describing fixed points, open circles stand for unstable fixed points with associated limit cycles, and mixed-color circles describe saddles. Bifurcations between different regimes are indicated as arrows.
circle" (SNIC) bifurcation point a c = ω is approached, the frequency of the oscillations vanishes, implying that the period becomes infinite, while the amplitude remains constant, as in the LG model (let us recall that type-I excitability can rely either on a "saddle-node into invariant circle" (SNIC) bifurcation, as in the equation above, or in a homoclinic bifurcation, as in the LG model [37]: both types share the common relevant feature of generating spike-like infinite-period oscillations at the bifurcation point [59]). Thus, finally, the full model readṡ where the sum runs over the (M j ) nearest neighbors of unit j ∈ 1, 2, ...N in a given network. We consider versions of the model embedded on fully-connected (FC) networks (M j = N, ∀j), which are useful for analytical approaches and on two-dimensional (2D) lattices (as, at large scales, the cortex can be treated as a 2D sheet [10,37]). Let us remark that Eq.(E1) is sometimes called Shinomoto-Kuramoto model or Winfree's ring model [64,65] and has been analyzed in diverse contexts [64,[66][67][68]. Its collective state can be quantified by the Kuramoto-Daido parameters: where k is any positive integer number; for k = 1 this coincides with the usual Kuramoto parameter. Analytical progress is achieved in a mean-field approximation (which, as usual, becomes exact for infinitely large FC networks) which consists in: (i) writing down a continuity (Fokker-Planck) equation for the probability distribution of phase values, P (ϕ, t) [53,64,[66][67][68], (ii) expanding P (ϕ, t) in power series, and (iii) writing an infinite hierarchy of coupled equations for its coefficients (which coincide with the Z k 's; see Appendix E and [64,67]): where the bar stands for complex conjugate. The associated phase diagram has been scrutinized in the literature by using different low-dimensional closures for this infinite hierarchy of coupled equations. For instance, using the Ott-Antonsen (OA) ansatz [69] or other more-refined closures [67,68], one can obtain the phase diagram, summarized in Figure 2 (a careful discussion can be found in Appendices 6 and 7). Inspection of the phase diagram reveals that there are two main types of collective dynamical regimes: oscillations (synchronous states) and stable fixed points (corresponding to either high-activity asynchronous states in the upper part of the diagram or low-activity states in the lower/right part). These are separated by different types of bifurcation lines. In particular, for low-noise amplitudes, σ, as the control parameter a is increased there is a collective SNIC bifurcation from the oscillatory regime to a phase characterized by a stable fixed point with very-low spiking activity, but susceptible to collectively react to external inputs, called the collective-excitability phase. In analogy with the classification of excitability types we refer to this as type-I synchronization transition. Conversely, for small values of a, by increasing σ one encounters a collective Hopf bifurcation, a type-II synchronization transition, to a highactivity asynchronous state (see Appendix D for a discussion of excitability types). Thus, there are two main lines of bifurcations to a collectively oscillatory phase in the two-dimensional phase diagram: one of type-I and one of type-II. Remarkably, the above two bifurcations lines cannot possibly intersect each other owing to topological reasons [59], so there is not such a thing as a "tricritical" point. Instead, in the region where they come close to each other, there exists a triangular-shaped region of bistability (green area in Figure 1a and b) [66,70] delimited by three bifurcation lines and three co-dimension-2 bifurcations (which are signaled by red circles in Figure  1b).  [64]. The location of the bistability region was established by numerically solving Eq.(E7) for the first k = 50 harmonics (Z51 = 0). The inset shows a zoom of the bistability region, where codimension-2 points are marked; a star indicates a point with scale-free avalanches. The three segments (red, green and blue) indicate three possible types of transition to synchronization.
For the sake of completeness and illustration, Figure 1b presents a sketch illustrating the (complex) order parameter in the different phases and the different bifurcation lines. In particular, there is a Bognadov-Takens (BT) point, where the Hopf-bifurcation line finishes, colliding tangentially with a line of saddle-node bifurcations; a saddle-node-loop (SNL) where the line of SNIC bifurcations ends, becoming a standard saddle-node line; and a cusp, where two saddle-node bifurcation lines collide. Observe that the bistability region is divided into two halves by the Hopf-bifurcation line, so the regime of collective excitability coexists with either oscillations below the Hopf line or the high-activity asynchronous state above the Hopf line. In other words, in this region, the Hopf bifurcation occurs in one of the branches of two co-existing solutions, i.e., in concomitance with bistability.
In order to obtain a highly-accurate phase diagram needed for forthcoming analyses, we complemented the above analytical approximations (i.e., closures) with extensive computational analyses of the complete stochastic system (size N = 10 4 ) as well as a direct numericalintegration of Eq.(E7) truncated at sufficiently large values of k (k = 50). Results of this combined approach are summarized in Figure 1a which shows that the overall shape of the phase diagram is qualitatively identical to the one predicted by low-dimensional closures, but with a slightly smaller bistability region. Remarkably, a very similar phase diagram, with the same type of phases and phase transitions is found in 2D lattices (see Appendix H).
The reader can gain insight into the dynamics in each phase and around the different transition points by inspecting Figure 3, as well as in Supplementary Material 1 [71] and contains a set of videos for diverse parameter choices.
The moral of these findings is that, to analyze general synchronization transitions, it is necessary to consider not only the standard type-I and type-II cases, but also more complex scenarios, including the case in which the transition to synchrony occurs in concomitance with bistability, i.e., when incipient oscillations coexist with low-activity asynchronous states; we call these "hybrid-type (HT) synchronization transitions". This is illustrated in Figure3 (obtained from simulations on FC networks), that shows representative raster plots around each of these possible transition types (colored segments in Figure 1a identify these three examples in the phase diagram): type-I (SNIC), type-II (Hopf), and HT synchronization transition, respectively. In particular, it shows results slightly within the synchronous phase (left), very close to criticality (central), and in the asynchronous phase (right). Surprisingly, even naked-eye inspection already reveals that raster plots nearby the HT transition exhibit a much larger dynamical richness (see below for a more quantitative analysis).

IV. SCALE-FREE AVALANCHES NEAR THE HYBRID-TYPE SYNCHRONIZATION TRANSITION
Following the protocol for avalanche detection (Appendix A) we determined the statistics of avalanche sizes and durations in these three scenarios. As shown in Figure 4 for FC networks, power-law distributed avalanches do not emerge at the type-II (Hopf) transition (red curves) nor at the type-I (SNIC) (green curves). For instance, for the Hopf-bifurcation case results resemble, not surprisingly, those for the annealed Kuramoto model which also exhibits a type-II collective bifurcation, while for the type-I case there are large collective events recruiting most of the units in the system into anomalously-large avalanches, separated by extremely long periods of qui- escence (as illustrated in the raster plot in Figure 3). In any case, we observe no signature whatsoever of scale-free avalanches at any point along the line of SNIC bifurcations (green dots in Figure 4; see also Appendix G).
On the other hand, when entering the synchronous phase through the HT transition, clean scale-invariant avalanche distributions are observed at criticality. These avalanches span across many decades and obey finite-size scaling (see Figure 4c, Figure 4d as well as Appendix G). More specifically, one can fit exponent values τ ≈ 2.1(1), α ≈ 2.5(1) and γ −1 ≈ 0.75 (5), respectively [72,73]. Observe also that there is always a peak, corresponding to anomalously system-size spanning avalanches, but these "bumps" scale with system size in a scale invariant way (much as in self-organized bistability [74]). Analogous results are found for 2D lattices -even if with slightly different exponent values: τ ≈ 1.7(1), α ≈ 1. 9(1) and γ −1 ≈ 0.75 (5), which are compatible with the values reported in experimental works [34]). In both cases, the scaling relation holds and the corresponding values of γ coincide with its seemingly robust value reported in [34]. Thus far we do not have a proper analytical understanding of these numerical values and their possible universality.
As a complementary measure of complexity at the different types of transition points, we also computed the probability distribution of the inter-spike intervals (ISI), and its associated coefficient of variation (CV). Only around the HT transition, within the bistability region, we found broad ISI distributions with significant coefficients of variation CV > 1 (see Appendix G). Thus, even if type-I and type-II synchronization transitions are well-known to exhibit signatures of criticality such as finite-size scaling (see e.g. [75] and [76], respectively) -unlike to HT transitions-they are not able to generate the higher levels of complexity required for scale-free avalanches and large dynamical repertoires as observed in the cortex.

V. DISCUSSION: ON HYBRID-TYPE SYNCHRONIZATION, CRITICALITY, AND UNIVERSALITY
One could wonder whether the rich phenomenology emerging around the HT transitions stems from any of the three special codimension-2 bifurcation points that appear in the phase diagram ( Figure 2). Remarkably, the saddle-node-loop (SNL) bifurcation has been previously argued to be necessary for the generation of high variability and dynamical richness in neural networks [77,78], and it has also been established that the crucial features of dynamically-rich neural networks stem from a phase diagram organized around a Bogdanov-Takens bifurcation [79,80]. As illustrated in the actual phase diagram (see Figure 2 and its 2D counterpart in Appendix H) the bistability region is rather small in the parameter space, so that all the discussed codimension-2 points are very close to each other. Thus, providing a clear-cut computational answer to the question above is a hard problem. Numerically, we can conclude that scale-free avalanches appear when entering the synchronous phase within the bistability region in the close vicinity of such codimension-2 bifurcations.
An aspect that needs further analysis is the relationship between criticality and bistability: these two features are usually opposed to each other, as they correspond to either continuous or discontinuous phase transitions, respectively (nevertheless, we should remind that scale-invariant avalanches can also appear at discontinuous phase transitions [74,81]). However, the onset of synchronization in the HT transition occurs within a region of bistability. The intertwining between criticality and bistability could well be at the basis of the very rich dynamical repertoires in this case. It is very plausible, that the above-described LG theory, and possibly other models [41], exhibit scale-free avalanches at the edge of synchronization, since some kind of bistability is also present around the synchronization transition.
Remarkably, the reported complex triangular-shaped structure -with its bifurcation lines and co-dimension-2 points-is rather universal and emerges in other models exhibiting a both type-I and type-II transitions (see e.g. [59,[82][83][84]); in particular, it appears in the paradigmatic and broadly used (Wilson-Cowan) model of excitatoryinhibitory networks [85]. Thus, the phenomenology discussed here is rather universal and not model-specific. It is noteworthy that other scenarios have been described to connect lines of type-I and type-II transitions -e.g., subcritical Hopf bifurcations followed by a fold of limit cycles-which also involve a regime of bistability (see [40]). Note that the phenomenology underlying these alternative scenarios is, at its core, similar to the one displayed at the HT synchronization transition, includ-ing collective excitability and some codimension-2 bifurcation points. The exciting possibility that scale-free avalanches can also emerge in such cases will be investigated elsewhere.
An aspect of our model that might need to be modified to better reproduce the results in the Landau-Ginzburg theory [37] is that here there are no true absorbing states, i.e. states from which the system cannot exit (neither as consequence of the deterministic dynamics nor as the result of stochasticity). In particular, the noise amplitude in [37] vanishes in the absence of activity, while in the present model such an ingredient cannot be easily implemented: large fluctuations can always excite individual nodes due to the absence of any bona-fide absorbing state. On the other hand, it is well established that absorbing states are needed to generate branchingprocess exponents (see e.g. [28]). Thus, further work is still needed to elucidate what happens in modelling approaches as ours if absorbing states are included. This problem will be tackled elsewhere.
Even if the minimal model studied here is exceedingly simple to be a realistic model of the cortex, it can provide us with insight on the basic dynamical mechanisms needed to generate its complex dynamical features. Furthermore, it is well-established that diverse control, self-organization (or "homeostatic") mechanisms are able to regulate a network to lie around some target regime or point [81,[86][87][88]. Thus, a cortical neural network with dynamics akin to that of the present simple model could be self-organized to the region near the HT synchronization transition [74,81,89], and by doing so, it could rapidly shift its behavior from synchronous to asynchronous, to collective excitability, or up-and-down transitions in a dynamical way, allowing for an extremely rich and flexible dynamical repertoire derived from operating at such an "edge of the edge". We hope this work opens the door for novel research lines, including renormalization group analyses [90]), paving the way to the long-term goal of constructing a statistical-mechanics of the cortex. The protocol to measure neuronal avalanches is based on the one first proposed by Beggs and Plenz [20], which has been widely employed in analyzing both experimental and theoretical data (e.g. [23,34,37,80,91]). This protocol allows one to study the structure of the spatiotemporal clusters of neuronal activity, e.g., spikes in individual neurons or peaks of the negative local field potentials. Here, we discuss how the protocol is adapted for the case considered in this paper of phase oscillators, illustrated in Fig. 5.

No. Oscillator
FIG. 5. Construction of raster plots. The individual noisy oscillators can display spike-like events; these can be detected by translating their phases into "activities" via the transformation y(t) = 1 + sin ϕ(t); then, by thresholding, the activity variable y(t) one can define "spikes" as done, e.g., in local-field-potential measurements. The top panel shows the activity of 3 randomly selected oscillators and the bottom one their corresponding spiking events as a function of time. The symbol size of each event accounts for the total time-integrated activity over the threshold. The left and right cases correspond to the synchronous phase and the bistability region, respectively.
1. For each unit j, one needs to define its "activity".
Once the activity is defined, a "spike" of a given oscillator j is a transient event occurring whenever its phase ϕ j (t) crosses a given threshold value which corresponds to a significant activity. In particular, it is possible to define the activity as y j = 1+sin ϕ j , that takes large (resp. small) values (close to 2, resp 0) when the phase crosses π/2 (resp. −π/2).

2.
A threshold y th is defined such that whenever y j crossed, a "spiking event" starts to be tracked; it finishes when y j goes again below threshold. The total integral of y j − y th along such a large-activity time window is the size s j of the local event at the initiation time t k j . k = 1, 2..., etc., label the sequence of spikes.
3. The complete set of spikes for all units j (at times t k j and sizes s k j ) defines a raster plot employed for the subsequent analyses.
4. The inter-event interval (IEI) between all couples of consecutive spikes (regardless of which unit generated them) is computed, and its probability distribution P (IEI) and its average value IEI are computed. A time scale ∆t = IEI is used to discretize the raster plots in time bins.
5. An avalanche is defined as a series of spikes in between two empty bins (with no spike) such that all consecutive time bins include some activity. The sum of all s k j in between such two empty bins is the avalanche size S, while avalanche duration T is defined as the time elapsed between the two limiting empty bins.
6. The probability distribution function (histogram) of avalanche sizes and durations is then computed.
The empirical detection of avalanches in noisy data and/or in continuous time series is often exposed to some potential pitfalls that are important to underline.
• First of all, it relies on an arbitrary choice of a threshold for detecting "spikes" of activity. In our case, selecting a threshold y th (e.g. y th = 1.6) ensures that background noise around the fixed point is filtered and only significant excursions in the phase value are considered.
• It also relies on the choice of the time bin size as the average inter-event interval; this choice has been made in accordance with the usual one in neuroscience analyses [20,23,34,92].
• Let us also remark that, instead of using the activity y j with its corresponding threshold, we could have alternatively defined a "spike" each time oscillator hits a predefined phase value, such as π/2. This method does not allow to integrate the "size" for each event, but it makes no difference with the previous one. However, in our particular case, we seek to mimic the dynamics by [37] and actual LFP recordings, where each event has a size. Therefore, we stick to the threshold y th = 1.6 for simulations.
• Let us finally stress that computational model analyses, like the ones reported here, do not suffer from severe subsampling effects that may strongly impair empirical measurements of avalanches [93][94][95].
Appendix B: Avalanches in the annealed Kuramoto model Without loss of generality, we fix J = 1 and ω = 1 -that is not relevant as a change of variables to a corotating reference frame can be used to set ω = 0-leaving σ as the only free parameter. Then, σ c = 1 indicates the critical point for infinitely large systems. Due to the finite size N of the considered networks, the precise location of the critical point needs to be computationally estimated; in particular, as usually done in finite-size scaling analyses, σ c is estimated as the value of σ such that the variance of the order parameter R is maximal [37,53]. As shown in Figure 6a and Figure 6b, showing results of our computer simulations for a system of N = 500 units, the critical point is located at σ c ≈ 0.98 (and shifts progressively towards 1 as N is increased). Observe that at the estimated critical point, owing to finite-size effects, the level of synchronization is R 0.2. This is illustrated in Figure 6c, which shows three characteristic raster plots within the synchronous phase σ = 0.9, the critical point σ c (N = 500) ≈ 0.98 and the asynchronous phase σ = 1.1, respectively.
As an important technical remark, we should emphasize that in numerical simulations of the annealed Kuramoto model, for any finite size N , the integration step δt needs to be small enough as to have sufficient time resolution, i.e. δt < ∆t = IEI , so that avalanches can be measured. Observe that this might depend on N ; in particular, in the asynchronous phase, as the number of neurons N grows, the IEI decreases, and thus one needs progressively smaller integration steps to measure avalanches. In the limit N → ∞, the asynchronous raster plot has an average interevent interval IEI → 0 and, consequently, one could say that avalanches are not well defined if a fixed time bin was considered. On the other hand by considering sufficiently small δt's for each case, our computational analyses -summarized in Figure  6 panels (d) and (e)-show that the avalanche statistics do not exhibit heavy tails.
Finally, a careful mathematical analysis of this model and its avalanches in the thermodynamics limit will be addressed elsewhere. The Landau-Ginzburg model [37] was pioneer in proposing that scale-free avalanches occur at the edge of a synchronization phase transition. It relies on a model for the dynamics of individual mesoscopic regions in the cortex. Each such region (or "unit") is characterized by two dynamical variables: its level of (mesoscopic) neuronal activity ρ(t) and the amount of available synaptic resources R(t). The dynamics (at a deterministic level, i.e. excluding fluctuations) is described by where a, b > 0 are constants, h is a very small external driving (that can be set to 0 in the absence of external stimuli), η(t) is a zero-mean Gaussian white noise, and ξ is the maximum amount of available synaptic resources, which serves as a control parameter which regulates the system state. In the second equation τ R and τ D represent the time scales of recovery and depletion of synaptic resources, respectively. Observe that the equation for the activity ρ -assuming the amount of synaptic resources R is fixed-is the minimal form of a first-order phase transition with hysteresis (or saddle-node bifurcation). It displays a quiescent (or "down") state ρ = 0 when R ≤ a, and an active or "up" state for R > a. On the other hand, the second equation accounts for the dynamics of the level of synaptic resources and includes a slow charge/recovery term (dominating when activity is low), and a fast discharge/consumption, which dominates the dynamics in the presence of activity, ρ = 0.
A simple analysis of Eqs.(C1) shows that the system behavior depends on the value of the maximum allowed synaptic resources, ξ (see Figure 8). If ξ < a, the only fixed point is a quiescent state of low activity. For larger values of ξ, the two nullclines of Eq.(C1) intersect at an unstable fixed point, giving rise to a limit cycle, i.e. to relaxation oscillations in which both ρ(t) and R(t) oscillate. Finally, if ξ is large enough, an up state (fixed point with non-vanishing activity) emerges. For more details we refer to [37] and [89]. Figure 8 illustrates the bifurcation diagram of Eq.(C1) as the control parameter ξ is varied. In agreement with what just described, values ξ < a leads to a "down" steady state with vanishing activity and non-depleted synaptic resources, i.e. R = ξ. At ξ c = a there is an infinite-period homoclinic bifurcation into a limit cycle. To determine if this bifurcation is homoclinic or rather a saddle-node into an invariant circle (SNIC) one, we have explicitly measured the average period between oscillations and plotted against the control parameter. The result displays a logarithmic decay of the period (see Figure 8b), a typical feature of homoclinic bifurcations [62]. Finally, as the control parameter is further increased, one encounters another homoclinic bifurcation at which the Points marked in red indicate that the number of spikes used for computing the average period were relatively low (because of costly statistics). Note that as ξ → ξc, larger simulations are required (increasing integration error). The data are much better fitted by a logarithm (characteristic of homoclinic bifurcations) than by a square-root fit (characteristic of SNIC bifurcations) [62]. All numerical solutions are found using Wolfram Mathematica.
limit cycle disappears, giving rise to an "up" fixed point.

Appendix D: A note on different types of excitability and bifurcations
A system is defined as excitable when it presents a single, stable equilibrium, but a sufficiently strong input can drive the system in a large excursion in the phase space before returning later to the stable fixed point [59,60,96]. Many physical and biological systems exhibit excitability [60]. Excitability is a concept of particular importance in the context of neuroscience, where neurons are at rest, and a super-threshold signal is able to evoke a significant response (e.g., a spike), returning at the end to the resting state. Different types of neurons may respond differently to the same input, leading to the so-called "excitability classes" [59,96]. The most usual classes are excitability classes I and II [59]. Type-I excitability is characterized by continuous growth of the spiking rate when the input current is continuously increased, while Type-II excitability involves a sudden jump in the spiking rate under the same circumstances. In bifurcation theory, type-I excitability corresponds to situations where the corresponding limit cycle appears with a vanishing frequency, i.e. infinite-period bifurcations, while in type-II excitability, limit cycles emerge with a finite, non-vanishing frequency [59].
Two representative examples of bifurcations corresponding to classes I and II are the SNIC and Hopf bifurcation, respectively, as sketched in Figure 9. The homoclinic bifurcation is common in neuronal models [78] and belongs to class I neuronal excitability, as the SNIC bifurcation. Both are differentiated only by the critical exponents of the firing rate. A comprehensive summary of the relationship between excitability classes and bifurcations can be found, for example, in [97].
It should be noted that the mean-field phase diagram described in the main text includes these two main types of excitability near the bistability region, where the system behaves as collectively excitable. The bistability region is surrounded by other types of bifurcations such as saddle nodes and codimension-2 bifurcations, where lines of standard (codimension-1) bifurcations intersect. Codimension-2 points, which in our case include Bogdanov-Takens, saddle-node-loop, and cusp bifurcations, can display richer dynamics [59,[77][78][79]).
We have classified synchronization transitions using a nomenclature that resembles that of excitability classes. In type-I synchronization, oscillations emerge at the transition point with zero frequency (infinite period) and finite amplitude, while in type-II synchronization, oscillations are born with a fixed non-vanishing frequency. The phenomenology becomes richer in "hybrid-type" synchronization transition, where the co-dimension 2 bifurcations and bistability are present.

Appendix E: Mathematical analysis of globally coupled oscillators
To assess the collective behavior of the system of coupled oscillators, here we reproduce with detail the derivation of the equations for the order parameter and discuss different closure methods. Most of these calculations can be found in the literature, but we reproduce them here in a self-consistent way for the sake of clarity. Our model is described by the set of stochastic equations, It is convenient to remove one parameter, fixing e.g. ω = 1. Here, we also set J = 1, leaving a and σ as the only free parameters. For completeness, we verify a posteriori that results are robust to changes in J (Appendix F).

Order-parameter equations
Equation (E1) is very similar to the previously discussed annealed Kuramoto model, except for the additional term a sin(ϕ) which induces an inhomogeneity in angular velocity across the unit circle of each oscillator. As discussed in Appendix G (in particular, in Section G 1), such an inhomogeneity makes the Kuramoto parameter inadequate to characterize the phase diagram of the present model as, for a = 0, there is a particular phase value around which each oscillator tends to spend most of the time. For uncoupled oscillators, this occurs for ϕ = ± arcsin(ω/a) where the angular velocity is minimal; this heterogeneity in angular velocity leads to a non-vanishing value of the Kuramoto order parameter, even when oscillators are uncoupled or, more in general, when they are asynchronous.
In order to circumvent this problem analytically, it is possible to consider the hierarchy of higher-order moments of the variable e iϕ , i.e., the so-called Kuramoto-Daido parameters, Z k : where k = 1, 2, ...∞ and of which the Kuramoto order parameter Z 1 is a particular case. For convenience, we will use either the notation in terms of amplitudes and phases (R k and ψ k ) to represent the complex-valued Kuramoto-Daido parameters Z k . Using standard trigonometric relations, Eq.(E1) can be rewritten as a function of Z 1 (t) = R 1 (t)e iψ1(t) , leading to the following set of Langevin equations, ϕ j (t) = ω + a sin ϕ(t) + JR 1 (t) sin (ψ 1 (t) − ϕ j (t)) + ση (t) (E3) where the mean-field nature of the coupling is evident.
In order to solve these equations, we employ a standard procedure to deal with coupled oscillators [50,64,67,98,99]. The first step is to consider a large number of oscillators, N → +∞, so that the system can be described in the continuum limit using the probability density to find an oscillator around any given phase value ϕ, i.e. P (ϕ)dϕ. The following Fokker-Planck equation gives the evolution of such a probability density, where the identity R 1 sin (ψ 1 − ϕ) = (Z 1 e −iϕ + Z 1 e iϕ )/(2i) has been used to simplify the forthcoming algebra. As the density P (ϕ, t) is periodic in the angle variable, it can be expanded in Fourier series: and p k =p −k , where the bar stands for complex conjugate. It turns out that the Kuramoto-Daido parameters coincide with these coefficients: Plugging the series expansion (E5) into the Fokker-Planck equation one obtains an infinite set of differential equations, one for each of the parameters Z k , i.e. for the Kuramoto-Daido parameters Z k (observe that, as Z −k =Z k , it suffices to analyze the order parameters with k ≥ 0). To obtain differential equations for each Z k , note that after performing the derivatives and doing some algebra, all the terms can be written as (2π) −1 k f (Z k , Z k+1 , Z k−1 , . . .)e ikϕ for some function f . Then, since the exponentials e ikϕ are the Fourierbasis elements, we can identify all parameters mode by mode, leading to an equation for the evolution of each Kuramoto-Daido order parameter [67,68]. The resulting set of equations, constitutes an exact description of the system. As reported below, we solved Eq.(E7) including up to k = 50 harmonics (i.e. imposing Z 51 = 0) and monitored the order parameters (see Section G 1). We found an excellent agreement with direct simulations, including a small bistable phase [67,98], as shown in Figure 2C of the main text. However, even if the bistable phase exists in the thermodynamic limit, its exact location is affected by finite-size effects in direct simulations.
In order to proceed analytically, given that all equations are coupled, one needs to find a dimensional sound reduction or "closure" to truncate the infinite hierarchy [99,100]. Different closures have been proposed in the literature of coupled oscillators. In what follows, we discuss three of them.

Approximate solutions or closures
In deterministic, noise-free systems, an exact solution to equation (E7) is provided by the Ott-Antonsen ansatz [69], which consists in writing Z k (t) = [Z(t)] k , so that the first moment already contains all the relevant information. On the other hand, the situation is more complicated for stochastic systems for which the Ott-Antonsen ansatz does not provide an exact solution [68,99]. In particular, using the Ott-Antonsen ansatz and writing Z(t) = R(t)e iψ(t) leads to the system of equationṡ Remarkably, these equations are the same as those obtained by Childs and Strogatz in the case of deterministic oscillators with heterogeneous (quenched) frequencies distributed as a Lorentzian [70]. However, only the deterministic system is exactly solved by the Ott-Antonsen ansatz, while in the stochastic case it gives an approximation of order O(σ 2 ) [68].
In order to increase the precision of the closure, one can assume that the global phase is distributed as a (wrapped) Gaussian with some mean ψ(t) and variance ∆(t) [67]: The Kuramoto-Daido parameters can be explicitly computed via direct integration, Plugging this ansatz into Eq.(E7) leads tȯ Observe that Eq.(E10) allows one to write Z k as a function of only the first mode [68], Z, giving a functional form very similar to the Ott-Antonsen ansatz: Let us remark that the Ott-Antonsen ansatz is equivalent to the assumption of a Lorentzian distribution for the angles. That is, changing Eq.(E9) to a Lorentzian distribution and following the same procedure, one recovers the Ott-Antonsen ansatz. The wrapped Gaussian approximation is slightly superior to the Ott-Antonsen one, but, as we will see shortly, it is not good enough as to generate a precise phase diagram (see Figure 11 below for a comparison). A way to go beyond the Ott-Antonsen and the Gaussian ansatzes is to consider more harmonics in the expansion. Tyulkina et al [68,99] proposed to use the circular cumulants of the order parameters to generate a better closure. The advantage of the cumulant expansion is that all cumulants, except χ 1 = Z, vanish when the Ott-Antonsen ansatz is selected, so choosing additional non-zero cumulants gives rise to systematic corrections to the Ott-Antonsen solution, order by order in such an expansion. In particular, the first three cumulants are given by χ 1 = Z, χ 2 = Z 2 −Z 2 , χ 3 = (Z 3 −3ZZ 2 +2Z 3 )/2. Selecting Z and χ 2 to be different from 0 but fixing χ 3 = 0 gives Z 3 = Z 3 + 3Zχ 2 , effectively closing the infinite system.
[101] Substituting this last in the equation for the first three harmonics, the resulting system reads: Note that imposing χ 2 = 0 leads to the Ott-Antonsen ansatz, as expected. This closure provides us with a noticeable quantitive improvement on where the Hopf and SNIC bifurcations are located, which coincide very well with the results of numerical integrations. However, simple parameter inspection did not render any bistability using the reduction by Tyulkina et al. Given the size of the bistable region in phase space, finding the bistability probably needs a more systematic study -such as the one we did solving the system for k = 50 harmonics, where bistability is clear as we showed in Fig. 2. Although losing the possibility of bistability in analytical approaches relying on some closures is a well-known problem in stochastic processes [100].

Appendix F: Bifurcation analysis of the Ott-Antonsen equations
To gain analytical insight into the structure and topological organization of the phase diagram, here we analyze the bifurcation diagram of the approximation provided by the simpler Ott-Antonsen closure, eqs. (E8), i.e.:Ṙ Observe first that, for a fixed value of R, the equation of the collective phase ψ is the normal form of a saddle-node into an invariant circle (SNIC) bifurcation. On the other hand, the equation forṘ is almost the same as in the annealed Kuramoto model, but adding a perturbation proportional to the "excitability parameter" a. The above set of equations is difficult to study analytically, but its bifurcations can be obtained following the same procedure of Childs and Strogatz who studied this system with a fixed value of σ = √ 2 [70]. The main idea is as follows: in the annealed Kuramoto limit (a = 0) the system undergoes a Hopf bifurcation; on the other hand, individual (uncoupled) oscillators, i.e., for J = 0, exhibit a SNIC bifurcation. Hence, by continuity in solutions, we expect two branches of these two types of bifurcations to be present in Eq.(E8).
Calling Q the Jacobian at a fixed point, at Hopf bifurcations tr Q = 0 while in a saddle-node det Q = 0 [62]. Thus, imposing one of these conditions, together with the fixed point equationsṘ =ψ = 0, leads to a set of equations for the parameters of the system as a function of the fixed point values, R * and ψ * . Since such values are bounded, one can use these equations as parametric equations of the bifurcation curve, without computing explicitly the values of the fixed points and their stability.
Let us start with the Hopf bifurcation. Between parameters and fixed points, there are 6 unknowns: R * , ψ * , ω, a, J and σ. Remember that in the simulations we fixed ω = J = 1 to leave a and σ as the only free parameters. In what follows, we obtain equations for the bifurcations, written in a parametric form a = a(R * , ψ * , ω, J) and σ = σ(R * , ψ * , ω, J). After some algebra, it turns out that not all the dependences are necessary. Solving foṙ R =ψ = tr Q = 0, there are 3 remaining parameters. We can choose any parameters to solve for, but it turns out (as shown in [70]) that solving for R, cos ψ and sin ψ is highly convenient. Note that these are only two parameters, since the sine and cosine are not independent function as sin 2 ψ + cos 2 ψ = 1. Of course, one could have tried to directly solve for ψ and a, but it is easier to obtain expressions for the trigonometric functions and then extract a: Note that, in this case, we obtained a parametric curve for the Hopf bifurcation, a H = a H (ω, J, σ), without requiring specific knowledge about the location of the fixed points. In particular, as J is kept fixed, it is possible to derive a curve a H = a H (ω, σ).
In what respects saddle-node bifurcations, the calculation is a bit more involved since solving for the same 3 variables gives high-degree polynomials for R that cannot be explicitly solved. For this reason, we choose to solve for ω, cos ψ and sin ψ, since in this way the parameters do not depend on ψ. After solving and simplifying the resulting equations, one obtains: Since 0 ≤ R ≤ 1, and one needs to explore all the possible fixed points, there are two free parameters to choose. Given that J is kept fixed in our simulations, we dismiss it and obtain ω and a as functions of σ. The saddle-node lines enclose a "collective excitable" phase, which is type-I excitable. This can be easily seen by realizing, as discussed above, that phase dynamics are again of the formψ = ω + c sin ψ for fixed Kuramoto order parameter R, with the angular speed taking the role of the external input usually employed to classify excitability classes. The manifolds of Hopf and SNIC bifurcations could be drawn together in a three-dimensional space, using ω, a and σ as coordinates. An alternative easy way to visualize them is to make projections into the (a, σ) twodimensional space for different values of J (see Figure 10). Since we set ω = 1, and in each such projection J is fixed, the Hopf bifurcation is then obtained as a H = a H (σ), while the saddle node is obtained by solving Eqs.(F4) as a parametric curve, depending on R and σ. Figure  10 shows that there is a bistability region delimited by a Hopf line and two saddle-node lines for different values of J. This region decreases in size as the coupling decreases until it disappears for sufficiently low values of J.
Appendix G: Computational analyses and results

Phase diagram of the full model
As we have seen, deriving all the phases and transitions between them analytically is a difficult task, and thus, one needs to resort to computational analyses. First of all, an adequate order parameter needs to be defined. As discussed above, the usual Kuramoto order parameter, R = |Z|, is not a good choice for inhomogeneous oscillators, because | Z t | = 0 even when they are uncoupled or asynchronous, due to the different amount of time they spend at diverse phase values.
We employed the so-called Shinomoto and Kuramoto (SK) parameter that solves this problem, being able to discriminate between synchronous and asynchronous regimes [64]. In particular, we consider thecomputationally more efficient-variant of such a parameter employed by Lima and Copelli [40]: to distinguish the synchronized phase (S = 0) from asynchronous or excitable states (S = 0). The main result of our computational analyses is the phase diagram reported in Figure 2 of the main text. The different phases are identified computationally as follows: • A non-vanishing S value characterizes the synchronous region.
• The asynchronous and "collective excitability" regions are both characterized by S = 0. To detect this second regime one needs to perturb the system and analyze its collective response (or lack of it).
• The regime of bistability is challenging to study numerically since in finite networks, fluctuations can drive the system to jump between the two coexisting states. To determine this region, the analytical solution was computed by solving for the first 50 terms in the series expansion (E7). It was then computationally verified that, for large enough network sizes, two alternative stable states exist within such a region.

Accuracy of different closures
For completeness, we checked the accuracy of the different considered closures by comparing them with the results of computational analyses. In order to do so, the different ansatzes or closures: (i) the Ott-Antonsen (eqs. (E8)), (ii) the Gaussian closure (eqs. (E12)), and (iii) the equations proposed by Tyulkina et al. (eqs. E14) were solved near the Hopf and SNIC bifurcation branches and compared the results -shown in Figure 11-with those of direct simulations of Eq.(E1) as reported above. We would like to remark that, up to our knowledge, it is the first time that the accuracy of the different closures to locate different kind of bifurcations has been formally checked. The conclusion is that the ansatz by Tyulkina et al. [99] fits more accurately both the Hopf and SNIC branches than the other ones. The Gaussian ansatz captures very well the phenomenology near the SNIC bifurcation but not near the Hopf one. On the other hand, the Ott-Antonsen solution does not predict the transition accurately at any of the bifurcations. However, the ansatz by Tyulkina et al. [99] is not able to predict the existence of the bistability phase, while the other two do so. Solving the complete system of equations for the Kuramoto-Daido parameters, at least 5 harmonics are needed in order to find the bistability region.

Avalanches at different types of bifurcations
Here, we further investigate whether scale-free avalanches emerged when we moved away from the bista-bility region. Results at the different type of bifurcations are reported in Figure 2 of the main text, which displays raster plots computed at either a (Type II) Hopf (upper panels) or at a (Type I) SNIC bifurcation point (lower panels), respectively, as well as within the synchronous and asynchronous phases surrounding them. Here, we show results for values slightly below, at, and above the bifurcation. As discussed at extent in the main text, scale-free avalanches (with power-law distributed sizes and durations) emerge only in the vicinity of the hybrid type synchronization transition; when the system is moved away from it, scaling behavior as well as scalefree behavior breaks down (see Figure 12), while near the bistability region scale-free distributed avalanches remain (see Figure 13).

Dynamical variability
As a complementary measure of complexity at the different transition points, we have also computed the probability distribution of the inter-spike intervals (ISI), along with its associated coefficient of variation (CV) where µ ISI and σ ISI are the mean and standard devia- 13. Avalanches close to the hybrid type synchronization transition Avalanche-size and avalanche-duration distributions for a network of size N = 5000 evaluated at the hybrid type synchronization transition (a = 1.07, σc = 0.496) and two other nearby points, slightly to the left of it. Power law behavior is observed only in the bistable region close to the hybrid-type synchronization.
tion of the inter-spike interval for each oscillator, respectively, and · indicates an average over units. It is important to distinguish between the interspike interval of each single unit, and the interval between any two consecutive spikes in a network. When looking at avalanches, the second option is preferred, since it gives a small timescale able to resolve the internal structure of avalanches. On the other hand, to measure the CV one focuses on interspike intervals of individual oscillators. For a Poisson process, one expects an exponential distribution of ISI's along with a CV = 1; as a rule of thumb, CV > 1 is the fingerprint of irregular spiking activity. Figure 14 shows the probability distribution of the ISIs, and the CVs for the different phases and transitions, as determined in computational analyses. Only two cases exhibit CV > 1: the hybrid type synchronization transition, as well as a small neighborhood around it in the bistable regime; they are also the only two cases characterized by a broad distribution of ISI values. Thus, a high level of variability -similar to that observed in the cortex-is only found in the region around the hybrid type synchronization transition, but not in the neighborhood of either standard Hopf or SNIC bifurcations. Computational analyses in 2D systems reveal a very similar phase diagram to the mean-field one but with a richer phenomenology (as graphically illustrated in Figure 15). There are three main phases, in a nutshell: a synchronous regime, an asynchronous one, and a collectively excitable phase, much as in the mean-field case. In the asynchronous regime, clusters of activity appear, propagate, and vanish dynamically, with an averaged constant network activation level (smaller/larger close to the SNIC/Hopf transition, respectively) but with no overall synchronization. On the other hand, within the synchronous phase, activity wane and washes and periods of overall quiescence are followed by bursts of overall activity that spreads quickly from different focuses. Such a collective propagation requires some level of synchronization (let us recall that a perfectly phase-synchronous state does not exist in 2D [98] as rotational symmetry cannot possibly be broken in low-dimensional systems; there are always some "topological defects" in the system preventing it to exhibit perfect synchronization, as it happens in well-known models of equilibrium statistical mechanics [102] and also in the Kuramoto model [53]). Finally, the excitable state has all units in the "down" regime, with almost no activity, but is susceptible to respond to external perturbations or inputs.
Let us discuss the observed phenomenology at the different bifurcation lines separating these phases: (i) Nearby the type-I synchronization transition, at the SNIC bifurcation (slightly within the synchronous phase), oscillating activity appears in the form of traveling waves. Visual inspection reveals, e.g. the presence of typical spiral patterns typical of two-dimensional excitable systems (see Supplementary Material 1 [71] for videos, and Figure 15 upper row).
(ii) Nearby the type-II synchronization transition, at the Hopf bifurcation (slightly within the synchronous phase) the overall level of activity oscillates in time (i.e., the system "breathes") and is spatially distributed in fragmented clusters (see Figure15, central row).
(iii) Near the hybrid type synchronization transition, the dynamical behavior is much more complex: somehow in between the overall oscillatory behavior along the Hopf Spatio-temporal dynamics on twodimensional systems The figure shows three rows, each one with six different frames corresponding to six different times on a running simulation, for the following cases: Upper row: near the different bifurcations SNIC (a = 1, σ = 0.08); central row: near a Hopf bifurcation (a = 0.5, σ = 0.65); and lower row: hybrid type synchronization transition (a = 0.98, σ = 0.185). Blue color indicates lack of activity, while red color stands for maximum activity levels (identified as 1 + sin φj). Near the SNIC transition, noise fluctuations generate wavefronts that propagate in the system. At the Hopf transition, there is some background activity whose level grows and shrinks periodically. At the hybrid type synchronization transition, the spatio-temporal patterns are more complex, being a mixture of the two aforementioned types; in particular, the system sometimes falls in the excitable state. Simulations performed for N = 64 2 with periodic boundary conditions. Videos showing the evolution for both cases and the other phases can be found in the Supplementary Material 1 [71].
line and the emergence of wavefronts at the SNIC line (see Figure 15 bottom row).
Thus, even if we are not attempting to quantify spatiotemporal complexity here, it is clear that more complex dynamics emerge around the HT synchronization transition. As shown in Figure 16 the computationallyobtained phase diagram in the case of two-dimensional lattices has a structure qualitatively identical to the mean-field case, including the same phases.
[1] W. R. Softky and C. Koch, The highly irregular firing of cortical cells is inconsistent with temporal integration of random epsps, J. Neurosci. 13, 334 (1993). . Reddish colors stand for asynchronous states, while blueish colors correspond to partially synchronized ones. As in the FC-network case, the system can de-synchronize either through a supercritical Hopf bifurcation or through a SNIC bifurcation. In an intermediate region it is possible to find a hybrid type synchronization transition (yellow star). This is characterized as the point where power-law distributed avalanches obeying finite-size scaling appear. Symbols describe the coordinates for which the videos of Supplementary Material 1 [71] were recorded: "Hopf" bifurcation (black triangles); "SNIC" bifurcation (green squares) and the HT transition point (white circles and yellow star).