††thanks<EMAIL_ADDRESS>
# Self-organized explosive synchronization in complex networks:
Emergence of synchronization bombs

Lluís Arola-Fernández  Departament d’Enginyeria Informàtica i Matemàtiques,
Universitat Rovira i Virgili, 43007 Tarragona, Spain  Sergio Faci-Lázaro
Department of Condensed Matter Physics, University of Zaragoza, 50009 Zaragoza
(Spain).  GOTHAM lab, Institute of Biocomputation and Physics of Complex
Systems, University of Zaragoza, 50018 Zaragoza, Spain  Per Sebastian Skardal
Department of Mathematics, Trinity College, Hartford, CT 06106, USA  Emanuel-
Cristian Boghiu  ICFO –- Institut de Ciències Fotòniques, The Barcelona
Institute of Science and Technology, 08860 Castelldefels (Barcelona), Spain
Jesús Gómez-Gardeñes  Department of Condensed Matter Physics, University of
Zaragoza, 50009 Zaragoza (Spain).  GOTHAM lab, Institute of Biocomputation and
Physics of Complex Systems, University of Zaragoza, 50018 Zaragoza, Spain
Alex Arenas<EMAIL_ADDRESS>Departament d’Enginyeria Informàtica i
Matemàtiques, Universitat Rovira i Virgili, 43007 Tarragona, Spain

###### Abstract

_We introduce the concept of synchronization bombs as large networks of
coupled heterogeneous oscillators that operate in a bistable regime and
abruptly transit from incoherence to phase-locking (or vice-versa) by adding
(or removing) one or a few links. Here we build a self-organized and
stochastic version of these bombs, by optimizing global synchrony with
decentralized information in a competitive link-percolation process driven by
a local rule. We find explosive fingerprints on the emerging network
structure, including frequency-degree correlations, disassortative patterns
and a delayed percolation threshold. We show that these bomb-like transitions
can be designed both in systems of Kuramoto –periodic– and Rössler –chaotic–
oscillators and in a model of cardiac pacemaker cells. We analytically
characterize the transitions in the Kuramoto case by combining a precise
collective coordinates approach and the Ott-Antonsen ansatz. Furthermore, we
study the robustness of the phenomena under changes in the main parameters and
the unexpected effect of optimal noise in our model. Our results propose a
minimal self-organized mechanism of network growth to understand and control
explosive synchronization in adaptive biological systems like the brain and
engineered ones like power-grids or electronic circuits. From a theoretical
standpoint, the emergence of synchronization explosions and bistability
induced by localized structural perturbations –without any fine-tuning of
global parameters– joins explosive synchronization and percolation under the
same mechanistic framework._

##  I Introduction

The emergence of abrupt, explosive transitions in the macroscopic behavior of
complex networked systems is a fascinating phenomenon, ubiquotuos 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 unconsciousness Joiner  _et
al._ (2013); Kim  _et al._ (2016, 2017), epileptic seizures Wang  _et al._
(2017a, b) and fibromyalgia Lee  _et al._ (2018) and it explains biological
switches displaying abrupt responses to external perturbations Chatterjee  _et
al._ (2008). Also, in infrastructural and power-grids networks, it is crucial
to detect and control small vulnerabilites that can lead to abrupt structural
damages and global desynchronization blackouts Dobson  _et al._ (2007); Newman
(2010).

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 delayed Achlioptas
_et al._ (2009). 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 disease
De Domenico  _et al._ (2016); Böttcher  _et al._ (2015); Matamalas  _et al._
(2020), opinion diffusion Gómez-Gardeñes  _et al._ (2016) or traffic flow
Echenique, P.  _et al._ (2005); Lampo  _et al._ (2021), to name a few D’Souza
_et al._ (2019); Boccaletti  _et al._ (2016).

A particularly suitable framework to model the birth of explosive transitions
is the synchronization process of coupled oscillators. The phenomena of
collective synchronization is widely spread in natural, social and
technological systems Pikovsky  _et al._ (2001); Arenas  _et al._ (2008). 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-all Arenas  _et al._ (2008); Pikovsky  _et
al._ (2001), abrupt transitions in synchrony as the coupling parameter is
increased were found to occur for a uniform distribution of frequencies Pazó
(2005) and the hysteresis cycle involving incoherence and partial synchrony
was exactly characterized for a bimodal distribution with a shallow dip
Martens  _et al._ (2009). 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 degrees
Gómez-Gardeñes  _et al._ (2011). 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 gaps Leyva  _et
al._ (2013) or adaptive anti-Hebbian rules for the weights Avalos-Gaytán  _et
al._ (2018), explosive transitions occur as the coupling constant is tuned.
Recently, it has been found that these degree and frequency correlations
associated to explosive behavior also optimize the global phase
synchronization in the system Arola-Fernández  _et al._ (2021); Wei Chen and
Xiao (2021). Furthermore, explosive transitions can also appear in multilayer
and dynamically coupled systems Zhang  _et al._ (2015); Soriano-Paños  _et
al._ (2019) and they can be enhanced by the presence of noise Skardal and
Arenas (2014) and higher-order -beyond pair-wise- interactions Skardal and
Arenas (2019).

Our fundamental understanding of explosive synchronization has significantly
increased in the last years, but due to the analytical challenges of
synchronization dynamics on arbitrary complex networks, a rigorous framework
analogous to explosive percolation is still missing D’Souza  _et al._ (2019).
Importantly, explosive synchronization and percolation focus on different
aspects of the system, namely in the abrupt changes on the macroscopic
dynamical and structural properties, respectively, when subject to small
variations of the control parameter (the coupling strength or the density of
links). Interestingly, in Zhang  _et al._ (2014), 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 unveil a deep connection between both
phenomena, the choice on the coupling dependence is heuristic and the system
produces the explosive behavior under changes in a global control parameter
(the coupling strength), unlike the explosive percolation which is induced
locally, by adding or removing single links.

Apart from the theoretical interest, there is an ongoing consensus that
explosive synchronization phenomena is behind the operation of biological
switches and neural systems Chatterjee  _et al._ (2008); Kim  _et al._ (2016,
2017); Wang  _et al._ (2017a, b); Lee  _et al._ (2018); Joiner  _et al._
(2013); Myung  _et al._ (2018). Biological units usually operate with limited,
decentralized information and are affected by noise Ishida (1997); Izhikevich
(2003); Orlandi  _et al._ (2013); Zhang  _et al._ (2016). This ubiquity of the
explosive transitions cannot be explained by means of global and deterministic
optimization routes, specific network and oscillator designs or global fine-
tuning of coupling parameters. Thus, it is still poorly understood how complex
biological systems like the brain can self-organize to display the observed
explosive behavior Chatterjee  _et al._ (2008); Wang  _et al._ (2017b);
Scarpetta  _et al._ (2018).

To tackle 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. 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 and stochastic route to explosive synchronization by invoking a
simple principle of synchrony maximization in a decentralized and noisy
environment.

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 dynamics Arola-
Fernández  _et al._ (2021) 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, dissasortative dynamical and structural
patterns, and a delayed percolation threshold. Third, we provide an analytical
characterization of the dynamics by means of the Collective Coordinates
Gottwald (2015); Hancock and Gottwald (2018) (CC) and Ott-Antonsen Ott and
Antonsen (2008) (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 (SI), we study the robustness of the presented
phenomenology under variations in system parameters, and we explore in depth
the effect of noisy sampling in our model, showing that the presence of noise
is beneficial because it improves the decentralized optimization of synchrony
driven by a local rule.

##  II Results

Figure 1: a) Illustrative network of $N=50$ oscillators where the size of the
node is proportional to its degree, the color is related to its natural
frequency (blue for $\omega=-1$, gray for $\omega=0$ and red for $\omega=1$)
and the black lines represent the links between the oscillators. Green lines
mark the $M=5$ potential links sampled in that $p$-step. The continuous line
represents the chosen link and the dashed ones are the discarded ones. b)
Histogram of the $\Delta r$ values for the existing links of the network,
where red lines correspond to the values of the five sampled links. c) Example
of the typical synchronization transition in our bomb-like model, with the
order parameter $r$ depending on the fraction of links, $p$ in a system of
$N=200$ Kuramoto oscillators for fixed coupling strength $\lambda=0.05$ and
three values of sampling $M$. Eq.(1) is numerically integrated using Heun’s
method, with $dt=0.05$ and $10^{4}$ time steps and temporal averages of $r$
are taken at each link change. In d), we represent the oscillators phases for
the $M=10$ case before (left) and after (right) the forward transition. It is
important to note that this jump in the Kuramoto order parameter from
incoherence to complete phase–locking takes place just with the addition of
one link.

