Document [original]

Fatemeh Parastesh, Karthikeyan Rajagopal, Sajad Jafari, Matjaž

Perc, Eckehard Schöll

Blinking coupling enhances network

synchronization

Open Access via institutional repository of Technische Universität Berlin

Document type

Journal article | Accepted version

(i. e. final author-created version that incorporates referee comments and is the version accepted for

publication; also known as: Author’s Accepted Manuscript (AAM), Final Draft, Postprint)

This version is available at

https://doi.org/10.14279/depositonce-17011

Citation details

Parastesh, F., Rajagopal, K., Jafari, S., Perc, M., & Schöll, E. (2022). Blinking coupling enhances network

synchronization. In Physical Review E (Vol. 105, Issue 5). American Physical Society (APS).

https://doi.org/10.1103/PhysRevE.105.054304.

This work is protected by copyright and/or related rights. You are free to use this work in any way permitted by

the copyright and related rights legislation that applies to your usage. For other uses, you must obtain

permission from the rights-holder(s).

Blinking coupling enhances network synchronization

Fatemeh Parastesh,1Karthikeyan Rajagopal,2Sajad Jafari,1, 3 Matjaˇz Perc,4, 5, 6 and Eckehard Sch¨oll7, 8, 9, ∗

1Department of Biomedical Engineering, Amirkabir University of Technology (Tehran polytechnic), Iran

2Center for Nonlinear Systems, Chennai Institute of Technology, India

3Health Technology Research Institute, Amirkabir University of Technology (Tehran polytechnic), Iran

4Faculty of Natural Sciences and Mathematics, University of Maribor, Koroˇska cesta 160, 2000 Maribor, Slovenia

5Department of Medical Research, China Medical University Hospital,

China Medical University, Taichung 404332, Taiwan

6Complexity Science Hub Vienna, Josefst¨adterstraße 39, 1080 Vienna, Austria

7Institut f¨ur Theoretische Physik, Technische Universit¨at Berlin, Hardenbergstrasse 36, 10623 Berlin, Germany

8Bernstein Center for Computational Neuroscience Berlin, Humboldt-Universit¨at, 10115 Berlin, Germany

9Potsdam Institute for Climate Impact Research, Telegrafenberg A 31, 14473 Potsdam, Germany

This paper studies the synchronization of a network with linear diffusive coupling, which blinks

between the variables periodically. The synchronization of the blinking network in the case of

sufficiently fast blinking is analyzed by showing that the stability of the synchronous solution depends

only on the averaged coupling and not on the instantaneous coupling. To illustrate the effect of

the blinking period on the network synchronization, the Hindmarsh-Rose model is used as the

dynamics of nodes. The synchronization is investigated by considering constant single-variable

coupling, averaged coupling, and blinking coupling through a linear stability analysis. It is observed

that by decreasing the blinking period, the required coupling strength for synchrony is reduced. It

equals that of the averaged coupling model times the number of variables. However, in the averaged

coupling, all variables participate in the coupling, while in the blinking model only one variable

is coupled at any time. Therefore, the blinking coupling leads to an enhanced synchronization in

comparison with the single-variable coupling. Numerical simulations of the average synchronization

error of the network confirm the results obtained from the linear stability analysis.

PACS numbers: 05.45.Xt

I. INTRODUCTION

Complex dynamical networks have attracted much at-

tention in recent years [1]. A universal phenomenon in

these networks is the synchronized behavior of the com-

ponents [2–4]. It has been shown that the structure of

the network plays a key role in synchronization [5]. For

this reason, many studies were focused on the influence

of the network topology on synchronization. Wang et al.

[6] proposed that the synchronizability of a homogeneous

network can be enhanced by considering weighted and

asymmetric couplings, similar to a scale-free network.

Nishikawa et al. [7] reported that networks with a ho-

mogeneous distribution of connectivity, although having

larger average path lengths, are more likely to synchro-

nize than those with heterogeneous connectivity.

Enhancing synchronization is of great importance in

many applications, including diverse brain functions

[8, 9]. For example, synchronization is essential in many

memory processes such as working memory and long-

term memory by enhancing neural plasticity [10]. In the

attention-related process, the neurons receiving the at-

tended stimuli exhibit enhanced synchrony in the gamma

band [11]. In contrast to these desirable functions, some

of the pathological brain states are caused by increased

synchrony in special regions. For instance, in patients

∗scho[email protected]

with Parkinson’s disease, abnormal synchronization is

observed in the cortico-basal ganglia circuits [12]. In the

past years, many efforts have been devoted to enhanc-

ing synchrony in dynamic networks. Lin and Chen [13]

demonstrated that two chaotic oscillators could become

synchronous by applying white-noise-based coupling. It

has been shown that assigning the direction to the links

of a network, for example, by the residual degree gradi-

ent (RDG) method or the residual edge-betweenness gra-

dient method [14, 15] can improve the synchronization.

Ramirez et al. [16] proposed a dynamic coupling for a

master-slave network and showed the enhanced synchro-

nization, which could not be obtained with static cou-

pling. Banerjee et al. [17] found that applying a param-

eter mismatch to well-defined oscillators of the network,

such that the identical oscillators interact indirectly, can

help achieve synchronization. Sevilla-Escoboza et al. [18]

investigated the stability of synchronization by consid-

ering multivariable couplings and extracted the optimal

scheme that resulted in maximum stability. Panahi et al.

[19] revealed that the optimal synchronization in circu-

lant oscillators is obtained by multivariable coupling with

equal coupling coefficients. Time delay is another factor

impacting synchronization [20, 21]. Kyrychko et al. [22]

examined the stability of different types of synchroniza-

tion in different topologies considering distributed time

delays in coupling. They considered uniform and gamma

delay distributions and found that the stability of the

synchronization is improved with increasing the width

of a uniform distribution or decreasing the mean of the

gamma distribution. Gjurchinovski et al. [23] studied

a network of coupled limit-cycle oscillators with time-

varying delays in the coupling and self-feedbacks. They

analyzed the stability of synchronization and found that

the time-varying delay leads to the formation of ampli-

