Hebbian plasticity rules abrupt desynchronization in pure simplicial complexes

This Letter investigates the upshots of adaptive development of pure 2- and 3- simplicial complexes (triad and tetrad) on the nature of the transition to desynchrony of the oscillator ensembles. The adaptation exercised in the pure simplicial coupling takes a cue from the Hebbian learning rule, i.e., the coupling weight of a triad (tetrad) is prone to increase if the oscillators forming it are in phase and decrease if they are out of phase. The coupling weights in these pure simplicial complexes experiencing such adaptation give rise to first-order routes to desynchronization, whose onsets are entirely characterized by respective Hebbian learning parameters. Mean-field analyses presented for the order parameters for the adaptive 2- and 3- simplicial complexes strongly corroborate with the respective numerical assessments.

For decades, consideration of the pairwise interaction between different networked dynamical units has been at the forefront to capture the underlying dynamics affecting distinct dynamical processes on various physical and biological complex systems. However, complex systems such as brains [1,2] and social interaction networks [3,4,5,6] have the underlying topology of higher-order connections. The higher-order interactions are often framed using simplicial complexes [7,8] that are topological structures. A simplicial complex may include simplices of different dimensions, namely, vertices (0-simplex), edges (1-simplex), triangles (2-simplex), tetrahedrons (3-simplex) and so on stuck to each other. A pure simplicial complex is where all facets have the same dimension. Simplicial complexes composed of geometrical objects of different dimensions offer a suitable framework to capture the underlying geometry in complex systems. For instance, simplicial complexes have been used to model the topological map encoded by the hippocampus to capture the environment's basic geometrical features [9,10]. Recently modeling many-body interactions using simplicial complexes is gaining momentum divulging the rich interplay between network geometry and dynamical processes [11,12,13,14,15,16,17,18,19,20,21,22,23,24].
The adaptation is at the backbone of the construction and functioning of many physical and biological complex systems. Experimental findings in neuroscience have led to an impression that the synaptic plasticity between the interacting neurons is the rationale of the learning process and long-term memory in the human brain [25,26]. Hebb [27] first conceptualized the fact that the synaptic strength between two interacting neurons is strengthened (weakened) if they are simultaneously firing (not firing), which was later supported by experimental evidence [26,28,29]. Further, a neural network can be represented by an ensemble of phase oscillators by encoding the relative spike timing of presynaptic and postsynaptic spikes of the interacting neurons in terms of phases. Thus, the plasticity of the neural network materializes through the adaptation of the synaptic strength between two interacting neurons. The systems of phase oscillators incorporating plasticity in the connection strength between interacting oscillators have revealed intriguing structures and phenomena, for instance, the existence of cluster states [30,31,32,33,34,35] or mesoscale structures [36,37], explosive synchronization arising from anti-Hebbian adaptation in monolayer networks [38] and interlayer Hebbian adaptation in multiplex networks [39]. Also an adaptive simplicial complex model provides working description of the spatial learning process in the Hippocampal [40]. The collective dynamics of man-made and natural complex systems represented by conventional graphs, for long, have been modeled using pairwise networked Kuramoto oscillators [41,42]. When subjected to distinct adaptation-based features, the pairwise networked Kuramoto oscillators have led to a phenomenon termed first-order or explosive synchronization. A first-order transition sports an abrupt jump to the coherence and then an irreversible abrupt collapse to the incoherence as the coupling strength between the interacting oscillators varies [38,43,44,45,46,47,48,49,50,51,52,53,54]. In this Letter, breaking away from the traditional approach of assuming adaptation in pairwise coupling, we consider adaptive higher-order coupling interactions among the oscillators. The dynamics of Kuramoto oscillators in the two models studied are respectively driven by pure triadic and pure tetradic couplings that adapt abiding by the Hebbian learning mechanism [39], i.e., a triadic (tetradic) coupling strengthens when the oscillators forming the triad (tetrad) are in phase and weakens when they are out of phase. Such dynamically adaptive pure triadic and tetradic couplings give birth to first-order routes to desynchronization, whose onsets are entirely manageable through respective Hebbian learning parameters. Mean-field analyses presented for the 2-and 3- when B ijk is static (a-b) and adaptive (c-d) with ε=1 and µ=1. The Hebbian plasticity negates the impact of initial asymmetry η on the outset of desynchronization, instead, makes its impact more profound at the degree of synchronization.
simplicial complexes successfully explain the dependence of the onset points of abrupt desynchronization on the respective Hebbian learning parameters.
To develop an extensive understanding of the repercussions of the Hebbian learninginspired adaptation in higher-order couplings on synchronization dynamics of oscillator ensembles, we numerically and analytically explore the synchronization dynamics on the pure 2-and 3-simplicial complexes.
Adaptive 2-simplex (triadic) interaction: We consider a phase ensemble of Kuramoto oscillators subjected to adaptive higher-order interactions among them. The higher-order coupling exercised among the oscillators is triadic (three-way) encoded by a 2-simplex (triangle) structure, whose weight dynamically adapt following a mechanism similar to the Hebbian learning rule for pair-wise interaction [39]. The dynamical evolution of a triadic weight assigned to a triplet of nodes {i, j, k} having instantaneous phases θ i,j,k is governed by the Hebbian-learning inspired adaptatioṅ where constants ε ∈ [0, 1] and µ > 0 are learning enhancement rate and saturation rate, respectively. The cosine term implies that the coupling weight B ijk tends to increase if triadic oscillators {i, j, k} are in phase, while it tends to decrease if they are out of phase. The linear saturating term µB ijk refrains the coupling weight from increasing or decreasing without any bound. The phase evolution of Kuramoto oscillators under the impression of such triadic interaction is ruled bẏ where natural frequency ω i (i=1, . . . , N ) of each oscillator i is drawn randomly from a uniform or unimodal distribution g(ω), and λ is a uniform coupling strength among Hebbian plasticity rules abrupt desynchronization in pure simplicial complexes  Figure 3. (Color online) Impact of (a) learning rate ε for d [2] =12 and (b) mean 2-simplex degree d [2] on transition in a random 2-simplicial complex of N =10 3 , µ=1 and η=1. Slower rate ε requires relatively higher critical λ for the system to desynchronize. A higher mean connectivity enhances the level of synchronization and requires relatively lower critical λ to desynchronize.
the interacting nodes. Here, B is an adjacency tensor of 2-simplex that encodes the topology of the network, i.e., B ijk =1 if there is a triadic coupling among the nodes {i, j, k}, otherwise B ijk =0. The number of triangles in the complex a node is part of, is defined as 2-simplex degree d [2] denotes the mean 2-simplex degree.
Eqs. (1) & (2) collectively determine the dynamics of the adaptive 2-simplicial complex. To capture the degree of synchronization of the oscillator ensemble, we define generalized order parameter as R m e iΨm = 1 N N j=1 e imθ j , where m=1, 2. R m and Ψ m are amplitude and argument, respectively, corresponding to m-cluster order parameter that captures the formation of m-clusters. R 1 quantifies one-cluster synchronization whereas R 2 quantifies one-cluster and two-cluster synchronization [32,55]. Hence, to measure only two-cluster (anti-phase) synchronization one must adjust R 2 for R 1 , i.e., |R 2 − R 1 |.
Next, we numerically compute the order parameters and discuss the microscopic dynamics of the triadic weights. First, we construct a pure 2-simplex structure by pinpointing each distinct triangle from the 1-simplex (Erdös-Rényi random graph) structures. The natural frequencies of the nodes are drawn uniform randomly, i.e., ω i ∼U (−∆, ∆), where ∆=0.5 or otherwise mentioned elsewhere. For a large λ, a uniformly selected fraction η of the nodes are assigned initial phases to θ i (0)=0 and the rest (1 − η) assigned to θ i (0)=π. The initial weights B ijk (0) are then determined by The initial phase asymmetry (η) does affect the degree of synchronization for R 1 . The maximum value to which R 1 plateaus (at large λ) decreases with the increase in the level of asymmetry, .i.e., η=0.5 (maximum asymmetry) yields R 1 0 while η=1 (symmetry) yields the largest possible R 1 . However, R 2 remains independent of the initial phase asymmetry. Moreover, both R 1 and R 2 do not synchronize with the forward increase in λ, hence the incoherent state persists for any λ > 0. These behavior of R 1 and R 2 also persist for unimodal g(ω) however at different critical coupling points (see Supplemental Material [56]). Note that the Hebbian plasticity of dyadic (1-simplex) interactions does not lead to an abrupt desynchronization (see Supplemental Material [56]).
Impact of learning rate and mean 2-simplex degree: Additionally, we demonstrate the effect of learning rate ε and mean 2-simplex degree d [2] on the nature of desynchronization in Fig. 3. As the rate ε is gradually decreased, the onset of firstorder desynchronization transition occurs at higher values of critical coupling point λ c [see Fig. 3(a)]. Thus, the rate ε provides control over the onset of desynchronization. Moreover, the mean triadic connectivity d [2] also affect the jump height of the abrupt desynchronization [see Fig. 3(b)]. As d [2] is increased, more and more triadic interactions are involved leading to an increment in the height of abrupt jump at rather lower value of critical coupling point.
Stationary triadic weights and phases: The microscopic dynamics of stationary weights and stationary phases during the backward sweep of λ are presented in Fig. 4 for different values of η. For η=1, stationary weights are clustered around B ijk →ε/µ at a large λ. As λ is gradually decreased, a gradually increasing number of B ijk starts departing away from ε/µ towards B ijk =0. At the critical point λ c , all B ijk get uniformly scattered between −ε/µ and ε/µ and remain scattered so in the incoherent state until λ=0 is reached. For η=1, the stationary phases corresponding to λ > λ c remain locked to form a single peak distribution. For the case of η=0.5, two almost equal-half populations of steady B ijk are flocked around −ε/µ and ε/µ, respectively. A gradual decrease in λ shifts a gradually increasing fraction of B ijk towards 0. At the critical point, the entire B ijk population now gets scattered between −ε/µ and ε/µ. The steady phases related to λ > λ c remain almost equally divided into two-clusters at a phase-difference of π.
Impact of saturation rate µ: The impact of parameter µ on stationary triadic weights is shown in Fig. 5. The ratio ε/µ determines the lower and upper bounds for the stationary weights B ijk , which can be corroborated from Eq. (1) in the steady state. Parameter 0<µ≤1 [ Fig. 5(a)] yields higher weight bounds which causes desynchronization at an early coupling strength. On the other hand, µ>1 [ Fig. 5(b)] yields lower weight bounds which leads to the desynchronization at rather higher coupling strength. Therefore, µ also serves as a control parameter in determining the onset of abrupt desynchronization.
Analytical discussion: To gain theoretical insight of the underlying dynamics, we turn our focus to an all-to-all connected 2-simplex structure modeled aṡ The system can achieve its stationary state only when both triadic weights and phases achieve their respective stationary states. The stationary triadic weights,Ḃ ijk =0, are given by B ijk =(ε/µ) cos(θ j + θ k − 2θ i ). Stationary B ijk are set apart into two weight populations, one of B ijk →ε/µ corresponding to (θ j +θ k −2θ i )=0 and other of B ijk →−ε/µ related to (θ j +θ k −2θ i )=π. For that matter, an initial asymmetric population comprising 0 or π phases with a probability η and 1 − η, respectively, is taken into account. The Eq. (3) can be re-expressed in terms of mean-field parameters by utilizing the expressions for stationary B ijk and R 2 aṡ Since phases in a frame rotating with the average intrinsic frequencyΨ m are given by θ i = θ i +Ψ m t. Nevertheless, considering a non-rotating frame withΨ m =0 and setting the average phase Ψ m =0 by selecting appropriate initial phases, the fixed points in the synchronous stateθ i = 0 (|ω i | ≤ bq) are determined by where η and (1−η) denote a constant probability [η(ω i ) = η] with which the initial phases are set to 0 and π, respectively. Hence, the order parameters for the locked oscillators in continuum limit can be expressed as Further mathematical simplifications supply the real part of the order parameters as After plugging in the expressions of cos(θ * ) and cos(2θ * ) from Eq. (5) and solving for a uniform distribution g(ω) = 1 2∆ , Eqs. (7) furnish us with (also see Supplemental Material [56]) It is apparent from Eq. (8) that only order parameter R 1 is impacted by the change in the initial phase asymmetry η while R 2 remains independent of it. The backward end points of R 1 − λ and R 2 − λ traces are obtained by setting ∆ = bq = bλ c R 2 2c : Hebbian plasticity rules abrupt desynchronization in pure simplicial complexes It is apparent that the critical coupling point λ c for the onset of desynchronization is the function of Hebbian plasticity rates ε and µ. The analytical predictions [Eq. (8)] for R 1 and R 2 are plotted in Fig. 6 for different values of ε and η. They are fairly in good agreement with their respective numerical estimations. The values of critical points λ c and R 1c and R 2c are also shown matching their respective numerical assessments.
Adaptive 3-simplex (tetradic) interaction: Now we consider the higher-order coupling interaction among the oscillators to be tetradic encoded by a pure 3-simplex (tetrahedron) structure. The weights of the tetradic coupling also dynamically adapt in a fashion similar to the Hebbian learning rule. The evolution of the phases and the tetradic weights can respectively be expressed aṡ where C is 3-simplex adjacency tensor encoding a network of tetrahedrons. C ijkl =1 if a tetrahedronic connection exists among the nodes {i, j, k, l}, otherwise C ijkl =0. 3simplex degree of a node d [3] i = 1 3! N j,k,l=1 C ijkl denotes the number of tetrahedrons it is part of. We consider the 1-simplex structure (Erdös-Rényi random graph) and then identify each distinct tetrahedron to construct a pure 3-simplicial complex.
The evolution of the system begins at a large λ with initial phases fixed to 0 and π with a constant probability η and 1 − η, respectively. The initial tetradic weights are assigned as C ijkl (0) fixed λ c for different values of η. R 2 remains independent of η whereas R 1 decreases with decrease in the value of η. Also, both R 1 and R 2 always exhibit incoherence for forward continuation in λ be it either static or adaptive weight C ijkl . Now we look at the nature of stationary weightsĊ ijkl =0, C ijkl = (σ/ν) cos(θ j + θ k + θ l − 3θ i ). Fig. 8 unveils that all C ijkl in the synchronous state are segregated into two clusters, one of C ijkl →σ/ν related to (θ j + θ k + θ l − 3θ i )→0 and other of C ijkl → − σ/ν related to (θ j +θ k +θ l −3θ i )→π. The size of two anti-phase clusters depends on parameter η. In the asynchronous state, all C ijkl are uniformly scattered between −σ/ν and σ/ν.
For the sake of analytical simplification we consider an all-to-all connected pure 3-simplicial complex, for which the interaction term in the phase evolution [Eq. (10)] is normalized by N 3 instead of 6 d [3] . Further, substituting C ijkl with the stationary C ijkl , the phase evolution then can be rewritten in terms of R 2 and Ψ 2 .
Hence the critical coupling points for the outset of desynchronization is controlled by the Hebbian learning parameters σ and ν. The theoretical predictions for R 1 and R 2 [Eq. (12)] and their numerical estimations fall in fine concurrence with each other for different values of σ and η (see Fig. 9), and so do the hypothesized critical points R 1c , R 2c and λ c .