Model. We consider a large system of heterogeneous coupled oscillators on top
of a network of interactions that evolves under a competitive link percolation
process Achlioptas  _et al._ (2009); D’Souza  _et al._ (2019). For the
dynamics, we begin with the classical Kuramoto model, a paradigmatic example
of the emergence of collective synchronization Kuramoto (2003); Pikovsky  _et
al._ (2001); Arenas  _et al._ (2008). An ensemble of $N$ heterogeneous
Kuramoto oscillators interacting on top of a network follows the equations of
motion

$\dot{\theta}_{i}=\omega_{i}+\lambda\sum_{j=1}^{N}a_{ij}\sin(\theta_{j}-\theta_{i}),\mbox{
}\forall\ i\in 1\dots N,$ (1)

where $\theta_{i}$ is the phase and $\omega_{i}$ is the intrinsic frequency of
the $i$-oscillator, $a_{ij}$ are the entries of the adjacency matrix $A$, that
capture the interactions among the units and $\lambda$ is a constant coupling
strength. As usual, the macroscopic behavior of the system is captured by the
modulus of the Kuramoto order parameter

$r(t)=\frac{1}{N}\left|\sum_{j=1}^{N}e^{i\theta_{j}(t)}\right|,$ (2)

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=\langle r(t)\rangle$. Although our results can be extended
to more general settings, in the following we restrict our study to the case
of unweighted ($a_{ij}=0,1$) and undirected networks $(a_{ij}=a_{ji})$, and
consider, for analytical convenience, a uniform frequency distribution
$g(\omega)\in[-\gamma,\gamma]$ with zero mean.

The growth of the synchronization bomb is made by keeping the coupling
strength $\lambda$ is 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 $(-\pi,\pi)$. 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

$\Delta
r_{ij}=\frac{1}{\lambda^{2}N}\left(\frac{\omega_{i}}{k_{i}}-\frac{\omega_{j}}{k_{j}}\right)\left(\frac{\omega_{i}}{k_{i}^{2}}-\frac{\omega_{j}}{k_{j}^{2}}\right).$
(3)

In practice, when connecting isolated nodes at the very initial steps of the
process, we add an infinitesimally small value to the degrees of the nodes
with $k_{i}=0$ to evaluate Eq. (3) in terms only of the natural frequencies.
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\rightarrow\infty$, and it reduces to
the random percolation case in $M=1$.

Eq. (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 IV 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 anti-correlations. These
correlations were imposed _ad hoc_ in previous models that induce explosive
synchronization D’Souza  _et al._ (2019), and they could indeed emerge from
applying a broader class of local percolation rules in the form
$p(\omega_{i},k_{i},\omega_{j},k_{j})$, 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 the left panels of Fig. 1, we illustrate the former basic mechanics for
assigning a link out of $M=5$ possible candidates. We take the forward process
(construction of the network by adding links) as an example. As shown in Fig.
1.a) 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) neighbours, 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 the panel Fig. 1.c) 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. 1.d) the transition from incoherence $(r\approx
0.05)$ to full phase-locking ($r\approx 0.9)$ when a unique link is added to
the system. This phenomenon motivates our choice for referring to these
growing networks as synchronization bombs.

Figure 2: Left column: a) Scatter plot of the tuples $(|\omega_{i}|,k_{i})$
and b) $(|\omega_{i}|,\langle\omega\rangle_{i})$ -where
$\langle\omega\rangle_{i}$ is the average frequency of the neighbours of the
$ith$-node- for each oscillator at $p=0.1$, in a single realization of the
forward process, for three values of noise and $N=200$, including in a) the
deterministic prediction of Eq. (4). Middle column: Scatter plot of (c) the
maximum and minimum eigenvalues of the Laplacian matrix, $\mu_{min}(L)$ and
$\mu_{max}(L)$ and (d) of the normalized Adjacency matrix,
${\mu}_{min}(\hat{A})$ and ${\mu}_{max}(\hat{A})$ at different $p$-steps of
the forward process (the size of the dots decreases for larger $p$) for
different values of $M$. In the inset, we plot the relative correlation
$\alpha_{i}=\langle\omega,v_{i}\rangle$ between the frequency vector and the
eigenvector of $L$ (or $\hat{A}$) associated with the maximum (or minimum)
eigenvalue. Right column: e) Evolution of the average size of giant component
(SGC), for noise/sampling ranging from the random scenario of $(M=1)$ to a
more deterministic one with $(M=20)$ in a network of size $N=5000$. Results
are averaged over $20$ realizations of the process. Inset in e) shows a single
realization in log scale. In f) we show the effect of network size on the
percolation transition for the $M=1$ (left) and $M=20$ (right) scenarios. It
is observed that larger and more deterministic networks under the rule of Eq.
(3) experiment sharper transitions.

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 transitions
D’Souza  _et al._ (2019); Arola-Fernández  _et al._ (2021).

_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 noise is reduced in the
process (for larger sampling $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 $\Delta r_{ij}$ 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 $\Delta r_{ij}$ among the potential links,
and frequencies and degrees become balanced in the precise way that makes
$\Delta r_{ij}$ more similar among these –still absent– links. Since the rule
predicts the scaling $\Delta r\sim\omega^{2}/k^{3}$, we find that, in the
deterministic (large $M$) regime (see Methods, section IV for details) the
relation is given by

$k_{i}\approx\dfrac{5}{3}pN\gamma^{-2/3}|\omega_{i}|^{2/3},$ (4)

where the scaling is controlled by the mean degree, expressed in terms of the
density of links $p$ and size $N$, since $\langle k\rangle\approx pN$. As
expected, Eq. (4) becomes more accurate as noise is reduced in the system, as
observed in Fig. 2.a). In Fig. 2.b), we see that frequency anti-correlations
among pairs of connected nodes are also present in the system and become
stronger for decreasing noise (large sampling $M$). Let us note that these
type of correlations are explicitly imposed in the majority of studied
mechanisms that induce explosive synchronization D’Souza  _et al._ (2019);
Boccaletti  _et al._ (2016) 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: towards optimal and bipartite networks._ We study the
evolution of the extreme eigenvalues $\mu_{max}$ and $\mu_{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}}AD^{-\frac{1}{2}}$ during
the percolation process, for different values of sampling $M$. In the central
panels of Fig. 2, we observe that the network evolves in a path that maximizes
both the largest positive eigenvalue of $L$, $\mu_{max}(L)$, ranging from zero
to $N$, and the largest negative eigenvalue of $\hat{A}$,
$\mu_{min}(\hat{A})$, ranging from minus one to zero, when noise 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 towards 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 $\mu_{max}(L)$ Skardal
_et al._ (2014). 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}_{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 towards
positive eigenvalues) or bipartite-like structures (shifting towards negative
eigenvalues) Arola-Fernández  _et al._ (2021). Thus, from Fig. 2.d) 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
towards 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 important 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. 2.e), 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. 2.f) 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 criterion Molloy and Reed (1995) 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(\omega)$, (see Methods IV)
that the threshold is estimated as

$p_{c}=\frac{42}{25N},$ (5)

which can be written in terms of the percolation threshold in a random network
Newman (2010) as $p_{c}\approx 1.68\cdot p_{c}^{rand}$. In Fig. 2.e) we
observe that Eq. (5) works quite well for sufficiently large $M$. More
sophisticated analytical tools, as the recently developed feature-enriched
percolation framework Artime and De Domenico (2021), 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, $\lambda$, and noise, $M$, values. It is important to 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 $\bm{\omega}$ Arenas  _et al._ (2008); D’Souza  _et al._ (2019).

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
systems Gottwald (2015) and extended to complex networks in Hancock and
Gottwald (2018). 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 IV 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. 3.a) for $\lambda=0.02$ and
$0.04$ ($M=10$).

