In nature and technology, complex networks serve as a ubiquitous paradigm with a broad range of applications from physics, chemistry, biology, neuroscience, socioeconomic, and other systems [1]. Dynamical networks consist of interacting dynamical units, such as neurons or lasers. Collective behavior in dynamical networks has attracted much attention in recent decades. Depending on the network and the specific dynamical system, various synchronization patterns with increasing complexity were explored [2–5]. Even in simple models of coupled oscillators, patterns such as complete synchronization [6, 7], cluster synchronization [8–11], and various forms of partial synchronization have been found, such as frequency clusters [12], solitary [13–15], or chimera states [16–20]. In particular, synchronization is believed to play a crucial role in brain networks, for example, under normal conditions in the context of cognition and learning [21, 22], and under pathological conditions, such as Parkinson’s disease [23–25], epilepsy [26–29], tinnitus [30, 31], schizophrenia, to name a few [32].

The powerful methodology of master stability function [33] has been a milestone for the analysis of synchronization phenomena. This method allows for the separation of dynamic and structural features in dynamical networks. It greatly simplifies the problem by reducing the dimension and unifying the synchronization study for different networks. Since its introduction, the master stability approach has been extended and refined for various complex systems [34–42], and methods beyond the local stability analysis have been developed [43–47]. More recently, the master stability approach has been extended to another class of oscillator networks with high application potential, namely adaptive networks [48].

Adaptive networks are commonly used models for various systems from nature and technology [49–57]. A prominent example are neuronal networks with spike-timing dependent plasticity, in which the synaptic coupling between neurons changes depending on their relative spiking times [58–61]. There are a large number of studies investigating the dynamic properties induced by this form of synaptic plasticity [62]. However, analysis is usually limited to only one or two forms of spike timing-dependent plasticity within a neuronal population. On the other hand, experimental studies indicate that different forms of spike timing-dependent plasticity may be present within a neuronal population, where the form depends on the connection structure between the axons and dendrites [63]. Among all structural aspects, an important factor for the specific form of the plasticity rule is the distance between neurons [64–66]. More specifically, it has been found that the plasticity rule between proximal or distal neurons, respectively, can change from Hebbian-like to anti-Hebbian-like [67, 68].

This work introduces a methodology to study synchronization in adaptive networks with heterogeneous plasticity (adaptation) rules. As a paradigmatic system, we consider an adaptively coupled phase oscillator network [69–75], which is proven to be useful for predicting and describing phenomena occurring in more realistic and detailed models [76–79]. More specifically, in the spirit of the master stability function approach, we consider the synchronization problem as the interplay between network structure and a heterogeneous adaptation rule arising from distance- (or location-)dependent synaptic plasticity. For a given heterogeneous adaptation rule, our master stability function provides synchronization criteria for any coupling configuration. As illustrative examples, we consider a nonlocally coupled ring with biologically motivated plasticity rule, and a network with a Gaussian distance-dependent coupling weights. We explained such intriguing effects as synchronization or desynchronization by enhancement of long-distance links.

We introduce the model in Section 2. Building on findings from [48], we develop a master stability approach in Section 3 that takes a heterogeneous adaptation rule in account. In Section 4.1, we provide an approximation of the structural eigenvalues that determine the stability of the synchronous state. We then consider two different setups: a nonlocally coupled ring in Section 4.2 and a weighted network with Gaussian distance distribution of coupling weights in Section 4.3. Both systems are equipped with a biologically motivated plasticity rule. In Section 5, we summarize the results.

