Higher-order simplicial synchronization of coupled topological signals
Reza Ghorbanchian, Juan G. Restrepo, Joaquín J. Torres, Ginestra Bianconi
HHigher-order simplicial synchronization of coupled topological signals
Reza Ghorbanchian, Juan G. Restrepo, Joaqu´ın J. Torres, and Ginestra Bianconi
1, 4 School of Mathematical Sciences, Queen Mary University of London, London, E1 4NS, United Kingdom Department of Applied Mathematics, University of Colorado at Boulder, Boulder, CO 80309, USA Departamento de Electromagnetismo y F´ısica de la Materia and Instituto Carlos I de F´ısica Te´orica y Computacional,Universidad de Granada, 18071, Granada, Spain The Alan Turing Institute, The British Library, London, United Kingdom
Simplicial complexes capture the underlying network topology and geometry of complex systemsranging from the brain to social networks. Representing complex systems as simplicial complexesopens the possibility to use the powerful tools of algebraic topology to characterize their structure.Here we show that algebraic topology can also be a fundamental tool to capture the higher-orderdynamics of simplicial complexes. In particular we consider topological signals, i.e., dynamicalsignals defined on simplices of different dimension, here taken to be nodes and links for simplicity.We use algebraic topology for characterizing the coupled dynamics of topological signals defined onnodes and links. This dynamical process is ideally suited for simplicial complexes but might alsobe used to treat topological dynamics on networks with exclusively pairwise interactions. We showthat this coupled synchronization leads to explosive topological synchronization in which phasesdefined on nodes synchronize simultaneously to phases defined on links at a discontinuous phasetransition. The model is here studied on real connectomes and on simplicial network models. Finally,we provide a comprehensive theoretical approach that captures this transition on random networkstreated within the annealed approximation.
I. INTRODUCTION
Higher-order networks [1–4] are attracting increasingattention as they are able to capture the many-body in-teractions of complex systems ranging from brain to so-cial networks. Simplicial complexes are higher-order net-works that encode the network geometry and topologyof real datasets. Using simplicial complexes allows thenetwork scientist to formulate new mathematical frame-works for mining data [5–10] and for understanding thesegeneralized network structures revealing the underlyingdeep physical mechanisms for emergent geometry [11–15]and for higher-order dynamics [16–33]. In particular, thisvery vibrant research activity is relevant in brain researchto analyse real brain data and its profound relation tobrain dynamics [1, 6, 15, 34–37] and in biology, e.g., forthe study of biological transport networks [10, 38].In networks, dynamical processes are typically definedover signals associated to the nodes of the network. Inparticular, the Kuramoto model [39–43] investigates thesynchronization transitions of phases associated to thenodes of the network. This scenario can change signif-icantly in the case of simplicial complexes [16, 17, 19].In fact, simplicial complexes can sustain dynamical sig-nals defined on simplices of different dimension, includ-ing nodes, links, triangles and so on, called topologicalsignals . For instance, topological signals defined on linkscan represent fluxes of interest in brain research and inbiological transportation networks as well. The intereston topological signals is rapidly growing with new resultsrelated to signal processing of topological signals [17, 19]and higher-order topological synchronization [16, 28]. Inparticular higher-order topological synchronization [16]demonstrates that topological signals (phases) associatedto higher dimensional simplices can undergo a synchro- nization phase transition. Interestingly this synchroniza-tion transition can go unnoticed if the dynamics is notprocessed with the suitable topological transformations.These results open new uncharted territory for the inves-tigation of higher-order synchronization.Here we build on the mathematical framework ofhigher-order topological synchronization proposed inRef. [16] and consider a synchronization model in whichtopological signals of different dimension are coupled.We focus in particular on the coupled synchronizationof topological signals defined on nodes and links, but wenote that the model can be easily extended to topologi-cal signals of higher dimension. The reason why we focuson topological signals defined on nodes and links is three-fold. First of all we can have a better physical intuition oftopological signals defined on nodes (traditionally stud-ied by the Kuramoto model) and links (like fluxes) thatis relevant in brain dynamics [44] and biological trans-port networks [10, 38]. Secondly, although the coupledsynchronization dynamics of nodes and links can be con-sidered as a special case of coupled synchronization dy-namics of higher-order topological signals on a genericsimplicial complex, this dynamics can be observed alsoon networks including only pairwise interactions. Indeednodes and links are the simplices that remain unchangedif we reduce a simplicial complex to its network skele-ton. Since currently there is more availability of networkdata than simplicial complex data, this fact implies thatthe coupled dynamics studied in this work has wide ap-plicability as it can be tested on any network data andnetwork model. Thirdly, defining the coupled dynam-ics of topological signals defined on nodes and links canopen new perspectives in exploiting the properties of theline graph of a given network which is the network whosenodes corresponds to the links or the original network a r X i v : . [ c ond - m a t . d i s - nn ] N ov FIG. 1:
Schematic representation of the Kuramoto and the topological Kuramoto model.
Panel (a) shows a networkin which nodes sustain a dynamical variable (a phase) whose synchronization is captured by the Kuramoto model. Panel (b)shows a simplicial complex in which not only nodes but also links sustain dynamical variables whose coupled synchronizationdynamics is captured by the higher-order topological Kuramoto model. [45].In this work, we show that by adopting a global adap-tive coupling of dynamics inspired by Refs. [46, 47] thecoupled synchronization dynamics of topological signalsdefined on nodes and links is explosive [48], i.e., it oc-curs at a discontinuous phase transition in which the twotopological signals of different dimension synchronize atthe same time. We also illustrate numerical evidenceof this discontinuity on real connectomes and on simpli-cial complex models including the configuration model ofsimplicial complexes [49] and the non-equilibrium simpli-cial complex model called Network Geometry with Fla-vor [12, 13]. We provide a comprehensive theory of thisphenomenon formulated within the annealed network ap-proximation of random networks offering a good analyt-ical understanding of the observed transition. The an-alytical results give solid ground to the claim that theinvestigated transition is discontinuous.
II. SIMPLICIAL COMPLEXES AND HIGHERORDER LAPLACIANSA. Definition of simplicial complexes
Simplicial complexes represent higher-order networkswhose interactions include two or more nodes. Thesemany-body interactions are captured by simplices . A n -dimensional simplex α is a set of n + 1 nodes α = [ i , i , . . . , i n ] . (1)For instance a node is a 0-dimensional simplex, a linkis a 1-dimensional simplex, a triangle is a 2-dimensionalsimplex, a tetrahedron is a 3-dimensional simplex, andso on. A face of a simplex is the simplex formed by aproper subset of the nodes of the original simplex. Forinstance the faces of a tetrahedron are 4 nodes, 6 linksand 4 triangles. A simplicial complex is a set of simplicesclosed under the inclusion of the faces of each simplex.Any simplicial complex can be reduced to its simplicialcomplex skeleton , which is the network formed by thesimplicial complex nodes and links. Simplices have a rel-evant topological and geometrical interpretation and con-stitute the topological structures studied by discrete al-gebraic topology. Therefore representing the many-body FIG. 2:
The boundary operators and their representation in terms of the incidence matrices.
Panel (a) and (b)describe the action of the boundary operator on an oriented link and on an oriented triangle respectively. Panel (c) shows atoy example of a simplicial complex and panel (d) indicates its incidence matrices B [1] and B [2] representing the boundaryoperators ∂ and ∂ respectively. interactions of a complex system with a simplicial com-plex opens the very fertile opportunity to use the toolsof algebraic topology [5, 50] to study the topology of thesystem under investigation. In this work we show thatalgebraic topology can also shed significant light on therole that topology has on higher-order synchronization. B. Oriented simplices and boundary map
In algebraic topology simplices are oriented. For in-stance a link α = [ i, j ] has the opposite sign of the link[ j, i ], i.e., [ i, j ] = − [ j, i ] . (2)Similarly to higher order simplices we associate an orien-tation such that[ i , i , . . . , i n ] = ( − σ ( π ) [ i π (0) , i π (1) , . . . , i π ( n ) ] , (3)where σ ( π ) indicates the parity of the permutation π .It is good practice to use as orientation of the simplicesthe orientation induced by the labelling of the nodes, i.e., giving, for example, a positive orientation to any simplex[ i , i , . . . , i n ] , (4)where i < i < i . . . < i n . (5)This will ensure that the spectral properties of the higher-order Laplacians that will be defined later are indepen-dent of the labelling of the nodes. Given a simplicial com-plex, a n -chain consists of the elements of a free abeliangroup C n with basis formed by the set of all oriented n -simplices. Therefore every element of C n can be uniquelyexpressed as a linear combination of the basis elements( n -simplices) with coefficients in Z . The boundary op-erator ∂ n is a linear operator ∂ n : C n → C n − whoseaction is determined by the action on each n -simplex ofthe simplicial complex given by ∂ n [ i , i . . . , i n ] = n (cid:88) p =0 ( − p [ i , i , . . . , i p − , i p +1 , . . . , i n ] . (6)As a concrete example, in Figure 2 we demonstrate theaction of the boundary operator on links and triangles.A celebrated property of the boundary operator is thatthe boundary of a boundary is null , i.e. ∂ n ∂ n +1 = 0 (7)for any n >
0. This relation can be directly proven byusing Eq. (6). Let us consider a simplicial complex K .Let us indicate with N [ n ] the number of simplices of thesimplicial complex with generic dimension n . Given a ba-sis for the linear space of n -chains C n and for the linearspace of ( n − C n − formed by an ordered listof the n simplices and ( n −
1) simplices of the simplicialcomplex, the boundary operator ∂ n can be representedas N [ n − × N [ n ] incidence matrix B [ n ] . In Figure 2 weshow a 2-dimensional simplicial complex and its corre-sponding incidence matrices B [1] and B [2] . Given thatthe boundary matrices obey Eq. (7) it follows that theincidence matrices obey B [ n ] B [ n +1] = , B (cid:62) [ n +1] B (cid:62) [ n ] = , (8)for any n > C. Higher order Laplacians
Using the incidence matrices it is natural to generalizethe definition of the graph Laplacian L [0] = B [1] B (cid:62) [1] (9)to the higher-order Laplacian L [ n ] (also called combina-torial Laplacians) [17, 19, 51] that can be represented asa N [ n ] × N [ n ] matrix given by L [ n ] = L down [ n ] + L up [ n ] (10)with L down [ n ] = B (cid:62) [ n ] B [ n ] , L up [ n ] = B [ n +1] B (cid:62) [ n +1] , (11)for n >
0. The higher-order Laplacian can be proven tobe independent of the orientation of the simplices as longas the simplicial complex has an orientation induced bya labelling of the nodes.The most celebrated property of higher-order Lapla-cian is that the degeneracy of the zero eigenvalue of the n Laplacian L [ n ] is equal to the Betti number β n and thattheir corresponding eigenvectors localize around the cor-responding n -dimensional cavities of the simplicial com-plex. The higher-order Laplacians can be used to definehigher-order diffusion [17] and can display a higher-orderspectral dimension on network geometries. Here we areparticularly interested in the use of incidence matricesand higher-order Laplacians to define higher-order topo-logical synchronization. III. HIGHER-ORDER TOPOLOGICALSYNCHRONIZATIONA. Higher-order topological Kuramoto model oftopological signals of a given dimension
The Kuramoto dynamics describes the synchonizationtransition for phases θ = ( θ , θ , . . . θ N [0] ) (12)associated to nodes, i.e., simplices of dimension n = 0.The Kuramoto model is typically defined on a networkbut it can treat also synchronization of the phases asso-ciated to the nodes of a simplicial complex. Each node i has associated an internal frequency ω i drawn froma given distribution, for instance a normal distribution ω i ∼ N (Ω , θ = ω − σ B [1] sin (cid:16) B (cid:62) [1] θ (cid:17) , (13)where here and in the following we use the notation sin( x )to indicate the column vector where the sine function istaken element wise. Note that here we have chosen towrite this system of equations in terms of the incidencematrix B [1] . However if we indicate with a the adjacencymatrix of the network and with a ij its matrix elements,this system of equations is equivalent to˙ θ i = ω i + σ N (cid:88) j =1 a ij sin( θ j − θ i ) , (14)valid for every node i of the network. For coupling con-stant σ = σ c the Kuramoto model [39–41] displays acontinuous phase transition above which the order pa-rameter R = 1 N [0] (cid:12)(cid:12)(cid:12)(cid:12)(cid:12)(cid:12) N [0] (cid:88) i =1 e i θ i (cid:12)(cid:12)(cid:12)(cid:12)(cid:12)(cid:12) (15)is non-zero also in the limit N [0] → ∞ . The higher-ordertopological Kuramoto model [16] describes synchroniza-tion of phases associated to the n dimensional simplicesof a simplicial complex. Although the definition of themodel applies directly to any value of n , here we considerspecifically the case in which the higher-order Kuramotomodel is defined on topological signals (phases) associ-ated to the links φ = ( φ (cid:96) , φ (cid:96) , . . . φ (cid:96) N [1] ) . (16)In particular here φ (cid:96) r indicates the phase associated tothe r -th link of the simplicial complex. The higher orderKuramoto dynamics defined on simplices of dimension n > ω theinternal frequencies associated to the links of the sim-plicial complex, sampled for example from a normal dis-tribution, ˜ ω (cid:96) ∼ N (Ω , φ = ˜ ω − σ B (cid:62) [1] sin( B [1] φ ) − σ B [2] sin( B (cid:62) [2] φ ) . (17)Once the synchronization dynamics is defined on higher-order topological signals of dimension n (here taken to be n = 1) an important question is whether this dynamicscan be projected on ( n + 1) and ( n −
1) simplices. In-terestingly, algebraic topology provides a clear solutionto this question. Indeed for n = 1, when the dynamicsdescribes the evolution of phases associated to the links,one can consider the projection φ [ − ] and φ [+] respectivelyon nodes and on triangles defined as φ [ − ] = B [1] φ , φ [+] = B (cid:62) [2] φ . (18)Note that in this case B [1] acts as a discrete divergenceand B (cid:62) [2] acts as a discrete curl. Interestingly, since theincidence matrices satisfy Eq. (8) these two projectedphases follow the uncoupled dynamics˙ φ [ − ] = B [1] ˜ ω − σ L [0] sin φ [ − ] , ˙ φ [+] = B (cid:62) [2] ˜ ω − σ L down [2] sin φ [+] . (19)These two projected dynamics undergo a continuous syn-chronization transition at σ c = 0 [16] with order param-eters R down = 1 N [0] (cid:12)(cid:12)(cid:12)(cid:12)(cid:12)(cid:12) N [0] (cid:88) i =1 e i φ [ − ] i (cid:12)(cid:12)(cid:12)(cid:12)(cid:12)(cid:12) ,R up = 1 N [2] (cid:12)(cid:12)(cid:12)(cid:12)(cid:12)(cid:12) N [2] (cid:88) i =1 e i φ [+] i (cid:12)(cid:12)(cid:12)(cid:12)(cid:12)(cid:12) . (20)In Ref. [16] an adaptive coupling between these two dy-namics is considered formulating the explosive higher-order topological Kuramoto model in which the topolog-ical signal follows the set of coupled equations˙ φ = ˜ ω − σR up B (cid:62) [1] sin( B [1] φ ) − σR down B [2] sin( B (cid:62) [2] φ ) . (21)The projected dynamics on nodes and triangles are nowcoupled by the modulation of the coupling constant σ with the order parameters R down and R up , i.e. the twoprojected phases follow the coupled dynamics˙ φ [ − ] = B [1] ˜ ω − σR up L [0] sin φ [ − ] , ˙ φ [+] = B (cid:62) [2] ˜ ω − σR down L down [2] sin φ [+] . (22) This explosive higher-order topological Kuramoto modelhas been shown in Ref. [16] to lead to discontinuoussynchronization model on different models of simplicialcomplexes and on clique complexes or real connectomesas well. B. Higher-order topological Kuramoto model ofcoupled topological signals of different dimension
In this work we want to discuss the case in whichtopological signals of different dimensions are coupledtogether leading to an explosive synchronization transi-tion. Specifically we focus on the coupling of the tradi-tional Kuramoto model (Eq.(13)) to higher-order topo-logical Kuramoto model defined for phases associated tothe links (Eq.(17)). The coupling between these two dy-namics is here performed considering the modulation ofthe coupling constant σ with the global order parametersof the node dynamics (defined in Eq. (15)) and the linkdynamics (defined in Eq. (20)). Specifically, we considertwo models denoted as Model 1 and Model 2. Model 1couples the dynamics of the phases of the nodes θ and ofthe links φ according to the following dynamical equa-tions˙ θ = ω − σR down B [1] sin( B (cid:62) [1] θ ) , (23)˙ φ = ˜ ω − σR B (cid:62) [1] sin( B [1] φ ) − σ B [2] sin( B (cid:62) [2] φ ) . (24)The projected dynamics for φ [ − ] and φ [+] then obey˙ φ [ − ] = B [1] ˜ ω − σR L [0] sin φ [ − ] , (25)˙ φ [+] = B (cid:62) [2] ˜ ω − σ L down [2] sin φ [+] . (26)Therefore the projection on the nodes φ [ − ] of the phases φ associated to the links (Eq. (25)) is coupled to the dy-namics of the phases θ (Eq. (23)) associated directly tonodes. However the projection on the triangles φ [+] ofthe phases φ associated to the links is independent of φ [ − ] and of θ as well. Model 2 also describes the cou-pled dynamics of topological signals defined on nodes andlinks but the adaptive coupling captured by the model isdifferent. In this case the dynamical equations are takento be ˙ θ = ω − σR down B [1] sin( B (cid:62) [1] θ ) , (27)˙ φ = ˜ ω − σR R up B (cid:62) [1] sin( B [1] φ ) − σR down B [2] sin( B (cid:62) [2] φ ) . (28)For Model 2 the projected dynamics for φ [ − ] and for φ [+] obey ˙ φ [ − ] = B [1] ˜ ω − σR R up L [0] sin φ [ − ] , (29)˙ φ [+] = B (cid:62) [2] ˜ ω − σR down L down [2] sin φ [+] . (30)Therefore, as in Model 1 the dynamics of the projection φ [ − ] of the phases φ associated to the links (Eq. (29))is coupled to the dynamics of the phases θ associateddirectly to nodes (Eq. (27)) and vice versa. Moreoverthe dynamics of the projection of the phases φ on thetriangles φ [+] (Eq. (30)) is now also coupled with thedynamics of φ [ − ] (Eq. (29)) and vice versa. Here and inthe following we will use the convenient notation (usingthe parameter m ) to indicate both models 1 and 2 withthe same set of dynamical equations given by˙ θ = ω − σR down B [1] sin( B (cid:62) [1] θ ) , (31)˙ φ = ˜ ω − σR ( R up ) m − B (cid:62) [1] sin( B [1] φ ) − σ (cid:0) R down (cid:1) m − B [2] sin( B (cid:62) [2] φ ) , (32)which reduce to Eqs. (24) for m = 1 and to Eqs. (28) for m = 2.We make two very relevant observations: • First, the proposed coupling between topologicalsignals of different dimension can be easily ex-tended to signals defined on higher-order simplicesproviding a very general scenario for coupled dy-namical processes on simplicial complexes. • Secondly, the considered coupled dynamics of topo-logical signals defined on nodes and links can bealso studied on networks with exclusively pairwiseinteractions where we assume that the number ofsimplices of dimension n > R and R down at a non-zero coupling constant σ = σ c , however R up follows an independent transition at zero coupling(see Figure 3 panels (a) and (c)). In Model 2, on thecontrary all order parameters R , R down and R up dis-play a discontinuous transition occurring for the samenon zero value of the coupling constant σ = σ c (see Fig-ure 3 panels (b) and (d)). This is a direct consequenceof the fact that in Model 1 the adaptive coupling leadingto discontinuous phase transition only couples the phases φ [ − ] and θ , while for Model 2 the coupling involves alsothe phases φ [+] .Additionally we studied both Model 1 and Model 2 ontwo real connectomes: the human connectome of Ref. [52]and the c. elegans connectome from Ref. [53] (see Figure4). Interestingly also for these real datasets we observethat in Model 1 the explosive synchronization involvesonly the phases θ and φ [ − ] while in Model 2 we ob-serve that also φ [+] undergoes an explosive synchroniza-tion transition at the same value of the coupling constant σ = σ c . Note that the connectome of homo sapiens dis-plays a very non-trivial dynamics and for observing thephase transition it is necessary to let the system equili-brate by considering very long time series. IV. ANALYTIC TREATMENT OF THECOUPLED DYNAMICS OF NODES AND LINKSA. Theoretical solution of the model
As mentioned earlier the higher-order topological Ku-ramoto model coupling the topological signals defined onnodes and on links can be defined on simplicial com-plexes and on networks as well. In this section we ex-ploit this property of the dynamics to provide an analyt-ical understanding of the synchronization transition onuncorrelated random network within the annealed net-work approximation. Using the Ott-Antonsen method[43] combined with the annealed approximation [42] wecharacterize the steady state of the coupled higher-orderKuramoto dynamics and compare the analytical resultsto simulations of the original model and of the annealedequations.
B. Annealed dynamics
The first step for providing an analytical understand-ing of the proposed coupled dynamics of nodes and linksis to determine the equations that capture the dynamicsin the annealed approximation. For the dynamics of thephases θ associated to the nodes - Eq. (31)- it is possi-ble to proceed as in the traditional Kuramoto model [42].However the annealed approximation for the dynamics ofthe phases φ defined in Eq. (32) needs to be discussed indetail as it is not directly reducible to previous results.For addressing this problem our aim is to directly definethe annealed approximation for the dynamics of the pro-jected variables φ [ − ] which, here and in the following areindicated as ψ = φ [ − ] , (33)in order to simplify the notation. Moreover we will indi-cate with N = N [0] the number of nodes in the networkor in the simplicial complex skeleton. Here we focus onthe model defined on networks, i.e. we assume that thereare no simplices of dimension two. Consequently R up ac-quires the trivial value R up = 1 and Model 2 reduces toModel 1 that we consider in detail here. We notice thatEq. (25), valid for Model 1 can be written as˙ ψ = B [1] ˜ ω − σR L [0] sin( ψ ) . (34)This equation can be also written elementwise as˙ ψ i = ˆ ω i + σR N (cid:88) j =1 a ij [sin( ψ j ) − sin( ψ i )] , (35) (a) R R R (b) R R R (c) R R R (d) R R R FIG. 3:
The Higher-order topological synchronization models (Models 1 and 2) coupling nodes and links onsimplicial complexes.
The order parameters R , R down and R up are plotted versus σ for the higher-order topologicalsynchronization Model 1 (panels (a) and (c)) and Model 2 (panels (b) and (d)) defined over the Network Geometry with Flavor[13] (panels (a) and (c)) and the configuration model of simplicial complexes [49] (panels (b) and (d)). The Network Geometrywith Flavor on which we run the numerical results shown in (a) and (b) includes N [0] = 500 nodes and has flavor s = − d = 3. The configuration model of simplicial complexes on which we run the numerical results shown in (c) and (d) includes N [0] = 500 nodes and has generalized degree distribution which is power-law with exponent γ = 2 . . In both Model 1 and inModel 2 we have set Ω = Ω = 2. where the vector ˆ ω is given byˆ ω = B [1] ˜ ω . (36)Let us now consider in detail these frequencies in thecase in which the generic internal frequency ˜ ω (cid:96) of a linkfollows a Gaussian distribution, specifically in the case inwhich ˜ ω (cid:96) ∼ N (Ω ,
1) for every link (cid:96) . Using the definitionof the boundary operator on a link it is easy to show thatthe expectation of ˆ ω i is given by (cid:104) ˆ ω i (cid:105) = (cid:88) ji a ij Ω . (37)Given that each node has degree k i , the covariancematrix C is given by the graph Laplacian L [0] of thenetwork, i.e. C ij = (cid:104) ˆ ω i ˆ ω j (cid:105) c = (cid:88) (cid:96),(cid:96) (cid:48) (cid:10) [ B [1] ˜ ω ] i [ B [1] ˜ ω ] j (cid:11) c = [ L [0] ] ij = k i δ ij − a ij , (38) where we have indicated with (cid:104) . . . (cid:105) c the connected corre-lation. Therefore the variance of ˆ ω in the MF approxi-mation is (cid:10) ˆ ω i (cid:11) c = (cid:104) ˆ ω i (cid:105) − (cid:104) ˆ ω i (cid:105) = k i . (39)Moreover, the projected frequencies are actually corre-lated and for i (cid:54) = j we have (cid:104) ˆ ω i ˆ ω j (cid:105) c = (cid:104) ˆ ω i ˆ ω j (cid:105) − (cid:104) ˆ ω i (cid:105) (cid:104) ˆ ω j (cid:105) = − a ij . (40)It follows that the frequencies ˆ ω are correlated Gaussianvariables with average given by Eq. (37) and correla-tion matrix given by the graph Laplacian. Here and inthe following we will denote by G ( ˆ ω ) the distribution ofthese frequencies. The fact that the frequencies ˆ ω i arecorrelated is an important feature of the dynamics of ψ and, with few exceptions (e.g., [54]), this feature has re-mained relatively unexplored in the case of the standardKuramoto model. Additionally let us note that the av-erage of ˆ ω over all the nodes of the network is zero. In (c) R R R (d) R R R (a) R R R (b) R R R FIG. 4:
The Higher-order topological synchronization models (Models 1 and 2) coupling nodes and links onreal connectomes.
The order parameters R , R down and R up are plotted versus σ on real connectomes. Panel (a) and (b)show the numerical results on the human connectome [52] for Model 1 and Model 2 respectively. Panel (c) and (d) show thenumerical results on the c. elegans connectome [53] for Model 1 and Model 2 respectively. In both Model 1 and in Model 2 wehave set Ω = Ω = 2. fact N (cid:88) i =1 ˆ ω i = T ˆ ω = T B [1] ω = 0 , (41)where with we indicate the N -dimensional column vec-tor of elements 1 i = 1. By using the symmetry of theadjacency matrix, i.e. the fact that a ij = a ji , Eq. (41)implies that the sum of ˙ ψ i over all the nodes of the net-work is zero, i.e. N (cid:88) i =1 ˙ ψ i = (cid:80) Ni =1 ˆ ω i + σR (cid:80) i,j a ij [sin( ψ j ) − sin( ψ i )] = 0 . We now consider the annealed approximation consist-ing in substituting the adjacency matrix element a ij withits expectation in an uncorrelated network ensemble a ij → k i k j (cid:104) k (cid:105) N , (42)where k i indicates the degree of node i and (cid:104) k (cid:105) is theaverage degree of the network. In the annealed approxi-mation we can put (cid:104) ˆ ω i (cid:105) (cid:39) k i Ω − (cid:88) j>i k j (cid:104) k (cid:105) N . (43) Also, in the annealed approximation the dynamicalEq. (31) and Eq. (34) reduce to˙ θ = ω − σR down ˆ R k · sin( θ − Θ) , (44)˙ ψ = ˆ ω + σR ˆ R down k sin Ψ − σR k (cid:12) sin ψ , (45)where (cid:12) indicates the Hadamard product (element by el-ement multiplication) and where two additional complexorder parameters are defined asˆ R e i Θ = N (cid:88) i =1 k i (cid:104) k (cid:105) N e i θ i , ˆ R down e i Ψ = N (cid:88) i =1 k i (cid:104) k (cid:105) N e i ψ i . (46) C. Solution of the dynamical equations in theannealed approximation
1. General framework for obtaining the solution of theannealed dynamical equations
In this section we will provide the analytic solutionsfor the order parameter of the higher-order topologicalsynchronization studied within the annealed approxima-tion, i.e., captured by Eqs. (44) and (45). In particularfirst we will find an expression of the order parameters R of the dynamics associated to the nodes (Eq. (44))and subsequently in the next paragraph we will derivethe expression for the order parameter R down associatedto the projection on the nodes of the topological signaldefined on the links (Eq. 45)). By combining the tworesults it is finally possible to uncover the discontinuousnature of the transition.
2. Dynamics of the phases of the nodes
When we investigate Eq. (44) we notice that thisequation can be easily reduced to the equation for thestandard Kuramoto model treated within the annealedapproximation [42] if one performs a rescaling of thecoupling constant σ R → σ . Therefore we can treatthis model similarly to the known treatment of the stan-dard Kuramoto model [40–42]. Specifically, starting fromEq. (44) and using a rescaling of the phases θ accordingto θ i → θ i − Ω t, (47)it is possible to show that we can set Θ = 0 and thereforeEq. (44) reduces to the well-known annealed expressionfor the standard order Kuramoto model given by˙ θ = ω − Ω − σR down ˆ R k · sin( θ ) . (48)Assuming that the system of equations reaches asteady state in which both R down and ˆ R become timeindependent, the order parameters of this system of equa-tions can be found to obey [40, 42, 48, 55]ˆ R = N (cid:88) i =1 k i (cid:104) k (cid:105) N (cid:90) | ˆ c i | < dωg ( ω ) (cid:118)(cid:117)(cid:117)(cid:116) − (cid:32) ω − Ω σk i ˆ R R down (cid:33) ,R = 1 N N (cid:88) i =1 (cid:90) | ˆ c i | < dωg ( ω ) (cid:118)(cid:117)(cid:117)(cid:116) − (cid:32) ω − Ω σk i ˆ R R down (cid:33) , (49)where c i indicates c i = ω − Ω σk j ˆ R R down . (50)and g ( ω ) is the Gaussian distribution with expectationΩ and standard deviation 1.
3. Dynamics of the phases of the links projected on thenodes
In this paragraph we will derive the expression of theorder parameters R down and ˆ R down which, together withEqs. (49), will provide the annealed solution of our model. To start with we assume that the frequencies ˆ ω areknown. In this case we can express the order parame-ters R down and ˆ R down as a function of the probabilitydensity function ρ ( i ) ( ψ, t | ˆ ω ) that node i is associated toa projected phase of the link equal to ψ . Since in the an-nealed approximation ψ i has a dynamical evolution dic-tated by Eq. (45) the probability density function obeysthe continuity equation ∂ t ρ ( i ) ( ψ, t | ˆ ω ) + ∂ ψ (cid:104) ρ ( i ) ( ψ, t | ˆ ω )) v i (cid:105) = 0 (51)with associated velocity v i given by v i = κ i − σR k i sin ψ (52)where we have defined κ i as κ i = ˆ ω i + σk i R ˆ R down sin Ψ . (53)In this case the complex order parameters are given byˆ R down e i Ψ = N (cid:88) i =1 k i (cid:104) k (cid:105) N (cid:90) dψρ ( i ) ( ψ, t | ˆ ω ) e i ψ R down e i ˜Ψ = N (cid:88) i =1 N (cid:90) dψρ ( i ) ( ψ, t | ˆ ω ) e i ψ (54)In order to solve the continuity equation we follow Ott-Antonsen [43] and we express ρ ( i ) ( ψ, t | ˆ ω ) in the Fourierbasis as ρ ( i ) ( ψ, t | ˆ ω ) = 12 π (cid:40) ∞ (cid:88) m =1 ˆ f ( i ) m (ˆ ω i , t ) e i mψ + c.c. (cid:41) . (55)Making the ansatzˆ f ( i ) m (ˆ ω i , t ) = [ b i (ˆ ω i , t )] m (56)we can derive the equation for the evolution of b i = b i (ˆ ω i , t ) given by ∂ t b i + i b i κ i + σk i R
12 ( b i −
1) = 0 . (57)Since we showed before that the average value of ˙ ψ i over nodes is zero, we look for non-rotating stationarysolutions of Eq. (57), ∂ t b i = 0, given by b i = − i d i ± (cid:113) − d i , (58)where d i is given by d i = ˆ ω i σk i R + ˆ R down sin Ψ . (59)0By plugging this expression into Eq. (54) we get the ex-pression of the order parameters given the projected fre-quencies ˆ ω ,ˆ R down cos Ψ = N (cid:88) i =1 k i (cid:104) k (cid:105) N (cid:113) − d i θ (1 − d i )ˆ R down sin Ψ = N (cid:88) i =1 k i (cid:104) k (cid:105) N (cid:26)(cid:113) d i − χ ( d i ) + d i (cid:27) R down = 1 N (cid:12)(cid:12)(cid:12)(cid:12)(cid:12) N (cid:88) i =1 d i − sign( d i + 1) (cid:113) d i − (cid:12)(cid:12)(cid:12)(cid:12)(cid:12) (60) where, indicating by θ ( x ) the Heaviside function, we havedefined χ ( d i ) = [ − θ ( d i −
1) + θ ( − − d i )] . (61)Finally, if the projected frequencies ˆ ω are not known wecan average the result over the frequency distribution G ( ˆ ω ) gettingˆ R down cos Ψ = N (cid:88) i =1 k i (cid:104) k (cid:105) N (cid:90) | d i |≤ d ˆ ω G ( ˆ ω ) (cid:115) − (cid:18) ˆ ω i σR k i + ˆ R down sin Ψ (cid:19) , ˆ R down sin Ψ = − N (cid:88) j =0 k i (cid:104) k (cid:105) N (cid:90) d i > d ˆ ω G ( ˆ ω ) (cid:115)(cid:18) ˆ ω i σR k i + ˆ R down sin Ψ (cid:19) − N (cid:88) i =1 k i (cid:104) k (cid:105) N (cid:90) d i < − d ˆ ω G ( ˆ ω ) (cid:115)(cid:18) ˆ ω i σR k i + ˆ R down sin Ψ (cid:19) − N (cid:88) i =1 k i (cid:104) k (cid:105) N (cid:90) ∞−∞ d ˆ ω G ( ˆ ω ) (cid:18) ˆ ω i σR k i + ˆ R down sin Ψ (cid:19) , (62)and R down = 1 N (cid:12)(cid:12)(cid:12)(cid:12)(cid:12)(cid:12) N (cid:88) i =1 (cid:90) ∞−∞ d ˆ ω G ( ˆ ω ) ˆ ω i σR k i + ˆ R down sin Ψ − sign( d i + 1) (cid:115)(cid:18) ˆ ω i σR k i + ˆ R down sin Ψ (cid:19) − (cid:12)(cid:12)(cid:12)(cid:12)(cid:12)(cid:12) . (63)By solving these equations together with Eqs. (49) wecan capture the critical behavior of the higher-order Ku-ramoto model coupling topological signals defined onnodes and links within the annealed approximation.In Figure 5 we compared the direct numerical integra-tion of the annealed Eqs. (44) and (45) with the obtainedtheoretical predictions finding very good agreement. Thisderivation therefore shows evidence of the discontinuousnature of the observed dynamics. However we have foundthat the annealed approximation of the considered modelcaptures the critical point σ c of the transition but is oth-erwise not a very controllable approximation of the orig-inal higher-order topological dynamics at least for themoderately small network sizes that can be considerednumerically. V. CONCLUSIONS
Until recently the synchronization phenomenon hasbeen explored only in the context of topological sig-nals associated to the nodes of a network. However thegrowing interest in simplicial complexes opens the per-spective of investigating synchronization of higher ordertopological signals, for instance associated to the links ofthe discrete networked structure. Here we uncover howtopological signals associated to nodes and links can becoupled to one another giving rise to an explosive syn-chronization phenomenon involving both signals at thesame time. The model has been tested on real connec-tomes and on major examples of simplicial complexes(the configuration model [49] of simplicial complex andthe Network Geometry with Flavor[13]). Moreover we1 R AAR
AAR APR AP FIG. 5:
The discontinuous phase transition obtained in the annealed approximation
The order parameters R (circles) and R down (squares) are plotted as a function of the coupling constant σ for the annealed approximation (AA)dynamics. The solid lines indicate the analytical predictions (AP) obtained in this approximation. Simulations (AA) are hereobtained by integrating numerically Eqs. (44) and (45) for a Poisson network of N = 1000 nodes and average degree c = 2 . provide an analytical solution of this model within theannealed approximation that provides a theoretical un-derstanding of the mechanism driving the emergence ofthis discontinuous phase transition. This work can beextended in different directions including the study ofthe de-synchronization dynamics of this coupled higher-order synchronization and the duality of this model withthe same model defined on the line graph of the samenetwork. Acknowledgements
This work is partially founded by SUPERSTRIPESOnlus. This research utilized Queen Mary’s Apoc- rita HPC facility, supported by QMUL Research-IT.http://doi.org/10.5281/zenodo.438045. J.J.T. acknowl-edges financial support from the Spanish Ministry of Sci-ence and Technology, and the Agencia Espa˜nola de In-vestigaci´on (AEI) under grant FIS2017-84256-P (FEDERfunds). [1] C. Giusti, R. Ghrist, and D. S. Bassett, Journal of com-putational neuroscience , 1 (2016).[2] F. Battiston, G. Cencetti, I. Iacopini, V. Latora, M. Lu-cas, A. Patania, J.-G. Young, and G. Petri, Physics Re-ports (2020). [3] L. Torres, A. S. Blevins, D. S. Bassett, and T. Eliassi-Rad, arXiv preprint arXiv:2006.02870 (2020).[4] V. Salnikov, D. Cassese, and R. Lambiotte, EuropeanJournal of Physics , 014001 (2018).[5] N. Otter, M. A. Porter, U. Tillmann, P. Grindrod, and H. A. Harrington, EPJ Data Science , 17 (2017).[6] G. Petri, M. Scolamiero, I. Donato, and F. Vaccarino,PloS one , e66506 (2013).[7] G. P. Massara, T. Di Matteo, and T. Aste, Journal ofcomplex Networks , 161 (2016).[8] R. Sreejith, K. Mohanraj, J. Jost, E. Saucan, andA. Samal, Journal of Statistical Mechanics: Theory andExperiment , 063206 (2016).[9] A. P. Kartun-Giles and G. Bianconi, Chaos, Solitons &Fractals: X , 100004 (2019).[10] J. W. Rocks, A. J. Liu, and E. Katifori, Physical ReviewResearch , 033234 (2020).[11] Z. Wu, G. Menichetti, C. Rahmede, and G. Bianconi,Scientific reports , 1 (2015).[12] G. Bianconi and C. Rahmede, Scientific reports , 41974(2017).[13] G. Bianconi and C. Rahmede, Physical Review E ,032315 (2016).[14] M. M. Dankulov, B. Tadi´c, and R. Melnik, PhysicalReview E , 012309 (2019).[15] B. Tadi´c, M. Andjelkovi´c, and R. Melnik, Scientific re-ports , 1 (2019).[16] A. P. Mill´an, J. J. Torres, and G. Bianconi, PhysicalReview Letters , 218301 (2020).[17] J. J. Torres and G. Bianconi, Journal of Physics: Com-plexity , 015002 (2020).[18] M. Reitz and G. Bianconi, Journal of Physics A: Mathe-matical and Theoretical (2020).[19] S. Barbarossa and S. Sardellitti, IEEE Transactions onSignal Processing (2020).[20] N. Landry and J. G. Restrepo, arXiv preprintarXiv:2006.15453 (2020).[21] P. S. Skardal and A. Arenas, Physical review letters ,248301 (2019).[22] P. S. Skardal and A. Arenas, arXiv preprintarXiv:1909.08057 (2019).[23] I. Iacopini, G. Petri, A. Barrat, and V. Latora, Naturecommunications , 1 (2019).[24] D. Taylor, F. Klimm, H. A. Harrington, M. Kram´ar,K. Mischaikow, M. A. Porter, and P. J. Mucha, Naturecommunications , 1 (2015).[25] M. Lucas, G. Cencetti, and F. Battiston, Physical Re-view Research , 033410 (2020).[26] Y. Zhang, V. Latora, and A. E. Motter, arXiv preprintarXiv:2010.00613 (2020).[27] P. S. Skardal and A. Arenas, Journal of Physics: Com-plexity (2020).[28] L. DeVille, arXiv preprint arXiv:2010.07421 (2020).[29] T. Carletti, D. Fanelli, and S. Nicoletti, arXiv preprintarXiv:2006.01243 (2020).[30] A. P. Mill´an, J. J. Torres, and G. Bianconi, ScientificReports , 1 (2018).[31] R. Mulas, C. Kuehn, and J. Jost, arXiv preprintarXiv:2003.13775 (2020).[32] L. Gambuzza, F. Di Patti, L. Gallo, S. Lepri, M. Ro-mance, R. Criado, M. Frasca, V. Latora, and S. Boc-caletti, arXiv preprint arXiv:2004.03913 (2020). [33] A. P. Mill´an, J. J. Torres, and G. Bianconi, Phys. Rev.E , 022307 (2019).[34] F. P. U. Severino, J. Ban, Q. Song, M. Tang, G. Bian-coni, G. Cheng, and V. Torre, Scientific reports , 29640(2016).[35] M. Saggar, O. Sporns, J. Gonzalez-Castillo, P. A. Ban-dettini, G. Carlsson, G. Glover, and A. L. Reiss, Naturecommunications , 1 (2018).[36] C. Giusti, E. Pastalkova, C. Curto, and V. Itskov, Pro-ceedings of the National Academy of Sciences , 13455(2015).[37] M. W. Reimann, M. Nolte, M. Scolamiero, K. Turner,R. Perin, G. Chindemi, P. D(cid:32)lotko, R. Levi, K. Hess, andH. Markram, Frontiers in computational neuroscience ,48 (2017).[38] M. Ruiz-Garcia and E. Katifori, arXiv preprintarXiv:2001.01811 (2020).[39] S. H. Strogatz, Physica D: Nonlinear Phenomena , 1(2000).[40] S. Boccaletti, A. N. Pisarchik, C. I. Del Genio, andA. Amann, Synchronization: from coupled systems tocomplex networks (Cambridge University Press, 2018).[41] F. A. Rodrigues, T. K. D. Peron, P. Ji, and J. Kurths,Physics Reports , 1 (2016).[42] J. G. Restrepo, E. Ott, and B. R. Hunt, Physical ReviewE , 036151 (2005).[43] E. Ott and T. M. Antonsen, Chaos: An InterdisciplinaryJournal of Nonlinear Science , 037113 (2008).[44] W. Huang, T. A. Bolton, J. D. Medaglia, D. S. Bassett,A. Ribeiro, and D. Van De Ville, Proceedings of theIEEE , 868 (2018).[45] T. S. Evans and R. Lambiotte, The European PhysicalJournal B , 265 (2010).[46] R. M. D’Souza, J. G´omez-Garde˜nes, J. Nagler, andA. Arenas, Advances in Physics , 123 (2019).[47] X. Zhang, S. Boccaletti, S. Guan, and Z. Liu, Physicalreview letters , 038701 (2015).[48] S. Boccaletti, J. Almendral, S. Guan, I. Leyva, Z. Liu,I. Sendi˜na-Nadal, Z. Wang, and Y. Zou, Physics Reports , 1 (2016).[49] O. T. Courtney and G. Bianconi, Physical Review E ,062311 (2016).[50] R. W. Ghrist, Elementary applied topology , Vol. 1 (Cre-atespace Seattle, 2014).[51] D. Horak and J. Jost, Advances in Mathematics , 303(2013).[52] P. Hagmann, L. Cammoun, X. Gigandet, R. Meuli, C. J.Honey, V. J. Wedeen, and O. Sporns, PLoS Biol , e159(2008).[53] L. R. Varshney, B. L. Chen, E. Paniagua, D. H. Hall,and D. B. Chklovskii, PLoS Comput Biol , e1001066(2011).[54] P. S. Skardal, J. G. Restrepo, and E. Ott, Physical Re-view E , 060902 (2015).[55] T. Ichinomiya, Physical Review E70