Figure 3: a) Two examples of synchronization curves, $r(p)$ in the forward and
backward processes, for a rule dependent $(M=10)$ case in a network of size
$N=200$ and fixed coupling to $\lambda=0.02$ and $\lambda=0.04$. Measurements
here are taken every 40 links and results are obtained in a single realization
of the process. b) Synchronization phase-space depending on $M$ and $p$ for a
fixed $\lambda=0.05$. In the three panels, dashed (solid) lines correspond to
the theoretical predictions of the forward (backward) synchronization
thresholds, and circle markers in a) give the analytical prediction of the
whole backward curve (see main text and IV for the derivations). Note that the
results displayed in Fig. 1.c) correspond to three $M-$slices here. c)
Synchronization phase-space depending on $\lambda$ and $p$ for a certain level
of noise $M=10$, where the colormap indicates whether the system is in the
incoherent, bistable or phase-locking regime. Results are averaged over $20$
realizations.

For the forward process we cannot use the CC approach and we rely on the
celebrated OA ansatz Ott and Antonsen (2008), which has been successfully used
to characterize systems in the presence of frequency and degree correlations
Restrepo and Ott (2014); Skardal  _et al._ (2015); Peron  _et al._ (2020).
Specifically, we benefit from a recent elegant development used to describe
the mean-field dynamics of Janus oscillators Peron  _et al._ (2020) 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 IV. For the particular
case of a uniform distribution $g(\omega)\in[-1,1]$ we obtain the closed form

$p_{c}^{f}\approx\frac{21}{25\lambda N}.$ (6)

The predicted value $p_{c}^{f}$ is plotted in Fig. 3.a) showing again a
remarkable agreeement. This analytical estimation allows addressing the
aforementioned issue about the relation between synchronization and
percolation onsets by making use of Eq. (5) and Eq. (6). Combining both
expressions we can write a simple relation for the percolation $p_{c}$ and
forward synchronization $p_{c}^{f}$ thresholds as

$p_{c}\approx 2\lambda p_{c}^{f},$ (7)

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. 3.b), and $(p,\lambda)$-plane, Fig. 3.c).
In Fig. 3.b), we observe that, fixing $\lambda=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 emerges Skardal
and Arenas (2019). This critical point for which explosive behavior shows up
takes place around $M\approx 5$. In Fig. 3.c), we focus on the coupling
strength $\lambda$, 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 $\lambda$, to occur much later
for smaller values of $\lambda$ and to finally disappear for sufficiently
small $\lambda$. Interestingly, the system transits more abruptly for large
$\lambda$ (low $p$), but has a larger region of hysteresis for low $\lambda$
(large $p$).

In the SI, 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 show that both phenomenology and
theory are robust to changes in the distribution of intrinsic frequencies,
$g(\omega)$. 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(\omega)$
increases the bistable regime and the abruptness of the transitions.
Furthermore, we show that the bomb-like transitions also occur for directed
networks, and we analyze in more detail the role of the noisy sampling in the
model, finding that an optimal amount of noise can enhance the explosive
performance of the synchronization bomb because it improves the self-organized
optimization process driven by a local rule.

Chaotic synchronization bombs. One of most relevant applications of
synchronization theory is its implementation when coupling chaotic systems
Boccaletti  _et al._ (2002), a counter–intuitive nonlinear 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 dynamics Rössler (1976).

Figure 4: a) Example of synchronization curves for several values of $M$ and
$\lambda$. It is observed that a hysteresis cycle appears when $M>1$ and that
lower (higher) values of $\lambda$ translate into wider (more narrow) cycles
and less (more) abrupt transitions. b) Synchronization phase-space, depending
on $p$ and $M$ for a fixed $\lambda=0.02$. Results are qualitatively similar
to the ones found in Fig. 3.b)), although here the transitions are less abrupt
and narrower than in the Kuramoto case, and the birth of hysteresis occurs for
higher noise (lower $M$). c) Evolution of the oscillators trajectories in the
Rössler attractor for $M=10$ and $\lambda=0.02$ at two different $p$-steps,
before (left) and after (right) the forward synchronization transition. The
color bar corresponds to the frequency of the oscillator relative to the mean.
Eq.(8) is numerically integrated using Heun’s method, with $dt=0.05$ and
$10^{3}$ time steps and temporal averages of $r$ are taken at every $5$ link
changes.

Here we use an ensemble of diffusively coupled heterogeneous chaotic
oscillators Rosenblum  _et al._ (1996); Leyva  _et al._ (2012); Skardal  _et
al._ (2017), a modified, piece-wise linear Rössler system Rössler (1976),
which evolves in a 3-dimensional space following

$\begin{array}[]{c}\dot{x}_{i}=-f_{i}\left[\tau\left(x_{i}-\lambda\sum_{j=1}^{N}a_{ij}(x_{j}-x_{i})\right)+\beta
y_{i}+\delta z_{i}\right],\\\ \dot{y}_{i}=-f_{i}\left(-x_{i}+\nu
y_{i}\right),\\\ \dot{z}_{i}-f_{i}\left(-g(x_{i})+z_{i}\right),\end{array}$
(8)

where the non-linear function that induces the chaotic behavior is defined as
$g(x)=0$ if $x\leq 3$ and $g(x)=\mu(x-3)$ if $x>3$. The remaining parameters
are set following Leyva  _et al._ (2012); Skardal  _et al._ (2017), with
$\tau=0.05$, $\beta=0.5$, $\delta=1$, $\nu=0.02-100/R$. $R=100$ ensures that
the system is in a phase-coherent regime Rosenblum  _et al._ (1996); Leyva
_et al._ (2012); Skardal  _et al._ (2017), where a phase can be defined after
projecting onto the xy-plane, i.e. $\theta_{i}=arctan(y_{i}/x_{i})$, such that
the synchronization order parameter $r$ can be measured by the standard Eq.
(2). See Fig. 4.c) 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), $\lambda$ is the fixed coupling strength and the entries
$a_{ij}$ 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 $f_{i}$,
which we assign proportional to the frequency, $f_{i}=10+0.2\omega_{i}$, drawn
again from a uniform distribution $g(\omega)$ in $[-1,1]$.

Figure 4.a) illustrates three examples of synchronization transitions $r(p)$
for a system of $N=200$ and different choices of $\lambda$ and $M$. Similarly
to the Kuramoto case, it is observed how, in the construction process, for
noise values of $M>1$ the order parameter experiments abrupt jumps from
dynamical incoherence of $r\sim 0.1$ to a more coherent state with $r\simeq
0.7$, that continues to continuously grow to stronger synchronization
($r\gtrsim 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
panel 4.c) 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 optimal Skardal  _et al._ (2017) and explosive Leyva
_et al._ (2012) synchronization properties of the Kuramoto Model on the
diffusively coupled Rössler system. However, as numerical results in Fig.
4.a)-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 circuits Leyva  _et al._ (2012).

Application to cardiac pacemaker cells.

Lastly, we demonstrate the existence of self-organized explosive
synchronization via 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 contraction Djabella  _et al._ (2007)–. 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
concentrations Djabella  _et al._ (2007). For a system of $N$ such pacemakers
the equations of motion are given by

$\displaystyle\dot{v}_{i}$
$\displaystyle=\tau_{i}^{-1}f(v_{i},h_{i})+K_{v}\sum_{j=1}^{N}a_{ij}(v_{j}-v_{i}),$
(9) $\displaystyle\dot{h}_{i}$
$\displaystyle=\tau_{i}^{-1}g(v_{i},h_{i})+K_{h}\sum_{j=1}^{N}a_{ij}(h_{j}-h_{i}),$
(10)

where the local dynamics of each $v_{i}$ and $h_{i}$ are described by