Conclusion:
We studied desynchronization transition in oscillator ensembles networked in triadic (pure 2-simplex) and tetradic (pure 3-simplex) structures with adaptive coupling. The adaptation exercised in triadic (tetradic) coupling is inspired by the Hebbian learning mechanism, i.e., a triadic (tetradic) weight tends to increase when the nodes forming the triad (tetrad) are in phase and decrease when they are out of phase. The dynamical evolution of such Hebbian learned triadic (tetradic) coupling leads to a first-order desynchronization in the network. For both pure 2-and 3-simplicial complexes, respective triggering points of the emergent abrupt desynchronization are manageable through respective Hebbian learning parameters, i.e., the learning rate and saturation rate. Besides, for both the complexes, the initial phase asymmetry characterizes the one-cluster order parameter while the two-clusters order parameter remains independent of it. Also, the triadic (tetradic) weights remain uniformly distributed in the incoherent state, whereas they form one-cluster or two-clusters depending on the initial phase asymmetry. The mean-field analyses presented for pure 2and 3-simplicial complexes too successfully explain the respective Hebbian adaptation governed nature of one-cluster and two-clusters order parameters. So far, all the investigations featuring the Hebbian learning mechanism have revolved around employing a pair-wise interaction between the oscillators. The present study would be advantageous in better understanding the role of adaptive higher-order neural synaptic connections in the learning process and procuring long term memory in the brain. The proposed models also open an avenue for investigating the imprint of concurrent adaptation in various combinations of dyadic, triadic and tetradic couplings on the synchronization dynamics of oscillator ensembles. This Supplemental Material includes auxiliary figures and mathematical steps, which support the results and arguments described in the main text.

