Synchronization in Networks With Heterogeneous Adaptation Rules and Applications to Distance-Dependent Synaptic Plasticity [original]

Synchronization in Networks With

Heterogeneous Adaptation Rules and

Applications to Distance-Dependent

Synaptic Plasticity

Rico Berner

* and Serhiy Yanchuk

Institut für Mathematik, Technische Universität Berlin, Berlin, Germany,

Institut für Theoretische Physik, Technische Universität

Berlin, Berlin, Germany

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.

Keywords: synaptic plasticity, adaptive networks, phase oscillator, synchronization, distance-dependent synaptic

plasticity, nonlocally coupled rings, master stability approach

1 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 speciﬁc 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 simpliﬁes 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 reﬁned 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

Edited by:

Jun Ma,

Lanzhou University of Technology,

China

Reviewed by:

Syamal Kumar Dana,

Jadavpur University, India

Gopal R.,

SASTRA University, India

*Correspondence:

Rico Berner

[email protected]rlin.de

Specialty section:

This article was submitted to

Dynamical Systems,

a section of the journal

Frontiers in Applied Mathematics and

Statistics

Received: 26 May 2021

Accepted: 21 June 2021

Published: 15 July 2021

Citation:

Berner R and Yanchuk S (2021)

Synchronization in Networks With

Heterogeneous Adaptation Rules and

Applications to Distance-Dependent

Synaptic Plasticity.

Front. Appl. Math. Stat. 7:714978.

doi: 10.3389/fams.2021.714978

Frontiers in Applied Mathematics and Statistics | www.frontiersin.org July 2021 | Volume 7 | Article 7149781

ORIGINAL RESEARCH

published: 15 July 2021

doi: 10.3389/fams.2021.714978

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 speciﬁc

form of the plasticity rule is the distance between neurons

[64–66]. More speciﬁcally, 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 speciﬁcally, 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

conﬁguration. 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 ﬁndings

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.

2 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:

dtϕiω+

j1

aijκijgϕi−ϕj,(1)

dtκij −ϵκij +hijϕi−ϕj,(2)

where ϕi∈S1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 jto i,ωdenotes the natural frequency of each

oscillator, and aij ∈[0,1]are the entries of the weighted

adjacency matrix Adescribing 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 gand hij denote the coupling and

the N2plasticity functions, respectively. For illustrative purposes,

the coupling function is set throughout the paper to g(ϕ)

−sin(ϕ+α)/Nwith 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,2is shift-symmetric,

i.e., invariant under the transformation ϕi1ϕi+ψfor any

ψ∈S1. 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,2from the models

considered previously in the literature [40,70,71,74,82], is that

the plasticity functions hij can be different for each network

connection j→i.

A solution to Eqs. 1,2is called phase-locked if, for all

i1,...,N, the phases evolve as ϕiΩt+ϑiwith some

collective frequency Ω∈Rand ϑi∈S1.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)−wg(0)t,(3)

κs

ij −hij(0),(4)

where we assume that the weighted row sum wN

j1aijhij(0)is

constant for all. Such an assumption of constant row sum is

necessary for the existence of the synchronous state. Moreover, it

is satisﬁed 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.

3 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

Frontiers in Applied Mathematics and Statistics | www.frontiersin.org July 2021 | Volume 7 | Article 7149782

Berner and Yanchuk Synchronization in Networks With Distance-Dependent Plasticity

master stability approach for homogeneous adaptations hij 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−ϕsand χij κij −κs

ij. The linearized equations for these

variations can be written in the following matrix form

dtξ

χJξ

χDg(0)Lhg(0)B

−ϵC−ϵIN2ξ

χ,(5)

where ξ(ξ1,...,ξN)Tis N-dimensional vector containing the

perturbations ξiδϕiof the phases and χ(χ11,χ12,...,χNN )T

are N2- dimensional vectorized perturbations of coupling weights

χvec[δκij], respectively. The N×Nweighted Laplacian matrix

Lhhas the following elements

ij ⎧

⎨

⎩−N

m1,m≠iaimhim(0),ij,

aijhij(0),i≠j.(6)

The time-independent matrices Band Care

B⎛

⎜

⎝a1

aN⎞

⎟

⎠,

C⎛

⎜

⎝(Dh)T

(Dh)T

N⎞

⎟

⎠−⎛

⎜

⎝diag (Dh)1

diag (Dh)N⎞

⎟

⎠,

where ai(ai1,...,aiN ),(Dh)i(Dhi1(0),...,DhiN (0)), and

diag (Dh(0))i⎛

⎜

⎝Dhi1(0)