$\displaystyle f(v,h)$
$\displaystyle=\frac{h(v+0.2)^{2}(1-v)}{0.3}-\frac{v}{6},$ (11) $\displaystyle
g(v,h)$ $\displaystyle=\frac{1}{150}+(8.333\times
10^{-4})[1-\text{sgn}(v-0.13)]$
$\displaystyle~{}~{}~{}~{}~{}~{}~{}~{}~{}~{}~{}~{}\times\\{0.5[1-\text{sgn}(v-0.13)]-h\\}.$
(12)

The timescales $\tau_{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
$K_{v}=0.009$ and $K_{h}=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. 5(a). Note that at roughly
$p\approx 0.02$ the system undergoes an explosive transition from (relatively)
large to small errors, indicating weak and strong synchronization. In Figs.
5(b) and (c) we present the actual voltage dynamics right before and after the
explosive transition, plotting each individual voltage time series $v_{i}(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. Lastly, note that, as with the
Kuramoto and Rössler dynamics, we observe a region of bistability where the
forward and backward percolation result in weak and strong synchronization
dynamics, respectively.

Figure 5: (a) Voltage error (quantified by the standard deviation) as a
function of the percolation parameter $p$, under forward and backward
percolation dynamics plotted in dot-dashed and solid lines, respectively. (b),
(c) Individual voltage time series $v_{i}(t)$ and the mean, plotted in light
and dark blue stokes, respectively, from right before and after the explosive
transition, as indicated in panel (a).

##  III 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 nonlinear 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 towards 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 and noise. 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 an
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. Importantly, 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.

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. 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, triggering the
bomb by means of a single local rule and imposing localized, instead of
global, perturbations in the system, deepens the explosive connection between
synchronization and competitive percolation Zhang  _et al._ (2014); D’Souza
_et al._ (2019), and it provides a missing explanation for the birth of abrupt
synchronization in pair-wise networks via a universal route Kuehn and Bick
(2021). By switching on a single additional parameter (the amount of sampling
in the percolation process), an oscillator network can self-organize towards a
high-dimensional correlated state where explosive behavior spontaneously
emerges.

##  IV Methods

Derivation of the local rule. We begin with two key assumptions: _i)_ the
system attempts to maximize the overall degree of synchronization, by adding
or removing undirected connections in a percolation process and _ii)_ 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 regime Dörfler  _et al._ (2013); Skardal  _et al._
(2014). Under the linearization, the resulting system reads in matrix form as

$\bm{\dot{\theta}}=\bm{\omega}-\lambda L\bm{\theta},$ (13)

where $L=D-A$ is the Laplacian of the network. The solution of Eq. (13) in the
stationary state is found by setting $\bm{\dot{\theta}}=0$. In the corotating
frame at speed $\langle w\rangle=0$, the solution reads as

$\bm{\theta}^{*}\approx\frac{1}{\lambda}{L}^{\dagger}\bm{\omega},$ (14)

where ${L}^{\dagger}$ is the Moore-Penrose pseudoinverse of the Laplacian
matrix, which can be constructed via the spectral decomposition of $L$ for
undirected networks (see Skardal  _et al._ (2014) for more details). Since we
are close to the synchronization attractor, the phases in Eq. (2) can also be
expressed in a Taylor expansion Skardal  _et al._ (2014). Invoking again
linearization, the order parameter is given by

$r\approx 1-||\bm{\theta}||^{2}/2N.$ (15)

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) Arola-Fernández  _et al._ (2021), 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
synchrony Arola-Fernández  _et al._ (2021) is obtained by truncating the
expansion at its second term (taking into account the effect of the nearest
neighbors of the nodes), leading to

$r\approx
1-\frac{1}{2\lambda^{2}N}\sum_{i=1}^{N}\left(\dfrac{\omega_{i}+z_{i}}{k_{i}}\right)^{2}.$
(16)

where $z_{i}=\sum_{j=1}^{N}a_{ij}(\omega_{j}/k_{j})$ is the contribution of
first neighbours. For more details on the accuracy of Eq. (16), see Arola-
Fernández  _et al._ (2021). From Eq. (16), we can estimate the local impact in
the synchrony of adding or removing a single link between oscillators $(p,q)$.
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 $z_{i}=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 $(\omega_{i},k_{i})$ of a given pair of
nodes. Explicitly

$\Delta r_{ij}=\frac{\pm
1}{\lambda^{2}N}\left(\frac{\omega_{i}}{k_{i}}-\frac{\omega_{j}}{k_{j}}\right)\left(\frac{\omega_{i}}{k_{i}^{2}}-\frac{\omega_{j}}{k_{j}^{2}}\right),$
(17)

where the $\pm$ sign accounts for the addition (removal) of a link. It is
important to 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 condition Arola-
Fernández  _et al._ (2021). 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-neighbours 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 criterion Molloy and Reed (1995), that
predicts the transition for random network without correlations for the value
of $p=p_{c}$ at which

$\langle k^{2}\rangle(p_{c})=2\langle k\rangle(p_{c}).$ (18)

To compute $\langle k\rangle$ and $\langle k^{2}\rangle$, we consider a
uniform distribution, such that $g(\omega)=1/(2\gamma)$ if
$\omega\in[-\gamma,\gamma]$ (the same analysis could be done for any other
frequency distribution) and also take into account that, explicitly, we have
the general correlation $k_{i}=c|\omega_{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(\omega)$. Using that $\langle
k\rangle=p(N-1)$, we find that $c=p(N-1)/\langle|\omega|^{2/3}\rangle$, with
$\langle|\omega|^{2/3}\rangle=\frac{3}{5}\gamma^{2/3}$. Thus we obtain the
correlation

$k_{i}\approx\frac{5}{3}pN\gamma^{-2/3}|\omega_{i}|^{2/3},$ (19)

and with it

$\langle
k^{2}\rangle\approx\dfrac{25}{9}p^{2}N^{2}\gamma^{-4/3}\langle|\omega|^{4/3}\rangle=\dfrac{25}{21}p^{2}N^{2}.$
(20)

Thus, substituting in Eq. (18) we obtain

$p_{c}\approx\dfrac{42}{25}\cdot\dfrac{1}{N}=1.68\cdot p_{c}^{rand},$ (21)

which corresponds to Eq. (5) in the results section.

Collective coordinates ansatz. We use the theory introduced in Gottwald
(2015); Hancock and Gottwald (2018) 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

$\theta_{i}=q(t)\psi_{i},$ (22)

where $\psi_{i}$ is the exact solution of the linearized dynamics of Eq. (1)
Skardal  _et al._ (2014), i.e.
$\psi_{i}=\frac{1}{\lambda}L^{\dagger}\bm{\omega}$. By minimizing the error
made by Eq. (3) in the full dynamics of Eq. (1) and after some manipulation
Hancock and Gottwald (2018), 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

$\dot{q}=1+\frac{1}{\psi^{T}L\psi}\sum_{i,j}\psi_{i}\sin(q(\psi_{j}-\psi_{i})).$
(23)

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,\lambda)$ 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 as
Hancock and Gottwald (2018)