A. Hebbian plasticity of dyadic interaction
The phase and weight evolutions of an all-to-all connected 1-simplex (dyadic) network of N Kuramoto oscillators, with its couplings subjected to the Hebbian plasticity (learning rate α and saturation rate κ), are given bẏ respectively. At a large λ, the oscillators are uniform randomly assigned phases θ i (0)=0 and π with probability η and 1 − η, respectively. The initial coupling weights are then determined by A ij (0) = α cos[θ j (0) − θ i (0)].   S1 illustrates that both R 1 and R 2 desynchronize espousing a second-order transition route. Nonetheless, the maximum value to which R 1 plateaus in the synchronous state is characterized by the asymmetry probability η. Moreover, R 1 does not see synchronization for the forward continuation of λ.

B. Robust desynchronization with unimodal frequency distributions g(ω)
A manifestation of the first-order desynchronization in a random pure 2-simplex and 3-simplex networks corresponding to a Gaussian g(ω)= 1 ∆ √ 2π exp(−ω 2 /2∆ 2 ) and a Lorentzian g(ω)= ∆ π(ω 2 +∆ 2 ) , where ∆ = half width at half maximum, are shown in Fig. S2 and Fig. S3, respectively. The abrupt desynchronization takes place at different critical coupling strength for the two symmetric distributions and forward transition to synchronization does not occur for any λ > 0.

C. Order parameters for the adaptive 2-simplicial complex
The expressions of the order parameters for the 2-simplicial complex treated for the uniform distribution g(ω)= 1 2∆ : where b = ε/2µ and q = λR 2 2 .