tude death. Liu et al. [24] considered a network with

time-varying delay and investigated local and global ex-

ponential synchronization. By using the average dwell

time method, delay-dependent sufficient conditions were

derived.

In realistic networks, the interactions between compo-

nents are dynamic and evolve in time [25]. The varying

communications between individuals in society or the in-

teractions between neurons in the brain are some exam-

ples. In this context, many researchers have focused on

synchronization in time-varying networks in the last two

decades. Belykh et al. [26] proposed a small-world net-

work with fixed nonlocal connections and time-varying

links between any pairs of oscillators with the chaotic

attractor. The blinking connections were switched on

and off with pre-defined probability. It was found that

the addition of random links reduces the effective path

length and also enhances synchronization with a lower

cost. Hasler et al. [27] considered blinking random con-

nections in non-neighboring nodes in a multi-stable net-

work. They showed that for sufficiently small blinking pe-

riods, i.e., fast switching, the behavior of the blinking net-

work and the time-averaged network is almost the same.

A comparison of the fast switching network and its aver-

aged network has also been performed for finite time and

infinite time intervals considering different cases accord-

ing to the invariance or non-invariance of the numbers

of attractors of the averaged system [28, 29]. Porfiri [30]

established a method for finding necessary and sufficient

conditions for the stochastic synchronization of a network

of chaotic maps with blinking links. Lu et al. [31] ana-

lyzed the synchronization of discrete-time networks with

time-varying topologies by using the Hajnal diameter,

which is equivalent to transverse Lyapunov exponents.

Besides blinking small-world networks, synchronization

of time-varying random networks has also been consid-

ered. Jeter and Belykh [32] studied the global stability of

the synchronization in networks with time-varying ran-

dom topology and intrinsic parameters. In another study,

they found an optimal window of frequency in which the

synchronization is stable [33]. Barabash and Belykh [34]

studied a random network with time-varying connections

and showed its equivalence with the averaged network for

fast blinking. They found that synchronization is main-

tained even when increasing the blinking period consid-

erably. Furthermore, with increasing number of oscilla-

tors in the network, the synchronization threshold be-

comes independent from the blinking period. The effect

of memory in the blinking links has also been studied. Liu

et al. [35] proposed a novel approach for deriving suffi-

cient conditions for the synchronization of networks with

time-varying coupling structure and weight. They con-

sidered two cases for stochastic switching processes where

in the first case, the sequences had an independent and

identical random distribution and in the second case, the

sequences created a Markov chain. Porfiri and Belykh

[36] investigated two one-dimensional coupled nonlinear

maps under Markovian switching with respect to memory

effects. Lanza et al. [37] analyzed synchronization in cou-

pled memristor-based oscillatory circuits since memris-

tors have represented a crucial function in the emergence

of synchronization in special cases [38]. They showed

the effects of blinking links on the invariant manifolds of

oscillators.

In the studies mentioned above, the topologies of the

networks were considered to be time-varying with fixed

coupling functions. However, there exists another con-

figuration where the network connections between the

nodes are fixed, but the coupling function is a time-

varying function. The coupling function determines

which dynamic variables of different nodes are connected

with each other, i.e., it determines how the interactions

among two (or more) dynamical systems evolve [39].

There are many applications in which the coupling func-

tion is time-varying such as the cardiorespiratory sys-

tem, transport grids and supply networks, and neural

cross-frequency coupling functions [40, 41]. In general,

a coupling function is defined by its strength and its

form. Some previous studies have focused on networks

with time-varying coupling strength [42–44], while some

have investigated the effects of time-varying coupling

forms [39, 45–47]. For example, Hagos et al. [45] reported

synchronization transitions induced by time-varying cou-

pling functions in phase oscillators. They indicated that

the collective behavior of the oscillators depends on the

shape of the coupling function, and the net coupling

strength has a negligible effect. As an application in

machine learning, Stelzer et al. [47] discussed deep neu-

ral networks with step-like switching functions between

multiple time-delayed feedbacks to achieve a better, more

efficient performance. These virtual networks consist of

a single node with multiple time-delayed feedback loops.

In this paper, we consider a network with fixed con-

nections and a special time-varying coupling function.

The variable through which the oscillators are coupled

is assumed to alternate between the system’s variables

periodically. Thus, the coupling function turns out to

be blinking. The synchronization of the network is ana-

lyzed under fast blinking. The Hindmarsh-Rose neuron

model is considered as the dynamics of each node of the

network. The synchronization stability is investigated by

computing the largest Lyapunov exponent of the varia-

tional equation for different blinking periods. It is found

that by decreasing the blinking period, synchronization is

achieved for lower coupling strengths. Numerical simula-

tions are also performed, and the average synchronization

error is calculated.

II. THE BLINKING NETWORK MODEL

A time-varying dynamical network composed of N

identical oscillators with linear diffusive coupling is con-

sidered. It is assumed that the links are constant, but

the variable used in the coupling is switched periodically

in time. Thus, the network can be described by

Xi=F(Xi) + σPN

j=1 GijH(t)Xj, i = 1, ..., N, (1)

where Xi= (xi1, xi2, ..., xid) is the d-dimensional state

variable of the oscillators, F(Xi) : Rd→Rdis the sys-

tem’s dynamics, and σis the coupling strength. The

network topology is determined by the adjacency ma-

trix GN×Nwith zero row sum (i.e., the negative of

the Laplacian matrix), where Gij = 1 if the ith and

jth oscillators are connected and Gij = 0 else, and

Gii =−PN

j=1,j=iGij . The time-varying matrix H(t) :

Rd→Rdis the internal-coupling matrix and determines

which variables are considered in the coupling at each

time t. Here, the matrix H(t) is assumed to blink be-

tween variables with equal time intervals with period

τ. Therefore, the matrix H(t) and its elements Hmn(t)

(m, n = 1, ..., d)can be described as follows