$\displaystyle J_{ij}=a_{ij}\cos(\hat{q}(\psi_{j}-\psi_{i})),\ i\neq j$
$\displaystyle J_{ij}=-\sum_{k}a_{ik}\cos(\hat{q}(\psi_{j}-\psi_{i})),\ i=j.$
(24)

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 network Hancock and Gottwald (2018). 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\approx 0$ remains stable
beyond the backward critical transition, thus creating a bistable region and a
delayed forward transtion. In order to analytically predict the forward
critical threshold, we consider the limit of large $N$ and also large $M$
(towards 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 it is
important to remark that the theory is derived in the infinite size and
deterministic limits of the model. Our approach is based on the celebrated OA
ansatz Ott and Antonsen (2008) and follows a very similar development to that
shown in Peron  _et al._ (2020).

We begin by defining the local order parameter

$R_{i}=\sum_{j}a_{ij}e^{i\theta_{j}},$ (25)

such that Eq. (1) can be written as

$\dot{\theta}_{i}=\omega_{i}+Im[e^{i\theta_{i}}R_{i}]\ \forall\ i\in
1,\dots,N.$ (26)

Following Peron  _et al._ (2020), we consider a large ensemble of systems,
described by the joint probability density $\rho(\theta,\omega,t)$, with
$\theta=(\theta_{i},\dots,\theta_{N})$ and
$\omega=(\omega_{i},\dots,\omega_{N})$. The evolution of the joint probability
has to satisfy the continuity equation Ott and Antonsen (2008)

$\frac{\partial\rho}{\partial
t}+\sum_{i=0}^{N}\frac{\partial}{\partial\theta_{i}}(\rho\dot{\theta}_{i})=0.$
(27)

where $\theta_{i}$ is given by Eq. (26). Multiplying the density function
$\rho$ by $\prod_{j\neq i}d\omega_{j}d\theta_{j}$ and integrating, one obtains
the evolution for the marginal oscillator density,
$\rho_{i}(\theta_{i},\omega_{i},t)$ which reads as Peron  _et al._ (2020)

$\frac{\partial\rho_{i}}{\partial
t}+\frac{\partial}{\partial\theta_{i}}(\rho\dot{\theta_{i}})=0.$ (28)

Now, the OA ansatz can be applied by expanding $\rho_{i}$ in a Fourier series
and setting the coefficients of the expansion as $b_{i,n}=\alpha_{i}^{n}$ Ott
and Antonsen (2008); Peron  _et al._ (2020). By inserting the Fourier series
with the ansatz in Eq. (28), one ends up with

$\displaystyle\dot{\alpha}_{i}+i\alpha_{i}\omega_{i}+\frac{\lambda}{2}(\alpha_{i}^{2}R_{i}-R_{i}^{*})=0,\
\forall\ i\in 1,\dots,N$ (29) $\displaystyle
R_{i}=\sum_{j=1}a_{ij}\int_{-\infty}^{\infty}\alpha_{j}^{*}(\omega_{j},t)g(\omega)d\omega_{j},\
\forall\ i\in 1,\dots,N.$ (30)

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. 2.c)-d)). Also, in this limit, the frequencies of the oscillators are
completely determined by their degrees. Then, we can look for solutions
$\alpha_{i}=\alpha_{k,\pm}$ Peron  _et al._ (2020), 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 as Peron  _et al._
(2020)

$R_{k,\pm}=\frac{k}{\langle
k\rangle}\sum_{k^{\prime}}k^{\prime}p_{k^{\prime}}\alpha_{k^{\prime},\pm}^{*}.$
(31)

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

$\omega_{k,\pm}=\pm\left(\frac{k}{c}\right)^{3/2},$ (32)

where $c$ is a scaling constant given in Eq. (19). After these considerations,
the resulting system can be written as

$\begin{split}\dot{\alpha}_{k,+}=&-i\left(\frac{k}{c}\right)^{3/2}\alpha_{k,+}+\\\
&\frac{\lambda k}{2\langle
k\rangle}\left[\sum_{k^{\prime}}k^{\prime}p_{k^{\prime}}\alpha_{k^{\prime},-}-\alpha_{k^{\prime},+}^{2}\sum_{k^{\prime}}k^{\prime}p_{k^{\prime}}\alpha_{k^{\prime},-}^{*}\right]\\\
\end{split}$ (33)
$\begin{split}\dot{\alpha}_{k,-}=&+i\left(\frac{k}{c}\right)^{3/2}\alpha_{k,-}+\\\
&\frac{\lambda k}{2\langle
k\rangle}\left[\sum_{k^{\prime}}k^{\prime}p_{k^{\prime}}\alpha_{k^{\prime},+}-\alpha_{k^{\prime},-}^{2}\sum_{k^{\prime}}k^{\prime}p_{k^{\prime}}\alpha_{k,+}^{*}\right].\end{split}$
(34)

Since we want to evaluate the stability of the incoherent state
$\alpha_{k,\pm}=0$, we linearize the system above and evaluate it around
$\alpha_{k,\pm}=\delta\alpha_{k,\pm}\ll 1$. After neglecting smaller terms of
order $\delta\alpha^{2}$, the dependence on the complex conjugates vanish and
we up with the following system for each degree class

$\displaystyle\delta\dot{\alpha}_{k,+}=-i\left(\frac{k}{c}\right)^{3/2}\delta\alpha_{k,+}+\frac{\lambda
k}{2\langle
k\rangle}\sum_{k^{\prime}}k^{\prime}p_{k^{\prime}}\delta\alpha_{k^{\prime},-}$
(35)
$\displaystyle\delta\dot{\alpha}_{k,-}=+i\left(\frac{k}{c}\right)^{3/2}\delta\alpha_{k,-}+\frac{\lambda
k}{2\langle
k\rangle}\sum_{k^{\prime}}k^{\prime}p_{k^{\prime}}\delta\alpha_{k^{\prime},+}.$
(36)

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

$\displaystyle\sum_{k}kp_{k}\delta\dot{\alpha}_{k,+}=-i\sum_{k}\left(\frac{k}{c}\right)^{3/2}kp_{k}\delta\alpha_{k,+}+\sum_{k}\frac{\lambda
k^{2}p_{k}}{2\langle k\rangle}\delta y$ (38)
$\displaystyle\sum_{k}kp_{k}\delta\dot{\alpha}_{k,-}=+i\sum_{k}\left(\frac{k}{c}\right)^{3/2}kp_{k}\delta\alpha_{k,-}+\sum_{k}\frac{\lambda
k^{2}p_{k}}{2\langle k\rangle}\delta x.$ (39)

With the approximation $\sum_{k}k^{5/2}p_{k}\delta\alpha_{k,+}\approx\langle
k^{3/2}\rangle\delta x$ and
$\sum_{k}k^{5/2}p_{k}\delta\alpha_{k,-}\approx\langle k^{3/2}\rangle\delta y$,
the set of equations reduces to a $2$-dimensional variational system for the
evolution of $\delta x$ and $\delta y$ that reads as

$\displaystyle\delta\dot{x}=-\frac{i\langle k^{3/2}\rangle}{c^{3/2}}\delta
x+\frac{\lambda\langle k^{2}\rangle}{2\langle k\rangle}\delta y$ (40)
$\displaystyle\delta\dot{y}=\frac{\lambda\langle k^{2}\rangle}{2\langle
k\rangle}\delta x++\frac{i\langle k^{3/2}\rangle}{c^{3/2}}\delta y$ (41)

It is straightforward to show that the critical condition for the stability of
the incoherent state is given by

$c^{3/2}\lambda\langle k^{2}\rangle=2\langle k^{3/2}\rangle\langle k\rangle.$
(42)

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 parameter Peron  _et al._ (2020).
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.(42), 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
$k_{i}=c|\omega_{i}|^{2/3}$, and for a general $g(\omega)$ the constant is
given by $c=\langle k\rangle/\langle|\omega|^{2/3}\rangle$, we obtain a
general closed form for the forward critical threshold $(p_{c},\lambda_{c})$
that is given by

$p_{c}^{f}=\frac{2\langle|\omega|^{2/3}\rangle^{2}\langle|\omega|\rangle}{\lambda
N\langle|\omega|^{4/3}\rangle}.$ (43)

For the particular case of a uniform distribution
$g(\omega)\in[-\gamma,\gamma]$ we can easily compute the expected moments and,
after plugging these results in Eq. (43), we end up with the simple formula

$p_{c}^{f}=\frac{21\gamma}{25\lambda N}.$ (44)

which corresponds to Eq. (6) in the results section.

###### Acknowledgements.

L.A.-F. and A.A. acknowledge the Spanish MINECO (Grant No.
PGC2018-094754-B-C2). JGG 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).

##  V Supplementary Information

