Abstract
Research on network percolation and synchronization has deepened our understanding of abrupt changes in the macroscopic properties of complex engineered and natural systems. While explosive percolation emerges from localized structural perturbations that delay the formation of a connected component, explosive synchronization is usually studied by fine-tuning of global parameters. Here, we introduce the concept of synchronization bombs as large networks of heterogeneous oscillators that abruptly transit from incoherence to phase-locking (or vice-versa) by adding (or removing) one or a few links. We build these bombs by optimizing global synchrony with decentralized information in a competitive percolation process driven by a local rule, and show their occurrence in systems of Kuramoto –periodic– and Rössler –chaotic– oscillators and in a model of cardiac pacemaker cells, providing an analytical characterization in the Kuramoto case. Our results propose a self-organized approach to design and control abrupt transitions in adaptive biological systems and electronic circuits, and place explosive synchronization and percolation under the same mechanistic framework.
Similar content being viewed by others
Introduction
The emergence of abrupt, explosive transitions in the macroscopic behavior of complex networked systems is a fascinating phenomenon, ubiquitous in fields ranging from neuroscience to biology and engineering. There is increasing empirical evidence that explosive synchronization in brain activity is associated with the onset of anesthesic-induced unconsciousness1,2,3, epileptic seizures4,5 and fibromyalgia6 and it explains biological switches displaying abrupt responses to external perturbations7. Also, in infrastructural and power-grids networks, it is crucial to detect and control small vulnerabilities that can lead to abrupt structural damages and global desynchronization blackouts8,9.
From a theoretical perspective, explosive percolation—an abrupt growth of the giant component of the network induced by the addition or removal of single links—was found to occur when competitive rules are applied on the choice of the links in a way that the formation of a giant cluster is delayed10. The discovery of this abrupt structural transition, which was shown to be continuous in the thermodynamic limit but with anomalous scaling properties, triggered further analyses to understand the mechanisms that can lead to the explosive behavior in the network growth. In parallel, abrupt transitions were explored in the collective dynamics of the system when considering a physical process among the units, as the spreading of a disease11,12,13, opinion diffusion14 or traffic flow15,16, to name a few17,18.
A particularly suitable framework to model the birth of explosive transitions is the synchronization process of coupled oscillators. The phenomena of collective synchronization are widely spread in natural, social and technological systems19,20. Its ubiquity has attracted the interest of the physics community, that have tackled its study through minimal models that capture the transition between a disordered phase and coherent dynamics. For populations of heterogeneous phase-oscillators coupled all-to-all19,20, abrupt transitions in synchrony as the coupling parameter is increased were found to occur for a uniform distribution of frequencies21 and the hysteresis cycle involving incoherence and partial synchrony was exactly characterized for a bimodal distribution with a shallow dip22. However, the bistable nature of explosive synchronization—an abrupt jump from incoherence to global synchrony induced by a change in the coupling parameter among the units, with an associated hysteresis cycle—was firstly discovered for scale-free networks (i.e. networks with very heterogeneous degree distributions) in the presence of positive correlations between the internal frequencies and the nodal degrees23. Further analyses showed that this is only one of the possible mechanisms that inhibit the emergence of a large synchronization cluster, and it was found that, by imposing frequency anti-correlations among connected units in the form of frequency gaps24 or adaptive anti-Hebbian rules for the weights25, explosive transitions occur as the coupling constant is tuned. Furthermore, explosive transitions can also appear in multilayer and dynamically coupled systems26,27 and they can be enhanced by the presence of noise28 and higher-order -beyond pair-wise- interactions29.
While our fundamental understanding of explosive synchronization has significantly increased in the last years, some key questions remain open. Due to the analytical challenges of studying synchronization dynamics on arbitrary complex networks, a rigorous framework analogous to explosive percolation is still missing17. Research on explosive synchronization and percolation has focused on different aspects of the system, namely in the abrupt changes on the macroscopic dynamical or structural properties, respectively, when subject to small variations of the control parameter (the coupling strength or the density of links). Interestingly, in ref. 30, the authors found that a particular choice of frequency-dependent coupling (which again induces anti-correlations) produces an explosive synchronization process where the formation of synchronized clusters is delayed analogously to its percolation counterpart. While these results unveiled a connection between both phenomena, the specific model used in ref. 30 is tailored ad hoc to show explosive phenomena and the system displays the transitions under changes in a global control parameter (the coupling strength), unlike the explosive percolation which is induced locally, by adding or removing single links. Also, it is well established that dynamical and structural correlations are necessary to induce the explosive behavior17,18,23 and recently, it has been suggested that these correlations can be related to maximal global synchrony31,32. However, it is not understood how empirical oscillator networks can self-organize toward these explosive regimes, when they usually evolve under uncertain or stochastic conditions and exploit only limited, decentralized information available in the immediate surroundings of the units5,7,17,33.
To tackle some of the aforementioned challenges from a theoretical perspective, here we present a self-organized dynamical network that act as a synchronization bomb, i.e. showing an abrupt synchronization transition in the course of a self-organized wiring process. In this way, our model attempts to bridge the conceptual gap between explosive synchronization and percolation by imposing local structural perturbations instead of global ones and proposes a self-organized route to explosive synchronization by invoking a simple principle of synchrony maximization in a decentralized and stochastic process.
The remainder of this paper is organized as follows. We first present our model of the synchronization bomb for an ensemble of Kuramoto oscillators. We introduce the optimal local rule for connecting or disconnecting units, derived from the truncated expansion of the linearized dynamics31 under the assumptions of maximizing global synchrony with local information, and explore the basic mechanisms and phenomenology of the synchrony-driven percolation process. Second, we analyze the explosive fingerprints that emerge on the underlying structure, in the form of degree-frequency correlations, disassortative dynamical and structural patterns, and a delayed percolation threshold. Third, we provide an analytical characterization of the dynamics by means of the Collective Coordinates34,35 (CC) and Ott-Antonsen36 (OA) model reduction techniques, unveiling the dependence of the main parameters and observing an excellent agreement with numerical simulations. Next, we extend the model to numerically produce synchronization bombs of coupled chaotic Rössler systems and cardiac pacemaker cells. We conclude with a discussion of our Results and a Methods’ section, including the mathematical machinery used to derive the local rule and the analytical predictions for the percolation and synchronization critical thresholds.
In the Supplementary Information, we study the robustness of the presented phenomenology under variations in system parameters. In particular, in Supplementary Note 1, we study the effects of size and noisy sampling, showing that the explosive transitions induced by single links are sustained for increasing size and that the presence of some randomness is beneficial because it improves the decentralized optimization of synchrony driven by a local rule. In Supplementary Note 2, we study the effect of different choices of frequency distributions. In Supplementary Note 3, we validate the extension of the model to directed networks, and in Supplementary Note 4, we explore in more detail the relation between percolation and synchronization onsets in the model.
Results
Model
We consider a large system of heterogeneous coupled oscillators on top of a network of interactions that evolves under a competitive link percolation process10,17. For the dynamics, we begin with the classical Kuramoto model, a paradigmatic example of the emergence of collective synchronization19,20,37. An ensemble of N heterogeneous Kuramoto oscillators interacting on top of a network follows the equations of motion:
where θi is the phase and ωi is the intrinsic frequency of the i-oscillator, aij are the entries of the adjacency matrix A, that capture the interactions among the units and λ is a constant coupling strength. As usual, the macroscopic behavior of the system is captured by the modulus of the Kuramoto order parameter:
which measures the degree of phase synchronization and is bounded between zero and one. In the following, we will make use of temporal averages of the order parameter, i.e. r = 〈r(t)〉. Although our results can be extended to more general settings, in the following we restrict our study to the case of unweighted (aij = 0, 1) and undirected networks (aij = aji), and consider, for analytical convenience, a uniform frequency distribution g(ω) ∈ [−γ, γ], shifting the resulting frequency vector to have zero mean.
The growth of the synchronization bomb is made by keeping the coupling strength λ constant and varying the density in the number of connections between the units, p, that acts as the control parameter and ranges from 0 (disconnected network) to 1 (fully-connected network). In the forward process we initialize our system from scratch, with a completely disconnected network (p = 0) of oscillators with assigned random phases drawn from (−π, π) (other choices are explored in Supplementary Note 2). We then run the percolation processes in which at each step one new link is added. This way the control parameter p changes sufficiently slow such that the system in Eq. (1) reaches the stationary state at each network step in the process. The addition of a new link at each p-step is made as follows: We uniformly sample M pairs of disconnected oscillators and select the connection (i, j) that maximizes the gain of synchrony given by:
In practice, when connecting isolated nodes at the very initial steps of the process, we add an infinitesimal ϵ to the degrees of the nodes with ki = 0 to evaluate Eq. (3) in terms only of the natural frequencies (ϵ = 1 does not alter the results). In the backward process, we just remove the links in the reversed order of the forward process. The proposed model is stochastic in nature but becomes completely deterministic in the limit M → ∞, and it reduces to the random percolation case in M = 1.
Equation (3) captures the actual change in the order parameter r in the strong phase-locking regime (i.e. after the transition) but it can be used to estimate the impact of each link in the whole synchronization process (see Methods section “Derivation of the local rule” for in depth derivation and discussion of this expression). Note that Eq. (3) only exploits local information of the considered nodes, and it is maximum when the ratios frequency-degree of the nodes are large and also when their difference is large as well, pinpointing a clear signature of frequency-degree correlations and frequency-frequency anti-correlations. These correlations were imposed ad hoc in previous models that induce explosive synchronization17, and they could indeed emerge from applying a broader class of local percolation rules in the form p(ωi, ki, ωj, kj), but we focus on Eq. (3) since it is the rule that is derived from a decentralized optimization of the phase-locking state, without other assumptions or guesses required.
In Fig. 1a, we illustrate the former basic mechanics for assigning a link out of M = 5 possible candidates and in Fig. 1b we show the histogram of the values of Eq. (3) for all the existing links at a given step of the process, highlighting the value of the five sampled, potential new links. We take the forward process (construction of the network by adding links) as an example. As shown in Fig. 1a the functional form of the basic rule, Eq. (3), induces some relevant features on the interplay between structural and dynamical patterns during the network growth. We observe that nodes with large (small) absolute frequencies accumulate more (less) neighbors, whereas links tend to be more present between nodes with alternate frequencies, producing bipartite-like structures, as we will explore in the following lines. In Fig. 1c we show the forward and backward explosive synchronization transitions by plotting the curves r(p) when different values of M are used. We observe that as M increases so it does the abruptness of the transition as well as the hysteresis region. To illustrate better the explosive nature of these transitions we show in Fig. 1d, e the transition from incoherence (r ≈ 0.05) to full phase-locking (r ≈ 0.9) when a unique link is added to the system (in Supplementary Note 1 we check that this effect holds for larger sizes). This phenomenon motivates our choice for referring to these growing networks as synchronization bombs.
Structural explosive fingerprints
Before characterizing the synchronization transition of explosive bombs in more depth, we now focus on the structural changes that the network undergoes during the percolation process governed by Eq. (3). In the following we analyze the emergence of several structural and dynamical patterns that are usually associated with explosive transitions17,31.
Degree-frequency correlations and frequency-frequency anticorrelations
During the network growth process, the system tends to a stationary degree distribution (which scales with system density) as sampling increases in the process (for larger M). To understand this effect, we recall that the rule of Eq. (3) tends to connect pairs of oscillators with large frequency differences and low degrees. When a link is chosen, the degrees of the adjacent nodes increase, reducing the value of Δrij for other potential links of these nodes. This constant competition between fixed frequencies and evolving degrees acts as a self-organized feedback that tends to homogenize the distribution of Δrij among the potential links, and frequencies and degrees become balanced in the precise way that makes Δrij more similar among these—still absent—links. Since the rule predicts the scaling Δr ~ ω2/k3, we find that, in the deterministic (large M) regime (see Methods section for details) the relation is given by:
where the scaling is controlled by the mean degree, expressed in terms of the density of links p and size N, since 〈k〉 ≈ pN. As expected, Eq. (4) becomes more accurate as sampling increases in the system, as observed in Fig. 2a. In Fig. 2b, we see that frequency anti-correlations among pairs of connected nodes are also present in the system and become stronger for large sampling M. Let us note that these type of correlations are explicitly imposed in the majority of studied mechanisms that induce explosive synchronization17,18 whereas here emerge from a decentralized optimization of the synchronized state. In the following, we explain how these dynamical anti-correlations translate into structural ones.
Spectral signatures: toward optimal and bipartite networks
We study the evolution of the extreme eigenvalues μmax and μmin of the Laplacian matrix L = D − A (D is the diagonal matrix of degrees) and of the normalized Adjacency matrix \(\hat{A}={D}^{-\frac{1}{2}}A{D}^{-\frac{1}{2}}\) during the percolation process, for different values of sampling M. In Fig. 2c, d, we observe that the network evolves in a path that maximizes both the largest positive eigenvalue of L, μmax(L), ranging from zero to N, and the largest negative eigenvalue of \(\hat{A}\), \({\mu }_{{{{{{\rm{min}}}}}}}(\hat{A})\), ranging from minus one to zero, when randomness is reduced in the process (larger sampling M). Also, the frequency of the oscillators tends to correlate with the entries of the associated extreme eigenvectors (see insets of both panels). These spectral signatures pinpoint that our model evolve toward optimal and bipartite configurations. First, it is well understood that optimal synchronization is achieved by the alignment of the frequencies with the largest eigenvector of the Laplacian matrix and by increasing the magnitude of the associated eigenvalue μmax(L)38. Second, note that the normalized Adjacency matrix, \(\hat{A}\), is a stochastic row sum and its spectra is bounded in \(\mu (\hat{A})\in [-1,1]\), with the largest eigenvalue \({\mu }_{{{{{{\rm{max}}}}}}}(\hat{A})=1\) if the network is connected. The remaining of the spectra follows Wigner’s semicircle law for random networks, becoming narrower as the link density increases, and it deviates from the random case in the presence of modules (shifting toward positive eigenvalues) or bipartite-like structures (shifting toward negative eigenvalues)31. Thus, from Fig. 2d we observe that bipartite patterns arise as determinism is increased (larger M) and the trajectory of the extreme eigenvalues tuple follows a clear asymmetric path toward the all-to-all (p = 1) limit. This effect shows that the rule derived in Eq. (3) induces negative structural correlations (bipartitivity) as a consequence of the negative dynamical correlations that emerge in terms of natural frequencies, and vice-versa.
Delayed percolation threshold
From the former results, it is clear that, as the percolation process evolves, the network self-organizes its architecture according to well-known explosive patterns. Two main issues are how this synchrony-driven percolation is related to the natural one, i.e. that observed when links are chosen at random, and, as we will cover below, how the emergence of a giant component (the proportion of the nodes connected in the largest cluster of the network) is related to the synchronization onset. To address these issues we study the emergence of the giant component as a function of the control parameter p when the rule of Eq. (3) is applied for different values of the sampling parameter M. In Fig. 2e, we observe that the proposed rule delays the percolation threshold with respect to the random case, and it produces more abrupt transitions. Looking more closely at the effect of the system parameters on the percolation threshold, in Fig. 2f, g we observe that, when increasing both the size of the system (large N) and the determinism in the rule (large M), percolation transitions become sharper and occur at higher p. Nevertheless, the nature of the transition appears to be continuous (i.e. second order) even for large system sizes. We can obtain a rough approximation for the average value of the percolation threshold by using the well-known Molloy and Reed criterion39 in the deterministic limit and neglecting the negative structural correlations that the rule induces. Using this criterion and leveraging the emergent degree-frequency correlation we obtain, for a uniform g(ω), (see Methods) that the threshold is estimated as:
which can be written in terms of the percolation threshold in a random network9 as \({p}_{c}\approx 1.68\cdot {p}_{c}^{rand}\). In Fig. 2e we observe that Eq. (5) works quite well for sufficiently large M. More sophisticated analytical tools, as the recently developed feature-enriched percolation framework40, could improve the predictions under local rules, such as Eq. (3), that exploit information both from the degrees and the frequencies of the units.
Analytical characterization of the Kuramoto bomb
Now we explore, by analytical and numerical means, the dynamical regimes of our system depending on the coupling, λ, and sampling, M, values. We remark that, despite the apparent simplicity of Eq. (1), the Kuramoto Model on complex networks does not have an analytical solution and approximations are required to predict the dynamical behavior using the information contained in A and ω17,20.
To the best of our knowledge, the current method that better captures the finite-size effects and the precise interplay between the structure and the oscillator dynamics in Eq. (1) is the model reduction technique based on Collective Coordinates, introduced first by Gottwald to globally coupled systems34 and extended to complex networks in35. We use this approach to estimate the value of the oscillator phases and the corresponding evolution of the order parameter r(p) in the backward branch and also to calculate numerically the backward synchronization threshold, \({p}_{c}^{b}\). See Methods section “Collective coordinates ansatz” for the precise details of this method. The agreement between CC theory and numerical simulations becomes evident in the backward synchronization diagrams shown in Fig. 3a for λ = 0.02 and 0.04 (M = 10).
For the forward process we cannot use the CC approach and we rely on the celebrated OA ansatz36, which has been successfully used to characterize systems in the presence of frequency and degree correlations41,42,43. Specifically, we benefit from a recent elegant development used to describe the mean-field dynamics of Janus oscillators43 and consider the limit of large N and M. Complete calculations to predict the loss of stability of the incoherent state, and therefore the forward synchronization threshold, are given in the Methods section “Ott-Antonsen ansatz”. For the particular case of a uniform distribution g(ω) ∈ [−1, 1] we obtain the closed form:
The predicted value \({p}_{c}^{f}\) is plotted in Fig. 3a showing again a remarkable agreement. This analytical estimation allows addressing the aforementioned issue about the relation between synchronization and percolation onsets by making use of Eqs. (5) and (6). Combining both expressions we can write a simple relation for the percolation pc and forward synchronization \({p}_{c}^{f}\) thresholds as:
which illustrates the natural connection between the structural and dynamical aspects of our model.
We extend our numerical and analytical characterization of the synchronization diagram in the (p, M)-plane, Fig. 3b), and (p, λ)-plane, Fig. 3c). In Fig. 3b, we observe that, fixing λ = 0.05, the collision of the theoretical backward curve and the approximated forward threshold successfully predicts the codimension-two point, where a saddle-node bifurcation collides/appears with a pitchfork bifurcation and bistability emerges29. This critical point for which explosive behavior shows up takes place around M ≈ 5. In Fig. 3c, we focus on the coupling strength λ, a parameter that does not play a role in the percolation process, but it is crucial to synchronization dynamics. The precise location of the synchronization thresholds can be controlled from occurring simultaneously with the percolation one for large values of λ (where the oscillators immediately synchronize when connected and the synchronization curves follow the percolation ones shown in Fig. 2e), to occur much later for smaller values of λ (where bistability and abruptness of transitions is maximized) and to finally disappear for sufficiently small λ. Interestingly, in the explosive range of λ, the emergence of a bipartite structure for sufficiently large M—as discussed in the previous section—hinders the formation and growth of microscopic small clusters in the incoherent regime (see for instance Fig. 1d), against what occurs in other models of explosive synchronization17,23,30. Furthermore, we note that the system transits more abruptly for a relatively large λ (low p), but has a larger region of hysteresis for low λ (large p).
In Supplementary Note 1 we explore the dynamics of the model for different system sizes, confirming that the abrupt jump in r occurring at single link changes remains large for increasing size, and we analyze in more detail the role of the noisy sampling in the model, finding that an optimal amount of sampling can enhance the explosive performance of the synchronization bomb because it improves the self-organized optimization process driven by a local rule. In Supplementary Note 2 we show that both phenomenology and theory are robust to changes in the distribution of intrinsic frequencies, g(ω). In particular, we show results for Gaussian (and bimodal) cases, exploring scenarios with less (and more) polarization than the uniform distribution, finding the expected result that polarization in g(ω) increases the bistable regime and the abruptness of the transitions. In Supplementary Note 3 we show that the bomb-like transitions also occur for directed networks. Lastly, in Supplementary Note 4 we also illustrate in more detail the relation between percolation and synchronization curves under the parametrization of Fig. 3c.
Chaotic synchronization bombs
One of the most relevant applications of synchronization theory is its implementation when coupling chaotic systems44, a counter–intuitive non-linear phenomenon as it achieves a perfect dynamical coherence between systems that, when isolated, display exponential divergence of nearby trajectories. Thus, to show the generality of our results, we round off by extending them beyond the Kuramoto framework and considering the Rössler system, a paradigmatic model for the emergence of chaotic dynamics45.
Here we use an ensemble of diffusively coupled heterogeneous chaotic oscillators46,47,48, a modified, piece-wise linear Rössler system45, which evolves in a 3-dimensional space following:
where the non-linear function that induces the chaotic behavior is defined as g(x) = 0 if x ≤ 3 and g(x) = μ(x − 3) if x > 3. The remaining parameters are set following47,48, with τ = 0.05, β = 0.5, δ = 1, ν = 0.02 − 100/R. R = 100 ensures that the system is in a phase-coherent regime46,47,48, where a phase can be defined after projecting onto the xy-plane, i.e. θi = arctan(yi/xi), such that the synchronization order parameter r can be measured by the standard Eq. (2). See Fig. 4c, d for a 3D representation of the trajectories of the chaotic, phase-coherent, oscillators at two different p-steps of the forward process. As in Eq. (1), λ is the fixed coupling strength and the entries aij of the adjacency matrix A capture the presence of undirected and symmetric interactions between the oscillators and evolve under the rule of Eq. (3). The instantaneous velocity of each unit is determined by fi, which we assign proportional to the frequency, fi = 10 + 0.2ωi, drawn again from a uniform distribution g(ω) in [−1, 1].
Figure 4a illustrates three examples of synchronization transitions r(p) for a system of N = 200 and different choices of λ and M. Similarly to the Kuramoto case, it is observed how, in the construction process, for sampling values of M > 1 the order parameter experiments abrupt jumps from dynamical incoherence of r ~ 0.1 to a more coherent state with r ≃ 0.7, that continues to continuously grow to stronger synchronization (r ≳ 0.9) as the link fraction, p, increases. For the backward transition the inverse process takes place but with the jump to incoherence happening for lower values of p, resulting in a small hysteresis cycle. In Fig. 4c we show the synchronization diagram in the (p, M)-plane, where it becomes clear that the bistable region shows up even for very small values of M.
The success of the chaotic synchronization bomb is grounded on previous research that exploits optimal48 and explosive47 synchronization properties of the Kuramoto Model on the diffusively coupled Rössler system. However, as numerical results in Fig. 4a, b manifest, the phenomenology is slightly noisier than in the Kuramoto case and the tuning of more parameters along with the chaotic behavior of the units may difficult its design and control. From a practical standpoint, these results show that synchronization bombs could be potentially implemented in the lab, at least by means of electronic circuits47.
Application to cardiac pacemaker cells
Lastly, we demonstrate the existence of self-organized explosive synchronization bombs in the biologically-plausible application of cardiac pacemaker cells—the collection of cells responsible for generating a strong, coherent pulse that propagates through the entire heart and initiates each contraction49. For simplicity, we consider a system of network-coupled pacemaker cells using, for each pacemaker, a two-variable system describing the dimensionless trans-membrane voltage v and gating variable h which summarizes ionic concentrations49. For a system of N such pacemakers the equations of motion are given by:
where the local dynamics of vi and hi are described by:
The timescales τi represent local heterogeneity between the different pacemakers, scaling the period of each isolated cell, ultimately resulting in an effective natural frequency for each pacemaker proportional to \({\tau }_{i}^{-1}\). Taking a system of N = 200 pacemakers with \({\tau }_{i}^{-1}\) uniformly distributed in [0.4, 1.6] and using coupling strengths Kv = 0.009 and Kh = 0.0044 (to indicate a stronger coupling via the voltage diffusion compared to ionic diffusion) we implement the coupled percolation and synchronization dynamics as presented previously in this work.
To measure the synchronization of the full system we consider the error in the voltage dynamics, quantified by the overall standard deviation. Taking temporal means of the error as the percolation dynamics are run forward and backwards, we plot the voltage error in Fig. 5a. In Fig. 5b, c we present the actual voltage dynamics right before and after the backward explosive transition, plotting each individual voltage time series vi(t) in a light blue stroke, and indicating the overall mean using a thick, dark blue stroke. Note here the physiological implications of the pacemakers ability or inability to produce a strong, coherent pulse for strongly and weakly synchronized behavior, respectively.
By inspecting Fig. 5a, we observe a region of bistability and explosive transitions in the voltage error, being more abrupt in the backward direction, where the transition occurs at a single link removal. As occurred with Rössler dynamics, results are noisier than in the Kuramoto case and the precise location of the thresholds—and the associated width of hysteresis—vary across realizations, presumably due to the small network size considered here. A detailed exploration of the synchronization bombs in these more nuanced oscillator models remains open for further research, but current results already confirm that the main phenomenology is sustained beyond systems of Kuramoto-like phase oscillators.
Discussion
Abrupt and explosive phenomena in the structure and dynamics of complex networks have been one of the most studied phenomena in non-equilibrium statistical physics and non-linear dynamics in recent years. Not only do they allow us to further our theoretical understanding of phase transitions, but also to develop models that are able to explain and reproduce the changes in the topology and behavior observed in natural and engineered systems, such as biological switches, brain activity and blackouts in power-grids. Motivated by the wide range of applications, network percolation and collective synchronization have become paradigmatic frameworks to understand the explosive changes in the structural and dynamical macroscopic properties of large complex systems. A crucial feature of explosive percolation is that it is induced by applying small localized structural perturbations to the system (addition or removal of a few links) by means of competitive rules that delay the formation of a connected component. This aspect was not explored in the synchronization counterpart, where explosive transitions were usually studied by fine-tuning of global coupling parameters in fixed or evolving structures. Furthermore, while the specific theoretical requirements for the explosive behavior become better understood, there is less knowledge on the actual routes that real systems may follow to self-organize toward these particular configurations.
In this work, we have attempted to bridge these gaps by deriving a local percolation rule for systems of heterogeneous phase-oscillators under the minimal assumption of maximizing global synchronization with decentralized information. We have shown that under this percolation rule the system behaves as a synchronization bomb. This way the network undergoes an explosive synchronization transition at some point of the wiring process, abruptly switching from incoherence to global phase-locking, and display a hysteresis cycle. We have also shown that as the network grows, it self-organizes in a way that several well-known explosive properties on the network structure show up. Crucially, this growth delays the percolation threshold as compared to the usual random case. We have provided an analytical characterization of the system using state-of-art model reduction techniques, obtaining a fair agreement with numerics and being able to reproduce the bistable region in the synchronization phase diagrams. As we show in the SI, all these results are robust under the variation of model assumptions and parameters, and also hold for directed networks. Interestingly, we find that a noisy, low sampling is beneficial in our model because it improves the decentralized optimization of synchrony driven by the proposed local rule. Finally, we have shown that synchronization bombs can be also obtained for systems of coupled chaotic units, paving the way to their implementation in the lab and in a model of cardiac pacemaker cells, proving potential applications in biological systems.
While the current results provide a self-organized and stochastic route to the emergence of these bombs, alternative, deterministic approaches could lead to a better optimization of the explosive behavior and control of the location of the transitions in empirical networked systems. From a theoretical perspective, our results propose a mechanistic origin of explosive synchronization in phase models, in terms of a local percolation rule that exploits both structural and dynamical information from the degrees and the frequencies of the oscillators. In this sense, our model can be seen as a synchrony-driven process in the recently established framework of feature-enriched percolation40. Furthermore, the analytical and numerical findings presented here align well with the novel hypothesis of self-organized bistability in neuronal dynamics50 and also provide a missing explanation to the emergence of explosive behavior in pair-wise networks via the activation of a single self-organized mechanism, somewhat analogous to the explosive transition described in globally coupled networks when a single parameter, the one controlling the strength of higher-order interactions, is switched on in the system51.
In a nutshell, our findings show that growing networks of heterogeneous dynamical units can develop to operate in a bistable regime, forming networked switches that display the dynamical-structural correlations that are observed when graphs are tuned to display explosive behavior. Also, engineered networks can be designed to be at the onset of total synchrony in which they show no dynamical coherence but, after a minimal wiring (just one or few links), experience synchronization explosions. This finding provides a justification for naming these systems as synchronization bombs.
Methods
Derivation of the local rule
We begin with two key assumptions: (1) the system attempts to maximize the overall degree of synchronization, by adding or removing undirected connections in a percolation process and (2) only limited information is available, making this percolation a decentralized process. This means that the units have access only to their immediate surroundings and they can exploit only local information to maximize synchronization, without having access to the overall network synchronization. In order to derive the rule under the previous assumptions, we invoke linearization arguments on the original system Eq. (1), which are shown to be valid when looking for optimal structural and dynamical properties even far away from the linearized regime38,52. Under the linearization, the resulting system reads in matrix form as:
where L = D �� A is the Laplacian of the network. The solution of Eq. (13) in the stationary state is found by setting \(\dot{{{{{{{{\boldsymbol{\theta }}}}}}}}}=0\). In the corotating frame at speed 〈w〉 = 0, the solution reads as:
where L† is the Moore-Penrose pseudoinverse of the Laplacian matrix, which can be constructed via the spectral decomposition of L for undirected networks (see ref. 38 for more details). Since we are close to the synchronization attractor, the phases in Eq. (2) can also be expressed in a Taylor expansion38. Invoking again linearization, the order parameter is given by:
In principle, one needs all the spectral information of the network to estimate the value of r. However, we can leverage recent results on the geometric expansion of Eq. (14)31, where it is shown that the linearized solution can be expressed as a sum of contributions from increasingly further neighborhoods. This way, the local approximation of synchrony31 is obtained by truncating the expansion at its second term (taking into account the effect of the nearest neighbors of the nodes), leading to:
where \({z}_{i}=\mathop{\sum }\nolimits_{j = 1}^{N}{a}_{ij}({\omega }_{j}/{k}_{j})\) is the contribution of first neighbors. For more details on the accuracy of Eq. (16), see ref. 31. From Eq. (16), we can estimate the local impact in the synchrony of adding or removing a single link between oscillators (i, j). Both discrete (considering single link perturbations) and continuous (using derivatives with respect to the degrees and the approximation \(\partial {z}_{n}/\partial {k}_{n^{\prime} }\approx \pm {\delta }_{nn^{\prime} }{\omega }_{n^{\prime} }/{k}_{n^{\prime} }\), evaluating the resulting expression at zi = 0) calculations, in the limit of large degree, lead to Eq. (3) in the results section, an expression that depends only on the local variables (ωi, ki) of a given pair of nodes. Explicitly:
where the ±sign accounts for the addition (removal) of a link. We remark that this result is derived assuming no bias in the coupling function of Eq. (1), symmetric and unweighted interactions and a frequency distribution of zero mean, meaning that the actual frequencies of the oscillators may need an appropriate shift to satisfy the condition31. Also, note that one could obtain more accurate rules for the maximization of r by using the exact result for the phases given by Eq. (14) or by including higher-order terms beyond the local approximation, although this increase of accuracy would require to use either spectral (thus global) information or to go beyond the local variables up to second-neighbors and so on. Furthermore, we note that a quadratic approximation of Eq. (17) as \(\Delta r \sim {({\omega }_{i}/{k}_{i}-{\omega }_{j}/{k}_{j})}^{2}\) also induces the explosive phenomena and may simplify the analytical treatment, but its study is left for further research.
Derivation of the percolation threshold
The percolation threshold is approximated by the Molloy and Reed criterion39, that predicts the transition for random network without correlations for the value of p = pc at which:
To compute 〈k〉 and 〈k2〉, we consider a uniform distribution, such that g(ω) = 1/(2γ) if ω ∈ [ − γ, γ] (the same analysis could be done for any other frequency distribution) and also take into account that, explicitly, we have the general correlation ki = c∣ωi∣2/3 where c is a normalization constant which depends on the network density, p, as well as the distribution of natural frequencies, g(ω). Using that 〈k〉 = p(N − 1), we find that c = p(N − 1)/〈∣ω∣2/3〉, with \(\langle | \omega {| }^{2/3}\rangle =\frac{3}{5}{\gamma }^{2/3}\). Thus we obtain the correlation:
and with it:
Thus, substituting in Eq. (18) we obtain:
which corresponds to Eq. (5) in the Results section.
Collective coordinates ansatz
We use the theory introduced in refs. 34,35 to predict the phases of the oscillators at any given p-step of the backward process and also the transition from phase-locking to incoherence. The main idea of the method is to reduce the dimensionality of the system by considering, as an ansatz, that the phases of the oscillators in the phase-locking regime are in the form:
where ψi is the exact solution of the linearized dynamics of Eq. (1)38, i.e. \({\psi }_{i}=\frac{1}{\lambda }{L}^{{{{\dagger}}} }{{{{{{{\boldsymbol{\omega }}}}}}}}\). By minimizing the error made by Eq. (3) in the full dynamics of Eq. (1) and after some manipulation35, one ends up with only one differential equation for the evolution of the q coefficient, thus drastically reducing the dimensionality from N coupled differential equations to a single one. The resulting equation reads as:
Solving the implicit Eq. (23) for \(\dot{q}=0\) allows estimating the phases of the oscillators in Eq. (1) beyond the linear regime of the system. Here, we use this theory to predict the phases of the oscillators and the corresponding curve for the order parameter in the full phase-locking regime of the system. Furthermore, to predict the appearance of the (backward) critical threshold \({p}_{c}^{b}\) within this theory, we use an explosive trick. We assume beforehand that in the explosive regime of our system, the backward process transits from full phase-locking to complete incoherence. With this idea in mind, we predict the backward threshold by looking at the last values (p, λ) for which Eq. (23) has a solution. Additionally, we check that the solution is linearly stable by numerically computing the eigenvalues of the Jacobian matrix of the full system in Eq. (1) around the equilibrium solution \(q=\hat{q}\). The Jacobian evaluated at the equilibrium point reads as35:
The system is stable if all the eigenvalues of J are negative. Thus, the backward critical threshold occurs at the last value of \({p}_{c}^{b}\) at which Eq. (23) admits a solution that is linearly stable. The explosive trick is particularly useful to simplify the calculation because, when considering transitions from full phase-locking to incoherence, we do not need to compute partial synchronized solution involving clusters of smaller size than the whole network35. In other words, we predict the loss of stability of the full phase-locking state, which in the explosive regime of our system corresponds to the desired backward synchronization threshold.
Ott-Antonsen ansatz
In the forward direction, we cannot use the Collective coordinates approach anymore because our system departs from the incoherent state where the ansatz Eq. (22) is not valid. Numerical simulations showed that, usually for M > 1, the incoherent state r ≈ 0 remains stable beyond the backward critical transition, thus creating a bistable region and a delayed forward transition. In order to analytically predict the forward critical threshold, we consider the limit of large N and also large M (toward a deterministic rule). In practice, the following results turn out to be valid even for quite small values such as N = 200 and M = 5, but we remark that the theory is derived in the infinite size and deterministic limits of the model. Our approach is based on the celebrated OA ansatz36 and follows a very similar development to that shown in ref. 43.
We begin by defining the local order parameter:
such that Eq. (1) can be written as:
Following43, we consider a large ensemble of systems, described by the joint probability density ρ(θ, ω, t), with θ = (θi, …, θN) and ω = (ωi, …, ωN). The evolution of the joint probability has to satisfy the continuity equation36:
where θi is given by Eq. (26). Multiplying the density function ρ by ∏j≠idωjdθj and integrating, one obtains the evolution for the marginal oscillator density, ρi(θi, ωi, t) which reads as43:
Now, the OA ansatz can be applied by expanding ρi in a Fourier series and setting the coefficients of the expansion as \({b}_{i,n}={\alpha }_{i}^{n}\)36,43. By inserting the Fourier series with the ansatz in Eq. (28), one ends up with:
where R* and \({\alpha }_{j}^{* }\) represent the complex conjugate and i the imaginary unit. Now we invoke the large M assumption. In this deterministic limit, the underlying network is purely bipartite, split between nodes with positive frequencies and nodes with negative ones (see the Results section and Fig. 2c, d). Also, in this limit, the frequencies of the oscillators are completely determined by their degrees. Then, we can look for solutions αi = αk,±43, reducing the problem to finding solutions for the coefficients of degree classes in the two groups. The local order parameter in this setting can be written as43:
The frequencies of the degree classes in the two groups are completely determined by the percolation rule for a wide range of p, leading to:
where c is a scaling constant given in Eq. (19). After these considerations, the resulting system can be written as:
Since we want to evaluate the stability of the incoherent state αk,± = 0, we linearize the system above and evaluate it around αk,± = δαk,± ≪ 1. After neglecting smaller terms of order δα2, the dependence on the complex conjugates vanish and we up with the following system for each degree class:
By defining the variables \(\delta x={\sum }_{k^{\prime} }k^{\prime} {p}_{k^{\prime} }\delta {\alpha }_{k^{\prime} ,+}\) and \(\delta y={\sum }_{k^{\prime} }k^{\prime} {p}_{k^{\prime} }\delta {\alpha }_{k^{\prime} ,-}\), and summing over degree classes (taking into account the degree distribution), we can write:
With the approximation ∑kk5/2pkδαk,+ ≈ 〈k3/2〉δx and ∑kk5/2pkδαk,− ≈ 〈k3/2〉δy, the set of equations reduces to a 2-dimensional variational system for the evolution of δx and δy that reads as:
It is straightforward to show that the critical condition for the stability of the incoherent state is given by:
In particular, the eigenvalues of the Jacobian matrix change from being both imaginary to become both real as density increases in the system. In fact, the fully imaginary spectrum predicts the existence of a center attractor, indicating a marginal stability of the incoherent state. Therefore, one might expect stationary oscillations of the order parameter43. Here we do not observe these oscillations in the forward process. The system is initialized with isolated units (in the incoherent state) and remains there as the network evolves in an adiabatic manner. Fortunately, the forward abrupt transition to phase-locking is well predicted by the critical condition given by Eq. (41), when the eigenvalues become real (one positive and one negative) indicating the appearance of an unstable saddle point. Accordingly, when the condition is achieved in the forward, growth process, the marginal stability of the incoherent state is lost and the system transits to phase-locking.
Using that in the deterministic limit we have that ki = c∣ωi∣2/3, and for a general g(ω) the constant is given by c = 〈k〉/〈∣ω∣2/3〉, we obtain a general closed form for the forward critical threshold (pc, λc) that is given by:
For the particular case of a uniform distribution g(ω) ∈ [−γ, γ] we can easily compute the expected moments and, after plugging these results in Eq. (42), we end up with the simple formula:
which corresponds to Eq. (6) in the Results section.
Data availability
Numerical data presented in this work are available from request to the authors.
Code availability
Computational and visualization codes to reproduce the results of this work are available from request to the authors.
References
Joiner, W. J. et al. Genetic and anatomical basis of the barrier separating wakefulness and anesthetic-induced unresponsiveness. PLoS Genet. 9, e1003605 (2013).
Kim, M. et al. Functional and topological conditions for explosive synchronization develop in human brain networks with the onset of anesthetic-induced unconsciousness. Front. Comput. Neurosci. 10, 1 (2016).
Kim, M., Kim, S., Mashour, G. A. & Lee, U. Relationship of topology, multiscale phase synchronization, and state transitions in human brain networks. Front. Comput. Neurosci. 11, 55 (2017).
Wang, C.-Q., Pumir, A., Garnier, N. B. & Liu, Z.-H. Explosive synchronization enhances selectivity: example of the cochlea. Front. Phys. 12, 1–9 (2017).
Wang, Z., Tian, C., Dhamala, M. & Liu, Z. A small change in neuronal network topology can induce explosive synchronization transition and activity propagation in the entire network. Sci. Rep. 7, 1–10 (2017).
Lee, U. et al. Functional brain network mechanism of hypersensitivity in chronic pain. Sci. Rep. 8, 1–11 (2018).
Chatterjee, A., Kaznessis, Y. N. & Hu, W.-S. Tweaking biological switches through a better understanding of bistability behavior. Curr. Opin. Biotechnol. 19, 475–481 (2008).
Dobson, I., Carreras, B., Lynch, V. & Newman, D. E. Complex systems analysis of series of blackouts: cascading failure, critical points, and self-organization. Chaos 17, 026103 (2007).
Newman, M. Networks: An Introduction. (Oxford University Press, Inc., 2010).
Achlioptas, D., D’Souza, R. M. & Spencer, J. Explosive percolation in random networks. Science 323, 1453–1455 (2009).
De Domenico, M., Granell, C., Porter, M. & Arenas, A. The physics of spreading processes in multilayer networks. Nat. Phys. 12 https://doi.org/10.1038/nphys3865 (2016).
Bttcher, L., Woolley Meza, O., Araujo, N., Herrmann, H. & Helbing, D. Disease-induced resource constraints can trigger explosive. Sci. Rep. 5 https://doi.org/10.1038/srep16571 (2015).
Matamalas, J. T., Gómez, S. & Arenas, A. Abrupt phase transition of epidemic spreading in simplicial complexes. Phys. Rev. Res. 2, 012049 (2020).
Gmez-Gardees, J., Lotero-Vlez, L., Taraskin, S. & Prez-Reche, F. Explosive contagion in networks. Sci. Rep. 6, 19767 (2016).
Echenique, P., Gómez-Gardeñes, J. & Moreno, Y. Dynamics of jamming transitions in complex networks. Europhys. Lett. 71, 325–331 (2005).
Lampo, A., Borge-Holthoefer, J., Gómez, S. & Solé-Ribalta, A. Multiple abrupt phase transitions in urban transport congestion. Phys. Rev. Res. 3, 013267 (2021).
D’Souza, R. M., Gmez-Gardees, J., Nagler, J. & Arenas, A. Explosive phenomena in complex networks. Adv. Phys. 68, 123–223 (2019).
Boccaletti, S. et al. Explosive transitions in complex networks? Structure and dynamics: percolation and synchronization. Phys. Rep. 660, 1–94 (2016).
Pikovsky, A., Rosenblum, M. G., & Kurths, J. Synchronization, A Universal Concept in Nonlinear Sciences (Cambridge University Press, 2001).
Arenas, A., Díaz-Guilera, A., Kurths, J., Moreno, Y. & Zhou, C. Synchronization in complex networks. Phys. Rep. 469, 93–153 (2008).
Pazó, D. Thermodynamic limit of the first-order phase transition in the Kuramoto model. Phys. Rev. E 72, 046211 (2005).
Martens, E. et al. Exact results for the Kuramoto model with a bimodal frequency distribution. Phys. Rev. E 79, 026204 (2009).
Gómez-Gardeñes, J., Gómez, S., Arenas, A. & Moreno, Y. Explosive synchronization transitions in scale-free networks. Phys. Rev. Lett. 106, 128701 (2011).
Leyva, I. et al. Explosive transitions to synchronization in networks of phase oscillators. Sci. Rep. 3, 1281 (2013).
Avalos-Gaytán, V. et al. Emergent explosive synchronization in adaptive complex networks. Phys. Rev. E 97, 042301 (2018).
Zhang, X., Boccaletti, S., Guan, S. & Liu, Z. Explosive synchronization in adaptive and multilayer networks. Phys. Rev. Lett. 114, 038701 (2015).
Soriano-Paños, D., Guo, Q., Latora, V. & Gómez-Gardeñes, J. Explosive transitions induced by interdependent contagion-consensus dynamics in multiplex networks. Phys. Rev. E 99, 062311 (2019).
Skardal, P. S. & Arenas, A. Disorder induces explosive synchronization. Phys. Rev. E 89, 062811 (2014).
Skardal, P. S. & Arenas, A. Abrupt desynchronization and extensive multistability in globally coupled oscillator simplexes. Phys. Rev. Lett. 122, 248301 (2019).
Zhang, X., Zou, Y., Boccaletti, S. & Liu, Z. Explosive synchronization as a process of explosive percolation in dynamical phase space. Sci. Rep. 4, 5200 (2014).
Arola-Fernández, L., Skardal, P. S. & Arenas, A. Geometric unfolding of synchronization dynamics on networks. Chaos 31, 061105 (2021).
Wei, C., Shengfeng, W., Yueheng, L., Weiqing, L. & Jinghua, X. Explosive synchronization caused by optimizing synchrony of coupled phase oscillators on complex networks. Eur. Phys. J. B 94, eabe3824 (2021).
Scarpetta, S., Apicella, I., Minati, L. & de Candia, A. Hysteresis, neural avalanches, and critical behavior near a first-order transition of a spiking neural network. Phys. Rev. E 97, 062305 (2018).
Gottwald, G. Model reduction for networks of coupled oscillators. Chaos 25 https://doi.org/10.1063/1.4921295 (2015).
Hancock, E. & Gottwald, G. Model reduction for kuramoto models with complex topologies. Phys. Rev. E 98 https://doi.org/10.1103/PhysRevE.98.012307 (2018).
Ott, E. & Antonsen, T. Low dimensional behavior of large systems of globally coupled oscillators. Chaos 18, 037113 (2008).
Kuramoto, Y. Chemical Oscillations, Waves, and Turbulence (Dover Publications, 2003).
Skardal, P. S., Taylor, D. & Sun, J. Optimal synchronization of complex networks. Phys. Rev. Lett. 113, 144101 (2014).
Molloy, M. & Reed, B. A critical point for random graphs with a given degree sequence. Random Struct. Algorithms 6, 161–180 (1995).
Artime, O. & De Domenico, M. Percolation on feature-enriched interconnected systems. Nat. Commun. 12 https://doi.org/10.1038/s41467-021-22721-z (2021).
Restrepo, J. & Ott, E. Mean field theory of assortative networks of phase oscillators. Europhys. Lett. 107 https://doi.org/10.1209/0295-5075/107/60006 (2014).
Skardal, P. S., Restrepo, J. G. & Ott, E. Frequency assortativity can induce chaos in oscillator networks. Phys. Rev. E 91, 060902 (2015).
Peron, T., Eroglu, D., Rodrigues, F. & Moreno, Y. Collective dynamics of random janus oscillator networks. Phys. Rev. Res. 2 https://doi.org/10.1103/PhysRevResearch.2.013255 (2020).
Boccaletti, S., Kurths, J., Osipov, G., Valladares, D. L. & Zhou, C. S. The synchronization of chaotic systems. Phys. Rep. 366, 1–101 (2002).
Rössler, O. E. An equation for continuos chaos. Phys. Lett. A 57, 397 (1976).
Rosenblum, M. G., Pikovsky, A. & Kurths, J. Phase synchronization of chaotic oscillators. Phys. Rev. Lett. 76, 1804–1807 (1996).
Leyva, I. et al. Explosive first-order transition to synchrony in networked chaotic oscillators. Phys. Rev. Lett. 108, 168702 (2012).
Skardal, P. S., Sevilla-Escoboza, V. P., Vera-Ávila, V. P. & Buldú, J. M. Optimal phase synchronization in networks of phase-coherent chaotic oscillators. Chaos 27, 013111 (2017).
Djabella, K., Landau, M. & Sorine, M. A two-variable model of cardiac action potential with controlled pacemaker activity and ionic current interpretation, in 2007 46th IEEE Conference on Decision and Control 5186–5191 https://doi.org/10.1109/CDC.2007.4434970 (2007).
Buendía, V., di Santo, S., Villegas, P., Burioni, R. & Muñoz, M. A. Self-organized bistability and its possible relevance for brain dynamics. Phys. Rev. Res. 2, 013318 (2020).
Kuehn, C. & Bick, C. A universal route to explosive phenomena. Sci. Adv. 7, eabe3824 (2021).
Dörfler, F., Chertkov, M. & Bullo, F. Synchronization in complex oscillator networks and smart grids. Proc. Natl. Acad. Sci. USA 110, 2005–2010 (2013).
Acknowledgements
L.A.-F. and A.A. acknowledge the Spanish MINECO (Grant No. PGC2018-094754-B-C2). J.G.-G. acknowledges the Spanish MINECO (Grant No. FIS2017-87519-P), the Departamento de Industria e Innovación del Gobierno de Aragón and Fondo Social Europeo through (Grant No. E36-17R FENOL), and Fundación Ibercaja and Universidad de Zaragoza (Grant No. 224220). A.A. acknowledges financial support from Generalitat de Catalunya (grant No. 2017SGR-896), Universitat Rovira i Virgili (grant No. 2019PFR-URV-B2-41), Generalitat de Catalunya ICREA Academia, and the James S. McDonnell Foundation (grant #220020325). This work was supported by MINECO and FEDER funds through Projects No. FIS2017-87519-P, No. FIS2017-90782-REDT (IBERSINC); from grant PID2020-113582GB-I00 funded by MCIN/AEI/10.13039/501100011033; and by the Departamento de Industria e Innovación del Gobierno de Aragón y Fondo Social Europeo through Grant No. E36-17R (FENOL). S.F.-L. acknowledges financial support by Gobierno de Aragón through the Grant defined in ORDEN IIU/1408/2018. E.-C.B. acknowledges support from the “Agencia Estatal de Investigación” (Ref. PRE2019-088482), Government of Spain (FIS2020-TRANQI; Severo Ochoa CEX2019-000910-S), Fundació Cellex, Fundació Mir-Puig, and Generalitat de Catalunya (CERCA, AGAUR).
Author information
Authors and Affiliations
Corresponding authors
Ethics declarations
Competing interests
The authors declare no competing interests.
Peer review
Peer review information
Communications Physics thanks the anonymous reviewers for their contribution to the peer review of this work.
Additional information
Publisher’s note Springer Nature remains neutral with regard to jurisdictional claims in published maps and institutional affiliations.
Supplementary information
Rights and permissions
Open Access This article is licensed under a Creative Commons Attribution 4.0 International License, which permits use, sharing, adaptation, distribution and reproduction in any medium or format, as long as you give appropriate credit to the original author(s) and the source, provide a link to the Creative Commons license, and indicate if changes were made. The images or other third party material in this article are included in the article’s Creative Commons license, unless indicated otherwise in a credit line to the material. If material is not included in the article’s Creative Commons license and your intended use is not permitted by statutory regulation or exceeds the permitted use, you will need to obtain permission directly from the copyright holder. To view a copy of this license, visit http://creativecommons.org/licenses/by/4.0/.
About this article
Cite this article
Arola-Fernández, L., Faci-Lázaro, S., Skardal, P.S. et al. Emergence of explosive synchronization bombs in networks of oscillators. Commun Phys 5, 264 (2022). https://doi.org/10.1038/s42005-022-01039-2
Received:
Accepted:
Published:
DOI: https://doi.org/10.1038/s42005-022-01039-2
This article is cited by
-
The spatial dynamics and phase transitions in non-identical swarmalators
Nonlinear Dynamics (2024)
-
Diffusion capacity of single and interconnected networks
Nature Communications (2023)