In this work, we study the synchronization on networks with adaptive coupling weights, where the adaptation (plasticity) rule depends on the distance between oscillators (neurons). We consider the model of adaptively coupled phase oscillators, which has proven to be useful for understanding dynamics in neuronal systems with spike timing-dependent plasticity [77, 79, 48]. The model reads as follows: d d t ϕ i = ω + ∑ j = 1 N a i j κ i j g ( ϕ i − ϕ j ) , (1) d d t κ i j = − ϵ ( κ i j + h i j ( ϕ i − ϕ j ) ) , (2) where ϕ i ∈ S 1 = ℝ / 2 π ℤ ( i = 1 , … , N ) is the phase of the ith oscillator, κ i j ( i , j = 1 , … , N ) is the dynamical coupling weight from oscillator j to i, ω denotes the natural frequency of each oscillator, and a i j ∈ [ 0,1 ] are the entries of the weighted adjacency matrix A describing the network connectivity. The time scales of the “fast” phase oscillators and “slow” coupling weights are separated by the parameter ϵ, which we assume to be small 0 < ϵ ≪ 1 . The functions g and h i j denote the coupling and the N 2 plasticity functions, respectively. For illustrative purposes, the coupling function is set throughout the paper to g ( ϕ ) = − sin ( ϕ + α ) / N with the phase lag parameter α [80]. Such a phase lag can account for a small synaptic propagation delay [81, 48]. For formal derivations, however, a generic coupling function is used. Note that the system Eqs. 1, 2 is shift-symmetric, i.e., invariant under the transformation ϕ i ↦ ϕ i + ψ for any ψ ∈ S 1 . This allows us to restrict our consideration to the case ω = 0 by introducing a new “co-rotating” coordinate system ϕ i , new = ψ i − ω t .

The main difference of system Eqs. 1, 2 from the models considered previously in the literature [40, 70, 71, 74, 82], is that the plasticity functions h i j can be different for each network connection j → i .

A solution to Eqs. 1, 2 is called phase-locked if, for all i = 1 , … , N , the phases evolve as ϕ i = Ω t + ϑ i with some collective frequency Ω ∈ ℝ and ϑ i ∈ S 1 . If ϑ i = ϑ for all i = 1 , … , N , the phase-locked state is called in-phase synchronous or, short, synchronous state.

In the case of in-phase synchronous state, we can set ϑ i = 0 for each oscillator due to the shift symmetry of Eqs. 1, 2. The in-phase synchronous state is given as ϕ s ( t ) = − w g ( 0 ) t , (3) κ i j s = − h i j ( 0 ) , (4) where we assume that the weighted row sum w = ∑ j = 1 N a i j h i j ( 0 ) is constant for all. Such an assumption of constant row sum is necessary for the existence of the synchronous state. Moreover, it is satisfied for commonly considered cases of global or nonlocal shift-invariant coupling.

In the following section, we show how the stability of the synchronous state is determined in a master-stability-like approach.

In Section 2, we have introduced a general class of models and the synchronous state, that are considered throughout this paper. In this section, we derive a framework for the local stability analysis of the synchronous states. We note that the master stability approach for homogeneous adaptations h i j = h was introduced in [48, 83]. Here we extend the methodology to heterogeneous adaptation rules.