In this supplementary section, we study in more detail the robustness of the
explosive phenomena in our model against changes in the main parameters. In
particular, focusing on the Kuramoto dynamics, we study the synchronization
transitions occurring at single links for different system sizes to validate
the abrupt nature of the process when localized perturbations are applied to
the system. We measure the maximum jump in the order parameter and the area of
hysteresis depending on both size and noise (in terms of the sampling
parameter $M$), finding that the explosive phenomena is maximized in the large
size limit and, unexpectedly, by an optimal amount of noise that compensates
the error made by constantly applying a local –decentralized– percolation
rule. We also check that our numerical and analytical results hold for a
Gaussian and a bimodal distribution of frequencies (note that all results in
the main text were presented for a uniform choice of $g(\omega)$, which
simplified calculations). We conclude by showing results for the Kuramoto bomb
in directed networks, where we use the modified percolation rule that
correctly accounts for directionality in the chosen links while preserving the
decentralized nature of the process. Interestingly, we find that the explosive
synchronization transitions also occur in the directed case, although with a
bistability window that is smaller than in its undirected counterpart.

###  V.1 Effect of size and noise

First, we study both the effect of size $N$ and noise (via sampling parameter
$M)$ in the synchronization transition emerging from the self-organized
network growth. We validate that the abruptness of the transitions occurring
at single link changes is sustained for increasing size, such that $\Delta r$
in a single step does not vanish as size grows D’Souza  _et al._ (2019). In
fact, we observe in Fig. 6.a) and Fig. 7.a) that the mean maximum jump value
increases monotonically with size (towards red colors) and leads to a
macroscopic jump in $r$ at single link changes even for large system sizes.
The same occurs in terms of the hysteresis area (normalized by size for proper
comparison), which increases for large system sizes, as shown in Fig. 6.b) and
Fig. 7.b). Interestingly, the dependence on noise, via the sampling parameter
$M$ is clearly non-monotonous, showing a peak around $M\approx 50$ regardless
of size.

Figure 6:  a) Maximum difference of the average order parameter between two
consecutive link changes in the forward (top) and backward (bottom) directions
for increasing values of $M$ in log scale and for different sizes (ranging
from blue to red). b) Normalized hysteresis area (the sum of differences
between the values in the backward and forward curves) for increasing values
of $M$ and for different sizes. We observe the monotonous dependence with $N$
and the non-monotonous one with $M$. Results are averaged over $10$
realizations of the process. Each percolation process is run in both
directions from $p=0$ to $p=20/N$, corresponding to a maximum mean degree
$\langle k\rangle=20$, integrating the KM equations using Heun’s method, with
$dt=0.05$ and $10^{3}$ steps, discarding the first half for averaging $r$.
Coupling strength is set to $\lambda=0.05$ and $g(\omega)\in[-1,1]$.

The counter–intuitive effect of the local rule better sustaining
synchronization as noise is widely present (in terms of a low sampling $M$)
can be explained by noting that the rule is derived with local information
(see the Methods section in the main text), such that higher-order effects are
neglected by assumption. However, applying the rule itself makes higher-order
effects more important (inducing structural and dynamical anti–correlations).
The local prediction of $\Delta r$ for the sampled links may deviate from the
exact one as we advance in the percolation process, producing negative
feedback that penalizes the maximization of $\Delta r$ as the mechanism
becomes more deterministic, and a precise amount of noise leads to the optimal
performance. Luckily, the optimal value of $M$ in a particular setting can be
estimated by leveraging the analytical results presented in the main text,
without running the dynamics. Due to the proven goodness of the CC ansatz in
the explosive regime of our system, finding the $M$ that maximizes the degree
of synchrony in the linearized solution (which can be directly computed via
the aforementioned pseudo-inverse Laplacian) will turn out to be the $M$ that
maximizes explosive behavior. However, the reader should note that the
proposed model is intrinsically noisy, and the location of the synchronization
transitions may vary between different realizations of the process. Note that
alternative –deterministic– methods to build synchronization bombs could
minimize this uncertainty, but the current mechanism is intentionally designed
in the presence of noise. This noisy aspect turns out to be crucial to the
optimal performance of the bombs in our model.

Figure 7: a) Maximum difference of the average order parameter between two
consecutive link changes in the forward (top) and backward (bottom) directions
for increasing values of size $N$ in log scale and for different sampling or
noise $M$ (ranging from blue to red). b) Normalized hysteresis area (the sum
of differences between the values in the backward and forward curves) for
increasing values of $N$ and for different samplings. We observe again the
monotonous dependence with $N$ and the non-monotonous one with $M$.
Parametrization is the same as in the previous figure.

###  V.2 Effect of the frequency distribution

In the main text we presented our results for a particular choice of the
frequency distribution $g(\omega)$ -the uniform one- in order to simplify the
analytical treatment. Here we show that our model is robust to different
choices of the intrinsic frequencies of the oscillators. In particular, we
consider a Normal distribution and a bounded bimodal one, generated with a
Beta(0.1,0.1) distribution, a family of continuous probability distributions
defined on the interval [0,1], fixing the mean to zero and the variance to
$\sigma^{2}=1/3$, in order to compare against the uniform case in [-1,1] used
in the main text, which has the aforementioned variance.

Figure 8: Synchronization phase-space depending on $M$ and $p$ for a fixed
$\lambda=0.05$ for a Gaussian (a) and bimodal (b) distribution of frequencies
with mean zero and variance $1/3$. Dashed (solid) lines correspond to the
theoretical predictions of the forward (backward) synchronization thresholds
(see main text for the derivations).

In Fig. 8 we observe that a clear bistable region in the plane $(p,M)$ also
emerges for these choices of $g(\omega)$. The bistable region is larger for
the bimodal distribution, which shows that having a more polarized
distribution of frequencies enhances the explosivity in the system. In the
Gaussian case (less polarized than the uniform one), the bistable region is
much narrower, as can be seen in Fig. 8.a). Furthermore, in this case, the
prediction of the backward synchronization threshold (solid line) is less
accurate than in the other scenarios for low $M$ (high noise). This inaccuracy
can be explained by noting that the CC method Gottwald (2015); Hancock and
Gottwald (2018) used to predict the threshold is based on an _explosive trick_
that assumes that the whole system is in the phase-locking state before the
backward transition. This assumption does not hold for a Gaussian distribution
of $g(\omega)$, where the global phase-locking state is not supported by the
overall network (just by a large fraction of the oscillators) in the backward
process, and the value at which the full phase-locking state loses the
stability does not coincide with the backward synchronization threshold.
Nevertheless, as discussed in the Methods section in the main text, this
theory could be improved by finding the largest synchronized cluster of a
given size smaller than $N$, although this improvement demands larger
computational costs Hancock and Gottwald (2018). The forward prediction
(dashed line), based on the OA ansatz Ott and Antonsen (2008), does not suffer
from this issue and captures well the critical threshold even for high values
of noise (low $M$).

###  V.3 Extension to directed networks

We close the supplementary information by considering an important extension
of our initial results. In particular, in the main text we restricted our
model to undirected networks (assuming symmetric interactions, meaning that if
$a_{ij}=1$ then $a_{ji}=1$). Now we extend these results to a more general
setting by allowing directed connections, which may not necessarily be
symmetric.

Convenient to our purpose here, we can leverage the results in Arola-Fernández
_et al._ (2021), which exploited the truncated expansion of the linearized
synchronization dynamics in directed networks to predict the existence of
links leading to the counterintuitive Braess’ paradox in synchronization (a
removal of link that increase the degree of synchrony). Following Arola-
Fernández  _et al._ (2021), a modified version of our percolation rule (Eq.
(3) in the main text) that accounts for the directionality of links and
predicts the change of synchrony with local information can be written as

$\Delta
r_{ij}=\pm\frac{1}{Nk_{i}}\left[\frac{\omega_{p}}{k_{p}}\left(\frac{\omega_{j}}{k_{j}}-\frac{\omega_{j}}{k_{j}}\right)\right],$
(45)

where $\Delta r_{ij}$ accounts for the change in $r$ after adding (or
removing) a directed link coming from $j$ to $i$, and $k_{i}=\sum_{j}a_{ij}$
is the in-degree of the $i$-node.