Hmm(t) = (1 if (m−1)τ

d<t<mτ

0 otherwise ,

Hmn(t)=0,

H(t+τ) = H(t).

(2)

This time-varying network can have a synchronous so-

lution s(t) = Xi(t), i= 1, ..., N. The stability of the

synchronous manifold is equivalent to the stability of the

error vector ηi(t) = Xi(t)−s(t) with respect to pertur-

bations around zero. Substituting ηi(t) into the network

equation (Eq. 1), one obtains the linearized variational

equation

˙ηi(t) = DFs(t)ηi(t) + σPN

j=1 GijDHs(t)ηi(t)(3)

where DFs:Rd→Rdand DHs:Rd→Rdare the

Jacobian matrices of Fand Hat s(t). Note that since

the coupling is considered to be linear diffusion, we have

DHs(t) = H(t). Thus, the system of linearized coupled

oscillators is given in terms of the Nd-dimensional vector

η(t) = (η1, η2, ..., ηN) by

η(t)=(IN⊗DFs(t) + σ(G⊗H(t)))η(t)(4)

where ⊗is the Kronecker product. The matrix Gcan be

diagonalized using Schur decomposition. Therefore, an

upper triangular matrix is formed with eigenvalues of G

appearing on its main diagonal. Finally, diagonalizing G,

and since the first term of (Eq. 4) is block diagonal, the

set of variational equations (Eq. 4) can be transformed

to the decoupled system

˙η(t)=[DFs(t) + σλkH(t)]η(t), k = 1, ..., N (5)

FIG. 1: a) Time series of the membrane potential x1of

the single Hindmarsh-Rose neuron model (Eq. 11) that

exhibits chaotic bursting. The parameters are set to

Iext = 3.2, r= 0.006, s= 4. (b) Synchronous time

series of two neurons coupled through x1variables and

coupling strength σ= 1. The blue and red colors

represent the time series of the membrane potential of

the first and second neurons x11 and x21, respectively.

The parameters are as in part (a).

where λkare the eigenvalues of the adjacency matrix (G)

and η(t) is a d-dimensional vector. For a network with a

connected graph, the first eigenvalue is always zero, and

thus, the system (Eq. 5) for k= 1 evolves along the

synchronous manifold. Consequently, the stability of the

system (Eq. 5) needs to be checked for k= 2, ..., N,

which corresponds to the directions transverse to the

synchronous manifold. This method is called the mas-

ter stability approach and has been proposed by Pecora

and Carroll [5]. The master stability function (MSF)

is defined as the largest Lyapunov exponent (Λ) of the

variational equation (5) as a function of the complex pa-

rameter µ=λk. If Λ(µ)<0 for all eigenvalues µ=λk,

k= 2, ..., N, the synchronized solution is stable. Huang

et al. [48] systematically studied the typical behavior of

the MSF with different coupling schemes for some well-

known chaotic systems and categorized the MSF behav-

ior into four classes.

III. SYNCHRONIZATION UNDER FAST

BLINKING

The stability of synchronization of the time-varying

network with a sufficiently fast blinking coupling function

can be estimated by the synchronization of the averaged

network. Stilwell et al. [49] proved this theory for the

network with time-varying topology. Here, we use the

same approach for the time-varying coupling function.

In this regard, the following theorem is given.

Theorem 1. It is supposed that the system of oscilla-

tors with linear diffusive coupling and the static internal

coupling function ( ¯

H) as,

Xi=F(Xi) + σPN

j=1 Gij ¯

HXj, i = 1, ..., N (6)

has a stable synchronization manifold. Then, there exists

a small ε∗, such that the set of oscillators with time-

varying internal coupling function as

Xi=F(Xi) + σPN

j=1 GijH(t/ε)Xj, i = 1, ..., N (7)

reaches a stable synchronization manifold for 0< ε < ε∗,

τZt+τ

H(α)dα =¯

H. (8)

The above theorem can be easily proved according to

the following lemma, which is given for the fast switching

systems.

Lemma 1. Supposing ˙x(t) = (A(t) + ¯

E)x(t)has a uni-

formly exponentially stable solution, where Eis a ma-

trix function satisfying ¯

E=1

τRt+τ

tE(α)dα for all t.

Then, there exists small ε∗such that for 0< ε < ε∗,

˙z(t)=(A(t) + E(t/ε))z(t)has a uniformly exponentially

stable solution.

The proof of the lemma is presented in the appendix

[49].

Proof of Theorem 1. In the previous section, it was de-

scribed that the stability of synchronization of the net-

work (Eq. 6) is equivalent to the stability of the linearized

variational equation,

˙η= [DFs+σλk¯

H]η, k = 2, ..., N. (9)

Supposing that the network (6) achieves synchroniza-

tion, then the system (9) is exponentially stable. Ac-

cording to the internal coupling matrix defined in the

previous section and Eq. (2), there exists an average ma-

trix ¯

H=1

τRt+τ

tH(α)dα =





1/d 0. . . 0

0 1/d ...0

.......0

0. . . 0 1/d







for all

t. Thus, according to the lemma, there is ε∗such that

the following system is exponentially stable,

˙η= [DFs+σλkH(t/ε)]η, k = 2, ..., N. (10)

The system (10) can be considered as the linearized vari-

ational equation of the network (7). Therefore, the time-

varying network (Eq. 7) can achieve stable synchroniza-

tion.

From this theory, it can be concluded that when the

blinking of the internal coupling function is sufficiently

fast, i.e., the period of blinking is considerably smaller

than the period of the oscillators, the stability of syn-

chronization of the blinking network is determined by the

synchronization stability of the averaged network. There-

fore, the synchronization stability is not dependent on the

coupling function H(t) at time tbut on the average of

H(t).

IV. BLINKING NEURONAL NETWORK

It has been shown that the interactions in the brain

are not static and evolve in time [50–54]. The synap-

tic connections and also the strength of the connections

change temporally to optimize the functionality of the

neurons [55]. Since the focus of our studies is the blink-

ing coupling function and not the network topology, we

use a very simple network consisting of two nodes in order

to gain insight into the mechanism of blinking coupling.

Therefore, as an illustration of blinking coupling, we

choose a well-known dynamic model from neuroscience,

and investigate synchronization in a simple network of