DhiN (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 Jare the Lyapunov exponents of the

synchronous state and hence determine its local stability. In the

following proposition, we exploit the fact that Jcontains a large

diagonal block −ϵIN2to 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

2N-dimensional linear system

dtvDg(0)Lhg(0)IN

ϵLDh−ϵINv,(7)

where Dg(0)and Lhare as in Eq. 5 and the N ×N weighted

Laplacian matrix LDhpossesses the following elements

lDh

ij ⎧

⎨

⎩−N

m1,m≠iaimDhim(0),ij,

aijDhij(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 N2×N2block given by the

simple diagonal matrix −ϵIN2. This implies that Eq. 5 possess

N2−Nstable directions with Lyapunov exponents −ϵ.Toﬁnd

these directions, we substitute (ξ,χ)e−ϵt(ξ0,χ0)into Eq. 5 and

obtain the linear system

Dg(0)Lh+ϵINg(0)B

−ϵC0ξ0

χ00.(9)

This system has at least N2−Nlinearly independent solutions, since

the matrix in Eq. 9 is degenerate due to the large N2×N2zero block.

The structure of the invariant subspaces in system Eq. 5 allows for

introducing new coordinates, which separate the N2−Nstable

directions (corresponding to the eigenvalues −ϵ)fromthe

remaining 2Ndirections. Explicitly, this transformation is given by

ξ

χRξ

χ,RIN00

0(1/r)BTK

with (N2+N)×(N2+N)matrix R. Here Kis an (N2−N)×

(N2−N)orthogonal matrix with BK 0. Applying this

transformation, we obtain the following system

dt⎛

⎜

⎝ξ

χN

χN2−N⎞

⎟

⎠⎛

⎜

⎝Dg(0)Lhg(0)IN0

ϵLDh−ϵIN0

−ϵKTC0−ϵIN2−N⎞

⎟

⎠⎛

⎜

⎝ξ

χN

χN2−N⎞

⎟

⎠,

(10)

where (ξ,χN,χN2−N)T(ξ,χ)T,with χNand χN2−Nare an Nand

N2−N-dimensional vectors, respectively, and the N×Nweighted

Laplacian matrix LDhas 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 χN2−N.Hence,separatingthe

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 signiﬁcantly

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 Eq. 7 into dynamically independent blocks of

dimension 2. For this, we restrict our consideration to the case

when Lhcan be diagonalized ShQ−1LhQby a nonsingular

complex-valued matrix Q. Note that the eigenvalues μiof Lhlie

on the diagonal of Sh.Ingeneral,thematricesLhand LDhdo not

commute. Therefore, Q−1LDhQis not necessarily of upper

triangular shape. Regardless of this fact, the following

proposition provides an explicit form for the eigenvalues of Jin

Eq. 5 in the limit of slow adaptation, i.e., ϵ≪1.

PROPOSITION 2. Assume that Lhis diagonalizable, with Sh

Q−1LhQ being the associated diagonal matrix and Q the

corresponding transformation. Let ϕiΩt be an in-phase

synchronous state of Eqs. 1,2Then, the local stability of this state

is determined by the solutions of N quadratic equations, which are

given up to the ﬁrst order in ϵas

λ2+ϵ−Dg(0)μiλ−ϵDg(0)μi+g(0)]i0,i1,...,N,

(11)

where μiare the eigenvalues of Lhlocated on the diagonal of Shand ]i

are the corresponding diagonal elements of Q−1LDhQ.If Lhand LDh

commute, then Eq. 11 is exact, and ]iare the eigenvalues of LDh.

Frontiers in Applied Mathematics and Statistics | www.frontiersin.org July 2021 | Volume 7 | Article 7149783

Berner and Yanchuk Synchronization in Networks With Distance-Dependent Plasticity

PROOF. Due to Proposition 1, the eigenvalues of the Jacobian in

Eq. 5 are given by

det Dg(0)Lh−λINg(0)IN

ϵLDh−(ϵ+λ)IN

det Dg(0)Sh−λINg(0)IN

ϵQ−1LDhQ−(ϵ+λ)IN0,

where we have used the transformation Qthat brings Lhto the

diagonal form ShQ−1LhQ. Making further use of the Schur

complement [84], we obtain

det Dg(0)Sh−λINg(0)IN

ϵQ−1LDhQ−(ϵ+λ)IN

det(λ+ϵ)λIN−Dg(0)Sh−ϵg(0)Q−1LDhQ0.(12)

The latter equation is almost diagonal. The only off-diagonal

components remain from Q−1LDhQand scale with ϵ. Let

us consider the Leibniz formula for the determinant of an

