Higher-order simplicial synchronization of coupled topological signals

ARTICLE Higher-order simplicial synchronization of coupled topological signals Reza Ghorbanchian1, Juan G. Restrepo 2✉, Joaquín J. Torres 3✉& Ginestra Bianconi 1,4✉ Simplicial complexes capture the underlying network topology and geometry of complex systems ranging from the brain to social networks. Here we show that algebraic topology is a fundamental tool to capture the higher-order dynamics of simplicial complexes. In particular we consider topological signals, i.e., dynamical signals defined on simplices of different dimension, here taken to be nodes and links for simplicity. We show that coupling between signals defined on nodes and links leads to explosive topological synchronization in which phases defined on nodes synchronize simultaneously to phases defined on links at a discontinuous phase transition. We study the model on real connectomes and on simplicial complexes and network models. Finally, we provide a comprehensive theoretical approach that captures this transition on fully connected networks and on random networks treated within the annealed approximation, establishing the conditions for observing a closed hysteresis loop in the large network limit. https://doi.org/10.1038/s42005-021-00605-4 OPEN 1School of Mathematical Sciences, Queen Mary University of London, London, UK. 2Department of Applied Mathematics, University of Colorado at Boulder, Boulder, CO, USA. 3Departamento de Electromagnetismo y Física de la Materia and Instituto Carlos I de Física Teórica y Computacional, Universidad de Granada, Granada, Spain. 4The Alan Turing Institute, The British Library, London, UK. ✉email: [email protected];[email protected];ginestra. [email protected] COMMUNICATIONS PHYSICS | (2021) 4:120 | https://doi.org/10.1038/s42005-021-00605-4 | www.nature.com/commsphys 1 1234567890():,; Higher-order networks1–4are attracting increasing attention as they are able to capture the many-body interactions of complex systems ranging from brain to social networks. Simplicial complexes are higher-order networks that encode the network geometry and topology of real datasets. Using simplicial complexes allows the network scientist to formulate new mathematical frameworks for mining data5–10 and for understanding these generalized network structures revealing the underlying deep physical mechanisms for emergent geometry11–15 and for higher-order dynamics16–33. In particular, this very vibrant research activity is relevant in neuroscience to analyze real brain data and its profound relation to dynamics1,6,15,34–37 and in the study of biological transport networks10,38. In networks, dynamical processes are typically defined over signals associated to the nodes of the network. In particular, the Kuramoto model39–43 investigates the synchronization of phases associated to the nodes of the network. This scenario can change significantly in thecaseofsimplicialcomplexes 16,17,19. In fact, simplicial complexes can sustain dynamical signals defined on simplices of different dimension, including nodes, links, triangles, and so on, called topological signals. For instance, topological signals defined on links can represent fluxes of interest in neuroscience and in biological transportation networks. The interest on topological signals is rapidly growing with new results related to signal processing17,19 and higherorder topological synchronization16,28. (Note that here higher-order refers to the higher-order interactions existing between topological signals and not to higher-order harmonics.) In particular, higherorder topological synchronization16 demonstrates that topological signals (phases) associated to higher dimensional simplices can undergo a synchronization phase transition. These results open a new uncharted territory for the investigation of higher-order synchronization. Higher-order topological signals defined on simplices of different dimension can interact with one another in non-trivial ways. For instance, in neuroscience the activity of the cell body of a neuron can interact with synaptic activity which can be directly affected by gliomes in the presence of brain tumors44. In order to shed light on the possible phase transitions that can occur when topological signals defined on nodes and links interact, here we build on the mathematical framework of higher-order topological synchronization proposed by Millán et al.16 and consider a synchronization model in which topological signals of different dimension are coupled. We focus in particular on the coupled synchronization of topological signals defined on nodes and links, but we note that the model can be easily extended to topological signals of higher dimension. The reason why we focus on topological signals defined on nodes and links is threefold. First of all we can have a better physical intuition of topological signals defined on nodes (traditionally studied by the Kuramoto model) and links (like fluxes) that is relevant in brain dynamics44,45 and biological transport networks10,38. Secondly, although the coupled synchronization dynamics of nodes and links can be considered as a special case of coupled synchronization dynamics of higher-order topological signals on a generic simplicial complex, this dynamics can be observed also on networks including only pairwise interactions. Indeed nodes and links are the simplices that remain unchanged if we reduce a simplicial complex to its network skeleton. Since currently there is more availability of network data than simplicial complex data, this fact implies that the coupled dynamics studied in this work has wide applicability as it can be tested on any network data and network model. Thirdly, defining the coupled dynamics of topological signals defined on nodes and links can open new perspectives in exploiting the properties of the line graph of a given network which is the network whose nodes corresponds to the links or the original network46. In this work, we show that by adopting a global adaptive coupling of dynamics47–49 the coupled synchronization dynamics of topological signals defined on nodes and links is explosive50, i.e., it occurs at a discontinuous phase transition in which the two topological signals of different dimension synchronize at the same time. We also illustrate numerical evidence of this discontinuity on real connectomes and on simplicial complex models, including the configuration model of simplicial complexes51 and the nonequilibrium simplicial complex model called Network Geometry with Flavor (NGF)12,13. We provide a comprehensive theory of this phenomenon on fully connected networks offering a complete analytical understanding of the observed transition. This approach can be extended to random networks treated within the annealed network approximation. The analytical results reveal that the investigated transition is discontinuous. Results and discussion Higher-order topological Kuramoto model of topological signals of a given dimension. Let us consider a simplicial complex Kformed by N [n] simplices of dimension n, i.e., N [0] nodes, N [1] links, N [2] triangles, and so on. In order to define the higher-order synchronization of topological signals we will make use of algebraic topology (see the Appendix for a brief introduction) and specifically we indicate with B [n] the nth incidence matrix representing the nth boundary operator. The higher-order Kuramoto model generalizes the classic Kuramoto model to treat synchronization of topological signals of higher-dimension. The classic Kuramoto model describes the synchonization transition for phases θ¼ðθ1;θ2;¼θN½0Þð1Þ associated to nodes, i.e., simplices of dimension n=0 (see Fig. 1). The Kuramoto model is typically defined on a network but it can treat also synchronization of the phases associated to the nodes of a simplicial complex. Each node ihas associated an internal frequency ω i drawn from a given distribution, for instance a Fig. 1 Schematic representation of the Kuramoto and the higher-order topological synchronization model. Panel a shows a network formed by nodes and links in which nodes (blue circles numbered from 1 to 8) sustain a dynamical variable (a phase θ i with i∈{1, 2, …, 8}) whose synchronization is captured by the Kuramoto model. Panel bshows a simplicial complex formed by nodes, links, and triangles (here shaded in orange) in which not only nodes but also links sustain dynamical variables (indicated θ i for the nodes i∈{1, 2, …, 8} and ϕ ij for the links [i,j] with i,j∈{1, 2, …, 8}) whose coupled synchronization dynamics is captured by the higher-order topological Kuramoto model. ARTICLE COMMUNICATIONS PHYSICS | https://doi.org/10.1038/s42005-021-00605-4 2COMMUNICATIONS PHYSICS | (2021) 4:120 | https://doi.org/10.1038/s42005-021-00605-4 | www.nature.com/commsphys normal distribution ωiNðΩ0;1=τ0Þ. In absence of any coupling, i.e., in absence of pairwise interactions, every node oscillates at its own frequency. However in a network or in a simplicial complex skeleton the phases associated to the nodes follow the dynamical evolution dictated by the equation: _ θ¼ωσB½1sin B> ½1θ  ;ð2Þ where here and in the following we use the notation sinðxÞto indicate the column vector where the sine function is taken elementwise. Note that here we have chosen to write this system of equations in terms of the incidence matrix B [1] . However if we indicate with athe adjacency matrix of the network and with a ij its matrix elements, this system of equations is equivalent to _ θi¼ωiþσ∑ N j¼1aij sinðθjθiÞ;ð3Þ valid for every node iof the network. For coupling constant σ =σ c the Kuramoto model39–41 displays a continuous phase transition above which the order parameter R0¼1 N½0 ∑ N½0 i¼1eiθi ð4Þ isnon-zeroalsointhelimitN [0] →∞. The higher-order topological Kuramoto model16 describes synchronization of phases associated to the ndimensional simplices of a simplicial complex. Although the definition of the model applies directly to any value of n,hereweconsider specifically the case in which the higher-order Kuramoto model is defined on topological signals (phases) associated to the links ϕ¼ðϕ‘1;ϕ‘2;¼ϕ‘N½1Þ;ð5Þ where ϕ‘rindicates the phase associated to the rth link ℓ r of the simplicial complex (see Fig. 1). The higher-order Kuramoto dynamics defined on simplices of dimension n> 0 is the natural extension of the standard Kuramoto model defined by Eq. (2). Let us indicate with ~ ωthe internal frequencies associated to the links of the simplicial complex, sampled for example from a normal distribution, ~ ω‘NðΩ1;1=τ1Þ. The higher-order topological Kuramoto model is defined as _ ϕ¼~ ωσB> ½1sinðB½1ϕÞσB½2sinðB> ½2ϕÞ:ð6Þ 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 dynamics can be projected on (n+1) and (n−1) simplices. Interestingly, algebraic topology provides a clear solution to this question. Indeed for n=1, when the dynamics describes the evolution of phases associated to the links, one can consider the projection ϕ[−]and ϕ[+], respectively, on nodes and on triangles defined as ϕ½ ¼B½1ϕ; ϕ½þ ¼B> ½2ϕ:ð7Þ Note that in this case B [1] acts as a discrete divergence and B> ½2 acts as a discrete curl. Interestingly, since the incidence matrices satisfy B [1] B [2] =0and B> ½2B> ½1¼0(see “Methods”) these two projected phases follow the uncoupled dynamics _ ϕ½ ¼B½1~ ωσL½0sin ϕ½; _ ϕ½þ ¼B> ½2~ ωσLdown ½2sin ϕ½þ;ð8Þ where L½0¼B½1B> ½1and Ldown ½2¼B> ½2B½2. These two projected dynamics undergo a continuous synchronization transition at σ c =016 with order parameters Rdown 1¼1 N½0 ∑ N½0 i¼1eiϕ½ i  ; Rup 1¼1 N½2 ∑ N½2 i¼1eiϕ½þ i  :ð9Þ In Millán et al.16 an adaptive coupling between these two dynamics is considered formulating the explosive higher-order topological Kuramoto model in which the topological signal follows the set of coupled equations _ ϕ¼~ ωσRup 1B> ½1sinðB½1ϕÞ σRdown 1B½2sinðB> ½2ϕÞ:ð10Þ The projected dynamics on nodes and triangles are now coupled by the modulation of the coupling constant σwith the order parameters Rdown 1and Rup 1, i.e. the two projected phases follow the coupled dynamics _ ϕ½ ¼B½1~ ωσRup 1L½0sin ϕ½; _ ϕ½þ ¼B> ½2~ ωσRdown 1Ldown ½2sin ϕ½þ:ð11Þ This explosive higher-order topological Kuramoto model hasbeenshowninMillánetal. 16 to lead to a discontinuous synchronization transition on different models of simplicial complexes and on clique complexes of real connectomes. Higher-order topological Kuramoto model of coupled topological signals of different dimension. Until now, we have captured synchronization occurring only among topological signals of the same dimension. However, signals of different dimension can be coupled to each other in non-trivial ways. In this work we will show how topological signals of different dimensions can be coupled together leading to an explosive synchronization transition. Specifically we focus on the coupling of the traditional Kuramoto model [Eq. (2)] to a higherorder topological Kuramoto model defined for phases associated to the links [Eq. (6)]. The coupling between these two dynamics is here performed considering the modulation of the coupling constant σwith the global order parameters of the node dynamics [defined in Eq. (4)] and the link dynamics [defined in Eq. (9)]. Specifically, we consider two models denoted as Model Nodes-Links (NL) and Model Nodes-LinksTriangles (NLT). Model NL couples the dynamics of the phases of the nodes θand of the links ϕaccordingtothefollowing dynamical equations _ θ¼ωσRdown 1B½1sinðB> ½1θÞ;ð12Þ _ ϕ¼~ ωσR0B> ½1sinðB½1ϕÞσB½2sinðB> ½2ϕÞ:ð13Þ The projected dynamics for ϕ[−]and ϕ[+]then obeys _ ϕ½ ¼B½1~ ωσR0L½0sin ϕ½;ð14Þ _ ϕ½þ ¼B> ½2~ ωσLdown ½2sin ϕ½þ:ð15Þ Therefore the projection on the nodes ϕ[−]of the phases ϕ associated to the links [Eq. (14)] is coupled to the dynamics of the phases θ[Eq. (12)] associated directly to nodes. However the projection on the triangles ϕ[+]of the phases ϕassociated to the links is independent of ϕ[−]and of θas well. Model NLT also describes the coupled dynamics of topological signals definedonnodesandlinksbuttheadaptive coupling captured COMMUNICATIONS PHYSICS | https://doi.org/10.1038/s42005-021-00605-4 ARTICLE COMMUNICATIONS PHYSICS | (2021) 4:120 | https://doi.org/10.1038/s42005-021-00605-4 | www.nature.com/commsphys 3 by the model is different. In this case the dynamical equations are taken to be _ θ¼ωσRdown 1B½1sinðB> ½1θÞ;ð16Þ _ ϕ¼~ ωσR0Rup 1B> ½1sinðB½1ϕÞ σRdown 1B½2sinðB> ½2ϕÞ:ð17Þ For Model NLT the projected dynamics for ϕ[−]and for ϕ[+] obeys _ ϕ½ ¼B½1~ ωσR0Rup 1L½0sin ϕ½;ð18Þ _ ϕ½þ ¼B> ½2~ ωσRdown 1Ldown ½2sin ϕ½þ:ð19Þ Therefore, as in Model NL, the dynamics of the projection ϕ[−] of the phases ϕassociated to the links [Eq. (18)] is coupled to the dynamics of the phases θassociated directly to nodes [Eq. (16)] and vice versa. Moreover, the dynamics of the projection of the phases ϕon the triangles ϕ[+][Eq. (19)] is now also coupled with the dynamics of ϕ[−][Eq. (18)] and vice versa. Here and in the following we will use the convenient notation (using the parameter r) to indicate both models NL and NLT with the same set of dynamical equations given by _ θ¼ωσRdown 1B½1sinðB> ½1θÞ;ð20Þ _ ϕ¼~ ωσR0Rup 1  r1B> ½1sinðB½1ϕÞ σRdown 1  r1B½2sinðB> ½2ϕÞ;ð21Þ which reduce to Eq. (13) for r=1 and to Eq. (17) for r=2. We make two relevant observations: ●First, the proposed coupling between topological signals of different dimension can be easily extended to signals defined on higher-order simplices providing a very general scenario for coupled dynamical processes on simplicial complexes. ●Second, the considered coupled dynamics of topological signals defined on nodes and links can be also studied on networks with exclusively pairwise interactions where we assume that the number of simplices of dimension n>1is zero. Therefore in this specific case this topological dynamics can have important effects also on simple networks. We have simulated Model NL and Model NLT on two main examples of simplicial complex models: the configuration model of simplicial complexes51 and the NGF12,13 (see Fig. 2). In the configuration model we have considered power-law distribution of the generalized degree with exponent γ< 3, and for the NGF model with have considered simplicial complexes of dimensions d=3 whose skeleton is a power-law network with exponent γ= 3. In both cases we observe an explosive synchronization of the topological signals associated to the nodes and to the links. On finite networks, the discontinuous transition emerge together with the hysteresis loop formed by the forward and backward synchronization transition. However the two models display a notable difference. In Model NL we observe a discontinuity for R 0 and Rdown 1at a non-zero coupling constant σ=σ c ; however, Rup 1 follows an independent transition at zero coupling (see Fig. 2, panels in the second and fourth column). In Model NLT, on the contrary, all order parameters R 0 ,Rdown 1, and Rup 1display a discontinuous transition occurring for the same non zero value of the coupling constant σ=σ c (see Fig. 2panels in the first and third column). This is a direct consequence of the fact that in Model NL the adaptive coupling leading to discontinuous phase transition only couples the phases ϕ[−]and θ, while for Model NLT the coupling involves also the phases ϕ[+]. Additionally we studied both Model NL and Model NTL on two real connectomes: the human connectome52 and the Caenorhabditis elegans (C. elegans) connectome53 (see Fig. 3). Interestingly also for these real datasets we observe that in Model NL the explosive synchronization involves only the phases θand ϕ[−]while in Model NLT we observe that also ϕ[+]undergoes an explosive synchronization transition at the same value of the coupling constant σ=σ c . Theoretical solution of the NL model. As mentioned earlier the higher-order topological Kuramoto model coupling the topological signals of nodes and links can be defined on simplicial complexes and on networks as well. In the following sections we exploit this property of the dynamics to provide an analytical understanding of the synchronization transition on uncorrelated random networks. It is well known that the Kuramoto model is challenging to study analytically. Indeed the full analytical understanding of the model is restricted to the fully connected case, while on a generic sparse network topology the analytical approximation needs to rely on some approximations. A powerful approximation is the annealed network approximation41 which consists in approximating the adjacency matrix of the network with its expectation in a random uncorrelated network ensemble. In order to unveil the fundamental theory that determines the coupled dynamics of topological signals described by the higher-order Kuramoto model here we combine the annealed approximation with the Ott–Antonsen method43. This approach is able to capture the coupled dynamics of topological signals defined on nodes and links. In particular, the solution found to describe the dynamics of topological signals defined on the links is highly non-trivial and it is not reducible to the equations valid for the standard Kuramoto model. Conveniently, the calculations performed in the annealead approximation can be easily recasted in the exact calculation valid in the fully connected case previous a rescaling of some of the parameters. The analysis of the fully connected network reveals that the discontinuous sychronization transition of the considered model is characterized by a non-trivial backward transition with a well defined large network limit. On the contrary, the forward transition is highly dependent on the network size and vanishes in the large network limit, indicating that the incoherent state remains stable for every value of the coupling constant σin the large network limit. This implies that on a fully connected network the NL model does not display a closed hysteresis loop as it occurs also for the model proposed in Skardal and Arenas21. This scenario is here shown to extend also to sparse networks with finite second moment of the degree distribution while scalefree networks display a well defined hysteresis loop in the large network limit. Annealed dynamics. For the dynamics of the phases θassociated to the nodes—Eq. (20)—it is possible to proceed as in the traditional Kuramoto model42,54,55. However, the annealed approximation for the dynamics of the phases ϕdefined in Eq. (21) needs to be discussed in detail as it is not directly reducible to previous results. To address this problem our aim is to directly define the annealed approximation for the dynamics of the projected variables ϕ[−]which, here and in the following, are indicated as ψ¼ϕ½;ð22Þ in order to simplify the notation. Moreover we will indicate with N=N [0] the number of nodes in the network or in the simplicial ARTICLE COMMUNICATIONS PHYSICS | https://doi.org/10.1038/s42005-021-00605-4 4COMMUNICATIONS PHYSICS | (2021) 4:120 | https://doi.org/10.1038/s42005-021-00605-4 | www.nature.com/commsphys complex skeleton. Here we focus on the NL Model defined on networks, i.e., we assume that there are no simplices of dimension two. We provide an analytical understanding of the coupled dynamics of nodes and links in the NL Model by determining the equations that capture the dynamics in the annealed approximation and predict the value of the complex order parameters R0eiΘ¼1 N∑ N i¼1eiθi; Rdown 1eiΨ¼1 N∑ N i¼1eiψi;ð23Þ (with R0;Rdown 1;Θ, and Ψreal) as a function of the coupling constant σ. We notice that Eq. (14), valid for Model NL, can be written as _ ψ¼B½1~ ωσR0L½0sinðψÞ:ð24Þ This equation can be also written elementwise as _ ψi¼^ ωiþσR0∑ N j¼1aij sinðψjÞsinðψiÞ hi ;ð25Þ where the vector ^ ωis given by ^ ω¼B½1~ ω:ð26Þ Let us now consider in detail these frequencies in the case in which the generic internal frequency ~ ω‘of a link follows a Gaussian distribution, specifically in the case in which ~ ω‘ NðΩ1;1=τ1Þfor every link ℓ. Using the definition of the boundary operator on a link it is easy to show that the expectation of ^ ωiis given by ^ ωi  ¼∑ j<i aij ∑ j>i aij  Ω1:ð27Þ Given that each node has degree k i , the covariance matrix Cis given by the graph Laplacian L [0] of the network, i.e. Cij ¼^ ωi^ ωj DE c¼∑ ‘;‘0½B½1~ ωi½B½1~ ωj DE c ¼½L½0ij τ2 1¼kiδij aij τ2 1 ;ð28Þ where we have indicated with ¼ hi cthe connected correlation. Therefore the variance of ^ ωin the annealed approximation is ^ ω2 i  c¼^ ω2 i  ^ ωi  2¼ki τ2 1 :ð29Þ Moreover, the projected frequencies are actually correlated and for i≠jwe have ^ ωi^ ωj DE c¼^ ωi^ ωj DE ^ ωi ^ ωj DE ¼aij τ2 1 :ð30Þ It follows that the frequencies ^ ωare correlated Gaussian variables with average given by Eq. (27) and correlation matrix given by the graph Laplacian. The fact that the frequencies ^ ωiare correlated is an important feature of the dynamics of ψand, with few exceptions56, this feature has remained relatively unexplored in the case of the standard Kuramoto model. Additionally let us note that the average of ^ ωover all the nodes of the network is 0246 0 0.5 1 0246 0 0.5 1 0246 0 0.5 1 0246 0 0.5 1 0246 0 0.5 1 0246 0 0.5 1 0246 0 0.5 1 0246 0 0.5 1 0246 0 0.5 1 0246 0 0.5 1 0246 0 0.5 1 0246 0 0.5 1 (d) (c) (b) (a) (f) (g) (h) (e) (i) (j) (k) (l) Fig. 2 The higher-order topological synchronization models coupling nodes and links on simplicial complexes. The hysteresis loop for the synchronization order parameters R 0 ,Rdown 1, and Rup 1are plotted versus σfor the higher-order topological synchronization Model Nodes-Links-Triangles (NLT) (panels a,e,iand c,g,k) and Model Nodes-Links (NL) (panels b,f,jand d,h,l)defined over the Network Geometry with Flavor13 (panels a,e,iand b,f,j) and the configuration model of simplicial complexes51 (panels c,g,kand d,h,l). The green lines indicate the backward transitions and the cyan lines indicate the forward transitions. The Network Geometry with Flavor on which we run the numerical results shown in a,b,e,f,i,jincludes N [0] =500 nodes and has flavor s=−1 and d=3. The configuration model of simplicial complexes on which we run the numerical results shown in c,d,g,h,k,l includes N [0] =500 nodes and has generalized degree distribution which is power-law with exponent γ=2.8. In both Model NL and in Model NLT we have set Ω 0 =Ω 1 =2 and τ 0 =τ 1 =1. COMMUNICATIONS PHYSICS | https://doi.org/10.1038/s42005-021-00605-4 ARTICLE COMMUNICATIONS PHYSICS | (2021) 4:120 | https://doi.org/10.1038/s42005-021-00605-4 | www.nature.com/commsphys 5 zero. In fact ∑ N i¼1 ^ ωi¼1T^ ω¼1TB½1ω¼0;ð31Þ where with 1we indicate the N-dimensional column vector of elements 1 i =1. By using the symmetry of the adjacency matrix, i.e. the fact that a ij =a ji , Eq. (31) implies that the sum of _ ψiover all the nodes of the network is zero, i.e. ∑ N i¼1 _ ψi¼∑ N i¼1 ^ ωiþσR0∑ i;jaij½sinðψjÞsinðψiÞ ¼ 0: We now consider the annealed approximation consisting in substituting the adjacency matrix element a ij with its expectation in an uncorrelated network ensemble aij !kikj k hi N;ð32Þ where k i indicates the degree of node iand k hiis the average degree of the network. Note that the considered random networks can be both sparse57 or dense58 as long as they display the structural cutoff, i.e. kiffiffiffiffiffiffiffiffiffiffi k hi N pfor every node iof the network. In the annealed approximation we can put ^ ωi  ’kiΩ112∑ j>i kj k hi N  :ð33Þ Also, in the annealed approximation the dynamical Eq. (20) and Eq. (24) reduce to _ θ¼ωσRdown 1^ R0ksinðθ^ ΘÞ;ð34Þ _ ψ¼^ ωþσR0^ Rdown 1ksin ^ ΨσR0ksin ψ;ð35Þ where ⊙indicates the Hadamard product (element by element multiplication) and where two auxiliary complex order parameters are defined as ^ R0ei^ Θ¼∑ N i¼1 ki k hi Neiθi; ^ Rdown 1ei^ Ψ¼∑ N i¼1 ki k hi Neiψi;ð36Þ with ^ R0;^ Θ;^ Rdown 1and ^ Ψreal. The dynamics on a fully connected network. On a fully connected network in which each node has degree k i =N−1 the dynamics of the NL Model is well defined provided its parameter are properly rescaled. In particular, we require a standard rescaling of the coupling constant with the network size, given by σ!σ=ðN1Þð37Þ which guarantees that the interaction term in the dynamical equations has a finite contribution to the velocity of the phases. The Model NL on fully connected networks requires also some specific model dependent rescalings associated to the dynamics on networks. Indeed, in order to have a finite expectation ^ ωi  of the projected frequencies ^ ωiand a finite of the covariance matrix C[given by Eqs. (27) and (28), respectively], we require that on a fully connected network both Ω 1 and τ 1 are rescaled according to Ω1!Ω1=N; τ1!τ1ffiffiffiffiffiffiffiffiffiffiffiffi N1 p:ð38Þ Considering these opportune rescalings and noticing that the order parameters obey ^ R0¼R0,^ Rdown 1¼Rdown 1,Θ¼^ Θ, and Ψ¼^ Ψ, we obtain that Model NL dictated by Eqs. (34)–(35) can 05 0 0.5 1 05 0 0.5 1 05 0 0.5 1 05 0 0.5 1 05 0 0.5 1 05 0 0.5 1 024 0 0.5 1 05 0 0.5 1 05 0 0.5 1 024 0 0.5 1 05 0 0.5 1 05 0 0.5 1 (b) (c) (d)(a) (f) (g) (h) (e) (i) (j) (k) (l) Fig. 3 The higher-order topological synchronization models coupling nodes and links on real connectomes. The hysteresis loop for the synchronization order parameters R 0 ,Rdown 1, and Rup 1are plotted versus σon real connectomes. The green lines indicate the backward transitions and the cyan lines indicate the forward transitions. Panels a,e,iand b,f,jshow the numerical results on the human connectome52 for Model Nodes-Links-Triangles (NLT) and Model Nodes-Links (NL) respectively. Panels c,g,kand d,h,lshow the numerical results on the C. elegans connectome53 for Model NLT and Model NL, respectively. In both Model NLT and in Model NL we have set Ω 0 =Ω 1 =2 and τ 0 =τ 1 =1. ARTICLE COMMUNICATIONS PHYSICS | https://doi.org/10.1038/s42005-021-00605-4 6COMMUNICATIONS PHYSICS | (2021) 4:120 | https://doi.org/10.1038/s42005-021-00605-4 | www.nature.com/commsphys be rewritten here as _ θ¼ωσRdown 1R0sinðθΘÞ;ð39Þ _ ψ¼^ ωþσR0Rdown 1sin ΨσR0sin ψ;ð40Þ with R0;Rdown 1;Θand Ψgiven by Eq. (23) and Cij ¼^ ωi^ ωj DE c¼δij 1 N1:ð41Þ Solution of the dynamical equations in the annealed approximation General framework for obtaining the solution of the annealed dynamical equations. In this section we will provide the analytic solutions for the order parameter of the higher-order topological synchronization studied within the annealed approximation, i.e., captured by Eqs. (34) and (35). In particular, first we will find an expression of the order parameters R 0 of the dynamics associated to the nodes (Eq. (34)) and subsequently in the next paragraph we will derive the expression for the order parameter Rdown 1associated to the projection on the nodes of the topological signal defined on the links (Eq. (35)). By combining the two results it is finally possible to uncover the discontinuous nature of the transition. Dynamics of the phases of the nodes. When we investigate Eq. (34) we notice that this equation can be easily reduced to the equation for the standard Kuramoto model treated within the annealed approximation42 if one performs a rescaling of the coupling constant σR 0 →σ. Therefore we can treat this model similarly to the known treatment of the standard Kuramoto model40–42. Specifically, starting from Eq. (34) and using a rescaling of the phases θaccording to θi!θiΩ0t;ð42Þ it is possible to show that we can set Θ=0 and therefore Eq. (34) reduces to the well-known annealed expression for the standard order Kuramoto model given by _ θ¼ωΩ01σRdown 1^ R0ksinðθÞ:ð43Þ Assuming that the system of equations reaches a steady state in which both Rdown 1and ^ R0become time independent, the order parameters of this system of equations in the coherent state ^ R0>0 and Rdown 1>0 can be found to obey40,42,50,54 ^ R0¼∑ N i¼1 ki k hi NZj^ cij<1 dωgðωÞffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffi 1ωΩ0 σki^ R0Rdown 1 ! 2 v u u t; R0¼1 N∑ N i¼1Zj^ cij<1 dωgðωÞffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffi 1ωΩ0 σki^ R0Rdown 1 ! 2 v u u t; ð44Þ where ^ ciindicates ^ ci¼ωΩ0 σkj^ R0Rdown 1 :ð45Þ and g(ω) is the Gaussian distribution with expectation Ω 0 and standard deviation 1. Dynamics of the phases of the links projected on the nodes. In this paragraph we will derive the expression of the order parameters Rdown 1and ^ Rdown 1which, together with Eq. (44), will provide the annealed solution of our model. To start with we assume that the frequencies ^ ωare known. In this case we can express the order parameters Rdown 1and ^ Rdown 1as a function of the probability density function ρðiÞðψ;tj^ ωÞthat node iis associated to a projected phase of the link equal to ψ. Since in the annealed approximation ψ i has a dynamical evolution dictated by Eq. (35) the probability density function obeys the continuity equation ∂tρðiÞðψ;tj^ ωÞþ∂ψρðiÞðψ;tj^ ωÞvi  ¼0ð46Þ with associated velocity v i given by vi¼κiσR0kisin ψi;ð47Þ where we have defined κ i as κi¼^ ωiþσkiR0^ Rdown 1sin ^ Ψ:ð48Þ In this case the complex order parameters are given by ^ Rdown 1eiΨ¼∑ N i¼1 ki k hi NZdψρðiÞðψ;tj^ ωÞeiψ; Rdown 1ei~ Ψ¼∑ N i¼1 1 NZdψρðiÞðψ;tj^ ωÞeiψ:ð49Þ In order to solve the continuity equation we follow Ott–Antonsen43 and we express ρðiÞðψ;tj^ ωÞin the Fourier basis as ρðiÞðψ;tj^ ωÞ¼ 1 2π1þ∑ 1 m¼1 ^ fðiÞ mð^ ωi;tÞeimψþc:c:  :ð50Þ Making the ansatz ^ fðiÞ mð^ ωi;tÞ¼½bið^ ωi;tÞmð51Þ we can derive the equation for the evolution of bi¼bið^ ωi;tÞgiven by ∂tbiþibiκiþσkiR0 1 2ðb2 i1Þ¼0:ð52Þ Since we showed before that the average value of _ ψiover nodes is zero, we look for non-rotating stationary solutions of Eq. (52), ∂ t b i =0. As long as R 0 > 0 these stationary solutions are given by bi¼idi±ffiffiffiffiffiffiffiffiffiffiffiffiffi 1d2 i q;ð53Þ where d i is given by di¼^ ωi σkiR0þ^ Rdown 1sin ^ Ψ:ð54Þ By inserting this expression into Eq. (49) we get the expression of the order parameters given the projected frequencies ^ ω, in the coherent phase in which R 0 >0 ^ Rdown 1cos ^ Ψ¼∑ N i¼1 ki k hi Nffiffiffiffiffiffiffiffiffiffiffiffiffi 1d2 i qθð1d2 iÞ; ^ Rdown 1sin ^ Ψ¼∑ N i¼1 ki k hi Nffiffiffiffiffiffiffiffiffiffiffiffiffi d2 i1 qχðdiÞþdi  ; Rdown 1cos Ψ¼∑ N i¼1 1 Nffiffiffiffiffiffiffiffiffiffiffiffiffi 1d2 i qθð1d2 iÞ; Rdown 1sin Ψ¼∑ N i¼1 1 Nffiffiffiffiffiffiffiffiffiffiffiffiffi d2 i1 qχðdiÞþdi  ; ð55Þ where, indicating by θ(x) the Heaviside function, we have defined χðdiÞ¼½θðdi1Þþθð1diÞ:ð56Þ Finally, if the projected frequencies ^ ωare not known we can average the result over the marginal frequency distribution of the COMMUNICATIONS PHYSICS | https://doi.org/10.1038/s42005-021-00605-4 ARTICLE COMMUNICATIONS PHYSICS | (2021) 4:120 | https://doi.org/10.1038/s42005-021-00605-4 | www.nature.com/commsphys 7 projected frequency ^ ωigiven by Gið^ ωÞ ^ Rdown 1cos ^ Ψ¼∑ N i¼1 ki k hi NZjdij≤1 d^ ωiGið^ ωiÞffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffi 1^ ωi σR0kiþ^ Rdown 1sin ^ Ψ  2 s; ^ Rdown 1sin ^ Ψ¼∑ N j¼0 ki k hi NZdi>1 d^ ωiGið^ ωiÞffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffi ^ ωi σR0kiþ^ Rdown 1sin ^ Ψ  2 1 s þ∑ N i¼1 ki khiNZdi<1 d^ ωiGið^ ωiÞffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffi ^ ωi σR0kiþ^ Rdown 1sin Ψ  2 1 s þ∑ N i¼1 ki k hi NZ1 1 d^ ωiGið^ ωiÞ^ ωi σR0kiþ^ Rdown 1sin Ψ  ; Rdown 1cos Ψ¼∑ N i¼1 1 NZjdij≤1 d^ ωiGið^ ωiÞffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffi 1^ ωi σR0kiþ^ Rdown 1sin ^ Ψ  2 s; ð57Þ and an analogous equations for Rdown 1sinðΨÞ(not shown). We note that in the case of distributions g(ω) and Gið^ ωÞthat are symmetric around their means the above equations always admit the solution Ψ¼^ Ψ¼0. Such values of the phases are also confirmed by direct numerical integration of the NL model. These equations together with Eq. (44) capture the steady-state behavior of the higher-order Kuramoto model coupling topological signals defined on nodes and links within the annealed approximation in the coherent synchronized phase. Note that by derivation, these equations cannot capture the asynchronous phase which is instead always a trivial solution of the dynamical equations corresponding to R0¼Rdown 1¼0. Finally we observe that for the NL Model as well as for the standard Kuramoto model on random networks, it is expected that the annealed approximation is more accurate for networks that are connected and are sufficiently dense. To illustrate the applicability of the theoretical analysis, we consider two examples of connected networks with N=1600 nodes: a Poisson network with average degree c=12 and an uncorrelated scale-free network with minimum degree m=6 and power-law exponent γ=2.5 In Fig. 4we compare the values of R 0 ,Rdown 1obtained from direct numerical integration of Eqs. (20) and (25) and the steady-state solutions obtained from the numerical solution of Eq. (55). The backward transition is fully captured by our theory, while the next paragraphs will clarify the theoretical expectations for the forward transition. Solution on the fully connected network. The integration of Eq. (57) requires the knowledge of the marginal distributions Gið^ ωÞ which does not have in general a simple analytical expression. However, in the fully connected networks with Gaussian distribution of the internal frequency of nodes and links this calculation simplifies significantly. Indeed, when the link frequencies are sampled from a Gaussian distribution with mean Ω 1 /Nand standard deviation 1=ðτ1ffiffiffiffiffiffiffiffiffiffiffiffi N1 pÞ, the marginal frequency distribution Gið^ ωÞof the internal frequency ^ ωiof a node iin a fully connected network is given by (see “Methods”for details) Gið^ ωÞ¼ τ1 ffiffiffiffiffiffiffiffiffiffi 2π= c pexp τ2 1 cð^ ωi^ ωi  Þ2 2 "# ;ð58Þ where  c¼N N1. By considering Ω0¼Ω1¼^ ωi  ¼0;and performing a direct integration of Eq. (57) we obtain (see “Methods” section for details) the closed system of equations for R 0 and Rdown 1 1¼σRdown 1hσ2R2 0ðRdown 1Þ2  ; Rdown 1¼σR0τ1ffiffi c phσ2τ2 1R2 0  ;ð59Þ where the scaling function h(x) is given by hðxÞ¼ ffiffiffi π 2 rex=4I0 x 4 þI1 x 4 hi ;ð60Þ with I 0 and I 1 indicating the modified Bessel functions. The 0.0 0.5 1.0 1.5 2.0 0.0 0.2 0.4 0.6 0.8 1.0 R0 0.0 0.5 1.0 1.5 2.0 0.0 0.2 0.4 0.6 0.8 1.0 R0 0.0 0.5 1.0 1.5 2.0 0 0.2 0.4 0.6 0.8 1 R1 down 0.0 0.5 1.0 1.5 2.0 0.0 0.1 0.2 0.3 0.4 0.5 R1 down (a) (b) (c) (d) Fig. 4 Comparison between the simulation results of the Nodes-Links (NL) Model and its solution in the annealed approximation. The hysteresis loop for the synchronization order parameters R 0 and Rdown 1of the NL Model are shown as a function of σfor a Poisson network with average degree c=12 (a,c) and for an uncorrelated scale-free network with minimum degree m=6 and power-law exponent γ=2.5 (b,d). Both networks have N=1600 nodes. The symbols indicate the simulation results for the forward (cyan diamonds) and the backward (green circles) synchronization transition. The solid black lines indicate the analytical solution for the backward transition obtained by integrating Eq. (55). ARTICLE COMMUNICATIONS PHYSICS | https://doi.org/10.1038/s42005-021-00605-4 8COMMUNICATIONS PHYSICS | (2021) 4:120 | https://doi.org/10.1038/s42005-021-00605-4 | www.nature.com/commsphys numerical solution of Eq. (59) reveals the following picture: for low values of σ, only the incoherent solution R0¼Rdown 1¼0 exists. At a positive value of σ, two solutions of Eq. (59) appear at a bifurcation point, with the upper solution corresponding to a stable synchronized state and the lower solution to an unstable synchronized solution. For larger values of σ, the values of R 0 and Rdown 1corresponding to the upper solution approach one (full phase synchronization), while those for the lower solution approach zero asymptotically, thus indicating that the incoherent state never loses stability. Indeed, it can be easily checked (see “Methods”for details) that for large σthe unstable solution of Eq. (59) has asymptotic behavior R0¼σ2J0; Rdown 1¼σ1J1;ð61Þ with J 0 and J 1 constants given by J0¼π 2 hi 2Gð0Þgð0Þ  1;ð62Þ J1¼gð0Þπ 2 hi 1:ð63Þ Therefore, the unstable branch approaches the trivial solution R0¼Rdown 1¼0 only asymptotically for σ→∞. This implies that the trivial solution remains stable for every possible value of σ although as σincreases it describes the stationary state of an increasingly smaller set of initial conditions. This scenario is confirmed by numerical simulations (see Fig. 5) showing that the backward transition is captured very well by our theory and does not display notable finite-size effects. The forward transition, instead, displays remarkable finite-size effects. Indeed, as σincreases, the system remains in the incoherent state until it explosively synchronizes at a positive value of σand reaches the stable synchronized branch. However the incoherent state is stable in the limit N→∞, and this forward transition is the result of finite-size fluctuations that push the system above the unstable branch, causing the observed explosive transition. This is consistent with the fact that for larger values of N, which have smaller finite-size fluctuations, the system remains in the incoherent state for larger values of σ. Therefore, while a closed hysteresis loop is not present in the NL model defined on fully connected networks, we observe fluctuation-driven hysteresis, in which finite-size fluctuations of the zero solution drive the system towards the synchronized solution, creating an effective hysteresis loop. Hysteresis on homogeneous and scale-free networks. In this section we discuss how the scenario found for the fully connected network can be extended to random networks with given degree distribution. We will start from the self-consistent Eq. (57) obtained within the annealed approximation model. These equations display a saddle point bifurcation with the emergence of two non-trivial solutions describing a stable and an unstable branch of these self-consistent equations. These solutions always exist in combination with the trivial solution R0¼Rdown 1¼0 describing the asynchronous state. Two scenarios are possible: either the unstable branch converges to the trivial solution only in the limit σ→∞or it converges to the trivial solution at a finite value of σ. In the first case, the scenario is the same as the one observed for the fully connected network, and the trivial solution remains stable for any finite value of σ. In this case the forward transition is not obtained in the limit N→∞and the transition observed on finite networks is only caused by finite-size effects. In the second case the trivial solution loses its stability at a finite value of σ. Therefore the forward transition is not subjected to strong finite-size effects and we expect to see a forward transition also in the N→∞limit. in order to determine which network 0123456789 0.0 0.2 0.4 0.6 0.8 1.0 R0 0123456789 0.0 0.2 0.4 0.6 0.8 1.0 R 1 down Fig. 5 The backward and the forward discontinuous phase transition on fully connected networks. The order parameters R 0 (circles) and Rdown 1 (squares) are plotted as a function of the coupling constant σon a fully connected network. The solid and the dashed lines indicate the stable branch and the unstable branch predicted by Eq. (59). Simulations (shown as data point) are here obtained by integrating numerically Eqs. (34) and (35) for a fully connected network of N=500 (cyan circles), N=1000 (green squares), and N=2000 (purple diamonds) with Ω 0 =Ω 1 =0 and (rescaled) τ 0 =τ 1 =1. The backward transition is perfectly captured by the theoretical prediction (solid black line) and is affected by finite-size effects very marginally. The forward transition is instead driven by stochastic fluctuations and moves to higher values of σas the network size increases. This is in agreement with the fact that the unstable branch of the self-consistent solution (black dashed line) does not cross the x-axis for any finite value of the coupling constant σ. COMMUNICATIONS PHYSICS | https://doi.org/10.1038/s42005-021-00605-4 ARTICLE COMMUNICATIONS PHYSICS | (2021) 4:120 | https://doi.org/10.1038/s42005-021-00605-4 | www.nature.com/commsphys 9