two coupled Hindmarsh-Rose neurons where the dynam-

ics of each node (F(X)) is described by

˙x1=x2+ 3x2

1−x3

1−x3+Iext

˙x2= 1 −5x2

1−x2

˙x3=r(s(x1+ 1.6) −x3)

(11)

where x1,x2and x3denote the membrane potential, the

fast and slow recovery variables. The parameters are set

to Iext = 3.2, r= 0.006, s= 4, where each node exhibits

chaotic bursting. The time series of the Hindmarsh-Rose

model with these parameters is shown in Fig. 1a. In our

simplest example, the network is assumed to be com-

posed of two nodes via G=−1 1

1−1with eigenvalues

λ1= 0 and λ2=−2. The synchronizability of the net-

work is obtained by linear stability analysis and numer-

ical simulations considering different cases for the inter-

nal coupling function as a) constant coupling function on

each state variable, b) averaged coupling function ( ¯

H),

c) blinking coupling function as Eq. 2. It is notable that

the synchronous manifold in all cases is also chaotic. An

example of the time series of the synchronous neurons is

illustrated in Fig. 1b.

FIG. 2: The largest Lyapunov exponent Λ of the variational equation (5) for µ=−2 with constant coupling

function (H) in dependence on the coupling strength σ. (a) Coupling in x1variables, (b) Coupling in x2variables,

A. Linear stability analysis

At first, the coupling function is considered to be con-

stant, with the coupling in one variable only. In this

case, the stability of synchronization can be obtained by

finding the stability of the variational equation (Eq. 5)

with constant Hwith respect to the zero solution. If

the master stability function, defined as the largest Lya-

punov exponent (Λ) of the variational equation (5) as a

function of the complex parameter µ, satisfies Λ(µ)<0

for all eigenvalues µ=λkof the adjacency matrix, the

synchronized solution is stable. For the considered cou-

pling matrix G=−1 1

1−1with eigenvalues λ1= 0 and

λ2=−2, we have to compute the largest Lyapunov ex-

ponent (Λ) of the variational equation for µ=λ2=−2.

Figure 2 shows the largest Lyapunov exponent Λ(µ=

−2) for constant Has a function of the coupling strength.

In Fig. 2a, the coupling of the oscillators is through the

x1variables, i.e., H=



100

000

000

. With this coupling,

synchronization can be attained for σ > 0.465. Figure

2b represents Λ(µ=−2) for coupling in the x2variables,

i.e., H=



000

010

000

. It is observed that for this coupling,

the network becomes synchronous for σ > 0.056. Finally,

the coupling is through the x3variables (Fig. 2c), i.e.,

FIG. 3: The largest Lyapunov exponent (Λ(µ=−2))

of Eq. 3 vs the coupling strength (σ) for different

blinking periods τ. Blue: τ=T= 30, Red:

τ=T/5 = 6, Orange: τ=T/10 = 3, Purple:

τ=T/100 = 0.3, Green: τ=T/1000 = 0.03. For faster

blinking, the largest Lyapunov exponent is the same as

the green line with threshold σ= 0.021. Other

parameters as in Fig. 1.

H=



000

001

. In this case, the synchronization is un-

stable for any value of σ. Next, H=



1/3 0 0

0 1/3 0

0 0 1/3

is

considered, which is the average of the blinking network.

The largest Lyapunov exponent, which is illustrated in

Fig. 2d, shows that the synchronization of the averaged

network is achieved for σ > 0.021. It can be clearly seen

that the largest Lyapunov exponent is a decreasing lin-

ear function of the coupling strength σ, which follows

from the diagonal coupling in all three variables [56, 57],

in contrast to the coupling in only one variable in Fig.

2a-c. The stability of the synchronized solution η= 0 is

governed by the linearized equation ˙η= [DFs+λ2σI/3]η

with λ2=−2 for perturbations of η= 0, where Iis the

unity matrix. In fact, for a constant matrix A0=DFs

the eigenvalues γ0of A0are related to the eigenvalues

γof A= [DFs+λ2σI/3] by γ=γ0+λ2σ/3 which fol-

lows from comparing the two eigenvalue equations. Since

the Lyapunov exponents correspond to the real part of

the eigenvalues averaged along the attractor’s orbits, the

eigenvalue equation γ=γ0+λ2σ/3 can be converted to

the equation Λ = Λ0+λ2σ/3, where Λ0is the maximum

Lyapunov exponent of DFs. Therefore, the maximum

Lyapunov exponent Λ as a function of σis given by a

straight line Λ = Λ0−2σ/3. Consequently, the threshold

for synchronization Λ = 0 can be calculated as σ= 3Λ0/2

which gives σ= 0.0207 for Λ0= 0.0138.

Next, we investigate the stability of the synchronous

solution of the blinking network. The largest Lyapunov

exponent of Eq. 5 is calculated for the time-varying cou-

pling function, i.e.,

H(t) =

















100

000





if 0 < t < τ/3







000

010

000





if τ/3<t<2τ/3







000

001





if 2τ/3< t < τ

H(t+τ) = H(t)

(12)

Figure 3 shows Λ(µ=−2) for different blinking periods

τas a function of the coupling strength σ. Although

the coupling function is discontinuous in time, the linear

stability analysis is straight-forward since the coupling is

linear. The master stability function has been extended

to nonsmooth (discontinuous) nonlinear coupling in [58].

The time-scale of the spikes in the chaotic bursting os-

cillations of the single Hindmarsh-Rose model is approx-

imately equal to T= 30. We have chosen the periods of

blinking as τ=T, T/5, T/10, T/100, T/1000, T/10000 =

30,6,3,0.3,0.03,0.003. The time step for solving the

equations is set at dt = 0.0002. It can be seen that

as the blinking occurs faster (with decreasing the period

of blinking τ), the threshold of the coupling strength at

which the synchronization is achieved is decreased. Fur-

FIG. 4: The synchronization error vs the coupling

strength σwith a constant coupling in (a) x1-variables,

(b) x2-variables, (c) x3-variables. Parameters as in Fig.

thermore, when the blinking period decreases sufficiently,

