Synchronization in networks with heterogeneous adaptation rules and applications to distance-dependent synaptic plasticity

This work introduces a methodology for studying synchronization in adaptive networks with heterogeneous plasticity (adaptation) rules. As a paradigmatic model, we consider a network of adaptively coupled phase oscillators with distance-dependent adaptations. For this system, we extend the master stability function approach to adaptive networks with heterogeneous adaptation. Our method allows for separating the contributions of network structure, local node dynamics, and heterogeneous adaptation in determining synchronization. Utilizing our proposed methodology, we explain mechanisms leading to synchronization or desynchronization by enhanced long-range connections in nonlocally coupled ring networks and networks with Gaussian distance-dependent coupling weights equipped with a biologically motivated plasticity rule.


INTRODUCTION
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,3,4,5]. Even in simple models of coupled oscillators, patterns such as complete synchronization [6,7], cluster synchronization [8,9,10,11], and various forms of partial synchronization have been found, such as frequency clusters [12], solitary [13,14,15], or chimera states [16,17,18,19,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,24,25], epilepsy [26,27,28,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,35,36,37,38,39,40,41,42], and methods beyond the local stability analysis have been developed [43,44,45,46,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,50,51,52,53,54,55,56,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,59,60,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,65,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,70,71,72,73,74,75], which is proven to be useful for predicting and describing phenomena occurring in more realistic and detailed models [76,77,78,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.

MODEL
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: where φ i ∈ S 1 = R/2πZ (i = 1, . . . , N ) is the phase of the ith oscillator, κ ij (i, j = 1, . . . , N ) is the dynamical coupling weight from oscillator j to i, ω denotes the natural frequency of each oscillator, and a ij ∈ [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 ij 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,82]. For formal derivations, however, a generic coupling function is used. Note that the system (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" The main difference of system (1)-(2) from the models considered previously in the literature [70,71,40,83,74], is that the plasticity functions h ij can be different for each network connection j → i.
A solution to (1)-(2) is called phase-locked if, for all i = 1, . . . , N , the phases evolve as φ i = Ωt + ϑ i with some collective frequency Ω ∈ R 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 (1)- (2). The in-phase synchronous state is given as where we assume that the weighted row sum w = N j=1 a ij h ij (0) is constant for all i = 1 . . . , N . 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 masterstability-like approach.

MASTER STABILITY 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 ij = h was introduced in [48,84]. Here we extend the methodology to heterogeneous adaptation rules.
To describe the local stability, we introduce the variations ξ i = φ i − φ s and χ ij = κ ij − κ s ij . The linearized equations for these variations can be written in the following matrix form d dt 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 [δκ ij ], respectively. The N × N weighted Laplacian matrix L h has the following elements The time-independent matrices B and C are Note that due to the shift symmetry of (1)-(2), the Jacobian J in (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 (1)- (2). Then its linear stability is determined by the 2N -dimensional linear system where Dg(0) and L h are as in (5) and the N × N weighted Laplacian matrix L Dh possesses the following elements PROOF. We remind that system (5) determines the spectrum (Lyapunov exponents) of the synchronous state. The Jacobian matrix in (5) is sparse with a large N 2 × N 2 block given by the simple diagonal matrix − I N 2 . This implies that (5) possess N 2 − N stable directions with Lyapunov exponents − . To find these directions, we substitute (ξ, χ) = e − t (ξ 0 , χ 0 ) into (5) and obtain the linear system This system has at least N 2 − N linearly independent solutions, since the matrix in (9) is degenerate due to the large N 2 × N 2 zero block. The structure of the invariant subspaces in system (5) allows for introducing new coordinates, which separate the N 2 − N stable directions (corresponding to the eigenvalues − ) from the remaining 2N directions. Explicitly, this transformation is given by Applying this transformation, we obtain the following system where respectively, and the N × N weighted Laplacian matrix L Dh as given in (8). For more details on the transformation, we refer the reader to [48,84]. 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 (7). This concludes the proof.
Proposition 1 reduces the problem's dimension significantly from N (N + 1) to 2N . In the spirit of the master stability approach [33], we aim for further decomposition of the 2N -dimensional coupled system (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 Dh do not commute. Therefore, Q −1 L Dh 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 (5) in the limit of slow adaptation, i.e., 1.
being the associated diagonal matrix and Q the corresponding transformation. Let φ i = Ωt be an in-phase synchronous state of (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 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 Dh Q.
If L h and L Dh commute, then Eq. (11) is exact, and ν i are the eigenvalues of L Dh .

Frontiers
PROOF. Due to Proposition 1, the eigenvalues of the Jacobian in (5) are given by Making further use of the Schur complement [85], we obtain The latter equation is almost diagonal. The only off-diagonal components remain from Q −1 L Dh Q and scale with . Let us consider the Leibniz formula for the determinant of an 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 (12) scale with , for any but the identical permutation (13) where µ i are the eigenvalues of L h , ν i are the diagonal elements of Q −1 L Dh Q and O( 2 ) denotes higher order terms ( m , m > 1). If L h and L Dh 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 Dh and the higher order terms O( 2 ) in (13) vanish.
The 2N solutions λ i of the N equations (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 equation (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 Dg(0)µ and g(0)ν. Therefore, we define the master stability function Λ : ) that maps each pair of parameters (Dg(0)µ, g(0)ν) to the corresponding Lyapunov exponent.
For an illustration, we consider a cross-section of (Dg(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 Dh 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 c 2 (α, µ, ν) = cos(α)µ + sin(α)ν > 0.
These conditions agree with the black dashed lines in figure 1 and are used subsequently to describe stability for certain network models.

SYNCHRONIZATION ON NETWORKS WITH DISTANCE DEPENDENT PLASTICITY
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 The plasticity function h depends on the phase difference φ i − φ j and on the distance d ij . In this work, we associate the distance to the difference of indices by d ij = |j − i|. For the plasticity function, we consider With this form of the adaptation function, we have a symmetric h ij (φ) = h ji (φ) and a circulant h i+l,j+l (φ) = h ij (φ) structure of the corresponding matrix with entries h ij . Particularly, for the numerical analysis, we useĥ where the distance dependence is encoded in the phase shift function In figure 2(a), we illustrate the distance-dependent plasticity function (18)-(20) for a network of N = 12 nodes. The illustration shows the different plasticity functions depending on the distance between the nodes d ij . 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 ij circulant.
If not indicated differently, we consider the coupling structure given by where a : [0, 1] → [0, 1] is a bounded and piece-wise continuous function. This corresponds to a distantdependent 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.

Approximation of the eigenvalues for large systems with circulant structure
In the previous part, we have defined the plasticity functions h ij in such a way that the structures of L h and L Dh inherit important properties from the underlying network structure a(d ij /N ). In particular, assuming that the adjacency matrix is circulant, renders L h and L Dh to be circulant, as well.
In this section, we briefly recall how one can derive the eigenvalues µ k and ν k (k = 0, . . . , N − 1) in case of a circulant structure. It is well-known that for a circulant matrix the eigenvalues are determined by applying a discrete Fourier approach [86]. More precisely, suppose L is a circulant N × N matrix where the elements of the first row are given by the entries l j with j = 1, . . . , N . Then the kth eigenvalue is explicitly given by For the case of L h as in (6), a ij and h ij as in (21) and (18), we obtain with x j = d 1j /N and Re(l h 11 ) = − 1 N N j=2 a(x j )h(0, x j ). Since the adjacency matrix A is assumed to be symmetric, the eigenvalues of L h are real. Therefore, we omit considering the imaginary part of µ k . Equation (22) provides exact expressions for the eigenvalues. However, the values depend on the total number of oscillators N that makes it harder to study the influence of other system properties, such as the coupling structure or the plasticity function. To remove this N -dependence, we consider the continuum limit N → ∞ (compare with [87]) and obtain Due to the definition of h and the symmetry of a(x), we find for any k. This explicit expression allows studying the distribution of the eigenvalues µ k for a given plasticity function h and coupling structure a. Note that a similar expression as (23) can be analogously derived for the eigenvalues of L Dh and reads We note that µ 0 = ν 0 = 0 due to the Laplacian structure of L h and L Dh .

Frontiers
The results from Eqs. (23) and (24) are applied in the next sections to analyze different networks.

Synchronization on nonlocally coupled ring networks
In this section, we analyze the effect of long distance connections on the stability of synchronous states in nonlocally coupled ring networks. We consider the coupling structure given by This means that any two oscillators are coupled if they are separated at most by the coupling range P . The coupling Eq. (25) defines a nonlocal ring structure with coupling range P to each side and two special limiting cases: local ring for P = 1 and globally coupled network for P = N/2 (if N is even, else P = (N + 1)/2). The matrix of the form (25) is circulant [86] and has constant row sum, i.e., N j=1 a ij = 2P for all i = 1, . . . , N . An illustration for N = 12 adn P = 3 is presented in figure 2(b). In order to study the influence of the coupling range, we use the approximations for the eigenvalues µ k and ν k derived in section 4.1. The nonlocally coupled ring structure is expressed by the piecewise continuous function a(x) = 0 for p < x < 1 − p and a(x) = 1 otherwise with relative coupling range p = P/N . Thus, for a nonlocally coupled ring (25) and plasticity function (18)-(20), we find for the eigenvalues µ k of L h . Analogously, we obtain for ν k of L Dh .
In figure 3(a), we provide an error analysis of the approximations (26) and (27) compared to the exact eigenvalues given by Eq. (22). As expected, the errors tend to zero as the number of oscillators increases. Additionally in figures 3(b,c), we display µ k and ν k for several values of k depending on the relative coupling range p. We observe that µ k ≥ 0 for all k. This is due to given plasticity function (18)- (20), for which the update is positive (or equal to zero) for all distances at φ = 0, i.e., h(0, d ij ) ≥ 0 for all d ij .
It is important to note, that our choice of the circulant adaptation functions imply that the matrices L h and L Dh are diagonalizable and commute. Hence, Proposition 2 holds with the master stability equation (11) being exact. Therefore, the stability criterium (14) is also exact.
Combining the fact µ k ≥ 0 with the stability criterium (14), we find cos(α) > 0 as a necessary condition for the stability of the synchronous state for → 0. This yields, that the synchronous state can be stable Figure 3. Panel (a) shows the errors e(µ) (black) and e(ν) (blue) with e(γ) = of the approximations (26) and (27), respectively, where γ exact k are the exact eigenvalues derived by a discrete Fourier transformation, see (22). The errors are displayed in dependence of the system size N (number of oscillators). The relative coupling range is set to p = 0.1. Panel (b) and (c) show the approximated eigenvalues given by (26) and (27), respectively, depending on the relative coupling range p for different values of k.
only for α ∈ (−π/2, π/2). In contrast to L h , the L Dh is in general neither positive nor negative definite, hence the eigenvalues ν k may take positive or negative values. This is due to the fact that the plasticity function may change sign at the origin, i.e., Dh ij may change signs depending on the distance d ij . In particular, we find that only the eigenvalue ν 1 changes the sign, see figure 3(c). This change may lead to a destabilization of the synchronous states as we show in the subsequent analysis. Finally, note that there exist µ ∞ = (1 − cos(2πp))/π and ν ∞ = − sin(2πp)/π to which the eigenvalues converge for large values of k. These limits are displayed in figures 3(b,c) as black lines.
In figure 4, we show different scenarios for the stability of the synchronous state depending on the phase lag parameter α and the coupling range p. Due to the necessary condition cos(α) > 0 as → 0, we consider α ∈ (−π/2, π, 2) only. Figures 4(a) and (b) show that for −π/2 < α < 0, the second stability condition (15) is only fulfilled for p larger than a critical value of the coupling range p c (α). In these cases, a higher coupling range stabilizes the synchronous state. Note that p c (α) → 0 as α → 0 with α < 0. The results seen in figure 4(a,b) are in agreement with the results for a network of N = 200 coupled phase oscillators. For this network, we calculate the Laplacian eigenvalues and plot them along with the master stability function in figure 4(e,f). The outcomes from numerical simulations are presented in figure 4(i,j).
The situation changes for 0 < α < π/2, as shown in figure 4(c,d). Here, for a large range of α, all nonlocally coupled networks lead to a stable synchronous state. However, closer to π/2, long distance connections destabilize the synchronous state. In particular, this destabilization can be traced back to the single negative eigenvalue ν 1 of the Laplacian L Dh , see figure 4(h). Hence, the unstable manifold of the synchronous state is only one-dimensional. This finding is in agreement with the example of N = 200 phase oscillators presented in figure 4(g,h,k,l). Particularly in figure 4(l), the low dimension of the unstable manifold manifests itself as follows: The black trajectory first tends to the synchronous state along the N (N + 1) − 1 stable directions before it is repelled along the direction corresponding to ν 1 .
of the corresponding distance. In the next section, we analyze a network with a more realistic structure with a distance-dependent distribution of weights.

Synchronization on isotropic and homogeneous network with Gaussian distance distribution
In the previous section, we used the prototypical example of a nonlocally coupled rings to study the effects of long-range interaction on synchronization. In this setup, however, all links are equally weighted. In realistic systems, in contrast, the number of links with a certain distance are distributed, see [67] for details. To incorporate this into our network model, we weight the links with respect to a distance distribution. Measurements suggest that the distance distribution can be estimated by a mean and a distribution width [67]. The Gaussian distributions is a paradigmatic distribution that allows for studying effects emanating from the mean and the distribution width. For the remainder of the section, we consider the link distance distribution given by a Gaussian distribution, and weight the links of the network connectivity structure A accordingly, i.e.
where ξ and σ are the mean value and the standard deviation, respectively. Note that the standard deviation characterizes the width of the distribution. For the numerical simulations, we normalize each row of A by N j=1 a ij . Here, we further make the assumption that the network is homogeneous and isotropic. This means that in any direction from a node and at each node the network looks the same. Hence, we obtain a circulant connectivity structure. An illustration of the weight distribution for N = 12 is presented in figure 2(c).
As we know from (14)-(15), for 1, the values of c 2 (α, µ k , ν k ) determine the stability of the synchronous state. In particular, the synchronous state is stable if c min = min k∈1,N −1 c 2 (α, µ k , ν k ) > 0 for a given N and unstable otherwise. In figure 5(a), we display c min for α = −0.4π and different mean values ξ and standard deviations σ of the weight distribution. In agreement with the finding in section 4.2, the synchronized state stabilizes due to an increase of long distance interaction expressed by an increase of σ. Complementing the finding in section 4.2, here, we note that the stability can be also achieved by distributions with peaks at long distance links alone. In this case, the width of the distribution is not important. Figure 5(b) shows how the boundary between regions corresponding to stable and unstable synchronization change for different values of α. As in the case of nonlocally coupled ring networks, with α → 0 (with α < 0) the boundary tends to the limiting point (σ, ξ) = (0, 0). On the contrary, if α → −π/2 (with α > −π/2), the width of the distribution has to increase to have stable synchronization for small values of the mean ξ.
An opposite scenario is shown in figure 5(c) for α = 0.4π. Here, an increase of the weights for long distance links destabilizes the synchronous state, as in figure 4(d,h,l). We also note that for small values α, the synchronous state is stable for almost all values of σ and ξ, see figure 5(d). Only in cases of distribution sharply peaked at long distances, i.e., ξ close to 1/2 and σ close to 0, the synchronous state is unstable. This effect could not be found in networks with nonlocally coupled rings, see section 4.2.

CONCLUSIONS
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,65,66] 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 Dh 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,84] 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 distancedependent plasticity. Here, we have used a ring-like network structure to study the impact of longand 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 towards systems with slightly different local dynamics at each node [88]. 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 [89,90]. 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 [91,92,93] or epilepsy [94,29] for the development of proper therapeutic treatments. We believe that our work provides an important step toward understanding synchronization under realistic conditions.