Figure 9: a) Synchronization curves depending on $p$ for $M=2$ (black) and
$M=10$ (red). Measurements are taken at every single link. Parametrization is
the same as in Fig. 1 of the main text, but coupling strength is set to
$\lambda=0.025$. b) Phase-space in the $(p,M)$. Colormap shows the
corresponding values of $r$ in the forward and backward directions. Results
are averaged over 25 realizations of the process.

We wonder to which extent the explosive synchronization transitions found in
our model remain present in the directed scenario. We numerically find that
the bomb-like transitions indeed occur at single directed link changes when
the rule of Eq. (45) is applied. In Fig. 9.a), we plot two examples of the
synchronization curves in the forward and backward directions for two values
of $M$, and we observe an abrupt synchronization diagram, with its associated
hysteresis, that occurs for a sampling parameter $M=10$, but it is completely
absent for a much lower value $M=2$ (close to random directed percolation). In
Fig. 9.b) we show the phase-space depending on both the density $p$ and
sampling or noise $M$. Interestingly, we see that, for the coupling value
$\lambda=0.025$, the hysteresis window is quite small, and hysteresis behavior
is only present for $M\in[5,20]$. Nevertheless, the phenomenology is
qualitatively similar to the undirected case (see Fig. 3 of main text for a
proper comparison).

These results confirm that the synchronization bomb can be extended to
directed networks, which may represent a more realistic scenario, at least in
biological systems as the brain. A more detailed theoretical study of the
directed synchronization bomb, including the analyses of the structural
properties such as the percolation threshold or frequency-degree correlations,
and its extension to other dynamical processes, is left for further work.
Furthermore, the already known appearance of the Braess’ Paradox in directed
networks Arola-Fernández  _et al._ (2021) points towards the counter–intuitive
possibility of designing _reversed_ synchronization bombs, where the
transition from incoherence to global synchrony (or vice-versa) is induced by
the removal (or addition) of a single directed link.