N×Nmatrix Fwith entries fij, that reads

det(F)σ∈Perm(N)sgn(σ)N

i1fiσ(i). In the latter expression

Perm(N)denotes the set of all permutations σof the

integer numbers 1,...,Nand 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 N

i1fiσ(i)scales with ϵ2

or higher. Hence, we are left with det(F)N

i1fii +O(ϵ2)

and ﬁnd

det (λ+ϵ)λIN−Dg(0)Sh−ϵg(0)Q−1LDhQ



i1λ2+ϵ−Dg(0)μiλ−ϵDg(0)μi+g(0)]i+Oϵ2

0,

(13)

where μiare the eigenvalues of Lh,]iare the diagonal elements

of Q−1LDhQand O(ϵ2)denotes higher order terms (ϵm,m>1).

If Lhand LDhcommute, both matrices share the same set of

eigenvectors and hence they can be brought to the diagonal

form with the same transformation Q.Inthiscase,the

diagonal elements ]iare the eigenvalues of LDhand the

higher order terms O(ϵ2)in Eq. 13 vanish.

The 2Nsolutions λiof the NEq. 11 determine the stability

of the synchronous state. More precisely, the real parts of

theses solutions determine the Lyapunov exponents. If

ΛmaxiRe(λ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 ﬁxed 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 Dg(0)μand g(0)].Therefore,wedeﬁne the master

stability function Λ:C2→Rwith Λ(Dg(0)μ,g(0)])

maxiRe(λi(Dg(0)μ,g(0)])) 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 Lhand LDhsince their eigenvalues are real. In Figure 1,we

present the master stability function for the coupling function

g(ϕ)−sin(ϕ+α)/Nand 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

c1α,μcos(α)μ>−ϵ,(14)

c2α,μ,]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.

4 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

dtϕiω−1

N

j1

aijκijsinϕi−ϕj+α,(16)

dtκij −ϵκij +hϕi−ϕj,dij.(17)

The plasticity function hdepends on the phase difference ϕi−ϕjand

on the distance dij . In this work, we associate the distance to the difference

of indices by dij j−i. For the plasticity function, we consider

hijϕhϕ,dij

N⎧

⎪

⎨

⎪

⎩hϕ,dij

Ndij ≤N2,

hϕ,1−dij

Ndij >N2.

(18)

With this form of the adaptation function, we have a

symmetric hij(ϕ)hji(ϕ)and a circulant hi+l,j+l(ϕ)hij(ϕ)

structure of the corresponding matrix with entries hij.

Particularly, for the numerical analysis, we use

hϕ,dijNsinϕ+βdijN,(19)

where the distance dependence is encoded in the phase shift function

βdij

N⎧

⎪

⎨

⎪

⎩2

Ndij −1π,Neven,

2

(N+1)dij −1π,Nodd.

(20)

Frontiers in Applied Mathematics and Statistics | www.frontiersin.org July 2021 | Volume 7 | Article 7149784

Berner and Yanchuk Synchronization in Networks With Distance-Dependent Plasticity

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 dij. 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 ﬁndings in [67]. Note the symmetry of the

plasticity function that renders the matrix with elements hij

circulant.

If not indicated differently, we consider the coupling structure

given by

aij adij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 Lhand LDh for large networks with circulant

FIGURE 1 | The master stability function Λ(Dg(0)μ,g(0)])for the coupling function g(ϕ)−sin(ϕ+α)/Nand real μand ν(Im(μ)0, Im(])0). The values of the

master stability function are color-coded in all panels (A–E). The dashed black line describes the border between regions corresponding to local stability and instability,

respectively. Parameters: ϵ0.01, (A) α−0.8π,(B) α−0.4π,(C) α0, (D) α0.4π, and (E) α0.8π.

FIGURE 2 | Panel (A) shows the plasticity function h1jgiven in Eqs. 18–20 depending on the distance d1jexempliﬁed for node i1 in a network with N12 nodes.

Note that the colors of the links in the network (left) correspond to the colors of the depicted plasticity function (right). Panel (B) displays the connectivity structure of a

nonlocally coupled ring network with N12 nodes and a coupling range P3. Panel (C) displays the weighted connectivity structure of a network with N12 nodes

(left) with distance-dependent Gaussian weight distribution (right). Note that the colors of the links in the network (left) correspond to the colors of the bars in the

weight distribution (right).

Frontiers in Applied Mathematics and Statistics | www.frontiersin.org July 2021 | Volume 7 | Article 7149785

Berner and Yanchuk Synchronization in Networks With Distance-Dependent Plasticity

Loading more pages...