Memory selection and information switching in oscillator networks with higher-order interactions
aa r X i v : . [ n li n . AO ] S e p Memory selection and information switching in oscillator networks with higher-order interactions
Per Sebastian Skardal ∗ and Alex Arenas Department of Mathematics, Trinity College, Hartford, CT 06106, USA Department d’Enginyeria Inform´atica i Matem´atiques, Universitat Rovira i Virgili, 43007 Tarragona, Spain
We study the dynamics of coupled oscillator networks with higher-order interactions and their ability to storeinformation. In particular, the fixed points of these oscillator systems consist of two clusters of oscillators thatbecome entrained at opposite phases, mapping easily to information more commonly represented by sequencesof 0’s and 1’s. While N such fixed point states exist in a system of N oscillators, we find that a relatively smallfraction of these are stable, as chosen by the network topology. To understand the memory selection of suchoscillator networks, we derive a stability criterion to identify precisely which states are stable, i.e., which piecesof information are supported by the network. We also investigate the process by which the system can switchbetween different stable states when a random perturbation is applied that may force the system into the basinof attraction of another stable state. PACS numbers: 05.45.Xt, 89.75.Hc
I. INTRODUCTION
Collective behavior in ensembles of network-coupled dy-namical units is an important area of research due to a widerange of both natural and engineered applications [1, 2]. Ex-amples include cardiac pacemakers [3], synthetic cell engi-neering [4], and power grid dynamics [5, 6]. Moving beyondsuch a system’s ability to either remain incoherent or syn-chronize into a single entrained group, various combinationsof dynamical and structural properties may give rise to morecomplicated synchronization patterns, for instance when os-cillators form coexisting synchronized and incoherent groupsknown as chimera states [7] or spontaneously partition them-selves in different clusters that remain entrained within but donot synchronize with one another [8]. Of particular interesthere are patterns where phase oscillators form multiple en-trained groups. In the illustrative case of two groups, whichwill be our focus in this work, oscillator systems can be usedas models for information storage, treating oscillators that en-train to one or the other group as a 0 or 1, i.e., a bit, in a pieceof information [9–12]. Given the mappings between phaseand integrate-and-fire oscillators [13], understanding the dy-namics that give rise to such patterns may shed light on cog-nitive function and memory [14–18].In addition to the collective behavior seen in network dy-namics that emerge from classical pair-wise interactions, agreat deal of interest has recently been paid to higher-orderinteractions, i.e., interactions that take place between threeor more units (and are fundamentally different than linearcombinations of pair-wise interactions) [19–21]. A large out-standing question in this direction is: What novel effects dosuch higher-order interactions have on macroscopic behavior?Phase reduction techniques have already shown that such in-teractions take place in generic oscillator systems [24, 25] andempirical data suggests that they may play an important rolein brain dynamics [26–28]. In particular, three-way interac-tions interactions may describe correlations in neuronal spik-ing activity in the brain that provides a missing link between ∗ [email protected] structure and function [29]. Some results have been developedfor treating synchronization of identical (possibly chaotic) os-cillators [22, 23] and higher-order interactions between phaseoscillators have been investigated in a handful of studies [30–34], but we demonstrated in recent work that such interactionsnaturally enable systems to store memory and information viathe kinds of patterns described above [35, 36].In this paper we study the process by which memory is gen-erated in oscillator systems with non-trivial network topolo-gies and higher-order interactions. We find that, comparedto the all-to-all coupled case studied in Refs. [35, 36], whenthere is an underlying network topology that constrains theinteractions the memory capacity of the system, quantified bythe number of stable fixed point states the dynamics supports,is also constrained. Specifically, we identify precisely whichoscillator states, i.e., pieces of information, are supported bya given network topology using a linear stability criterion. Wethen study the process by which the system switches from onepiece of information to another as perturbations are applied tothe different stable states and may move the system into thebasin of attraction of a different state.The remainder of this paper is organized as follows. InSec. II we present the model and discuss the connection be-tween stable states and memory. In Sec. III we present a sta-bility analysis and explore memory capacity in oscillator net-works with higher-order interactions with an illustrative ex-ample. In Sec. IV we investigate the process by which thesystem can move between different stable states under randomperturbations, thereby switching between pieces of informa-tion. In Sec. V we explore the effects of network structureand investigate the asymmetry of stable states. In Sec. VI weconclude with a discussion of our results. II. THE 2-SIMPLEX OSCILLATOR MODEL
We begin by consider an extension of the Kuramoto modelwith higher-order interactions, specifically three-way, or 2-simplex, interactions [37]. The all-to-all coupled case wastreated in detail in Refs. [35, 36], but when placed on a non-trivial network topology the system is governed by the follow-ing equations of motion for N oscillators, ˙ θ i = ω i + K N X j =1 N X l =1 B ijl sin( θ j + θ l − θ i ) , (1)where θ i and ω i are the phase and natural frequency, respec-tively, of oscillator i and K is the global coupling strength.Importantly, the coupling function is the sine of a linear com-bination of three oscillators, rather than two, and interactionsare determined by the adjacency tensor B . We note that thiscoupling function is one of the two possibilities for 2-simplexinteractions that arise from phase reduction. The other pos-sibility, sin(2 θ j − θ l − θ i ) , does not admit stable solutionswith multiple clusters, and therefore does not support mem-ory storage, so we focus on the coupling in Eq. (1). For sim-plicity we assume that the entries of B are unweighted andundirected, so that B ijl = 1 if a three-way interaction existsbetween oscillators i , j , and l (and other wise B ijl = 0 ) and B ijl = B ilj = B jil = B jli = B lij = B lji . This last propertyessentially states that interactions are not directed and, withthe factor of in the denominator of Eq. (1), each three-wayinteractions is counted exactly once. In principle, the higher-order structure encoded in B may be defined in different waysdepending on the application. In other words, the entries of B may or may not depend on the entries of an underlying net-work consisting of nodes and links (i.e., 0- and 1-simplexes)encoded in an adjacency matrix A . However, here we considerthe case where higher-order interactions are in fact dictatedby the adjacency matrix A (which we assume is unweightedand undirected) so that B ijl = A ij A jl A li , i.e., a higher-orderinteraction exists between a group of three nodes if they areconnected in a triangle.The dynamics of Eq. (1) naturally give rise to synchronizedstate patterns that entrain oscillators in two groups at oppositesides of the torus. To see this we introduce the set of localorder parameters z i = r i e iψ i = 12 N X j =1 N X l =1 B ijl e iθ j e iθ l , (2)which allows us to rewrite Eq. (1) as ˙ θ i = ω i + Kr i sin( ψ i − θ i ) . (3)First, by entering the rotating reference frame θ θ + ωt ,where ω = N P Ni =1 ω i is the mean natural frequency, we mayset the mean frequency to zero. Next, assuming that the dy-namics relax to a stationary state, oscillator i becomes phase-locked if | ω i | ≤ Kr i (we will return to this condition later)and converges to one of two stable fixed points defined by θ ∗ i = φ i or φ i + π. (4) where φ ∗ i = ψ i + arcsin (cid:16) ω i Kr i (cid:17) . (5)A suitable shift initial conditions allows us to center the meanphases ψ i about zero, and under typical conditions we ex-pect that phase-locking oscillators will have ψ i ≈ , which,combined with Eq. (4), yields two entrained groups of phase-locked oscillators: one about θ = 0 and the other about θ = π .While this demonstrates that phase-locked oscillators tendto fall into one of two possible groups, the presence of driftingoscillators poses an issue in terms of interpretation. To thisend, we consider the dynamics of Eq. (1) in the regime wherethe coupling strength is sufficiently strong compared to thespread of frequencies, so that K ≫ | ω i | for all i = 1 , . . . , N .By considering the rescaled time τ = Kt and approximating ω i /K ≈ , Eq. (1) simplifies to ˙ θ i = 12 N X j =1 N X l =1 B ijl sin( θ j + θ l − θ i ) , (6)in which case all oscillators converge to either or π .Before we move on to the analysis of these states, it’sworthwhile to discuss their representation of information. Inparticular, a fixed point of a system of N oscillators will begiven by some θ ∗ ∈ R N , where each θ ∗ i = 0 or π . Thisis easily interpretable as a string of bits: for instance, in asystem of N = 5 oscillators we may converge to the state θ ∗ = [0 , , π, , π ] T , which one can be mapped to the se-quence of bits (0 , , , , . This manner of memory in oscil-lator systems is not dissimilar to that observed in neural mod-els consisting of Ising spin particles, best exemplified in theseminal 1988 papers by Gardner [38] and Gardner and Der-rida [39]. In these models information bits, i.e., 0’s and 1’s,are equivalent not to oscillator phases, but Ising spins: plusesand minuses. The oscillator systems studied here can then beinterpreted as a different model for capturing and studying asimilar phenomenon. III. STABILITY ANALYSIS
To understand the nature of fixed point states of the dynam-ics given by Eq. (6), we begin by analyzing their asymptoticstability. Recall that we are interested in fixed points θ ∗ ∈ R N with entries θ ∗ i = 0 or π . First, N such states exist. However,the dynamics are invariant under rotations, including rotationby π , i.e., θ θ + π , so without any loss of generality we“anchor” θ ∗ = 0 , leaving N − possible fixed point solutions.Then, given a fixed point θ ∗ , its stability is governed by theeigenvalues of the Jacobian matrix DF ( θ ∗ ) whose entries aregiven by -10 -8 -6 -4 -2 0 2 Re( ) -3-2-10123 I m () (a)(b) (a)(b) (c)(d)stableunstable FIG. 1.
Example: Stable and unstable states. (a), (b) Two possible fixed points of the dynamics that are stable and unstable, respectively, ona network of size N = 8 . Oscillator values θ = 0 and π are colored blue and yellow, respectively. (The top right oscillator is constrained tobe θ ∗ = 0 .) (c) Eigenvalue spectrum of the Jacobian DF ( θ ∗ ) for states depicted in (a) and (b), plotted as circles and crosses, respectively.(d) The time series of the norm of perturbations δθ ( t ) = θ ( t ) − θ ∗ , initially set to k δθ (0) k = 10 − , for all N − = 128 possible states.Perturbations to states depicted in (a) and (b) are plotted as thick solid and dashed curves, respectively. Perturbations to all other possible statesare plotted as thin curves, solid or dashed if they are predicted to be stable or unstable, respectively. DF ij ( θ ∗ ) = ( − P Nk =1 P Nl =1 B ikl cos( θ ∗ k + θ ∗ l − θ ∗ i ) if i = j, P Nl =1 B ijl cos( θ ∗ j + θ ∗ l − θ ∗ i ) if i = j. (7)Due to the rotational invariance of the dynamics, DF ( θ ∗ ) isguaranteed to have one trivial eigenvalue λ = 0 , which canalso be seen by noting that the rows of DF ( θ ∗ ) all sum tozero, so v ∝ is an eigenvector with eigenvalue λ = 0 .Therefore, if all other eigenvalue have strictly negative realpart, i.e., are contained in the left-half complex plane, then θ ∗ is stable, and otherwise it is unstable. This gives us themachinery to evaluate the stability of each of the N − fixedpoints described above.Before a more in-depth exploration of the stability of pos-sible states we consider a small, illustrative example. Specif-ically, consider the network and two states depicted in Fig. 1(a) and (b), where blue and yellow filled nodes correspond tostates θ i = 0 and π , respectively. (Here we index the nodes i = 1 , , . . . , N = 8 in the clockwise direction, starting fromthe top, so we constrain θ = 0 .) For the given network struc-ture (where all nodes are part of at least one connected triangleand therefore are part of at least one higher-order interaction)the state depicted in (a) is stable, while the state in (b) is un-stable. To see this we plot in Fig. 1(c) the eigenvalue spectraof the Jacobian matrices DF ( θ ∗ ) for the two states (a) and (b)using circles and crosses, respectively. Not that in with excep-tion to the trivial eigenvalue λ = 0 for each, the full spectrumfor (a) is contained in the left-half complex plane, whereas four eigenvalues for (b) are located in the right-half complexplane. To test this stability criterion we directly simulate theevolution of a small perturbation δθ ( t ) = θ ( t ) − θ ∗ initiallyset at k δθ (0) k = 10 − and plot the results in Fig. 1(d). Timeseries for perturbations to states (a) and (b) are plotted as thicksolid and dashed curves, respectively, showing that the pertur-bations exponentially decay and grow until saturation, respec-tive. We also plot the time series for perturbations to all otherstates as solid or dashed curves based on the prediction fromthe stability criterion, getting agreement.Continuing with the example above, only a relatively smallfraction (10 out of a possible N − = 128 ) of states are in factstable. In Fig. 2 we plot these stable states. In terms of stringsof bits (where θ = 0 and π correspond to and ) these statesare precisely (0 , , , , , , , , (0 , , , , , , , , (0 , , , , , , , , (0 , , , , , , , , (0 , , , , , , , , (0 , , , , , , , , (8) (0 , , , , , , , , (0 , , , , , , , , (0 , , , , , , , , (0 , , , , , , , Inspecting these ten states more closely, we find two important (a)( f) (b) (c) (d) (e)(g) (h) (i) (j)
FIG. 2.
Stable states.
For the example network depicted in Fig. 1, all stable states (out of a possible N − = 128 possible fixed points). properties. First, a number of the stable states appear to becombinations of one another in the sense that oscillators with θ ∗ i = π in different stable states often combined to create an-other stable state. For example, the state depicted in Fig. 2(d)is a combination of those in Figs. 2(b) and (c). This does notalways happen, however, as states depicted in Figs. 2(c) and(e) do not combine to form a new stable state. More specifi-cally, it also points to the role of node i = 8 (the location of θ i = π in the state at the top of the right column) which hasno effect on the stability of a given state. To see this, note thatany stable state with θ = 0 is still stable when θ is changedto , and vice versa. It is worth noting that node i = 8 is themost poorly connected node in the network is the sense that ishas only one three-way interaction.The second important property that this example illustratesis a certain asymmetry in the collection of stable states. Notefirst, that of the ten total stable states, one has zero θ ∗ = π entries, four have one π entry, four have two π entry, and onehas three π entries. However, of all N − = 128 possiblefixed points, one has zero π entries, seven have one π entry,21 have two π entries, 35 have three π entries, 35 have four π entries, and so on. In particular, the majority of all fixed pointstates have a roughly even distribution of and π entries, butwhen we restrict our attention to only stable states we find astrong asymmetry as they tend to have an uneven distributionof and π entries. IV. SWITCHING
Lastly, with the stability analysis above as a tool for iden-tifying stable states of the system, i.e., supported pieces ofinformation, we investigate the process by which the systemmay switch between different pieces of memory, i.e., stablestates. What we describe here is not quite a homoclinic net-work [40, 41] but is in many senses similar. Sticking with ourexample network illustrated above, none of the 10 stable statesdepicted in Fig. 2 are connected with a heteroclinic orbit since (a)(e) (c)(b)(d)(h) (i) (j)(f)(g)
FIG. 3.
Switching.
Illustration of the transition network betweenstates (a)–(j) illustrated in Fig. 2. Each directed link is drawn withwidth proportional to the rate at which the system switches from onestate to another under a random perturbation of size k δθ k = 1 . .(Self-links, describing rates at which perturbations do not cause atransition, are not illustrated.) they are all asymptotically stable. Instead, we again considerperturbations δθ to these stable states θ ∗ , but now we allowperturbations to not be arbitrarily small, letting k δθ k = O (1) ,so that the perturbed state θ ∗ + δθ may fall outside of the basinof attraction of θ ∗ and eventually converge to another stablestate.More precisely, we consider the following setup. For eachof the 10 stable states (a)–(j) depicted in Fig. 2, we introduce random perturbations of size k δθ k = 1 . . Then, foreach of these perturbations to the states j we find the frac-tion that eventually converge to each other state i , and letthis fraction populate the entry π ij of the transition matrix Π .Thus, each entry π ij describes the rate at which perturbationsof this size cause a switch in states j → i . This transition st. dev. k (2)
Effect of network structure.
Illustration of the transitionnetwork between states (a)–(j) illustrated in Fig. 2. Each directedlink is drawn with width proportional to the rate at which the systemswitches from one state to another under a random perturbation ofsize k δθ k = 1 . . (Self-links, describing rates at which perturba-tions do not cause a transition, are not illustrated.) network is by nature directed and weighted, and is illustratedin Fig. 3, letting more strongly (weakly) weighted links bethicker (thinner). We have neglected to illustrate self-links,i.e., entries π ii that describe rates at which the perturbationdoes not result in a switch in state. It should be noted that thechoice k δθ k = 1 . is important; choosing perturbations toosmall eventually leads to no transitions from one state to an-other, while choosing perturbations too large eventually leadsto much denser networks whose links reflect only the size ofthe basins of attraction of each state. The choice used here liesin between these two extremes. In this example, for instance,we see that a few state pairs have particularly strong rates oftransitions, i.e., (d) → (c) and (j) → (i). This stems from the sim-ilarity of the states themselves (see Fig. 2 where it can be seenthat these states differ only by a single oscillator). Switch-ing between such similar states appears to exist between otherpairs as well, e.g., (b) → (a), (f) → (e), and (h) → (g), but thestrength of the transition is variable. Moreover, the uniformstate (a) (where θ ∗ = ) is a sink where the system cannotescape from. This property is not particularly surprising forreasonable perturbations since state (a) is the most linearlystable in the sense that the largest nontrivial eigenvalue is themost negative compared to those for other states. V. NETWORK STRUCTURE AND STATE ASYMMETRY
Before concluding we investigate two aspects of the sta-bility properties of the system. We begin by probing the ef-fect of network structure on the overall stability of differentstates. Keeping in mind the combinatorial complexity of thenumber of possible states to probe (recall that for a networkof N oscillators we have N − states to examine) we main-tain a relatively small network size, specifically N = 12 , butstudy an ensemble of soft spherical networks built us-ing the hyperbolic graph generator [42]. In particular, we usea small but finite temperature of T = 1 guaranteeing that,while some long-range links exist, the networks are signif-icantly clustered. (We throw out any networks that are notconnected by its triangles.) Moreover, while we use a target asymmetry, ( ) s t ab l e f r a c t i on FIG. 5.
Asymmetry of stable states.
Illustration of the transitionnetwork between states (a)–(j) illustrated in Fig. 2. Each directedlink is drawn with width proportional to the rate at which the systemswitches from one state to another under a random perturbation ofsize k δθ k = 1 . . (Self-links, describing rates at which perturba-tions do not cause a transition, are not illustrated.) mean 1-simplex degree of h k (1) i = 4 we note that the mean 2-simplex degree, h k (2) i , where k (2) i = P Nj,l =1 A ij A jl A li , aswell as its standard deviation, vary significantly from networkto network.We begin by investigating the overall stability of states by,for each network, calculating the fraction of all states that arestable. In Fig. 4(a) and (b) we plot the fraction of stable statesfor each network vs the corresponding mean and standard de-viation of the 2-simplex degree. It’s clear to see that the frac-tion of stable states that a network admits has a positive cor-relation with both the mean and standard deviation of the 2-simplex degree.We also investigate more closely the nature of the states thatend up being stable. Recall that from the example depicted inFig. 2 we observed a noticeable asymmetry. We now examinethis by defining an asymmetry identifier for each state θ ∗ as η ( θ ∗ ) = P Ni =1 χ θ i =0 − P Ni =1 χ θ i = π N , (9)where χ θ ∗ i =0 = 1 if θ ∗ i = 0 and χ θ ∗ i =0 = 0 otherwise (andsimilarly, χ θ ∗ i = π = 1 if θ ∗ i = π and χ θ ∗ i = π = 0 otherwise).In particular, η ( θ ∗ ) = 1 or in the cases that all θ ∗ i sharethe same value or are even split between and π . Using thisasymmetry identifier to organize the entire ensemble of pos-sible states, we then evaluate, for all states θ ∗ that share anasymmetry identifier value η ( θ ∗ ) the fraction of states thatare stable. In Fig. 5 we plot, as a function of the asymmetryidentifier η ( θ ∗ ) the fraction of stable states for states sharingthe value η ( θ ∗ ) averaged over the ensemble of networks usedabove. (Error bars indicate one standard deviation above andbelow the mean.) We observe that the trend observed in ourexample depicted in Fig. 2 persists, i.e., states that are moreasymmetric are more likely to be stable. VI. DISCUSSION
In this paper we have studied the process by which oscilla-tor networks with higher-order interactions can select memoryand information, often represented as sequences of bits, bystabilizing fixed points consisting of two groups of clusteredoscillators at opposite phases. While oscillator networks havelong been used as models for encoding and storing such infor-mation [16], such systems have often required fine-tuning andengineering of, e.g., individual coupling strengths, to attain agiven state. However, the discovery that higher-order interac-tions naturally leads to a a wide variety of such states providesa more natural formalism for information storage without such fine-tuning [35, 36].Here we have presented a stability analysis for determin-ing precisely which of the N potential states are stable. Thisselection depends on the underlying network topology that en-codes the higher-order interactions between oscillators. More-over, we have also used our stability analysis to study the pro-cess by which the oscillator network can switch between dif-ferent stable states, as a perturbation can be applied to forcethe system into a basin of attraction for another stable state. [1] S. H. Strogatz, Sync: the Emerging Science of Spontaneous Or-der (Hypernion, 2003).[2] A. Pikovsky, M. Rosenblum, and J. Kurths,
Synchronization: AUniversal Concept in Nonlinear Sciences (Cambridge Univer-sity Press, 2003).[3] L. Glass and M. C. Mackey,
From Clocks to Chaos: TheRhythms of Life (Princeton University Press, Princeton, 1988).[4] A. Prindle, P. Samayoa, I. Razinkov, T. Danino, L. S. Tsim-ring, and J. Hasty, A sensing array of radically coupled genetic’biopixels’, Nature , 39 (2012).[5] M. Rohen, A. Sorge, M. Timme, and D. Witthaut, Self-organized synchronization in decentralized power grids, Phys.Rev. Lett. , 064101 (2012).[6] P. S. Skardal and A. Arenas, Control of coupled oscillator net-works with application to microgrid technologies, Sci. Adv. ,e1500339 (2015).[7] M. J. Panaggio and D. M. Abrams, Chimera states; coexistenceof coherence and incoherence in networks of coupled oscilla-tors, Nonlinearity , R67 (2015).[8] L. M. Pecora, F. Sorrentino, A. M. Hagerstrom, T. E. Murphy,and R. Roy, Cluster synchronization and isolated desynchro-nization in complex networks with symmetries, Nat. Commun. , 4079 (2014).[9] H. Daido, Multibranch entrainment and scaling in large popu-lations of coupled oscillators, Phys. Rev. Lett. , 1406 (1996).[10] P. Ashwin, G. Orosz, J. Wordsworth, and S. Townley, Dynamicson networks of cluster states for globally coupled phase oscil-lators, SIAM J. Applied Dynamical Systems , 728 (2007).[11] P. S. Skardal, E. Ott, and J. G. Restrepo, Cluster synchronyin systems of coupled phase oscillators with higher-order cou-pling, Phys. Rev. E , 036208 (2011).[12] M. Komarov and A. Pikovsky, Multiplicity of singular syn-chronous states in the Kuramoto model of coupled oscillators,Phys. Rev. Lett. , 204101 (2013).[13] A. Politi, and M. Rosenblum, Equivalence of phase-oscillatorand integrate-and-fire models. Phys. Rev. E, , 4, 042916(2015).[14] A. Arenas and C. J. P´erez Vicente, Phase locking in a networkof neural oscillators, Europhys. Lett. , 79 (1994).[15] C. J. P´erez Vicente, A. Arenas, and L. L. Bonilla, On the short-term dynamics of networks of Hebbian coupled oscillators, J.Phys. A: Math. Gen. , L9 (1996).[16] F. C. Hoppensteadt and E. M. Izhikevich, Oscillatory neuro-computers with dynamic connectivity, Phys. Rev. Lett. , 2983(1999).[17] J. Fell and N. Axmacher, The role of phase synchronization inmemory processes, Nat. Rev. Neurosci. , 105 (2011).[18] J. F. Hipp, D. J. Hawallek, M. Corbetta, M. Siegel, and A. K.Engel, Large-scale cortical correlation structure of spontaneousoscillatory activity, Nat. Neurosci. , 884 (2012). [19] V. Salnikov, D. Cassese, and R. Lambiotte, Simplicial com-plexes and complex systems, Eur. J. Phys. , 014001 (2019).[20] F. Battiston, G. Cencetti, I. Iacopini, V. Latora, M. Lucas, A.Patania, J.-G. Yong, and G. Petri, Networks beyond pairwiseinteractions: Structure and dynamics, Phys. Rep. , 1 (2020).[21] T. Carletti, D. Fanelli, and S. Nicoletti, Dynamical systems onhypergraphs, J. Phys. Complex. , 035006 (2020).[22] R. Mulas, C. Kuehn, and J. Jost, Coupled dynamics on hyper-graphs: Master stability of steady states and sycnhronization,Phys. Rev. E , 062313 (2020).[23] L. V. Gambuzza, F. Di Patti, L. Gallo, S. Lepri, M. Romance,R. Criado, M. Frasca, V. Latora, and S. Boccaletti, The MasterStability Function for synchronization in simplicial complexes,arXiv:2004.03913.[24] P. Ashwin and Ana Rodrigues, Hopf normal form with S N sym-metry and reduction to systems of nonlinearly coupled phaseoscillators, Physica D , 14 (2016).[25] I. Le´on and D. Paz´o, Phase reduction beyond the first order:The case of the mean-field complex Ginzburg-Landau equation,Phys. Rev. E , 012211 (2019).[26] G. Petri, P. Expert, F. Turkheimer, R. Carhart-Harris, D. Nutt,P.J. Hellyer, and F. Vaccarino, Homological scaffolds of brainfunctional networks, J. R. Soc. Interface , 20140873 (2014).[27] C, Giusti, R. Ghrist, and D. S. Bassett, Two’s company, three(or more) is a simplex, J. Comput. Neurosci. , 1 (2016).[28] A. E. Sizemore, C. Giusti, A. Kahn, J. M. Vettel, R. Betzel, andD. S. Bassett, Cliques and cavities in the human connectome, J.Comput. Neurosci. , 115 (2018).[29] M.W. Reimann, M. Nolte, M. Scolamiero, K. Turner, R. Perin,G. Chindemi, P. Dlotko, R. Levi, K. Hess, and H. Markram,Cliques of neurons bound into cavities provide a missing linkbetween structure and function, Frontiers in Comp. Neuro., ,48 (2017).[30] T. Tanaka and T. Aoyagi, Multistable attractors in a network ofphase oscillators with three-body interactions, Phys. Rev. Lett. , 224101 (2011).[31] M. Komarov and A. Pikovsky, Finite-size-induced transitionsto synchrony in oscillator ensembles with nonlinear global cou-pling, Phys. Rev. E , 020901(R) (2015).[32] C. Bick, P. Ashwin, and Ana Rodrigues, Chaos in genericallycoupled phase oscillator networks with nonpairwise interac-tions, Chaos , 094814 (2016).[33] A. P. Mill´an, J. J. Torres, and G. Bianconi, Explosive higher-order Kuramoto dynamics on simplicial complexes, Phys. Rev.Lett. , 218301 (2020).[34] M. Lucas, G. Cencetti, and F. Battiston, A multi-orderLaplacian for synchronization in higher-order networks,arXiv:2003.09734.[35] P. S. Skardal and A. Arenas, Abrupt desynchronization and ex-tensive multistability in globally coupled oscillator simplexes, Phys. Rev. Lett. , 248301 (2019).[36] C. Xu, X. Wang, and P. S. Skardal, Bifurcation analysis andstructural stability of simplicial oscillator populations, Phys.Rev. Res. , 023281 (2020).[37] Y. Kuramoto, Chemical Oscillations, Waves, and Turbulence (Springer, New York, 1984).[38] E. Gardner, The space of interactions in neural network models,J. Phys. A: Math. Gen. , 257 (1988).[39] E. Gardner and B. Derrida, Optimal storage properties of neural network models, J. Phys. A: Math. Gen. , 271 (1988).[40] P. Ashwin and C. Postlethwaite, On designing heteroclinic net-works from graphs, Phys. D , 26 (2013).[41] C. Bick, Heteroclinic switching between chimeras, Phys. Rev.E , 050201(R) (2018).[42] R. Aldecoa, C. Orsini, and D. Krioukov, Hyperbolic graph gen-erator, Cump. Phys. Commun.196