Periodic coupling suppresses synchronization in coupled phase oscillators

The responses of complex dynamical systems to external perturbations are of both theoretical importance and practical significance. Different from the conventional studies in which perturbations are added on to the system variables, here we consider the scenario of perturbed couplings. Specifically, by the scheme of periodic coupling in which the coupling strength varies periodically with time, we investigate how the synchronization behaviors of an ensemble of phase oscillators are affected by the properties of periodic coupling, including the coupling frequency and amplitude. It is found that, by comparison with the constant coupling scenario, the presence of period coupling will suppress synchronization in general and, with a decrease in the frequency (or the increase of the amplitude) of the periodic coupling, the synchronization performance gradually deteriorates. The influences of periodic coupling on synchronization are demonstrated by numerical simulations in different network models, and the underlying mechanism is analyzed by the method of dimension reduction.


Introduction
The fact that realistic systems are inevitably influenced by their surrounding environments and the finding that the dynamics of nonlinear systems is sensitive to small perturbations make the dynamical responses of complex nonlinear systems to external perturbations an outstanding problem in science and engineering. In biological and neural systems, a question of particular interest is how the system collective behaviors are influenced by externally added periodic signals, e.g. the effects of seasonal variations on ecological processes [1], the influences of circadian rhythm on neural activities [2], the impacts of diurnal cycles on gene expressions [3], and the influences of electrical and magnetic stimulations on brain functions [4]. In particular, with the development of modern techniques in medicine and neuroscience (e.g. the transcranial electrical/magnetic stimulations), in recent years tremendous efforts have been paid to the treatment of neurological disorders, or the improvement of brain functions, by noninvasive electrical and magnetic simulations [5,6]. Despite the significant progresses achieved in experiment, the mechanisms behind the interaction between external simulations and neural activities remains mysterious [4][5][6].
A typical collective behavior observed in complex systems is synchronization, in which an ensemble of oscillators are observed to be moving in a coherent fashion when the coupling between them exceeds certain critical strength [7,8]. As the dynamical basis for the function and operation of many realistic systems, synchronization behaviors have been extensively studied by researchers from different fields over the past decades [9][10][11]. In exploring oscillator synchronization, one of the focusing issues is about the robustness of synchronization to external perturbations, including environmental noises and periodic drivings [12][13][14][15][16][17][18][19][20][21]. Due to the nonlinear feature of the system dynamics, intriguing phenomena could be generated in the presence of external perturbations, e.g. the enhancement of synchronization by random noise [19,20], the control of synchronization by periodic signals [14,15], the manipulation of synchronization by periodic drivings [16], and the inducing of synchronization by configuring initial conditions [21]. It is noted that in these studies, the perturbation signals are normally added on the system variables. While this setting is appropriate for many Any further distribution of this work must maintain attribution to the author(s) and the title of the work, journal citation and DOI. realistic situations, there do exist situations where external perturbations are acting on the couplings. For instance, in brain stimulations, the electromagnetic signals affect also the conductance of the synapses, leading to the plasticity of neural connections [5,6]. Similar concerns exist also in ecosystems (where the interactions between species are varying with seasons) and metabolic systems (where gene expressions are modulated by diurnal cycles) [1,3]. With these concerns, a question raises naturally is: how the synchronization behaviors of coupled oscillators are affected by periodic couplings? Inspired by the above question, we investigate in the present work the impacts of periodic coupling on the synchronization of coupled phase oscillators. More specifically, by adopting periodic coupling in the classical Kuramoto model [22], we investigate how the degree of synchronization is affected by the frequency and amplitude of the periodic coupling. We are able to demonstrate numerically and argue theoretically that the adoption of periodic coupling will suppress synchronization in general, with the suppression degree dependent of the coupling frequency and amplitude. Specifically, the smaller (larger) the coupling frequency (amplitude) is, the stronger the synchronization suppression will be. In addition, it is found that the central frequency of the oscillators, which characterizes the time scale of the synchronous motion, has no influence on synchronization. These findings shed lights on the responses of complex dynamical systems to external drivings, and provides a timely contribution to the study of stimulation-related system functions, e.g. brain stimulations [4][5][6].
In section 2, we will present the model of coupled phase oscillators and investigate numerically the influences of periodic coupling on synchronization. By the method of dimension reduction, in section 3 we will conduct a theoretical analysis on the impacts of periodic coupling on synchronization, and explore the underlying mechanism by the scheme of on-off coupling. Generalizations to other network models, including the random complex network and the neural network of nematode C. elegans, will be presented in section 4. Discussions and conclusion will be given in section 5.