## References

  * Joiner  _et al._ (2013) William J Joiner, Eliot B Friedman, Hsiao-Tung Hung, Kyunghee Koh, Mallory Sowcik, Amita Sehgal, and Max B Kelz, “Genetic and anatomical basis of the barrier separating wakefulness and anesthetic-induced unresponsiveness,” PLoS genetics 9, e1003605 (2013).
  * Kim  _et al._ (2016) Minkyung Kim, George A Mashour, Stefanie-Blain Moraes, Giancarlo Vanini, Vijay Tarnal, Ellen Janke, Anthony G Hudetz, and Uncheol Lee, “Functional and topological conditions for explosive synchronization develop in human brain networks with the onset of anesthetic-induced unconsciousness,” Frontiers in computational neuroscience 10, 1 (2016).
  * Kim  _et al._ (2017) Minkyung Kim, Seunghwan Kim, George A Mashour, and UnCheol Lee, “Relationship of topology, multiscale phase synchronization, and state transitions in human brain networks,” Frontiers in computational neuroscience 11, 55 (2017).
  * Wang  _et al._ (2017a) Chao-Qing Wang, Alain Pumir, Nicolas B Garnier, and Zong-Hua Liu, “Explosive synchronization enhances selectivity: Example of the cochlea,” Frontiers of Physics 12, 1–9 (2017a).
  * Wang  _et al._ (2017b) Zhenhua Wang, Changhai Tian, Mukesh Dhamala, and Zonghua Liu, “A small change in neuronal network topology can induce explosive synchronization transition and activity propagation in the entire network,” Scientific reports 7, 1–10 (2017b).
  * Lee  _et al._ (2018) UnCheol Lee, Minkyung Kim, KyoungEun Lee, Chelsea M Kaplan, Daniel J Clauw, Seunghwan Kim, George A Mashour, and Richard E Harris, “Functional brain network mechanism of hypersensitivity in chronic pain,” Scientific reports 8, 1–11 (2018).
  * Chatterjee  _et al._ (2008) Anushree Chatterjee, Yiannis N. Kaznessis, and Wei-Shou Hu, “Tweaking biological switches through a better understanding of bistability behavior,” Current Opinion in Biotechnology 19, 475–481 (2008), tissue, cell and pathway engineering.
  * Dobson  _et al._ (2007) Ian Dobson, Benjamin Carreras, V. Lynch, and D.E. Newman, “Complex systems analysis of series of blackouts: Cascading failure, critical points, and self-organization,” Chaos (Woodbury, N.Y.) 17, 026103 (2007).
  * Newman (2010) M. Newman,  _Networks: An Introduction_ (Oxford University Press, Inc., New York, NY, USA, 2010).
  * Achlioptas  _et al._ (2009) Dimitris Achlioptas, Raissa M. D’Souza, and Joel Spencer, “Explosive percolation in random networks,” Science 323, 1453–1455 (2009).
  * De Domenico  _et al._ (2016) Manlio De Domenico, Clara Granell, Mason Porter, and Alex Arenas, “The physics of spreading processes in multilayer networks,” Nature Physics 12 (2016), 10.1038/nphys3865.
  * Böttcher  _et al._ (2015) Lucas Böttcher, Olivia Woolley Meza, Nuno Araujo, Hans Herrmann, and Dirk Helbing, “Disease-induced resource constraints can trigger explosive,” Scientific reports 5 (2015), 10.1038/srep16571.
  * Matamalas  _et al._ (2020) Joan T. Matamalas, Sergio Gómez, and Alex Arenas, “Abrupt phase transition of epidemic spreading in simplicial complexes,” Phys. Rev. Research 2, 012049 (2020).
  * Gómez-Gardeñes  _et al._ (2016) Jesus Gómez-Gardeñes, Laura Lotero-Vélez, S. Taraskin, and Francisco Pérez-Reche, “Explosive contagion in networks,” Scientific Reports 6, 19767 (2016).
  * Echenique, P.  _et al._ (2005) Echenique, P., Gómez-Gardeñes, J., and Moreno, Y., “Dynamics of jamming transitions in complex networks,” Europhys. Lett. 71, 325–331 (2005).
  * Lampo  _et al._ (2021) Aniello Lampo, Javier Borge-Holthoefer, Sergio Gómez, and Albert Solé-Ribalta, “Multiple abrupt phase transitions in urban transport congestion,” Phys. Rev. Research 3, 013267 (2021).
  * D’Souza  _et al._ (2019) Raissa M. D’Souza, Jesus Gómez-Gardeñes, Jan Nagler, and Alex Arenas, “Explosive phenomena in complex networks,” Advances in Physics 68, 123–223 (2019), https://doi.org/10.1080/00018732.2019.1650450 .
  * Boccaletti  _et al._ (2016) S. Boccaletti, J.A. Almendral, S. Guan, I. Leyva, Z. Liu, I. Sendiña-Nadal, Z. Wang, and Y. Zou, “Explosive transitions in complex networks’ structure and dynamics: Percolation and synchronization,” Physics Reports 660, 1–94 (2016), explosive transitions in complex networks’ structure and dynamics: Percolation and synchronization.
  * Pikovsky  _et al._ (2001) A. Pikovsky, M. G. Rosenblum, and J. Kurths,  _Synchronization, A Universal Concept in Nonlinear Sciences_ (Cambridge University Press, 2001).
  * Arenas  _et al._ (2008) A Arenas, A. Díaz-Guilera, J. Kurths, Y Moreno, and C. Zhou, “Synchronization in complex networks,” Phys. Rep. 469, 93–153 (2008).
  * Pazó (2005) Diego Pazó, “Thermodynamic limit of the first-order phase transition in the kuramoto model,” Phys. Rev. E 72, 046211 (2005).
  * Martens  _et al._ (2009) Erik Martens, Ernest Barreto, S Strogatz, E Ott, Paul So, and T Antonsen, “Exact results for the kuramoto model with a bimodal frequency distribution,” Physical review. E, Statistical, nonlinear, and soft matter physics 79, 026204 (2009).
  * Gómez-Gardeñes  _et al._ (2011) Jesús Gómez-Gardeñes, Sergio Gómez, Alex Arenas, and Yamir Moreno, “Explosive synchronization transitions in scale-free networks,” Phys. Rev. Lett. 106, 128701 (2011).
  * Leyva  _et al._ (2013) I. Leyva, Adrián Navas, Irene Sendiña-Nadal, J Almendral, Javier Buldu, Massimiliano Zanin, David Papo, and Stefano Boccaletti, “Explosive transitions to synchronization in networks of phase oscillators,” Scientific reports 3, 1281 (2013).
  * Avalos-Gaytán  _et al._ (2018) Vanesa Avalos-Gaytán, Juan A. Almendral, I. Leyva, F. Battiston, V. Nicosia, V. Latora, and S. Boccaletti, “Emergent explosive synchronization in adaptive complex networks,” Phys. Rev. E 97, 042301 (2018).
  * Arola-Fernández  _et al._ (2021) Lluís Arola-Fernández, Per Sebastian Skardal, and Alex Arenas, “Geometric unfolding of synchronization dynamics on networks,” Chaos 31, 061105 (2021).
  * Wei Chen and Xiao (2021) Yueheng Lan Weiqing Liu Wei Chen, Shengfeng Wang and Jinghua Xiao, “Explosive synchronization caused by optimizing synchrony of coupled phase oscillators on complex networks,” Eur. Phys J. B 94, eabe3824 (2021).
  * Zhang  _et al._ (2015) Xiyun Zhang, Stefano Boccaletti, Shuguang Guan, and Zonghua Liu, “Explosive synchronization in adaptive and multilayer networks,” Phys. Rev. Lett. 114, 038701 (2015).
  * Soriano-Paños  _et al._ (2019) D. Soriano-Paños, Q. Guo, V. Latora, and J. Gómez-Gardeñes, “Explosive transitions induced by interdependent contagion-consensus dynamics in multiplex networks,” Phys. Rev. E 99, 062311 (2019).
  * Skardal and Arenas (2014) Per Sebastian Skardal and Alex Arenas, “Disorder induces explosive synchronization,” Physical Review E 89, 062811 (2014).
  * Skardal and Arenas (2019) Per Sebastian Skardal and Alex Arenas, “Abrupt desynchronization and extensive multistability in globally coupled oscillator simplexes,” Phys. Rev. Lett. 122, 248301 (2019).
  * Zhang  _et al._ (2014) Xiyun Zhang, Yong Zou, Stefano Boccaletti, and Zonghua Liu, “Explosive synchronization as a process of explosive percolation in dynamical phase space,” Scientific reports 4, 5200 (2014).
  * Myung  _et al._ (2018) Jihwan Myung, Christoph Schmal, Sungho Hong, Yoshiaki Tsukizawa, Pia Rose, Yong Zhang, Michael J Holtzman, Erik De Schutter, Hanspeter Herzel, Grigory Bordyugov, _et al._ , “The choroid plexus is an important circadian clock component,” Nature communications 9, 1–13 (2018).
  * Ishida (1997) Yoshiteru Ishida, “The immune system as a prototype of autonomous decentralized systems: an overview,” in  _Proceedings of the Third International Symposium on Autonomous Decentralized Systems. ISADS 97_ (IEEE, 1997) pp. 85–92.
  * Izhikevich (2003) Eugene M Izhikevich, “Simple model of spiking neurons,” IEEE Transactions on neural networks 14, 1569–1572 (2003).
  * Orlandi  _et al._ (2013) Javier G Orlandi, Jordi Soriano, Enrique Alvarez-Lacalle, Sara Teller, and Jaume Casademunt, “Noise focusing and the emergence of coherent activity in neuronal cultures,” Nature Physics 9, 582–590 (2013).
  * Zhang  _et al._ (2016) Wen-Hao Zhang, Aihua Chen, Malte J Rasch, and Si Wu, “Decentralized multisensory information integration in neural systems,” Journal of Neuroscience 36, 532–547 (2016).
  * Scarpetta  _et al._ (2018) Silvia Scarpetta, Ilenia Apicella, Ludovico Minati, and Antonio de Candia, “Hysteresis, neural avalanches, and critical behavior near a first-order transition of a spiking neural network,” Physical Review E 97, 062305 (2018).
  * Gottwald (2015) Georg Gottwald, “Model reduction for networks of coupled oscillators,” Chaos (Woodbury, N.Y.) 25 (2015), 10.1063/1.4921295.
  * Hancock and Gottwald (2018) Edward Hancock and Georg Gottwald, “Model reduction for kuramoto models with complex topologies,” Physical Review E 98 (2018), 10.1103/PhysRevE.98.012307.
  * Ott and Antonsen (2008) Edward Ott and Thomas Antonsen, “Low dimensional behavior of large systems of globally coupled oscillators,” Chaos (Woodbury, N.Y.) 18, 037113 (2008).
  * Kuramoto (2003) Y. Kuramoto,  _Chemical Oscillations, Waves, and Turbulence_ (Dover Publications, 2003).
  * Skardal  _et al._ (2014) P. S. Skardal, D. Taylor, and J. Sun, “Optimal Synchronization of Complex Networks,” Phys. Rev. Lett. 113, 144101 (2014).
  * Molloy and Reed (1995) Michael Molloy and Bruce Reed, “A critical point for random graphs with a given degree sequence,” Random Structures & Algorithms 6, 161–180 (1995).
  * Artime and De Domenico (2021) Oriol Artime and Manlio De Domenico, “Percolation on feature-enriched interconnected systems,” Nature Communications 12 (2021), 10.1038/s41467-021-22721-z.
  * Restrepo and Ott (2014) Juan Restrepo and Edward Ott, “Mean field theory of assortative networks of phase oscillators,” EPL (Europhysics Letters) 107 (2014), 10.1209/0295-5075/107/60006.
  * Skardal  _et al._ (2015) Per Sebastian Skardal, Juan G. Restrepo, and Edward Ott, “Frequency assortativity can induce chaos in oscillator networks,” Phys. Rev. E 91, 060902 (2015).
  * Peron  _et al._ (2020) Thomas Peron, Deniz Eroglu, Francisco Rodrigues, and Yamir Moreno, “Collective dynamics of random janus oscillator networks,” Physical Review Research 2 (2020), 10.1103/PhysRevResearch.2.013255.
  * Boccaletti  _et al._ (2002) S. Boccaletti, J. Kurths, G. Osipov, D.L. Valladares, and C.S. Zhou, “The synchronization of chaotic systems,” Physics Reports 366, 1–101 (2002).
  * Rössler (1976) O.E. Rössler, “An equation for continuos chaos,” Phys. Lett. A 57, 397 (1976).
  * Rosenblum  _et al._ (1996) M.G. Rosenblum, Arkady Pikovsky, and Juergen Kurths, “Phase synchronization of chaotic oscillators,” Physical review letters 76, 1804–1807 (1996).
  * Leyva  _et al._ (2012) I. Leyva, Ricardo Sevilla-Escoboza, Javier Buldu, Irene Sendiña-Nadal, Jesus Gómez-Gardeñes, Alex Arenas, Yamir Moreno, Sergio Gomez, Rider Jaimes-Reategui, and Stefano Boccaletti, “Explosive first-order transition to synchrony in networked chaotic oscillators,” Physical review letters 108, 168702 (2012).
  * Skardal  _et al._ (2017) P.S. Skardal, V.P. Sevilla-Escoboza, V.P. Vera-Ávila, and J.M. Buldú, “Optimal phase synchronization in networks of phase-coherent chaotic oscillators,” Chaos 27, 013111 (2017).
  * Djabella  _et al._ (2007) Karima Djabella, Mayer Landau, and Michel Sorine, “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_ (2007) pp. 5186–5191.
  * Kuehn and Bick (2021) Christian Kuehn and Christian Bick, “A universal route to explosive phenomena,” Science Advances 7, eabe3824 (2021).
  * Dörfler  _et al._ (2013) Florian Dörfler, Michael Chertkov, and Francesco Bullo, “Synchronization in complex oscillator networks and smart grids,” Proceedings of the National Academy of Sciences 110, 2005–2010 (2013), https://www.pnas.org/content/110/6/2005.full.pdf .