i.e., for τ < 0.03, synchronization can be obtained for

σ > 0.021. For fast blinking the stability is the same as

for the averaged coupling, and hence the maximum Lya-

punov exponent decreases linearly with coupling strength

as in Fig. 2d. One can see that the function becomes

more and more like a straight line with decreasing blink-

ing period τ. For the smallest value of τ(green line)

the graph becomes identical to Fig. 2d, with the same

critical coupling strength.

B. Numerical results

We have also solved the network dynamics of two cou-

pled Hindmarsh-Rose neurons numerically. To evaluate

the synchronization between two neurons, the temporally

averaged synchronization error [59] is computed as

Error =⟨∥X1(t)−X2(t)∥⟩t(13)

FIG. 5: The synchronization error vs coupling strength σwith blinking coupling for different blinking periods τ.

(a) τ=T= 30, (b) τ=T/5 = 6, (c) τ=T/10 = 3, (d) τ=T/100 = 0.3, (e) τ=T/1000 = 0.03, (f)

τ=T/10000 = 0.003. Other parameters as in Fig. 1.

where ⟨.⟩tdenotes the time average. At first, the coupling

is assumed to be constant through one of the variables.

The synchronization error for the x1-variable coupling is

shown in Fig. 4a. It can be seen that the neurons are

synchronized for σ > 0.48. When the coupling is only

in the x2-variables, the neurons become synchronous for

σ > 0.05. The synchronization error for this coupling is

shown in Fig. 4b. Finally, when the neurons are coupled

through the x3-variables, synchronization is not achieved

at all by increasing the coupling strength (Fig. 4c). Com-

paring Fig. 2 and Fig. 4 shows that the critical coupling

strengths found by the numerical solution are almost the

same as those obtained from the linear stability analysis.

Now it is assumed that the coupling blinks between the

three variables. Therefore, in the interval 0 < τ < τ/3,

the coupling is in the x1-variables, in τ/3< t < 2τ/3,

the coupling is in the x2-variables, and in 2τ/3< t < τ,

the coupling is in the x3-variables. The synchronization

errors for various blinking periods as in Fig. 3 are shown

in Fig. 5. It is observed that when the period of blink-

ing is long, the synchronization threshold is large. With

decreasing blinking period τ, synchronization is achieved

for smaller coupling strengths, such that for τ= 0.03

the synchronization is stable approximately for σ > 0.02.

The numerically found thresholds for stable synchroniza-

tion are consistent with the results derived from the linear

stability analysis.

To obtain an overview of the effect of the blinking

period on the synchronization threshold, the parameter

plane is presented in Fig. 6. The blinking period is con-

sidered to vary in the range τ∈[0.01,30] in a logarithmic

scale to represent the short periods better. The values of

the largest Lyapunov exponent of the linearized system

(Eq. 5) for µ=−2 are shown in Fig. 6a. The black

curve in this figure marks Λ = 0 and separates the re-

gions of asynchronous and synchronous behavior. Figure

6b represents the numerically obtained synchronization

error (Eq. 13) of the network. Note the nonmonotonic

behavior of the synchronization threshold as a function

of τwhich is visible both in Fig. 6a and 6b.

FIG. 6: The region of stable synchronization in the parameter plane of coupling strength σand blinking period τ.

(a) The largest Lyapunov exponent of the linearized system (Eq. 5) for µ=−2; the asynchronous and synchronous

states regions are separated by the black curve (Λ = 0). (b) Numerically calculated synchronization error.

Parameters as in Fig. 1.

V. CONCLUSIONS

In this paper, the synchronization of a time-varying

network with blinking coupling variables was investi-

gated. It was assumed that the links of the network

are fixed and time-independent, but the variable that

participates in the coupling switches between the sys-

tem’s variables with a well-defined periodicity. For the

case of fast blinking, a theorem was presented which

shows that the synchronization of the network under fast

blinking is equivalent to the synchronization of the av-

eraged network. Therefore, the instantaneous coupling

configuration does not affect the synchrony of the net-

work. Thus, in general, the blinking coupling function

can enhance synchronization of systems in which the

coupling strength needed for synchronization in single-

variable coupling is higher than that for averaged diag-

onal coupling. This is a significant result and can be

applied for enhancing the synchronization of any phys-

ical systems such as chaotic electronic circuits in which

the coupling in all variables can be readily implemented

[18], if those systems satisfy the above condition.

As an illustration, the Hindmarsh-Rose model was

used to describe the dynamics of the nodes. At first, the

synchronization of the network with constant coupling

was studied by calculating the master stability function.

It was observed that the best case for synchronization

is the coupling through the x2-variables, which provides

synchrony for σ > 0.056, but synchronization is never

achieved for coupling in the x3-variables. Then, the nec-

essary conditions for the stability of the synchronization

of the blinking network were obtained for different blink-

ing periods. It was found that a smaller period of switch-

ing leads the network towards synchronization at smaller

coupling strengths. For sufficiently fast blinking, the

blinking network with coupling strength σis dynamically

equivalent to fixed all-variable diagonal coupling with

coupling strength σ/3, and the coupling strength needed

for synchrony is considerably smaller than that needed in

constant single-variable coupling, namely σ > 0.021. It

is intriguing that the blinking coupling provides a better

synchronization for the systems, although being periodi-

cally coupled in the x3variable, which as a stand-alone

coupling cannot lead to synchronization. The network

was also studied numerically, and the time-averaged syn-

chronization error was computed. The numerical results

agreed well with the conditions found from the linear sta-

bility approach.

ACKNOWLEDGMENTS

E.S. acknowledges support by DFG project Nos.

429685422 and 440145547. M.P. is supported by the

Slovenian Research Agency (Grant Nos. P1-0403 and

J1-2457).

Appendix: Proof of Lemma 1

Here, the proof of Lemma 1 is presented [49].

In this lemma, it is supposed that the system ˙x(t) =

(A(t) + ¯

E)x(t) is uniformly exponentially stable. There-

fore, there exists a Lyapunov function v(x(t), t) =

xT(t)Q(t)x(t) with symmetric matrix Q(t) such that