Model and numerical results
Our model of coupled phase oscillators reads Here, i,j=1, K, N are the oscillator (node) indices, θ i (t) is the instant phase of the ith oscillator at time t, and ω i is the natural frequency of the ith oscillator. The natural frequency follows the distribution g(ω). The coupling relationship of the oscillators is captured by the adjacency matrix = { } a A ij , with = = a a 1 ij ji if oscillators i and j are directly connected by a link, otherwise a ij =0. The quantity = å d a i j ij denotes the number of connections associated to oscillator i, i.e. its degree. Different from existing models, we make the coupling strength vary periodically with time, as described by equation (2). The frequency and amplitude of the periodic coupling are represented by f and K 0 , respectively. If f=0, we have K=K 0 and the model is identical to the generalized Kuramoto model with constant coupling [9,11]. Furthermore, if the oscillators are globally connected (i.e. a ij =1 for all the non-diagonal elements in matrix A), equation (1) describes the classical Kuramoto model. The main task of our present work is to investigate how the variations of f and K 0 will affect the synchronization degree of the oscillators.
We start by demonstrating the impacts of periodic coupling on synchronization based on numerical simulations. To make the model theoretical tractable, we employ the structure of global coupling (the results for complex network structures will be discussed later). Following [23,24], we choose the natural frequency ω i of the oscillators randomly from the Lorentzian distribution with ω 0 the central frequency and Δ the scale parameter of the distribution. (We note that the main phenomena to be reported are independent of the specific form of the frequency distribution, which has been checked numerically for the Gaussian, uniform and triangle frequency distributions.) Without loss of generality, we set ω 0 =0 and Δ=1. The initial phases of the oscillators are randomly chosen within the range (0, 2π). The degree of synchronization is evaluated by the order parameter where |·| and á ñ · denote, respectively, the module and time averaged functions, and R(t) is the instant order parameter. We have Î [ ] R 0, 1 , with R=0 and 1 corresponding to the completely incoherent and global synchronization states, respectively. In numerical simulations, we fix the system size as N=1×10 3 , and update equations (1) and (2) by the fourth-order Runge-Kutta algorithm, with the time step being chosen as δt=5×10 −2 . In obtaining the order parameter, R, we first evolve the system for a transient period of T tr =6×10 3 , and then average the instant order parameter, R(t), over a period of T=200.
By numerical simulations, we plot in figure 1(a) the variation of R with respect to K 0 for the special case of constant coupling ( f = 0) and for two typical cases of periodic coupling ( f=0.2 and 2.0). We see that when K 0 is small ( , the values of R are almost identical for all the cases. However, as K 0 increases from K c , the value of R for the case of low-frequency periodic coupling ( f=0.2) is clearly smaller to that of constant coupling, signifying that synchronization is suppressed by low-frequency couplings. Furthermore, as K 0 increases from K c , the gap between the two curves is gradually enlarged. For instance, we have However, for the case of high-frequency periodic coupling ( f=2.0), we see that the value of R is almost identical to that of constant coupling for the whole region, indicating that synchronization is less affected by high-frequency couplings.
To investigate systematically the influence of coupling frequency on synchronization, we plot in figure 1(b) the variation of R with respect to f for different K 0 . We see that for small coupling amplitude (K 0 =1), despite the variation of f, the value of R is staying always around 0. This observation is consistent with the finding in figure 1(a) that R is independent of f in the region of K 0 <K c . For larger coupling amplitudes (K 0 =3 and 5), we see that as f increases from 0, the value of R is suddenly decreased first, then increases gradually and approaches finally the value generated by constant coupling = [ ( ) R f 0 ]. This observation is also consistent with the finding in figure 1(a) that synchronization is more suppressed by low-frequency than high-frequency couplings.
In previous studies on the synchronization of periodically coupled chaotic oscillators, an interesting finding is that when the time scale of the nodal dynamics is comparable to that of periodic coupling, the synchronizability of the system could be significantly affected [25,26]. In particular, it is shown that at some characteristic frequencies of the periodic coupling, the upper boundary of the synchronization region in the parameter space of the coupling strength can be significantly extended. For our model of periodically coupled phase oscillators, the averaged time scale of the nodal dynamics is characterized by the central frequency, ω 0 , in the Lorentzian distribution. This arouses our interest of the influence of ω 0 on the synchronization of periodically coupled phase oscillators. Changing ω 0 to 1.0 (ω 0 ≈2πf ) and 0.4π (ω 0 =2πf ), we plot in figure 1(a) again the variation of R with respect to K 0 for the case of low-frequency coupling ( f=0.2). We see that the results of ω 0 =1.0 and 0.4π are identical to that for ω 0 =0. Numerical results thus suggest that, unlike periodically coupled chaotic oscillators, the synchronization performance of periodically coupled phase oscillators is independent of the time scale of the nodal dynamics (the central frequency of the Lorentzian distribution). The independence of R on ω 0 is also observed in the variation of R with respect to f, as depicted in figure 1(b) (the case of K = 5 and ω 0 =1.0). (As predicted in [26], the characteristic frequencies observed in periodically coupled chaotic oscillators lies in the nonlinear response of the temporal stability of the local dynamics, which are absent in coupled one-dimensional oscillators. Our present study confirms this prediction.) In numerical studies of Kuramoto-type models, an important concern is whether the observed phenomena is dependent of the system size, namely the finite-size effect [27,28]. To check the size-effect, we plot in figure 2 the variation of R with respect to N for different values of K 0 and f. We see that for small coupling amplitudes , as N increases, the value of R is decreased to 0 by roughly the power-law scalings (figures 2(a)-(d)); while for large coupling amplitude (K 0 >K c ), the value of R is size-independent (figures 2(e) and (f)). In particular, comparing with the results of f=0 (figure 2(e)) and f=0.2 (figure 2(f)), we see that the value of R is decreased by nearly the same amount (ΔR≈0.1) for different values of N. Numerical simulations thus suggest that the phenomenon we have observed, i.e. synchronization is suppressed by periodic coupling, is independent of the system size.
To gain insights on the influence of periodic coupling on synchronization, we continue to study the temporal behavior of the order parameter. By the parameter K 0 =1, we plot in figure 3(a) the time evolution of R(t) for the cases of constant, low-frequency ( f=0.2) and high-frequency ( f = 2) couplings. We see that as time increases, the value of R(t) is quickly decreased and is staying around 0 with small fluctuations for all three cases. Increasing K 0 to 5, we plot in figure 3(b) the time evolution of R(t) again. We see that, different from the cases of small coupling amplitude (figure 3(a)), the time evolution of R(t) is clearly affected by the coupling frequency. Specifically, for the case of constant coupling, the value of R(t) is increased with time gradually and saturated at about 0.75 with almost no fluctuation; for the case of low-frequency coupling ( f=0.2), R(t) is oscillating periodically with a large amplitude, with the oscillating frequency identical to that of periodic coupling; for the case of high-frequency coupling ( f=2.0), R(t) is also varying periodically, but with a smaller amplitude. We thus infer from figure 3 that under large-amplitudes periodic coupling (K>K c ) the instant order parameter, R (t), becomes oscillatory, with the oscillating frequency identical to that of periodic coupling and the oscillating amplitude dependent of the coupling frequency.

Theoretical analysis
We move on to explore the dynamical mechanism underlying the observed phenomena, based on the method of dimension reduction proposed by Ott and Antonsen [23,24]. In the case of thermodynamic limit, i.e.  ¥ N in equation (1), the state of the system at time t can be described by a continuous probability density function ρ (θ, ω, t), with r q w q w ( ) t , , d d the fraction of oscillators with phases between θ and θ+dθ and natural frequencies between ω and ω+dω. The evolution of ρ(θ, ω, t) is governed by the continuity equation with * u q w w = + - 2i e e , 5 Utilizing the Ott-Antonsen ansatz [23,24]: , , n n and substituting equation (7) into equations (5) and (6), we obtain To evaluate the integral that gives *( ) Z t , we analytically continue ω to the complex ω-plane and carry out contour integration [23,24], where the analyticity of a holds in the lower half plane of the complex variable ω. For large negative values of ( ) a Im , equation (8) can be approximated as . Following the setting of numerical simulations, we have  We continue to analyze the case of periodic coupling. When periodic coupling is adopted, equation (10) can not be solved analytically. We thus turn to numerical simulations. By solving equation (10) numerically, we plot in figure 1 the variation of R with respect to K 0 for different values of f. We see that the theoretical results given by equation (10) are in good agreement with the results of direct simulations. Furthermore, equation (10) indicates that the evolution of R(t) is independent of the central frequency, ω 0 , which is also in consistent with the numerical observation in figure 1 (i.e. the changes of ω 0 from 0 to 1 and 0.4π do not affect the values of R). Equation (10) can be also used to trace the time evolution of R(t), as depicted in figure 3. We see that the theoretical results calculated from equation (10) fit the numerical results very well for both the cases of lowfrequency and high-frequency couplings.
Having justified the efficiency of the method of dimension reduction in characterizing the dynamics of the synchronization order parameter, we next exploit this method to investigate the mechanism behind the observed phenomenon, i.e. synchronization is suppressed by periodic coupling. To facilitate the analysis, here we consider the simple case of square wave coupling, the sign function. That is, the coupling strength is alternating between 0 and 2K 0 periodically with the same duration. Square wave coupling is also known as on-off coupling [25], which has been widely used in literature as approximations of pulsed couplings (in neural systems) and diurnal-cycle-based interactions (in ecosystems). To be consistent with numerical studies, we set Δ=1 in the following analysis.
We first analyze the situation of weak coupling amplitude, K 0 =1. In this situation, the coupling strength, K (t), is alternating between 2 (the 'on' state) and 0 (the 'off' state), and equation (10) is simplified to dR/dt=−R and dR/dt=−R 3 for K(t)=0 and 2, respectively. Figure 4(a1) shows the phase portraits of the two states. As dR/dt<0 for both states, we see that R(t) approaches the stable point R s =0 monotonically with time for both the cases of low-frequency and high-frequency couplings ( figure 4(a2)). We next analyze the situation of critical coupling amplitude, K 0 =2. In this situation, the 'on' state is governed by the equation = -R t R R d d 2 3 which has the stable point = R 1 2 s (figure 4(b1)). As time increases, the value of R(t) is approaching 1 2 and 0 during the episodes of 'on' and 'off' states, respectively. This leads to the oscillation of R(t), as depicted in figure 4(b2). Since dR/dt<R during the whole process (i.e. the origin is always more stable than R s ), we see that R(t) is finally attracted to 0. Like the situation of weak coupling ( figure 4(a2)), here the role of f is also reflected in affecting the oscillating amplitude of R(t): the smaller the coupling frequency is, the larger the oscillating amplitude will be.
We go on to analyze the situation of strong coupling amplitude, K 0 =5, with which the coupling is switching between 0 (the 'off' state) and 10 (the 'on' state). The dynamics of R(t) in the 'on' state is governed by the equation = -R t R R d d 4 5 3 , which has the stable point = R 0.8 s ( figure 4(c1)). Setting dR/dt=R, we have the equilibrium point = R 0.6 e where the attracting of the origin and R s is balanced with each other. Due to the equilibrium point R e , the value of R(t) is no longer attracted to the origin like the case of weak coupling (figures 4(a1) and (a2)), but will be oscillating around R e ( figure 4(c2)). Please note that R e is also the stable point associated with the constant coupling K 0 =5. (As a matter of fact, the relation = ( ) R R K e s 0 holds for any coupling strength larger than 2K c .) So, under the periodic coupling, the value of R(t) is oscillating around R s (K 0 )-the order parameter generated by the constant coupling of strength K 0 . If the coupling frequency is high enough, e.g. f=2, the value of R(t) will be trapped in the vicinity of R e . As a result of this, R(t) is oscillating around R e roughly in a symmetric fashion, resulting in R( f=2)≈R( f=0). This explains why synchronization is less affected by high-frequency couplings. The symmetric oscillation, however, is broken when the coupling frequency is low, e.g. f=0.2. As depicted in figure 4(c2), in this case the oscillation of R(t) is truncated at R s during the 'on' episode, but is approaching 0 continuously during the 'off' episode. It is just the asymmetric oscillation of R(t) that results in the suppressed synchronization under low-frequency couplings.
With the decrease of f (or the increase of K 0 ), the asymmetric feature of the oscillation becomes more prominent, making the synchronization order parameter gradually decreased. This understanding explains the numerical results shown in figure 1.

Generalizations
We finally check the generality of the observed phenomena in systems of complex network structures. The first system we investigate is the Erdös-Renyi (ER) random network [29]. The ER network is generated by connecting nodes in the system randomly with the probability p=0.6. Besides the structure, the other settings of the ER network are the same to that of the globally connected network studied in figures 1 and 3, including the network size, the frequency distribution, and the coupling forms figure 5(a) shows the variation of R with respect to K 0 for the cases of constant coupling ( f=0) and low-frequency coupling ( f=0.2) with on-off and sinusoidal coupling functions. We see that before the onset of synchronization (which occurs at about K c =3), the values of R for all three cases are very close; after the onset point, the value of R for the case of f=0.2 is clearly below to that of constant coupling and, with the increase of K 0 , the difference between the two cases is gradually enlarged. These results are consistent with the results obtained in globally connected network (figures 1 and 3). It should be note that due to the sparse connectivity of the ER network, the dynamics of the synchronization order parameter is not precisely described by equation (10) [23,24]. However, the dynamical mechanism we have revealed ( figure 4) is general, which is independent of the network structure. Figure 5(a) shows also the variation of R with respect to K 0 for the case of on-off periodic coupling ( f=0.2). We see that comparing to the case of sinusoidal coupling, network synchronization is more suppressed by on-off coupling.
The second system we investigate is the neural network of the nematode C. elegans, which consists of 297 nodes and totally 2148 links [30]. The natural frequencies of the oscillators are still chosen randomly from the Lorentzian distribution, with the parameters identical to those used in the model of globally connected network (figures 1 and 2). Figure 5(b) shows the variation of R with respect to K 0 for the cases of constant, sinusoidal, and on-off couplings. We see that, similar to the results of ER network ( figure 5(a)), the values of R are diverged from each other after the onset of synchronization (which occurs at about K c =30); and, comparing to the case of constant coupling, synchronization is suppressed by periodic couplings.

Discussions and conclusion
As the dynamical basis for the normal functioning of brain, the synchronization activities of distributed neurons has long interested neuroscientists [31]. In particular, aiming to improve the brain functions or recover brain from disorders, various stimulating techniques have been proposed in the past decades [4,5], saying, for instance, transcranial direct-current stimulations, transcranial magnetic stimulations, and deep-brain stimulations. Whereas the efficacy of these techniques has been testified by neurological experiments and in clinical treatments, the underlying mechanism remains not clear. In recent years, stimulated by the rapid progress in complex dynamical systems, there have been theoretical efforts to investigate the responses of neural system to stimulations, in which different models consisting of coupled phase oscillators have been proposed and studied [31]. In conventional studies, the external stimulation is normally introduced to the system variables, with the concern that the activities of ion channels are directly affected by electrical/magnetic signals.
Recent evidences from pharmacological and neurological studies, however, show that low-frequency stimuli could also increase the expression of certain proteins (e.g. NMDA receptors) at the synapses [5], which in turn could change the molecular and cellular pathology associated with some neurological disorders (e.g. Alzheimer diseases) [6]. Inspired by these experimental findings, we have proposed in the present work the scheme of timedependent coupling and investigated, numerically and theoretically, the influences of periodic coupling on synchronization performance. The finding that synchronization is suppressed by periodic coupling sheds new lights on the dynamical responses of neural systems to external stimuli, and might have implications to the therapeutic interventions of neurological disorders such as epileptic seizures, Alzheimer, and Parkinson diseases.
Summarizing up, motivated by the fact that the interaction strength between elements in some realistic systems are subjected to periodic perturbations, we have studied the synchronization behavior of coupled phase oscillators under periodic couplings. It is found that comparing to the situation of constant coupling, the synchronization propensity of the system is deteriorated under periodic coupling and, with the decrease of the frequency of the coupling, the synchronization degree is gradually decreased. By numerical simulations, we have demonstrated this phenomenon in different network models, and justified its independence to the system size. Furthermore, by the method of dimension reduction, we have conducted a detailed analysis on the influence of periodic coupling on synchronization order parameter, with the theoretical predictions in good agreement with the numerical results. In particular, by the scheme of on-off periodic coupling, we have analyzed the mechanism underlying the observed phenomenon, and found that the suppressed synchronization by periodic coupling is attributed to the asymmetric oscillation of the temporal synchronization order parameter. The finding extends our knowledge on the dynamical responses of complex nonlinear systems to external perturbations, and