Higher-order interactions in complex networks of phase oscillators promote abrupt synchronization switching

Synchronization processes play critical roles in the functionality of a wide range of both natural and man-made systems. Recent work in physics and neuroscience highlights the importance of higher-order interactions between dynamical units, i.e., three- and four-way interactions in addition to pairwise interactions, and their role in shaping collective behavior. Here we show that higher-order interactions between coupled phase oscillators, encoded microscopically in a simplicial complex, give rise to added nonlinearity in the macroscopic system dynamics that induces abrupt synchronization transitions via hysteresis and bistability of synchronized and incoherent states. Moreover, these higher-order interactions can stabilize strongly synchronized states even when the pairwise coupling is repulsive. These findings reveal a self-organized phenomenon that may be responsible for the rapid switching to synchronization in many biological and other systems that exhibit synchronization without the need of particular correlation mechanisms between the oscillators and the topological structure.

The collective dynamics of network-coupled dynamical systems has been a major subject of research in the physics community during the last decades [1][2][3][4] due to a wide range of applications applications including cardiac rhythms [5], power grid dynamics [6], and proper cell circuit behavior [7]. In particular, our understanding of both natural and man-made systems has significantly improved by studying how network structures and dynamical processes combine to shape overall system behaviors. This interplay gives rise to novel nonlinear phenomena like switch-like abrupt transitions to synchronization [8][9][10] and cluster states [11,12]. Recent work in physics and neuroscience have specifically highlighted the importance of higher-order interactions between dynamical units, i.e., three-and four-way interactions in addition to pairwise interactions, and their role in shaping collective behavior [13][14][15][16][17][18][19][20], prompting the network science community to turn its attention to higher-order structures to better represent the kinds of interactions that one can find beyond typical pairwise interactions [21]. These higher-order interactions are often encoded in simplicial complexes [22] that describe the different kinds of simplex structure present in the network: a filled clique of m + 1 nodes is known as an m-simplex, and together a set of 1-simplexes (links), 2-simplexes (filled triangles), etc. comprise the simplicial complex. While simplicial complexes have been proven to be very useful for analysis and computation in high dimensional data sets, e.g., using persistent homologies [17], little is understood about their role in shaping dynamical processes, save for a handful of examples [23][24][25].
In parallel to the previous developments, there has been also a lot of attention on another phenomena related to the collective dynamics of network-coupled oscillators namely the explosive synchronization phenomenon, see [10] and references therein. Explosive synchronization consists of an abrupt switch between incoherent and synchronized states, that can be achieved by the interplay between the network structure and the oscillators dynamics, being the most simple prescription that of each oscillator having a natural frequency proportional to the number of connections in the network. This mathematical finding is becoming particularly important in neuroscience, where bistability and fast switching of states are very relevant to understand, bistable perception [26], epileptic seizures in the brain [27,28], or hypersensitivity in chronic pain of Fibromyalgia patients [29]. However, the mechanisms for this abrupt switching to happen are still unclear.
The specificities of the networks should not be the most relevant parameter, given that human wiring is not equivalent between individuals [30], and then we rely on another aspect, the higher-order (beyond pairwise) interactions in the network. The collective dynamics of sources and loads in large-scale power grids provides another important application where abrupt synchronization transitions play an important role [31].
Motivated by the above mentioned dynamical processes, here we study the dynamics of heterogeneous phase oscillators with higher-order interactions on simplicial complexes with 1-, 2-, and 3-simplex interactions. Similarly to our recent approach in [9], where we investigate the desynchronization phenomena taking in to account 2-simplexes, here we aim to understand the effect of higher-order interactions that combine 1-, 2-, and 3-simplex interactions in the emergence of synchronization. In this previous study we showed that although 2-simplex interactions alone did not lead to any synchronization transition (i.e., they do not destabilize the incoherent state) the synchronization, they do give rise to abrupt desynchronization transitions. Here we show that the combination of multiple higher-order interactions gives rise to both abrupt synchronization and desynchronization transitions, allowing the system to easily switch between incoherent and synchronized states with relatively small changes to system parameters. We use the celebrated Kuramoto model [32] to scrutinize the higher order dynamics in complex networks. Previous studies already revealed a rich phase diagram where nonpairwise interactions are considered, showing multi-stability [33], quasiperiodicity [34], and even chaos [35]. Our contribution aligns with these previous works and demonstrates that higher-order interactions provide a natural mechanism for the emergence of explosive synchronization.
For a simplicial complex of N nodes we propose an extension of the Kuramoto-Sakaguchi phase rotator model [36] on networks to the higher-order Kuramoto model whose equations of motion are given bẏ where θ i is the phase of oscillator i, ω i is its natural frequency (typically assumed to be drawn from a distribution g(ω)), and K 1 , K 2 , and K 3 are the coupling strengths of 1-, 2-, and 3-simplex interactions, respectively. Importantly, these addition forms of coupling (i.e., those with K 2 and K 3 coefficients) come directly from higher-order terms that emerge from phase-reductions of limit-cycle oscillators [15,20]. The network structure (assumed to be undirected and unweighted) is encoded in the 1-simplex adjacency matrix A, 2-simplex adjacency tensor B, and 3-simplex adjacency tensor C, where A ij = 1 if nodes i and j are connected by a link (and otherwise A ij = 0), B ijl = 1 if nodes i, j, and l belong to a common 2-simplex (and otherwise B ijl = 0), and C ijlm = 1 if nodes i, j, l, and m belong to a common 3-simplex (and otherwise C ijlm = 0). For each node i we denote the q-simplex degree k q i as the number of distinct q-simplexes node i is a part of, and k q is the mean q-simplex degree across the network. (Note that each division by k q in equation (1) amounts to a rescaling of the respective coupling strength).
Taking inspiration from the importance of simplicial complexes in the brain, which displays rich synchronization dynamics [37], we consider as a motivating example the dynamics of equation (1) on the Macaque brain dataset which consists of 242 interconnected regions of the brain [38]. The adjacency matrix A is taken to be undirected and 2-and 3simplex structures are constructed by identifying each distinct triangle and tetrahedron from the 1-simplex structures. The 2-and 3-simplex coupling strengths are fixed to K 2 = 1.6 and K 3 = 1.1 as the 1-simplex coupling strength is varied and natural frequencies are drawn identically and independently from the standard normal distribution. In Fig. 1(a) we plot the amplitude r of the complex order parameter z = re iψ = N −1 N j=1 e iθj as K 1 is first increased adiabatically from K 1 = −0.6 to 0.4, then decreased. These simulations reveal that the presence of higher-order interactions in simplicial complexes give rise to abrupt (a.k.a. explosive) synchronization transitions [8], as the system quickly transitions Abrupt synchronization in simplicial complexes: Macaque brain and UK power grid networks. (a) The synchronization profile describing the macroscopic system state by the order parameter r as a function of 1-simplex coupling K1 for higher-order coupling strengths K2 = 1.6 and K3 = 1.1 using the Macaque brain network. Results are obtained by adiabatically increasing K1 from −0.6 to 0.4, then subsequently decreasing K1 from 0.4 back to −0.6. This protocol reveals a hysteresis loop with abrupt synchronization and desynchronization transitions at K sync from the incoherent state (r ≈ 0) to a partially synchronized state (r ∼ 1) at K sync 1 ≈ 0.25 as K 1 is increased, then another abrupt transition from synchronization to incoherence occurs at K desync the system admits a bistability where both incoherent and synchronized states are stable. In Figs. 1(b) and (c) we highlight this bistabiliy by showing the incoherent and synchronized states, respectively, for K 1 = 0.1, illustrating for 40% of the oscillators (chosen randomly) placed appropriately on the unit circle with their respective order parameter values r ≈ 0.07 and 0.46.
The results presented above illustrate two new critical findings using a real brain dataset. First, the presence of higherorder interactions, i.e., 2-and 3-simplexes, can induce abrupt synchronization transitions without any additional dynamical or structural ingredients. Incoherent and synchronized states have been mapped to resting and active states of the   (5), respectively, and circles denote results taken from direct simulations of equation (2) with N = 10 4 oscillators. (c) The full stability diagram describing incoherent, synchronized, and bistable states as a function of 1-simplex coupling K1 and higher-order coupling K2+3. Blue and red curves correspond to pitchfork and saddle-node bifurcations, which collide at a codimension-two point (black circle) at (K1, K2+3) = (2, 2). For K2+3 < 2 and K2+3 > 2 the pitchfork bifurcation is supercritical and subcritical, respectively. brain [39], respectively, with abrupt transitions representing quick and efficient mechanisms for switching cognitive tasks. However, previous work has shown that in the presence of only 1-simplex coupling, properties such as time-delays [40] or degree-frequency correlations [8] are needed to induce such transitions. Second, the presence of higher-order interactions can create and stabilize a synchronized state even when 1simplex coupling is negative, i.e., repulsive. Thus, higherorder interactions nonlinear effects that support synchronization on the macroscopic scale. To emphasize the broader implications of this finding, we plot in Fig. 1(d) the synchronization profile of the order parameter r vs 1-simplex coupling K 1 (again both increasing and decreasing K 1 to highlight the explosive transitions and bistability) for higher-order coupling strengths K 2 = 2.2 and K 3 = 3.3 on the UK power grid network [6]. Here, since the network is strongly geometric and adjacent nodes are geographically close to one another, and therefore likely similarly affected by local events, we identify 2-simplexes as 3-paths (i.e., paths of three connected nodes, a.k.a., wedges) and 3-simplexes as 4-paths and 4-stars (i.e., three nodes all connected to a fourth central node). The qualitatively similar behavior displayed here demonstrates a wide range of important synchronization applications where higher-order interactions may significantly affect the dynamics.
To better understand the dynamics that emerge in the system above, we turn our focus to a population of all-to-all coupled oscillators. The governing equations, which also serves as the mean-field approximation for equation (1), is given bẏ In the all-to-all case given by equation (2) the system can be treated using the dimensionality reduction of Ott and Antonsen [41], yielding a low dimensional system that governs the macroscopic dynamics via the order parameter z = re iψ . In particular, by considering the continuum limit of infinitelymany oscillators and applying the Ott-Antonsen ansatz (see Methods for details), we obtain for the amplitude r and angle ψ the simple differential equationṡ where we have assumed that the natural frequency distribution g(ω) is Lorentzian with mean ω 0 and the new coupling strength is given by the sum of the 2-and 3-simplex coupling strengths, i.e., K 2+3 = K 2 + K 3 . Note first that the amplitude and angle dynamics of r and ψ completely decouple and that the angle dynamics evolve with a constant angular velocity equal to the mean of the frequency distribution. Thus, by entering an appropriate rotating frame and shifting initial conditions we may set ψ = 0 without any loss of generality. Moreover, the higher-order interactions, i.e., 2-and 3-simplexes mediated by the coupling strength K 2+3 , surface in the form of cubic and quintic nonlinear terms. This implies that the stability of the incoherent state, given by r = 0, (which is always an equilibrium) is not affected by the higherorder interactions. However, these nonlinear terms that originate from the higher-order interactions mediate the possibility of synchronized states. In particular, one or two synchronized states also exists, given by where the plus and minus signs correspond to stable and unstable solutions when they exist.
We now show that the all-to-all case illustrates, in an analytically tractable setting, all the novel dynamics observed in the Macaque example (see Fig. 1). First, in Fig. 2(a) we plot steady-state solutions of the order parameter r as a function of the 1-simplex coupling strength K 1 for a variety of higherorder coupling strengths K 2+3 = 0, 2, 5, 8, and 10 (blue to red). Analytical predictions given by equation 5 are plotted as solid and dashed curves (for stable and unstable branches, respectively), and circles represent results from direct simulation of equation (1) with N = 10 4 oscillators. For sufficiently small higher-order coupling (e.g., K 2+3 = 0) the transition to synchronization is second-order, occurring via a supercritical pitchfork bifurcation. However, as K 2+3 is increased through a critical value of K sync 1 = 2 the synchronized branch folds over itself, giving rise to hysteresis and abrupt transitions between incoherence and synchronization for larger values of higher-order coupling (e.g., K 2+3 = 5, 8, and 10). In this regime the pitchfork bifurcation at K sync 1 = 2 becomes subcritical and a saddle-node bifurcation emerges at a lower value of K 1 , denoted K desync 1 , where the synchronized branch first appears. These two bifurcations correspond to the abrupt transitions observed in Fig. 1. We also observe that for K 2+3 ≥ 8 the synchronized branch stretches into the negative region K 1 < 0 (e.g., K 2+3 = 10), again demonstrating that higher-order interactions can stabilize synchronized states even when pairwise interactions are repulsive. In Fig. 2(b) we plot similar results as the higher-order coupling strength K 2+3 is varied for a variety of 1-simplex coupling strengths, K 1 = −0.5, 1, 1.8, 2, and 2.2 (blue to red). These curves highlight the existence and absence of bistability for K 1 < 2 and K 1 > 2, respectively. In Fig. 2(c) we provide the full stability diagram for the system, denoting the pitchfork bifurcations at K sync 1 = 2 (supercritical and subcritical for K 2+3 < 2 and K 2+3 > 3) in blue and the saddle-node bifurcation, given by K desync The region bounded by these curves corresponds to bistability between synchronization and incoherence, and is born at the intersection between the two bifurcations at the codimensiontwo point (K 1 , K 2+3 ) = (2,2).
Having demonstrated the novel synchronization dynamics that arise from higher-order interactions in simplicial complexes in a real brain dataset and the all-to-all scenario, we lastly turn to a synthetic network example, constructing a sim- plicial complex via a three-layer multiplex, where the q th layer consists of q-simplexes. In particular, aiming for such a multiplex with mean degrees k 1 , k 2 , and k 3 , we construct each layer randomly, placing M 1 = N k 1 /2 1-simplexes (i.e., links) in the first layer, M 2 = N k 2 /3 2-simplexes (i.e., filled triangles) in the second layer, and M 3 = N k 3 /4 3-simplexes (i.e., filled tetrahedra) in the third layer. (Note that the first layer is a classical Erdős-Rényi network [42] and the second and third layers are the generic extensions using 2-and 3-simplexes instead of typical links.) In Figs. 3(a) and (b) we plot the the order parameter r vs 1-simplex coupling K 1 and higher-order coupling K 2+3 , resecptively, for a multiplex network of N = 10 4 oscillators with mean degrees k 1 = k 2 = k 3 = 30 in circles. Similar to Figs. 2(a) and (b), solid and dashed curves represent the analytical results for the mean-field approximation from the all-to-all case. These results illustrate that the mean-field approximation accurately describes the dynamics of such randomly generated simplicial complexes.
The results presented above demonstrate that higher-order interactions in networks of coupled oscillators, which are encoded on the microscopic scale of by a simplicial complex, give rise to added nonlinearities in the macroscopic system dynamics. These nonlinearities give rise to two new phenomena that are not present in the absence of higher-order interactions, i.e., when interactions are solely pairwise. First, these nonlinearities induce abrupt transitions between incoherent and synchronized states without additional characteris-tics such as time delays or network-dynamics correlations. In the context of brain dynamics, incoherent and synchronized states correspond to resting and active states, with abrupt transitions facilitating efficient switching between cognitive tasks [39]. Second, when nonlinearities are sufficiently strong they create and stabilize synchronized states even when pairwise coupling is repulsive. Thus, even as certain kinds of coupling may degrade over time due to synaptic plasticity, the presence of other kinds of coupling may be enough to sustain bistability regimes between incoherence and synchronization. We note that in this paper we have taken brain dynamics as our primary motivating example due to the existence of direct evidence of higher-order interactions in a system with synchronization properties [13,14,16,19]. However, more general results suggest that higher-order interactions may be important in broader classes of physical systems [15,20] including large-scale power grids [46], indicating that the nonlinear phenomena observed in this context may point to other novel behaviors that arise from such interactions in different contexts.