η∥x(t)∥2≤v(x(t), t)≤ρ∥x(t)∥2(A.1)

and

dtv(x(t), t)≤ −µ∥x(t)∥2(A.2)

where η > 0, ρ > 0, and µ > 0.

To prove the exponential stability of ˙z(t) = (A(t) +

E(t/ε))z(t), it is shown that v(z(t), t) is its Lyapunov

function and the following equation is negative definite.

∆v(z(t+εT), t)≡v(z(t+εT, t +εT)−v(z(t), t) (A.3)

Eq. A.3 can be expanded to

∆v(z(t+εT), t) = zT(t+εT)Q(t+εT)z(t+εT)

−zT(t)Q(t)z(t) = zT(t)(ΦT

E(t+εT, t)Q(t+εT)

ΦE(t+εT, t)−Q(t))z(t)

(A.4)

where ΦE(t, τ) denotes the transition matrix of A(t) +

E(t/ε). Then, H(t+εT, t) = ΦE(t+εT, t)−Φ¯

E(t+εT, t)

is defined, where Φ ¯

E(t+εT, t) denotes the transition ma-

trix of A(t) + ¯

H(t+εT, t) = Φ(t+εT, t)−Φ¯

E(t+εT, t) =

I+Rt+εT

tA(σ1) + E(σ

ε)dσ+

P∞

i=2 Rt+εT

tA(σ1) + E(σ1

ε)Rσ1

t... Rσi−1

tA(σi) + E(σi

ε)dσi...dσ1

−I−Rt+εT

tA(σ1) + ¯

E dσ−

P∞

i=2 Rt+εT

tA(σ1) + ¯

ERσ1

t... Rσi−1

tA(σi) + ¯

E dσi...dσ1

(A.5)

With considering

Zt+εT

E(σ

ε)dσ =εT ¯

E(A.6)

Eq. A.5 turns into Eq. A.7

H(t+εT, t) = P∞

i=2 Rt+εT

tA(σ1) + E(σ1

ε)

Rσ1

t... Rσi−1

tA(σi) + E(σi

ε)dσi...dσ1

−P∞

i=2 Rt+εT

tA(σ1) + ¯

ERσ1

t... Rσi−1

tA(σi) + ¯

E dσi...dσ1

(A.7)

Therefore, a bound for H(t+εT, t) can be determined as

∥H(t+εT)∥ ≤ 2(eεT α −1−εTα) (A.8)

where

α≡sup

t≥0

(max(∥A(t) + ¯

E∥,∥A(t) + E(t/ε)∥)) (A.9)

Substituting ΦE= Φ ¯

E+Hin Eq. A.4 yields

∆v(z(t+εT), t) = zT(t)(ΦT

E(t+εT, t)Q(t+εT)

Φ¯

E(t+εT, t)−Q(t))z(t) + zT(t)(ΦT

E(t+εT, t)Q(t+εT)

H(t+εT, t) + HT(t+εT, t)Q(t+εT )ΦT

E(t+εT, t)+

HT(t+εT, t)Q(t+εT)H(t+εT, t))z(t)

(A.10)

To obtain an upper bound for ∆v(z(t+εT), t), we

should use some relationships which are derived from

Eqs. A.1 and A.2, as

∥Q(t)∥ ≤ ρ(A.11)

∥ΦT

E(t, t0)∥ ≤ rρ

ηe−µ

2ρ(t−t0)(A.12)

v(x(t), t)≤e−µ

2ρ(t−t0)v(x(t0), t0), t ≥t0(A.13)

The first term of Eq. A.10, with considering x(t) = z(t)

as the initial condition, can be written as

zT(t)(ΦT

E(t+εT, t)Q(t+εT)Φ ¯

E(t+εT, t)

−Q(t))z(t) = v(x(t+εT), t +εT)−v(x(t), t)(A.14)

and, using Eqs. A.13 and A.1 yields

v(x(t+εT), t +εT)−v(x(t), t)≤(e−µεT

ρ−1)v(x(t), t)

≤ρ(e−µεT

ρ−1)∥x(t)∥2

(A.15)

Therefore,

zT(t)(ΦT

E(t+εT, t)Q(t+εT)Φ ¯

E(t+εT, t)

−Q(t))z(t)≤ρ(e−µεT

ρ−1)∥z(t)∥2(A.16)

With combining Eqs. A.8, A.11, A.12, and A.16, the

upper bound can be obtained as

∆v(z(t+εT), t)≤

(ρ(e−µεT

ρ−1) + 4ρ(√ρηe−µεT

2ρ)(eεT α −1−εTα)

+4ρ(eεT α −1−εTα)2)∥z(t)∥2

(A.17)

By describing the right-hand side of A.17 by a continu-

ously differentiable function g(ε, x), we have g(0, z) = 0,

and ∂

∂ε g(0, z) = −µT∥z∥2≤0. Then, since for ε→ ∞,

g(ε, z)→ ∞, there is an ε∗such that for 0 < ε < ε∗

and z= 0, g(ε∗, z) = 0, and g(ε, z)<0. Therefore,

∆v(z(t+εT), t) is negative definite.

Now it is shown that the negative definiteness of ∆v

results in the stability of ˙z(t)=(A(t) + E(t/ε))z(t).

With finding εand γ > 0 to satisfy

∆v(z(t0+εT), t0) =

v(z(t0+εT), t0+εT)−v(z(t0), t0)≤ −γ∥z(t0)∥2

(A.18)

and rewriting A.1 as v(z(t0), t0)≤ρ∥z(t0)∥2, we obtain

v(z(t0+εT), t0+εT)−v(z(t0), t0)≤ −(γ/ρ)v(z(t0), t0)

(A.19)

or v(z(t0+εT), t0+εT)≤(1 −γ/ρ)v(z(t0), t0) which

yields v(z(t0+kεT), t0+kεT )≤(1 −γ/ρ)kv(z(t0), t0)

for positive integer k.

Therefore, for k→ ∞,v(z(t0+kεT), t0+kεT)→0 and

then z(t0+kεT)→0. Thus, the exponential stability of

˙z(t)=(A(t) + E(t/ε))z(t) is proved.

[1] S. Boccaletti, V. Latora, Y. Moreno, M. Chavez, and D.-

U. Hwang, Complex networks: Structure and dynamics,

Phys. Rep. 424, 175 (2006).

[2] A. Arenas, A. D´ıaz-Guilera, J. Kurths, Y. Moreno, and

C. Zhou, Synchronization in complex networks, Phys.

Rep. 469, 93 (2008).

[3] S. Boccaletti, A. N. Pisarchik, C. I. Del Genio, and

A. Amann, Synchronization: from coupled systems to

complex networks (Cambridge University Press, 2018).

[4] R. Berner, Patterns of synchrony in complex networks of

adaptively coupled oscillators (Springer Nature, 2021).

[5] L. M. Pecora and T. L. Carroll, Master stability functions

for synchronized coupled systems, Phys. Rev. Lett. 80,

2109 (1998).

[6] X. Wang, Y.-C. Lai, and C. H. Lai, Enhancing synchro-

nization based on complex gradient networks, Phys. Rev.

E75, 056205 (2007).

[7] T. Nishikawa, A. E. Motter, Y.-C. Lai, and F. C. Hop-

pensteadt, Heterogeneity in oscillator networks: Are

smaller worlds easier to synchronize?, Phys. Rev. Lett.

91, 014101 (2003).

[8] Y. Gong, B. Wang, and H. Xie, Spike-timing-dependent

plasticity enhanced synchronization transitions induced

by autapses in adaptive newman-watts neuronal net-

works, Biosyst. 150, 132 (2016).

[9] S. Rakshit, A. Ray, B. K. Bera, and D. Ghosh, Synchro-

nization and firing patterns of coupled rulkov neuronal

map, Nonlinear Dyn. 94, 785 (2018).

[10] J. Fell and N. Axmacher, The role of phase synchroniza-

tion in memory processes, Nature Rev. Neurosci. 12, 105

(2011).

[11] P. Fries, J. H. Reynolds, A. E. Rorie, and R. Desimone,

Modulation of oscillatory neuronal synchronization by se-

lective visual attention, Science 291, 1560 (2001).

[12] C. Hammond, H. Bergman, and P. Brown, Pathological

synchronization in parkinson’s disease: networks, models

and treatments, Trends Neurosci. 30, 357 (2007).

[13] W. Lin and G. Chen, Using white noise to enhance

synchronization of coupled chaotic systems, Chaos 16,

013134 (2006).

[14] A. Zeng, S.-W. Son, C. H. Yeung, Y. Fan, and Z. Di,

Enhancing synchronization by directionality in complex

networks, Phys. Rev. E 83, 045101 (2011).

[15] S.-W. Son, B. J. Kim, H. Hong, and H. Jeong, Dynamics

and directionality in complex networks, Phys. Rev. Lett.

103, 228702 (2009).

[16] J. P. Ramirez, A. Arellano-Delgado, and H. Nijmeijer,

Enhancing master-slave synchronization: The effect of

using a dynamic coupling, Phys. Rev. E 98, 012208

(2018).

[17] R. Banerjee, B. K. Bera, D. Ghosh, and S. K. Dana, En-

hancing synchronization in chaotic oscillators by induced

heterogeneity, Eur. Phys. J. Spec. Top. 226, 1893 (2017).

[18] R. Sevilla-Escoboza, R. Guti´errez, G. Huerta-Cuellar,

S. Boccaletti, J. G´omez-Garde˜nes, A. Arenas, and

J. Buld´u, Enhancing the stability of the synchroniza-

tion of multivariable coupled oscillators, Phys. Rev. E

92, 032804 (2015).

[19] S. Panahi, F. Nazarimehr, S. Jafari, J. C. Sprott,

M. Perc, and R. Repnik, Optimal synchronization of cir-

culant and non-circulant oscillators, Appl. Math. Com-

put. 394, 125830 (2021).

[20] E. Sch¨oll, Synchronization in delay-coupled complex net-

works, in Advances in analysis and control of time-

delayed dynamical systems (World Scientific, 2013) pp.

57–84.

[21] J. Sawicki, Delay controlled partial synchronization in

complex networks (Springer Nature, 2019).

[22] Y. Kyrychko, K. Blyuss, and E. Sch¨oll, Synchronization

of networks of oscillators with distributed delay coupling,

Chaos 24, 043117 (2014).

[23] A. Gjurchinovski, A. Zakharova, and E. Sch¨oll, Ampli-

tude death in oscillator networks with variable-delay cou-

pling, Phys. Rev. E 89, 032915 (2014).

[24] T. Liu, J. Zhao, and D. J. Hill, Exponential synchro-

nization of complex delayed dynamical networks with

switching topology, IEEE Trans. Circuit. Syst. I 57, 2967

(2010).

[25] S. Rakshit, S. Majhi, B. K. Bera, S. Sinha, and D. Ghosh,

Time-varying multiplex network: Intralayer and inter-

layer synchronization, Phys. Rev. E 96, 062308 (2017).

[26] I. V. Belykh, V. N. Belykh, and M. Hasler, Blinking

model and synchronization in small-world networks with

a time-varying coupling, Physica D 195, 188 (2004).

[27] M. Hasler and I. Belykh, Blinking long-range connections

increase the functionality of locally connected networks,

IEICE T. on Fund. Electr. 88, 2647 (2005).

[28] M. Hasler, V. Belykh, and I. Belykh, Dynamics of

stochastically blinking systems. part i: Finite time prop-

erties, SIAM J. Appl. Dyn. Syst. 12, 1007 (2013).

[29] M. Hasler, V. Belykh, and I. Belykh, Dynamics of

stochastically blinking systems. part ii: Asymptotic

properties, SIAM J. Appl. Dyn. Syst. 12, 1031 (2013).

[30] M. Porfiri, Stochastic synchronization in blinking net-

works of chaotic maps, Phys. Rev. E 85, 056114 (2012).

[31] W. Lu, F. M. Atay, and J. Jost, Synchronization of

discrete-time dynamical networks with time-varying cou-

plings, SIAM J. Math. Anal. 39, 1231 (2008).

[32] R. Jeter and I. Belykh, Synchronization in on-off stochas-

tic networks: Windows of opportunity, IEEE Trans. Cir-

cuits Syst. I 62, 1260 (2015).

[33] R. Jeter and I. Belykh, Dynamical networks with on-off

stochastic connections: Beyond fast switching, in Proc.

IEEE Int. Symp. Circuits Syst. (ISCAS-2014) (IEEE,

2014) pp. 1788–1791.

[34] N. Barabash and V. Belykh, Synchronization thresholds

in an ensemble o kuramoto phase oscillators with ran-

domly blinking couplings, Radiophys. Quantum El. 60,

761 (2018).

[35] B. Liu, W. Lu, and T. Chen, Synchronization in complex

networks with stochastically switching coupling struc-

tures, IEEE Trans. Autom. Control 57, 754 (2012).

[36] M. Porfiri and I. Belykh, Memory matters in synchro-

nization of stochastically coupled maps, SIAM J. Appl.

Dyn. Syst. 16, 1372 (2017).

[37] V. Lanza, J. Secco, and F. Corinto, Blinking networks of

memristor oscillatory circuits in the flux-dharge domain,

Front. Neurosci. 15, 618607 (2021).

[38] A. Ascoli, V. Lanza, F. Corinto, and R. Tetzlaff, Synchro-

nization conditions in simple memristor neural networks,

J. Franklin Inst. 352, 3196 (2015).

[39] T. Stankovski, Time-varying coupling functions: dynam-

ical inference and cause of synchronization transitions,

Phys. Rev. E 95, 022206 (2017).

[40] T. Stankovski, T. Pereira, P. V. McClintock, and A. Ste-

fanovska, Coupling functions: universal insights into dy-

namical interaction mechanisms, Rev. Mod. Phys. 89,

045001 (2017).

[41] B. Kralemann, M. Fr¨uhwirth, A. Pikovsky, M. Rosen-

blum, T. Kenner, J. Schaefer, and M. Moser, In vivo car-

diac phase response curve elucidates human respiratory

heart rate variability, Nature Commun. 4, 2418 (2013).

[42] Z. Li, L. Jiao, and J.-J. Lee, Robust adaptive global syn-

chronization of complex dynamical networks by adjust-

ing time-varying coupling strength, Physica A 387, 1369

(2008).

[43] V. Chandrasekar, J. H. Sheeba, B. Subash, M. Laksh-

manan, and J. Kurths, Adaptive coupling induced multi-

stable states in complex networks, Physica D 267, 36

(2014).

[44] Q. Ren and J. Zhao, Adaptive coupling and enhanced

synchronization in coupled phase oscillators, Phys. Rev.

E76, 016207 (2007).

[45] Z. Hagos, T. Stankovski, J. Newman, T. Pereira, P. V. E.

McClintock, and A. Stefanovska, Synchronization transi-

tions caused by time-varying coupling functions, Philos.

Trans. R. Soc. A 377, 20190275 (2019).

[46] T. Stankovski, A. Duggento, P. V. McClintock, and

A. Stefanovska, Inference of time-evolving coupled dy-

namical systems in the presence of noise, Phys. Rev. Lett.

109, 024101 (2012).

[47] F. Stelzer, A. R¨ohm, R. Vicente, I. Fischer, and

S. Yanchuk, Deep neural networks using a single neuron:

folded-in-time architecture using feedback-modulated de-

lay loops, Nature Commun. 12, 1 (2021).

[48] L. Huang, Q. Chen, Y.-C. Lai, and L. M. Pecora, Generic

behavior of master-stability functions in coupled nonlin-

ear dynamical systems, Phys. Rev. E 80, 036204 (2009).

[49] D. J. Stilwell, E. M. Bollt, and D. G. Roberson, Suf-

ficient conditions for fast switching synchronization in

time-varying network topologies, SIAM J. Appl. Dyn.

Syst. 5, 140 (2006).

[50] S. Rakshit, B. K. Bera, D. Ghosh, and S. Sinha, Emer-

gence of synchronization and regularity in firing patterns

in time-varying neural hypernetworks, Phys. Rev. E 97,

052304 (2018).

[51] L. Zhang, Y. Wang, and Q. Wang, Synchronization for

time-varying complex dynamical networks with different-

dimensional nodes and non-dissipative coupling, Com-

mun. Nonlinear Sci. Numer. Simul. 24, 64 (2015).

[52] S. Rakshit, S. Majhi, J. Kurths, and D. Ghosh, Neu-

ronal synchronization in long-range time-varying net-

works, Chaos 31, 073129 (2021).

[53] S. Chakravartula, P. Indic, B. Sundaram, and T. Killing-

back, Emergence of local synchronization in neuronal net-

works with adaptive couplings, PLoS One 12, e0178975

(2017).

[54] D. Cumin and C. Unsworth, Generalising the kuramoto

model for the study of neuronal synchronisation in the

brain, Physica D 226, 181 (2007).

[55] A. Destexhe and E. Marder, Plasticity in single neuron

and circuit computations, Nature 431, 789 (2004).

[56] A. Stefa´nski, J. Wojewoda, T. Kapitaniak, and

S. Yanchuk, Simple estimation of synchronization thresh-

old in ensembles of diffusively coupled chaotic systems,

Phys. Rev. E 70, 026217 (2004).

[57] A. Pikovsky, On the interaction of strange attractors, Z

Phys B: Condens Matter 55, 149 (1984).

[58] J. Ladenbauer, J. Lehnert, H. Rankoohi, T. Dahms,

E. Sch¨oll, and K. Obermayer, Adaptation controls syn-

chrony and cluster states of coupled threshold-model neu-

rons, Phys. Rev. E 88, 042713 (2013).

[59] J. Sawicki, I. Omelchenko, A. Zakharova, and E. Sch¨oll,

Delay controls chimera relay synchronization in multiplex

networks, Phys. Rev. E 98, 062224 (2018).