wp header logo 660

Emergent hypernetworks in weakly coupled oscillators – Nature.com

Local News MP News

Thank you for visiting nature.com. You are using a browser version with limited support for CSS. To obtain the best experience, we recommend you use a more up to date browser (or turn off compatibility mode in Internet Explorer). In the meantime, to ensure continued support, we are displaying the site without styles and JavaScript.
Carousel with three slides shown at a time. Use the Previous and Next buttons to navigate three slides at a time, or the slide dot buttons at the end to jump three slides at a time.
27 November 2020
Ricardo Gutiérrez, Massimo Materassi, … Stefano Boccaletti
29 October 2022
Lluís Arola-Fernández, Sergio Faci-Lázaro, … Alex Arenas
30 November 2020
Per Sebastian Skardal & Alex Arenas
04 January 2021
Ralf Tönjes, Carlos E. Fiore & Tiago Pereira
19 March 2021
David S. Glass, Xiaofan Jin & Ingmar H. Riedel-Kruse
19 April 2022
Ying Tang, Dinghua Shi & Linyuan Lü
28 October 2022
Luca Gallo, Riccardo Muolo, … Timoteo Carletti
19 January 2023
Jannik Luboeinski, Luis Claro, … Eduardo Mizraji
26 August 2021
Yuanzhao Zhang, Vito Latora & Adilson E. Motter
Nature Communications volume 13, Article number: 4849 (2022)
2584 Accesses
1 Citations
26 Altmetric
Metrics details
Networks of weakly coupled oscillators had a profound impact on our understanding of complex systems. Studies on model reconstruction from data have shown prevalent contributions from hypernetworks with triplet and higher interactions among oscillators, in spite that such models were originally defined as oscillator networks with pairwise interactions. Here, we show that hypernetworks can spontaneously emerge even in the presence of pairwise albeit nonlinear coupling given certain triplet frequency resonance conditions. The results are demonstrated in experiments with electrochemical oscillators and in simulations with integrate-and-fire neurons. By developing a comprehensive theory, we uncover the mechanism for emergent hypernetworks by identifying appearing and forbidden frequency resonant conditions. Furthermore, it is shown that microscopic linear (difference) coupling among units results in coupled mean fields, which have sufficient nonlinearity to facilitate hypernetworks. Our findings shed light on the apparent abundance of hypernetworks and provide a constructive way to predict and engineer their emergence.
Networks of weakly coupled oscillators are prolific models for a variety of natural systems ranging from biology1,2 and chemistry3,4 to neuroscience5,6 via ecology7 to engineering8. Such networks serve as stepping stones to understand collective dynamics9,10,11,12 and other emergent phenomena in networks13,14. In these models, the interactions are described in a pairwise manner and the collective dynamics of a network can be predicted by the superposition of such pairwise interactions.
Recent work, however, suggests that many networks described as pairwise interactions can be better described in terms of hypernetworks with triplet and quadruplet interactions among nodes15,16,17,18. In fact, hypernetworks appear as suitable representations of certain dynamical processes found in physics19,20, chemistry21 and neuroscience22,23. This has ignited research aimed at understanding the impact of higher-order interactions on the dynamical behavior of complex systems24,25,26,27. Moreover, besides considering hypernetworks as a good description of such models, we observed that hypernetworks could be revealed in data-driven model reconstructions when the original model is a network. Therefore, a major puzzle is why hypernetworks emerge as the fitting description of actual network data.
Here, we show that hypernetworks can describe experimental data of networks of electrochemical oscillators with nonlinear coupling. We uncover a mechanism that generates higher-order interactions as a model to describe oscillator networks from data. First, we show that sparse model recovery from data reveals higher-order interactions. We then develop a theory for the emergence of such higher-order interactions when the isolated system is close to a Hopf bifurcation. We provide an algorithm to reveal emergent hypernetwork and its emergent coupling functions for any network in disciplines ranging from neuroscience to chemistry. The emergent hypernetworks provide a dimension reduction that allows the characterization of critical transitions.
We designed an experimental system with four oscillatory chemical reactions coupled with nonlinear feedback and delay arranged in a ring network (see Fig. 1a). The set-up consists of a multichannel potentiostat interfaced with a real-time controller and connected to a Pt counter, a Hg/Hg2SO4 sat K2SO4 reference, and four Ni working electrodes in 3.0 M sulfuric acid electrolyte. At a constant circuit potential (V0 = 1100 mV with respect to the reference electrode) and with an external resistance (Rind = 1.0 kohm) attached to each nickel wire, the electrochemical dissolution of nickel exhibits periodic current and electrode potential oscillations with a natural frequency of 0.385 Hz.
a Experimental setup. b Schematic illustration of the electrochemical experiment with the nonlinear feedback. The blue, orange, yellow, and green lines represent the elements 1 to 4, respectively. The electrode potential signals (Ek) of the four (nearly) isolated electrodes are nonlinearly modulated and fed back with a delay τ to the corresponding circuit potential (Vk), which drives the metal dissolution. (The delay is implemented by storing the past data in the memory of the computer.) c Representation of the in a ring network topology used in the experiment. d Electrode potential time series. e Filtered and fitted (dark red line) instantaneous frequency using LASSO for hypernetwork reconstruction corresponding from top to bottom to oscillators 1 to 4, respectively. f Experimental recovery of the phase interactions given by a hypernetwork.
Without coupling, we adjusted the natural frequency of each oscillator to have a ratio with respect to oscillator 1 as ω2/ω1 = 2.53 (≈2.5), ω3/ω1 = 1.56 (≈1.5) and ω4/ω1 = 2.53 (≈2.5) with a set of resistors and capacitors (Cind), see Supplementary Note 1.) The natural frequencies create opportunities for triplet resonances, as there are small detunings for ω1 − ω2 + ω3 and ω1 − ω4 + ω3, as well as pairwise resonances ω2 ≈ ω4.
The individual electrode potentials (Ek) were recorded and rescaled and offset corrected
where ok and Ok are the time-averaged electrode potential and amplitude rescaling factor, respectively. (The rescaling factors, Ok = 0.5, 1, 0.5, 1 were applied to counter the different amplitudes of the slow oscillators.) A ring-coupling can be introduced with external feedback (see Fig. 1b, c) according to
where Vk(t) and V0,k are the applied and the offset circuit potential of the kth electrode, respectively, K is the coupling strength, Ak is the adjacency matrix, τ is a time delay, and
This delayed nonlinear feedback modulates the impact of the coupled units with a bias towards positive values (similar to a diode operation in the (−1, 1) interval). Note that this form of feedback is fundamentally different from previously applied nonlinear schemes4 in that it does not produce obvious synchronization patterns, for example, one and multi-cluster states.
Figure 1d shows the time series of the electrode potential for K = 5.2 and τ = 1.65 s. The slow oscillators (1 and 3) have larger amplitudes and the time series exhibit nonlinear waveform modulations without any obvious synchronization pattern (one-cluster state).
From the potentials ({tilde{E}}_{k}) we extract the frequencies ({dot{theta }}_{k}) and apply a first-order Savitzky-Golay filter with a time window of 45 s to remove the in-cycle and short-range phase fluctuation, as shown in Fig. 1e (solid line). For each oscillator, a slow variation is seen as the oscillators slow down and speed up on a timescale of about 100 s (or 40 cycles); notably, the elements 1 and 3 exhibit similar ({dot{theta }}_{k}) oscillations, which are different from those in elements 2 and 4.
To describe the nature of the phase dynamics, we consider the slow triplet phase differences
which correspond to the triplet frequency detunings.
The impact of triplet interactions on the dynamics can be extracted with a LASSO fit to
where ({hat{omega }}_{k}(t)={hat{omega }}_{k}^{0}+{hat{omega }}_{k}^{1}t+{hat{omega }}_{k}^{2}{t}^{2}) is the fitted, slowly drifting (up to quadratic variation in time) natural frequency, and ({C}_{j}^{k}) and ({D}_{j}^{k}) are the amplitudes of the sin and cos phase coupling functions corresponding to the appropriate triplet phase differences. The strength of the triplet interactions j = 1, 2 (for ϕj) on oscillator k is given by the amplitudes ({H}_{j}^{k}=sqrt{{({C}_{j}^{k})}^{2}+{({D}_{j}^{k})}^{2}}).
The dynamics of oscillators 1 and 3 are impacted by both triplet interactions; ϕ1 impacts oscillators 1 and 3 with amplitudes 4.9 × 10−3 and 4.4 × 10−3, and ϕ2 with 2.3 × 10−3 and 3.2 × 10−3, respectively. However, the dynamics of oscillators 2 and 4 are only impacted by triplet interactions ϕ1 (with amplitude 1.33 × 10−2) and ϕ2 (1.7 × 10−2), respectively. These triplet interactions describe phase fluctuations over the long time scale (red curves in Fig. 1e). Therefore, we can conclude that the phase dynamics of the oscillators coupled in a ring can be described by a hypernetwork shown in Fig. 1f.
The fact that model recovery provides triplets as the best description is rather puzzling. Also given that the resonant behavior ω2 ≈ ω4 did not appear in the model recovery from data. This suggests an interplay between the resonant frequencies and the network topology. The question arises, which resonances/triplet interactions emerge from a large number of possibilities in a given network, natural frequencies, and nonlinear coupling? An outstanding question is what is the origin of these triplet interactions that were generated by pairwise physical coupling?
To answer these questions, we develop a theory that captures the important characteristics of the experiments: nonlinear coupling and triplet resonance conditions. We consider the networks
where ({z}_{k}in {mathbb{C}}) is the state of the kth oscillator, ({h}_{k}:{mathbb{C}}times {mathbb{C}}to {mathbb{C}}) is the pairwise coupling function, Ak is the adjacency matrix, and α > 0 is the coupling strength. When the isolated system is close to a Hopf bifurcation, the dynamics is described by fk(zk) = γkzk + βkzkzk228. The Hopf bifurcation is a common route to oscillations in nonlinear systems and describes the appearance of oscillations in applications2,3,5,6,8. Our proofs are valid for γk = λ + iωk with small λ and ωk satisfying resonance conditions. We fix βk = − 1, but this value is immaterial. We develop a normal form theory to eliminate unnecessary terms of h(zk, z) and to expose higher-order ones that predict the dynamics. To a network of the form of Eq. (6) we associate non-resonance conditions that allow us to get rid of the leading interaction terms in α.
Since h(zk, z) is a linear combination of monomials and the theory can be applied to each monomial independently, we assume first that h(zk,z) is a single monomial of the form
for non-negative numbers d1,…,d4. Our major theoretical result is a formulation of a non-resonance condition given by
This condition shows up naturally in our approach, as a monomial Eq. (7) can only be eliminated by a transformation that divides by the left-hand side of Eq. (8). Hence, an interaction term in the coupling function h given by Eq. (7) can only be removed if the non-resonance condition is satisfied. The non-resonance condition is defined as the union over all non-resonance conditions of its monomial terms. The network non-resonance conditions are given by the union over all non-resonance conditions of h(zk, z) for which Ak ≠ 0. Our result is the following:
In Methods, we show that given Eq. (6) with (h:{mathbb{C}}times {mathbb{C}}to {mathbb{C}}) a smooth map with vanishing constant terms, under the network non-resonance conditions, there is a coordinate transformation that eliminates pairwise interaction terms and reveals the higher-order interactions. The proof consists of two main steps:
(i) Existence of a polynomial change of variables. Consider
for some polynomials Pk. The goal is to design Pk such that in the variables uk interaction terms linear in α vanish. We obtain higher-order interactions of order α2. For Eq. (6) we use
where ({tilde{h}}_{kell }(z,w)) is the function obtained from h(z, w) by transforming each monomial according to the following replacement rule:
Note that the imaginary part of the denominator in Eq. (11) is precisely the left-hand side of Eq. (8). While bringing the equations to the new form, we face a major challenge to understand the combinatorial behavior of the Taylor coefficients during the transformation. We define a bracket on the space of polynomials to track these coefficients.
(ii) Dealing with transformed isolated dynamics. The second major challenge lies in the fact that another coordinate transformation is needed to eliminate terms coming from the isolated dynamics fk. Indeed, as we eliminate coupling terms linear in α, other terms linear in α appear due to the isolated dynamics. A remarkable fact is that the same non-resonance conditions also ensure that the second transformation exists.
Our theorem is applicable to a much broader class of coupling functions and network formalisms than what is described by Eq. (6). A rich variety of new interaction rules can emerge, depending on the specifics of the set-up (see Supplementary Note 2).
Applying the replacement rule Eq. (11) we obtain
up to higher-order terms in α and u. In Methods, we discuss the new coupling functions 1Gk and 2Gk some their properties. The coupling is now α2 explaining anomalous synchronization transitions that appears in networks (see Supplementary Note 3).
Similar to the experiments we consider a ring of four oscillators with coupling function
Instead of delay, the oscillators are coupled through a conjugate variable that enables a streamlined theoretical treatment. Close to a Hopf bifurcation, the delay would have an effect of advancing the oscillations over half a period. As before, we consider ω1 − ω2 + ω3 and ω1 − ω4 + ω3 to be close to zero, so, capturing the triplet resonance in the experiments. We apply our theory to this case to unravel how higher-order interactions appear in the data.
The coupling function is a combination of (zbar{w}) and ({z}^{2}bar{w}), providing d1 = 1 and d4 = 1 for the first monomial and d1 = 2 and d4 = 1 for the latter. The resonance condition Eq. (8) is satisfied for both. Using the replacement rule Eq. (11), we find
Each node equation contains 16 interaction terms as in Eq. (12). We discuss some of these terms for the first node. ({}^{2}{G}_{1}^{23}) appears as node 1 is connected to node 2 and 2 to 3. This interaction is resonant, see Fig. 2a. ({}^{2}{G}_{1}^{43}) appears because node 1 is connected to 4 and node 4 to 3. This term is also resonant, see Fig. 2b. ({}^{1}{G}_{1}^{24}) is nonzero and nonresonant. This term appear as 1 is directed connected to 2 and 4, see Fig. 2c. Finally, the term ({}^{2}{G}_{1}^{24}) is a forbidden, the term would appear from an interaction of 1 to 2 and from 2 to 4, however, in the original network the later interaction is absent, see Fig. 2d. Remarkably, not all interactions are relevant when the goal is to describe slow oscillations in the phases.
Coupling functions appearing in Eq. (12) of node 1. Colors correspond to signs in the phase combination with blue standing for positive and orange for negative. a Resonant interaction term appearing as ({}^{2}{G}_{1}^{23}). b Resonant interaction term appearing as ({}^{2}{G}_{1}^{43}). Finally, c is a nonresonant term and d ({}^{2}{G}_{1}^{24}) is a forbidden term (it does not appear). These new interaction terms can be predicted from the combinatorics of the original network and coupling function.
Indeed, once we analyse the phases in the new equations, the coupling term coming from ({}^{2}{G}_{1}^{23}) will lead to oscillations with frequency close to ω1 − ω2 + ω3 while the term coming from ({}^{2}{G}_{1}^{43}) leads to a frequency close to ω1 − ω4 + ω3. This implies that both terms are slowly varying. In contrast, the term coming from ({}^{2}{G}_{1}^{24}) leads to oscillations with frequency ω1 − ω2 + ω4 ≈ ω1 and is fast oscillating in comparison to the slow terms with small frequencies. In virtue of the averaging theory, such fast oscillating terms can be neglected. In fact, only resonant terms connected by local trees in the original graph will survive such as the resonant ones involving ω1 − ω2 + ω3 and ω1 − ω4 + ω3. This yields
where ({eta }_{pq}=frac{1}{{gamma }_{p}+{bar{gamma }}_{q}}) and ({zeta }_{pqr}=frac{2}{{gamma }_{p}+{bar{gamma }}_{q}}+frac{2}{{gamma }_{p}+{bar{gamma }}_{r}}+frac{1}{bar{{gamma }_{q}}}+frac{1}{bar{{gamma }_{r}}}). Writing u = reiθ we obtain equations for the phases θ. The averaging theorem gives
where the phases ϕ1 and ϕ2 are given in Eq. (4). The functions ρ and σ are provided in the Supplementary Note 4. The emergent hypernetwork explains the experimental fitting found in Eq. (5). These functions represent hyperlinks as shown in Fig. 1f.
The phase triplets ϕ1 and ϕ2 are revealed from phase reduction in the normal form and they are not obvious from the original Eq. (6). We confirm these predictions by direct simulations of Eq. (6) (Supplementary Note 5). We present examples for a three-node path in Supplementary Note 6 and a six-node network in Supplementary Note 7.
In Supplementary Note 3, we show that the experimental recovery of a hypernetwork is not an artifact. Rather, we prove that imposing sparsity unavoidably leads to the recovery of the normal form instead. Indeed, as the recovery allows for a small least square deviation between the data and the model, the recovery finds the hypernetwork as a simpler description of the system. So, by measuring the original variables and attempting a model recovery while imposing sparsity, model recovery learns only the higher-order interactions. We now use the emergent network prediction for the ring network with the corresponding resonance conditions as in the experiment to explain the slow phase dynamics.
From the data we extract the slow phases ϕ1 and ϕ2 as shown in Fig. 3 in solid lines. Using our theory, from Eq. (16), we obtain that
where a’s and b’s are given in terms of the functions σ and ρ in Eq. (16) see Supplementary Note 5. We treat a’s and b’s as fitting parameters from the vector field in Eq. (17) obtained from first principles, since the corresponding coupling parameter and amplitudes are unknown. The resulting solutions agree with the experimental data as seen in Fig. 3. Our findings are not strictly limited to electrochemical oscillators. As shown in Supplementary Note 9, we detected the same hypernetworks in nonlinearly coupled integrate-and-fire neuron models.
We show the time series of the slow phase ϕ1 and ϕ2 from experimental data (solid) and the prediction of the emergent hypernetwork (dashed) capturing higher-order interactions. The vector field describing the phase interaction is obtained from first principles. The coefficients of the vector field are obtained by least-square minimization.
The requirement of a nonlinear coupling, at first sight, seems to be a limitation for practical applications. However, here we analyze how hypernetworks emerge in modular networks with microscopic pairwise coupling through phase differences.
We consider four subpopulations of N interacting Kuramoto oscillators13. Nodes in each subpopulation interact strongly among themselves with coupling strength μ and weakly between subgroups with coupling strength α, see Fig. 4. As we will show at the macroscopic mean-field level, the interaction is nonlinear. According to our theory, although the mean-fields have a pairwise interaction, their model recovery will be in terms of hypernetworks. We first consider the microscopic description; each oscillator is described by
or in terms of mean-fields ({dot{psi }}_{km}={omega }_{km}+{{{{{{{rm{Im}}}}}}}}left(mu {z}_{k}+alpha {sum }_{}{A}_{kl}{z}_{l}right){e}^{-i{psi }_{km}}) where
is the mean-field of the subpopulation k. The frequencies ωkm are distributed according to a Lorenzian ρ(ω, Ωk, σk) where Ωk is the mean subpopulation frequency and σk is the frequency dispersion. Applying the Ott-Antonsen ansatz15, we obtain the macroscopic equations describing the mean-fields in the limit N →  as
where fk is the Hopf normal form with constants γk = (iΩk + μ − σk) and βk = − μ and
thus, in the macroscopic description the coupling is nonlinear. We interpret α as a bifurcation parameter and deal with αzl as a nonlinear term as in bifurcation theory. We consider the ensemble frequencies to satisfy the resonance conditions Ω1 + Ω3 ≈ 2Ω2 and Ω2 + Ω4 ≈ 2Ω1. At α = 0 each subpopulation will have an order parameter behaving as ({z}_{k}(t)={r}_{k}{e}^{i{theta }_{k}(t)}) where ({r}_{k}=sqrt{frac{mu -{sigma }_{k}}{mu }}) and ({dot{theta }}_{k}={{{Omega }}}_{k}). To obtain the phase model, we bring the network to its normal form and apply the phase reduction. In Supplementary Note 10, we perform the calculations of such resonance conditions to obtain the new normal form equations. After discarding nonresonant terms the phase equations of the mean-fields read as
where Fi is a linear combination of sine and cosine.
a The original network of coupled subpopulations (with four distinct colours, namely, red, yellow, blue and orange). Oscillators are interacting by an internal coupling constant μ and inter-subpopulations coupling constant α. b Higher order phase interaction of the mean-fields represented with the same colors as in a (red, yellow, blue and orange). Applying our approach we uncover that the phase interaction between the mean-fields is described by a hypernetwork. c The mean-field slow phase variables φ1 (green) and φ2 (purple) were computed from the data collected from the simulations of mean fields on the associated network. The dashed curve is the simulation of the vector field of the slow phases φ1,2 reconstructed from data using the Lasso method.
Next, we fix the ensemble frequencies as Ω1 = 2, Ω2 = 3, Ω3 = 4 and Ω4 = 1 as well as the coupling strengths μ = 0.5, σk = 0.48 yielding rk = 0.15 and α = 0.1 for all subpopulations. We numerically integrate the mean-field equations and obtain the complex fields z1(t), z2(t), z3(t) and z4(t) which enables us to extract the phase dynamics θ1(t), θ2(t), θ3(t) and θ4(t). Performing a Lasso regression we recover the vector fields of Eq. (22) confirming the theoretical prediction of higher order interactions, see Supplementary Note 10.
As before, we introduce the slow phases
The theory predicts the higher order interaction between the slow phases as ({dot{varphi }}_{k}={varepsilon }_{k}+{G}_{k}({varphi }_{1},{varphi }_{2})), as shown in Supplementary Note 10. The fitting the predicted vector field of φ to the data is excellent as can be observed in Fig. 4c.
For these four subpopulation on a ring, the condition on the frequencies is close to the subspace Vres = {Ω1 + Ω3 = 2Ω2,  Ω2 + Ω4 = 2Ω1} , forming a co-dimension 2 resonance surface. That is, the emergence of hypernetworks is generic in a two parameter family of frequencies.
We have uncovered a mechanism by which nonlinear pairwise interactions with triplet resonance conditions result in nontrivial phase dynamics on a hypernetwork. Such interactions traditionally were attributed in brain dynamics to synaptic transmission between two neurons mediated by chemical messengers from a third neuron (heterosynaptic plasticity)29. Our findings provide an alternative mechanism. On one hand, this finding shows that phase dynamics can be mediated through ‘virtual’ interactions not physically present in the system. On the other hand, such a mechanism could be leveraged to design interactions between remote components not directly connected but instead having correlations in natural frequencies.
The experimental system with a generic network motif with a ring of four electrochemical oscillators presented here was an example, where a relatively simple nonlinear modulation of the coupling induced a hypernetwork driven phase dynamics. Networks with a ring topology are selected for the experiment since they are common for many network based complex systems, e.g., in lasers, biological systems, neuronal dynamics and many disciplines30,31. Such nonlinear modulation of the coupling can be quite general in gene expressions; for example, it was used to describe the coupling among circadian cells through Michaelis-Menten mechanism where coupling from one cell modulated the maximum gene expression rate in the other32.
Strikingly, we showed that the coupling resulting in mean-field coupling among network modules has sufficient nonlinearity to facilitate hypernetwork interactions. In particular, event related modulation of spectral responses of magnetoencephalogram (MEG) recordings (i.e., modulation of frequency-specific oscillations in the motor network established by a handgrip task) have shown very strong evidence for nonlinear, between-frequency coupling of remote brain regions33. Our results strongly suggest that in these MEG recordings, given the appropriate resonances and nonlinearities, hypernetwork description could facilitate the long-range modulation of frequencies. In conclusion, the findings open new avenues for hypernetwork based description and engineering of complex systems with heterogeneous frequencies and nonlinear interactions.
Our results give an algorithmic procedure for obtaining a hypernetwork that accurately describes the observed behavior of the original system. This emergent higher order system depends on details of the given network, the original coupling function and the resonance relations among the phases.
In Supplementary Note 2, we consider ODEs of the general form
with ({z}_{k}in {mathbb{C}}) and (alpha in {mathbb{R}}). The numbers ({beta }_{k},{gamma }_{k}in {mathbb{C}}) are assumed non-zero, and we furthermore write γk = λ + iωk. Here (lambda in {mathbb{R}}) is seen as the bifurcation parameter for a Hopf bifurcation, and we assume the interaction functions ({H}_{k}:{{mathbb{C}}}^{n}to {mathbb{C}}) to be smooth (i.e. C) for convenience. Moreover, we initially assume each Hk satisfies Hk(0) = 0 and DHk(0) = 0, though the condition on its derivative is later dropped.
Our main result shows that the ODE (24) can be put in a normal form that allows us to predict the phase dynamics of the oscillators. We do this by using two successive transformations:
for some appropriately chosen polynomials Pk and Qk. The first of these coordinate transformations is used to remove the term αHk(z) from the Eq. (24). This will generate additional terms in α2 that may be expressed in the coefficients of Hk and Pk following certain combinatorial rules. We manage this combinatorial behavior by introducing a special bracket [••] on the space of polynomials. In addition to these new interaction terms, the transformation will also produce terms in α involving Pk and βkzkzk2, which obscure an interpretation of the system as a (hyper) network. We therefore remove these additional terms using the second coordinate transformation. A crucial observation here is that the non-resonance conditions needed for the first transformation are sufficient to ensure the second. We are able to prove this using the precise bookkeeping enabled by the aforementioned bracket.
When dealing with the case where DHk(0) ≠ 0, we instead remove only the non-linear terms in Hk using the transformations (25) and (26). This reveals higher order terms as before. Even though DHk(0) accounts only for nonresonant terms by assumption, this linear term will nevertheless cause an overall frequency shift that has to be accounted for. More precisely, if we denote by Ω the diagonal matrix with entries the frequencies ω1, …, ωn, then the natural frequencies in the coupled case will be given by the imaginary part of the eigenvalues of iΩ + αDH(0). Here we have set H = (H1, …, Hn). These new frequencies can be approximated by standard eigenvalue perturbation techniques.
Applying the transformation of the theorem to Eq. (6) yields a new system of the form Eq. (12). In Supplementary Note 2, we show that
In Eq. (27) a term of degree d in h and a term of degree (tilde{d}) in ({tilde{h}}_{kell }) combine to form a term of degree (d+tilde{d}-1) in ({}^{1}{G}_{k}^{ell p}). As both h and ({tilde{h}}_{kell }) have terms of degree 2 and higher, we see that ({}^{1}{G}_{k}^{ell p}) only has terms of degree 3 and higher. The same holds true for ({}^{2}{G}_{k}^{ell p}), which means that a classical network description involving directed edges is no longer possible.
The third order terms are moreover easily found by replacing h and ({tilde{h}}_{kell }) in Eq. (27) by their quadratic terms. Likewise, the fourth order terms are found by replacing h by its quadratic terms and ({tilde{h}}_{kell }) by its cubic terms and vice versa in Eq. (27). We may also argue that these higher order terms in ({}^{1}{G}_{k}^{ell p}) and ({}^{2}{G}_{k}^{ell p}) are non-vanishing in general. Indeed, the coefficients in front of these terms are rational functions of γk and the coefficients of h. Such functions are either identical to the zero function (which Eq. (27) excludes) or non-vanishing on an open dense set.
New terms emerge that have an interpretation as higher-order interactions. The two double sums in Eq. (12) have a combinatorial interpretation. The first double sum counts all pairs of nodes (, p) that both influenced node k in the original network. The second double sum counts all pairs (, p) where influenced k and p influenced and p need not influence k directly in the old network, so that new node-dependency is formed.
We present an algorithm for obtaining an emergent hypernetwork from a given network system. Its input consists of the adjacency matrix A, the function h and the phases ω1 through ωn, and we assume the nonresonance conditions of the theorem to hold. The algorithm is as follows:
Emergent Hypernetworks
Input: Adjacency matrix A, coupling function h, frequencies and amplitudes γi’s
Output: Hypernetwork and Coupling functions
1:   for each (kin {{{{{{{mathscr{S}}}}}}}}) do
2:   for each (ell in {{{{{{{mathscr{S}}}}}}}}) do
3:        if Ak ≠ 0 then
4:         form the polynomials ({tilde{h}}_{kell }({u}_{k},{u}_{ell })) by the replacement rule
(hskip8pc{z}^{{d}_{1}}{bar{z}}^{{d}_{2}}{w}^{{d}_{3}}{bar{w}}^{{d}_{4}}mapsto frac{{z}^{{d}_{1}}{bar{z}}^{{d}_{2}}{w}^{{d}_{3}}{bar{w}}^{{d}_{4}}}{({d}_{1}-1){gamma }_{k}+{d}_{2}{bar{gamma }}_{k}+{d}_{3}{gamma }_{ell }+{d}_{4}{bar{gamma }}_{ell }})
5:         for each (pin {{{{{{{mathscr{S}}}}}}}}) do
6:         if AkAkp ≠ 0 then
7:         Compute ({}^{1}{G}_{k}^{ell p})
8:         if AkAp ≠ 0 then
9:         Compute ({}^{2}{G}_{k}^{ell p})
10:   procedure Resonant terms in the coupling functionsG
11:    for each ({u}_{k}^{{d}_{1}}{bar{u}}_{k}^{{d}_{2}}{u}_{ell }^{{d}_{3}}{bar{u}}_{ell }^{{d}_{4}}{u}_{p}^{{d}_{5}}{bar{u}}_{p}^{{d}_{6}}) monomial of ({}^{1}{G}_{k}^{ell p}) and ({}^{2}{G}_{k}^{ell p}) do
12:         if (d1 − d2 − 1)ωk + (d3 − d4)ω + (d5 − d6)ωp ≠ 0 then
13:       discard term
14:   procedure Remaining monomials are the couplings of node k
We provide the experimental time-series and the extracted phases of the oscillations (Fig. 1) at ref. 34Source data are provided with this paper.
The source code for reconstructing the functions representing hypernetwork dynamics from oscillatory networks dynamics is available ref. 35.
Watts, M., Tabak, J., Zimliki, C., Sherman, A. & Bertram, R. Slow variable dominance and phase resetting in phantom bursting. J. Theor. Biol. 276, 218–228 (2011).
Article  ADS  Google Scholar 
Kralemann, B. et al. In vivo cardiac phase response curve elucidates human respiratory heart rate variability. Nat. Commun. 4, 1–9 (2013).
Article  Google Scholar 
Sebek, M., Tönjes, R. & Kiss, I. Z. Complex rotating waves and long transients in a ring network of electrochemical oscillators with sparse random cross-connections. Phys. Rev. Lett. 116, 068701 (2016).
Article  ADS  Google Scholar 
Bick, C., Sebek, M. & Kiss, I. Z. Robust weak chimeras in oscillator networks with delayed linear and quadratic interactions. Phys. Rev. Lett. 119, 168301 (2017).
Article  ADS  Google Scholar 
Schneidman, E., Berry, M. J., Segev, R. & Bialek, W. Weak pairwise correlations imply strongly correlated network states in a neural population. Nature 440, 1007–1012 (2006).
Article  ADS  CAS  Google Scholar 
Ermentrout, G. B. & Terman, D. H. Mathematical foundations of neuroscience, vol. 35 (Springer Science & Business Media, 2010).
Blasius, B., Huppert, A. & Stone, L. Complex dynamics and phase synchronization in spatially extended ecological systems. Nature 399, 354–359 (1999).
Article  ADS  CAS  Google Scholar 
Matheny, M. H. et al. Exotic states in a simple network of nanoelectromechanical oscillators. Science 363, eaav7932 (2019).
Smeal, R. M., Ermentrout, G. B. & White, J. A. Phase-response curves and synchronized neural networks. Philos. Trans. R. Soc. Lond., B, Biol. Sci. 365, 2407–2422 (2010).
Article  Google Scholar 
Omel’chenko, E. & Wolfrum, M. Nonuniversal transitions to synchrony in the sakaguchi-kuramoto model. Phys. Rev. Lett. 109, 164101 (2012).
Article  ADS  Google Scholar 
Hong, H. & Strogatz, S. H. Kuramoto Model of Coupled Oscillators with Positive and Negative Coupling Parameters: An Example of Conformist and Contrarian Oscillators. Phys. Rev. Lett. 106, 054102 (2011).
Article  ADS  Google Scholar 
Kuramoto, Y. Chemical oscillations, waves, and turbulence (Courier Corporation, 2003).
Stankovski, T., Pereira, T., McClintock, P. V. E. & Stefanovska, A. Coupling functions: Universal insights into dynamical interaction mechanisms. Rev. Mod. Phys. 89, 045001 (2017).
Article  ADS  MathSciNet  Google Scholar 
Rodrigues, F. A., Peron, T. K. D., Ji, P. & Kurths, J. The kuramoto model in complex networks. Phys. Rep. 610, 1–98 (2016).
Article  ADS  MathSciNet  Google Scholar 
Tönjes, R., Fiore, C. E. & Pereira, T. Coherence resonance in influencer networks. Nature Communi. 12, 1–8 (2021).
Article  Google Scholar 
Giusti, C., Pastalkova, E., Curto, C. & Itskov, V. Clique topology reveals intrinsic geometric structure in neural correlations. Proc. Natl. Acad. Sci. U.S.A. 112, 13455–13460 (2015).
Article  ADS  MathSciNet  CAS  Google Scholar 
Reimann, M. W. et al. Cliques of neurons bound into cavities provide a missing link between structure and function. Front. Comput. Neurosci. 11, 48 (2017).
Article  Google Scholar 
Bassett, D. S., Zurn, P. & Gold, J. I. On the nature and use of models in network neuroscience. Nat. Rev. Neurosci. 19, 566–578 (2018).
Article  CAS  Google Scholar 
Millán, A. P., Torres, J. J. & Bianconi, G. Explosive higher-order kuramoto dynamics on simplicial complexes. Phys. Rev. Lett. 124, 218301 (2020).
Article  ADS  MathSciNet  Google Scholar 
Bick, C., Ashwin, P. & Rodrigues, A. Chaos in generically coupled phase oscillator networks with nonpairwise interactions. Chaos 26, 094814 (2016).
Article  ADS  MathSciNet  Google Scholar 
Kori, H., Kuramoto, Y., Jain, S., Kiss, I. Z. & Hudson, J. L. Clustering in globally coupled oscillators near a hopf bifurcation: theory and experiments. Phys. Rev. E 89, 062906 (2014).
Article  ADS  Google Scholar 
Giusti, C., Ghrist, R. & Bassett, D. S. Two’s company, three (or more) is a simplex. J. Comput. Neurosci. 41, 1–14 (2016).
Article  MathSciNet  Google Scholar 
Bassett, D. S. & Sporns, O. Network neuroscience. Nat. Neurosci. 20, 353–364 (2017).
Article  CAS  Google Scholar 
Grilli, J., Barabás, G., Michalska-Smith, M. J. & Allesina, S. Higher-order interactions stabilize dynamics in competitive network models. Nature 548, 210–213 (2017).
Article  ADS  CAS  Google Scholar 
Skardal, P. S. & Arenas, A. Abrupt desynchronization and extensive multistability in globally coupled oscillator simplexes. Phys. Rev. Lett. 122, 248301 (2019).
Article  ADS  CAS  Google Scholar 
Mulas, R., Kuehn, C. & Jost, J. Coupled dynamics on hypergraphs: Master stability of steady states and synchronization. Phys. Rev. E 101, 062313 (2020).
Article  ADS  MathSciNet  CAS  Google Scholar 
Bilal, S. & Ramaswamy, R. Synchronization and amplitude death in hypernetworks. Phys. Rev. E 89, 062923 (2014).
Article  ADS  Google Scholar 
Shil’nikov, L., Shil’nikov, A., Turaev, D. & Chua, L. Methods of qualitative theory in nonlinear dynamics, vol. 5 (World Scientific, 2001).
Chistiakova, M., Bannon, N. M., Bazhenov, M. & Volgushev, M. Heterosynaptic plasticity. Neuroscientist 20, 483–498 (2014).
Article  Google Scholar 
Popovych, O. V., Yanchuk, S. & Tass, P. A. Delay-and coupling-induced firing patterns in oscillatory neural loops. Phys. Rev. Lett. 107, 228102 (2011).
Article  ADS  Google Scholar 
Takamatsu, A. et al. Spatiotemporal symmetry in rings of coupled biological oscillators of physarum plasmodial slime mold. Phys. Rev. Lett. 87, 078102 (2001).
Article  ADS  CAS  Google Scholar 
Schroder, S., Herzog, E. D. & Kiss, I. Z. Transcription-Based Oscillator Model for Light-Induced Splitting as Antiphase Circadian Gene Expression in the Suprachiasmatic Nuclei. J. Biol. Rhytms 27, 79–90 (2012).
Article  Google Scholar 
Chen, C. C. et al. Nonlinear Coupling in the Human Motor System. J. Neurosci. 30, 8393–8399 (2010).
Article  CAS  Google Scholar 
Nijholt, E., Ocampo-Espindola, J. L., Eroglu, D., Kiss, I. Z. & Pereira, T. Emergent hypernetworks in weakly coupled oscillators (this paper). GitHub repository which includes the data measured from the experiments. https://github.com/jloespindola/Hypernetwork_data (2021).
Nijholt, E., Ocampo-Espindola, J. L., Eroglu, D., Kiss, I. Z. & Pereira, T. Emergent hypernetworks in weakly coupled oscillators (this paper). Zenodo repository: code for computing the functions representing hypernetwork dynamics from the networks dynamics. https://doi.org/10.5281/zenodo.5749164 (2021).
Download references
We thank Sajjad Bakrani, Zachary G. Nicolaou, Marcel Novaes, Edmilson Roque, Robert Ronge and Jeroen Lamb for enlightening discussions. T.P. was supported in part by FAPESP Cemeai Grant No. 2013/07375-0 and is a Newton Advanced Fellow of the Royal Society NAFR1180236. T.P. and E.N. were partially supported by Serrapilheira Institute (Grant No. Serra-1709-16124). D.E. was supported by TUBITAK Grant No. 118C236 and the BAGEP Award of the Science Academy. JLO-E acknowledges financial support from CONACYT. I.Z.K. acknowledges support from National Science Foundation (grant CHE-1900011).
Instituto de Ciências Matemáticas e Computação, Universidade de São Paulo, São Carlos, Brazil
Eddie Nijholt & Tiago Pereira
Department of Chemistry, Saint Louis University, St. Louis, MO, USA
Jorge Luis Ocampo-Espindola & István Z. Kiss
Faculty of Engineering and Natural Sciences, Kadir Has University, Istanbul, Turkey
Deniz Eroglu
Department of Mathematics, Imperial College London, London, UK
Tiago Pereira
You can also search for this author in PubMed Google Scholar
You can also search for this author in PubMed Google Scholar
You can also search for this author in PubMed Google Scholar
You can also search for this author in PubMed Google Scholar
You can also search for this author in PubMed Google Scholar
E.N. and T.P. designed the overall study and formulated the theory. J.L.O.-E. and I.Z.K. designed and performed the experiments. D.E. implemented the numerical simulations and analyses. All authors contributed to the writing of the manuscript. All authors reviewed and approved the final manuscript.
Correspondence to Tiago Pereira.
The authors declare no competing interests.
Nature Communications thanks the other anonymous reviewer(s) for their contribution to the peer review of this work. Peer review reports are available.
Publisher’s note Springer Nature remains neutral with regard to jurisdictional claims in published maps and institutional affiliations.
Open Access This article is licensed under a Creative Commons Attribution 4.0 International License, which permits use, sharing, adaptation, distribution and reproduction in any medium or format, as long as you give appropriate credit to the original author(s) and the source, provide a link to the Creative Commons license, and indicate if changes were made. The images or other third party material in this article are included in the article’s Creative Commons license, unless indicated otherwise in a credit line to the material. If material is not included in the article’s Creative Commons license and your intended use is not permitted by statutory regulation or exceeds the permitted use, you will need to obtain permission directly from the copyright holder. To view a copy of this license, visit http://creativecommons.org/licenses/by/4.0/.
Reprints and Permissions
Nijholt, E., Ocampo-Espindola, J.L., Eroglu, D. et al. Emergent hypernetworks in weakly coupled oscillators. Nat Commun 13, 4849 (2022). https://doi.org/10.1038/s41467-022-32282-4
Download citation
DOI: https://doi.org/10.1038/s41467-022-32282-4
Anyone you share the following link with will be able to read this content:
Sorry, a shareable link is not currently available for this article.

Provided by the Springer Nature SharedIt content-sharing initiative
Nature Communications (2022)
By submitting a comment you agree to abide by our Terms and Community Guidelines. If you find something abusive or that does not comply with our terms or guidelines please flag it as inappropriate.
Nature Communications (Nat Commun) ISSN 2041-1723 (online)
© 2023 Springer Nature Limited
Sign up for the Nature Briefing newsletter — what matters in science, free to your inbox daily.


Leave a Reply

Your email address will not be published. Required fields are marked *