To describe the local stability, we introduce the variations ξ i = ϕ i − ϕ s and χ i j = κ i j − κ i j s . The linearized equations for these variations can be written in the following matrix form d d t ( ξ χ ) = J ( ξ χ ) = ( D g ( 0 ) L h g ( 0 ) B − ϵ C − ϵ I N 2 ) ( ξ χ ) , (5) where ξ = ( ξ 1 , … , ξ N ) T is N-dimensional vector containing the perturbations ξ i = δ ϕ i of the phases and χ = ( χ 11 , χ 12 , … , χ N N ) T are N 2 - dimensional vectorized perturbations of coupling weights χ = vec [ δ κ i j ] , respectively. The N × N weighted Laplacian matrix L h has the following elements l i j h = { − ∑ m = 1 , m ≠ i N a i m h i m ( 0 ) , i = j , a i j h i j ( 0 ) , i ≠ j . (6)

The time-independent matrices B and C are B = ( a 1 ⋱ a N ) , C = ( ( D h ) 1 T ⋱ ( D h ) N T ) − ( diag   ( D h ) 1 ⋮ diag   ( D h ) N ) , where a i = ( a i 1 , … , a i N ) , ( D h ) i = ( D h i 1 ( 0 ) , … , D h i N ( 0 ) ) , and diag   ( D h ( 0 ) ) i = ( D h i 1 ( 0 ) ⋱ D h i N ( 0 ) ) .

Note that due to the shift symmetry of Eqs. 1, 2, the Jacobian J in Eq. 5 is time independent. Therefore, the real parts of the N ( N + 1 ) eigenvalues λ of J are the Lyapunov exponents of the synchronous state and hence determine its local stability. In the following proposition, we exploit the fact that J contains a large diagonal block − ϵ I N 2 to reduce the dimension of the eigenvalue problem for J.

Proposition 1. Suppose ϕ i = Ω t is an in-phase synchronous state of Eqs. 1, 2 . Then its linear stability is determined by the 2 N -dimensional linear system d d t v = ( D g ( 0 ) L h g ( 0 ) I N ϵ L D h − ϵ I N ) v , (7) where D g ( 0 ) and L h are as in Eq. 5 and the N × N weighted Laplacian matrix L D h possesses the following elements l i j D h = { − ∑ m = 1 , m ≠ i N a i m D h i m ( 0 ) , i = j , a i j D h i j ( 0 ) , i ≠ j . (8)

Proof.

We remind that system Eq. 5 determines the spectrum (Lyapunov exponents) of the synchronous state. The Jacobian matrix in Eq. 5 is sparse with a large N 2 × N 2 block given by the simple diagonal matrix − ϵ I N 2 . This implies that Eq. 5 possess N 2 − N stable directions with Lyapunov exponents − ϵ . To find these directions, we substitute ( ξ , χ ) = e − ϵ t ( ξ 0 , χ 0 ) into Eq. 5 and obtain the linear system ( D g ( 0 ) L h + ϵ I N g ( 0 ) B − ϵ C 0 ) ( ξ 0 χ 0 ) = 0. (9)

This system has at least N 2 − N linearly independent solutions, since the matrix in Eq. 9 is degenerate due to the large N 2 × N 2 zero block. The structure of the invariant subspaces in system Eq. 5 allows for introducing new coordinates, which separate the N 2 − N stable directions (corresponding to the eigenvalues − ϵ ) from the remaining 2 N directions. Explicitly, this transformation is given by ( ξ χ ) = R ( ξ χ ^ ) ,   R = ( I N 0 0 0 ( 1 / r ) B T K ) with ( N 2 + N ) × ( N 2 + N ) matrix R. Here K is an ( N 2 − N ) × ( N 2 − N ) orthogonal matrix with B K = 0 . Applying this transformation, we obtain the following system d d t ( ξ χ N χ N 2 − N ) = ( D g ( 0 ) L h g ( 0 ) I N 0 ϵ L D h − ϵ I N 0 − ϵ K T C 0 − ϵ I N 2 − N ) ( ξ χ N χ N 2 − N ) , (10) where ( ξ , χ N , χ N 2 − N ) T = ( ξ , χ ^ ) T , with χ N and χ N 2 − N are an N and N 2 − N -dimensional vectors, respectively, and the N × N weighted Laplacian matrix L D h as given in Eq. 8. For more details on the transformation, we refer the reader to [48, 83]. We observe that the variables ( ξ , χ N ) are independent on χ N 2 − N . Hence, separating the master from the slave system, the resulting coupled differential equations that determine the stability of the synchronous state are given by system Eq. 7. This concludes the proof.

 Proposition 1 reduces the problem’s dimension significantly from N ( N + 1 ) to 2 N . In the spirit of the master stability approach [33], we aim for further decomposition of the 2 N - dimensional coupled system Eq. 7 into dynamically independent blocks of dimension 2. For this, we restrict our consideration to the case when L h can be diagonalized S h = Q − 1 L h Q by a nonsingular complex-valued matrix Q. Note that the eigenvalues μ i of L h lie on the diagonal of S h . In general, the matrices L h and L D h do not commute. Therefore, Q − 1 L D h Q is not necessarily of upper triangular shape. Regardless of this fact, the following proposition provides an explicit form for the eigenvalues of J in Eq. 5 in the limit of slow adaptation, i.e., ϵ ≪ 1 .

Proposition 2. Assume that L h is diagonalizable, with S h = Q − 1 L h Q being the associated diagonal matrix and Q the corresponding transformation. Let ϕ i = Ω t be an in-phase synchronous state of Eqs. 1, 2 Then, the local stability of this state is determined by the solutions of N quadratic equations, which are given up to the first order in ϵ as λ 2 + ( ϵ − D g ( 0 ) μ i ) λ − ϵ ( D g ( 0 ) μ i + g ( 0 ) ν i ) = 0 ,   i = 1 , … , N , (11) where μ i are the eigenvalues of L h located on the diagonal of S h and ν i are the corresponding diagonal elements of Q − 1 L D h Q . If L h and L D h commute, then Eq. 11 is exact, and ν i are the eigenvalues of L D h .

Proof. Due to Proposition 1, the eigenvalues of the Jacobian in Eq. 5 are given by det ( D g ( 0 ) L h − λ I N g ( 0 ) I N ϵ L D h − ( ϵ + λ ) I N ) = det ( D g ( 0 ) S h − λ I N g ( 0 ) I N ϵ Q − 1 L D h Q − ( ϵ + λ ) I N ) = 0 , where we have used the transformation Q that brings L h to the diagonal form S h = Q − 1 L h Q . Making further use of the Schur complement [84], we obtain det ( D g ( 0 ) S h − λ I N g ( 0 ) I N ϵ Q − 1 L D h Q − ( ϵ + λ ) I N ) = det ( ( λ + ϵ ) ( λ I N − D g ( 0 ) S h ) − ϵ g ( 0 ) Q − 1 L D h Q ) = 0. (12)

The latter equation is almost diagonal. The only off-diagonal components remain from Q − 1 L D h Q and scale with ϵ. Let us consider the Leibniz formula for the determinant of an N × N matrix F with entries f i j , that reads det ( F ) = ∑ σ ∈ Perm ( N ) sgn ( σ ) ∏ i = 1 N f i σ ( i ) . In the latter expression Perm ( N ) denotes the set of all permutations σ of the integer numbers 1 , … , N and sgn ( σ ) ∈ { − 1,1 } is the sign of the permutation. Since all off-diagonal terms of the matrix considered in Eq. 12 scale with ϵ, for any but the identical permutation each term ∏ i = 1 N f i σ ( i ) scales with ϵ 2 or higher. Hence, we are left with det ( F ) = ∏ i = 1 N f i i + ∏ ( ϵ 2 ) and find det ( ( λ + ϵ ) ( λ I N − D g ( 0 ) S h ) − ϵ g ( 0 ) Q − 1 L D h Q ) = ∏ i = 1 ( λ 2 + ( ϵ − D g ( 0 ) μ i ) λ − ϵ ( D g ( 0 ) μ i + g ( 0 ) ν i ) ) + ∏ ( ϵ 2 ) = 0 , (13) where μ i are the eigenvalues of L h , ν i are the diagonal elements of Q − 1 L D h Q and ∏ ( ϵ 2 ) denotes higher order terms ( ϵ m , m > 1 ). If L h and L D h commute, both matrices share the same set of eigenvectors and hence they can be brought to the diagonal form with the same transformation Q. In this case, the diagonal elements ν i are the eigenvalues of L D h and the higher order terms ∏ ( ϵ 2 ) in Eq. 13 vanish.

The 2 N solutions λ i of the N Eq. 11 determine the stability of the synchronous state. More precisely, the real parts of theses solutions determine the Lyapunov exponents. If Λ = max i Re ( λ i ) < 0 , then the synchronous state is locally stable, while for Λ > 0 it is locally unstable. The case Λ = 0 provides the stability boundary.

Note that for a fixed time scale parameter ϵ ≪ 1 , the Eq. 11 and hence its solutions depend on the coupling function g, the connectivity, and the adaptation structure. This dependence, however, is only encoded in the two complex parameters D g ( 0 ) μ and g ( 0 ) ν . Therefore, we define the master stability function Λ : ℂ 2 → ℝ with Λ ( D g ( 0 ) μ , g ( 0 ) ν ) = max i Re ( λ i ( D g ( 0 ) μ , g ( 0 ) ν ) ) that maps each pair of parameters ( D g ( 0 ) μ , g ( 0 ) ν ) to the corresponding Lyapunov exponent.

For an illustration, we consider a cross-section of ( D g ( 0 ) μ , g ( 0 ) ν ) - space by setting Im ( μ ) = 0 and Im ( ν ) = 0 . This cross-section is of particular interest in cases of symmetric matrices L h and L D h since their eigenvalues are real. In Figure 1, we present the master stability function for the coupling function g ( ϕ ) = − sin ( ϕ + α ) / N and different values of the parameter α. In case of real μ and ν, we obtain two explicit stability conditions from Eq. 11: The synchronous state is locally stable ( Λ < 0 ) if c 1 ( α , μ ) = cos ( α ) μ > − ϵ , (14) c 2 ( α , μ , ν ) = cos ( α ) μ + sin ( α ) ν > 0. (15) These conditions agree with the black dashed lines in Figure 1 and are used subsequently to describe stability for certain network models.

In the previous section, we established a generic analytic tool for studying stability of synchronous states. In this section, we focus on the application of the tool to certain network models. For the rest of the work, we restrict our attention to the following generalization of the Kuramoto-Sakaguchi system with distance-dependent synaptic plasticity d d t ϕ i = ω − 1 N ∑ j = 1 N a i j κ i j sin ( ϕ i − ϕ j + α ) , (16) d d t κ i j = − ϵ ( κ i j + h ( ϕ i − ϕ j , d i j ) ) . (17)

The plasticity function h depends on the phase difference ϕ i − ϕ j and on the distance d i j . In this work, we associate the distance to the difference of indices by d i j = | j − i | . For the plasticity function, we consider h i j ( ϕ ) = h ( ϕ , d i j N ) = { h ^ ( ϕ , d i j N ) d i j ≤ N / 2 , h ^ ( ϕ , 1 − d i j N ) d i j > N / 2. (18)

With this form of the adaptation function, we have a symmetric h i j ( ϕ ) = h j i ( ϕ ) and a circulant h i + l , j + l ( ϕ ) = h i j ( ϕ ) structure of the corresponding matrix with entries h i j . Particularly, for the numerical analysis, we use h ^ ( ϕ , d i j / N ) = sin ( ϕ + β ( d i j / N ) ) , (19) where the distance dependence is encoded in the phase shift function β ( d i j N ) = { ( 2 N d i j − 1 ) π , N  even , ( 2 ( N + 1 ) d i j − 1 ) π , N  odd . (20)

In Figure 2A, we illustrate the distance-dependent plasticity function Eqs. 18–20 for a network of N = 12 nodes. The illustration shows the different plasticity functions depending on the distance between the nodes d i j . The plasticity function changes from a Hebbian to anti-Hebbian rule for proximal and distal node, respectively. This change, particularly in the proximity of ϕ = 0 , is in qualitative agreement with the experimental findings in [67]. Note the symmetry of the plasticity function that renders the matrix with elements h i j circulant.

If not indicated differently, we consider the coupling structure given by a i j = a ( d i j / N ) , (21) where a : [ 0,1 ] → [ 0,1 ] is a bounded and piece-wise continuous function. This corresponds to a distant-dependent coupling, and it results to a dihedral symmetry in the coupling structure (ring-like).

In the following section, we provide an approximation for the eigenvalues of L h and L Dh for large networks with circulant connectivity and plasticity structure. Using this approximation, we subsequently analyze the stability of the synchronous state on nonlocally coupled networks and on isotropic networks with Gaussian weight distribution.

In summary, we have investigated the phenomenon of synchronization on adaptive networks with heterogeneous plasticity rules. In particular, we have modeled systems with distance-dependent plasticity as they have been found in neuronal networks experimentally [64–67] as well as computational models [68]. For the realization, we have used a ring-like network architecture and associated the distance of two nodes with the distance of their placement on the ring.

In Section 3, we have developed a generalized master stability approach for phase oscillator models that are adaptively coupled and where each link has its own adaptation rule (plasticity). By using an explicit splitting of the time scales between fast dynamics of the phase oscillators and slow dynamics of the link weights, we have established an explicit stability condition for the synchronous state. More precisely, we found that the stability is governed by the coupling function and the eigenvalues of two structure matrices. These structure matrices L h and L D h are determined by the connectivity of the network and the plasticity rules of the link weights. Note that for the structural matrices, the plasticity rule needs only to be known in the vicinity of 0, which greatly facilitates the application of the approach to realistic forms of synaptic plasticity. Thus, we have extended previous work on the master stability function of adaptive networks [48, 83] and broaden the scope of potential future applications for this methodology.

In Section 4, we applied the novel technique to a system of adaptively coupled oscillators with distance-dependent plasticity. Here, we have used a ring-like network structure to study the impact of long- and short-distance connections on the stability of synchronization. For this purpose we introduced an approximation of the eigenvalues for the structure matrices in Section 4.1. This approximation allows for a comprehensive analysis of the stability as a function of various system parameters. Moreover, it enables us to identify critical eigenvalues that govern the stability of the synchronous state. In Sections 4.2 and 4.3, we have brought together all methodological findings and applied them to systems with a nonlocally coupled ring structure and with a Gaussian distribution of link weights. The latter structure accounts for the fact that in realistic neuronal populations the number of links with different distances are not uniformly distributed [67]. We found that long-distance connections can stabilize or destabilize the synchronous state, depending on the coupling function between the oscillators. A remarkable fact with respect to neuronal applications relates to the destabilization scenario. Here we observed that the destabilization can be attributed to the pronounced change of the plasticity rule from Hebbian to anti-Hebbian. For more realistic connectivity structures, we found that weight distributions of the connectivity structure with sharp peaks at long distances lead to destabilization for a wide range of the coupling function.

All in all, in this article, we have provided a general framework to study the emergence of synchronization in neuronal system with a heterogeneous plasticity rule. The developed methodology is not limited to distance-dependent types of plasticity and can also be used for non-symmetric setups. For the latter case, we have provided the necessary analytical result. In this work, we have restricted our attention to the case of phase oscillators, but the methods can be extended to more realistic neuron models by using techniques established, for example, in [48]. Moreover, techniques are available that allow for further generalization toward systems with slightly different local dynamics at each node [87]. On the one hand, the master stability approach offers a great tool to study the stability of the synchronous state depending on the networks structure. On the other hand, this approach allows for characterizing the network structures that are, in some sense, optimal for synchronization [88, 89]. In this regard, it remains an open question as to how plasticity optimizes the synchronizability of the network in a self-organized way. In addition, recent studies have shown that there is a great interest in synchronization phenomena to understand diseases such as Parkinson’s disease [90–92] or epilepsy [29, 93] for the development of proper therapeutic treatments. We believe that our work provides an important step toward understanding synchronization under realistic conditions.