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.

Advertisement

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

1

26

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 biology^{1,2} and chemistry^{3,4} to neuroscience^{5,6} via ecology^{7} to engineering^{8}. Such networks serve as stepping stones to understand collective dynamics^{9,10,11,12} and other emergent phenomena in networks^{13,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 nodes^{15,16,17,18}. In fact, hypernetworks appear as suitable representations of certain dynamical processes found in physics^{19,20}, chemistry^{21} and neuroscience^{22,23}. This has ignited research aimed at understanding the impact of higher-order interactions on the dynamical behavior of complex systems^{24,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/Hg_{2}SO_{4} sat K_{2}SO_{4} reference, and four Ni working electrodes in 3.0 M sulfuric acid electrolyte. At a constant circuit potential (*V*_{0} = 1100 mV with respect to the reference electrode) and with an external resistance (*R*_{ind} = 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 (*E*_{k}) of the four (nearly) isolated electrodes are nonlinearly modulated and fed back with a delay *τ* to the corresponding circuit potential (*V*_{k}), 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 (*C*_{ind}), 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 (*E*_{k}) were recorded and rescaled and offset corrected

where *o*_{k} and *O*_{k} are the time-averaged electrode potential and amplitude rescaling factor, respectively. (The rescaling factors, *O*_{k} = 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 *V*_{k}(*t*) and *V*_{0,k} are the applied and the offset circuit potential of the *k*th electrode, respectively, *K* is the coupling strength, *A*_{kℓ} 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 schemes^{4} 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 *k*th oscillator, ({h}_{k}:{mathbb{C}}times {mathbb{C}}to {mathbb{C}}) is the pairwise coupling function, *A*_{kℓ} 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 *f*_{k}(*z*_{k}) = *γ*_{k}*z*_{k} + *β*_{k}*z*_{k}∣*z*_{k}∣^{2}^{28}. The Hopf bifurcation is a common route to oscillations in nonlinear systems and describes the appearance of oscillations in applications^{2,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*(*z*_{k}, *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*(*z*_{k}, *z*_{ℓ}) is a linear combination of monomials and the theory can be applied to each monomial independently, we assume first that *h*(*z*_{k},*z*_{ℓ}) is a single monomial of the form

for non-negative numbers *d*_{1},…,*d*_{4}. 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*(*z*_{k}, *z*_{ℓ}) for which *A*_{kℓ} ≠ 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 *P*_{k}. The goal is to design *P*_{k} such that in the variables *u*_{k} 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 *f*_{k}. 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 ^{1}*G*_{k} and ^{2}*G*_{k} 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 *d*_{1} = 1 and *d*_{4} = 1 for the first monomial and *d*_{1} = 2 and *d*_{4} = 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* = *r**e*^{iθ} 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 oscillators^{13}. 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 ansatz^{15}, we obtain the macroscopic equations describing the mean-fields in the limit *N* → *∞* as

where *f*_{k} 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 *α**z*_{l} 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 *F*_{i} 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 *r*_{k} = 0.15 and *α* = 0.1 for all subpopulations. We numerically integrate the mean-field equations and obtain the complex fields *z*_{1}(*t*), *z*_{2}(*t*), *z*_{3}(*t*) and *z*_{4}(*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 *V*_{res} = {Ω_{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 disciplines^{30,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 other^{32}.

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 regions^{33}. 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 *H*_{k} satisfies *H*_{k}(0) = 0 and *D**H*_{k}(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 *P*_{k} and *Q*_{k}. The first of these coordinate transformations is used to remove the term *α**H*_{k}(*z*) from the Eq. (24). This will generate additional terms in *α*^{2} that may be expressed in the coefficients of *H*_{k} and *P*_{k} 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 *P*_{k} and *β*_{k}*z*_{k}∣*z*_{k}∣^{2}, 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 *D**H*_{k}(0) ≠ 0, we instead remove only the non-linear terms in *H*_{k} using the transformations (25) and (26). This reveals higher order terms as before. Even though *D**H*_{k}(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*Ω + *α**D**H*(0). Here we have set *H* = (*H*_{1}, …, *H*_{n}). 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** *A*_{kℓ} ≠ 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** *A*_{kℓ}*A*_{kp} ≠ 0 **then**

7: Compute ({}^{1}{G}_{k}^{ell p})

8: **if** *A*_{kℓ}*A*_{ℓp} ≠ 0 **then**

9: Compute ({}^{2}{G}_{k}^{ell p})

10: **procedure** Resonant terms in the coupling functions *G*

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** (*d*_{1} − *d*_{2} − 1)*ω*_{k} + (*d*_{3} − *d*_{4})*ω*_{ℓ} + (*d*_{5} − *d*_{6})*ω*_{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. 34. Source 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

Received:

Accepted:

Published:

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.

Focus

Advertisement

© 2023 Springer Nature Limited

Sign up for the *Nature Briefing* newsletter — what matters in science, free to your inbox daily.