Synchronous bursts on scale-free neuronal networks with attractive and repulsive coupling

This paper investigates the dependence of synchronization transitions of bursting oscillations on the information transmission delay over scale-free neuronal networks with attractive and repulsive coupling. It is shown that for both types of coupling, the delay always plays a subtle role in either promoting or impairing synchronization. In particular, depending on the inherent oscillation period of individual neurons, regions of irregular and regular propagating excitatory fronts appear intermittently as the delay increases. These delay-induced synchronization transitions are manifested as well-expressed minima in the measure for spatiotemporal synchrony. For attractive coupling, the minima appear at every integer multiple of the average oscillation period, while for the repulsive coupling, they appear at every odd multiple of the half of the average oscillation period. The obtained results are robust to the variations of the dynamics of individual neurons, the system size, and the neuronal firing type. Hence, they can be used to characterize attractively or repulsively coupled scale-free neuronal networks with delays.


Introduction
It is well known that synchronization in neuronal networks is particularly relevant for the efficient processing and transmission of information (see e.g. [1,2]). Experiments have shown that synchronized states can occur in many special areas of the brain, such as the olfactory system or the hippocampal region [3][4][5]. By using functional magnetic resonance imaging (fMRI) to record brain activity from both speakers and listeners during natural verbal communication, a recent study has shown that speakerlistener neural coupling underlies successful communication by means of synchronization [6]. Theoretically, neuronal synchronization on complex networks has been explored in detail [7][8][9][10][11][12][13], leading to several insights that have applicability on real problems in neuroscience. For example, synchronization of gap-junctioncoupled neurons has been investigated [11], and by means of the phase resetting curve, phase locking and synchronization in neuronal networks have been investigated as well [14,15]. Moreover, noise-induced and noise-enhanced synchronization have also been reported in realistic neuronal systems [16,17]. Interestingly, it was reported that chemical and electrical synapses perform complementary roles in the synchronization of interneuronal networks [18]. Indeed, synchronization, information transmission and signal sensitivity on complex networks are currently hot topics in theoretical neuroscience [19,20], as evidenced by several recent studies that are devoted to the explorations of this subject [21][22][23][24][25][26][27][28][29][30][31][32].
Previous research highlighted that information transmission delays are inherent to the nervous system because of the finite speed at which action potentials propagate across neuron axons, as well as due to time lapses occurring by both dendritic and synaptic processing [33]. It has been reported, for example, that the beta frequency is able to synchronize over long conduction delays, which corresponds to signals traveling a significant distance in the brain [34]. Thus far, it has also been reported that different time delay lengths can change both qualitative as well as quantitative properties of the dynamics [35]. For example, delays can introduce or destroy stable oscillations, enhance or suppress synchronization, as well as generate complex spatiotemporal patterns on regular neuronal networks. It has also been suggested that time delays can facilitate neural synchronization and lead to many interesting and even unexpected phenomena [36,37], including zigzag fronts of excitations, clustering antiphase synchronization and in-phase synchronization [38]. Most recently, the synchronizability threshold for an arbitrary network incorporating delays and noise has been derived, and additionally, by means of the scaling theory of the underlying fluctuations, the absolute limit of synchronization efficiency in a noisy environment with uniform time delays has been established [39].
Both phase-attractive (which can be related to excitatory synapses) and phase-repulsive (which can be related to inhibitory synapses) coupling exists in realistic neuronal systems. Hence, it is important to take this explicitly into account in theoretical studies. Effects of phase-repulsive coupling on neuronal dynamics have also been investigated in the past [40][41][42], where such coupling was considered to be related to inhibitory synapses. For example, it has been shown that a pair of excitable FitzHugh-Nagumo neurons can exhibit various firing patterns including multistability and chaotic firing when elements interact phase-repulsively [40]. Moreover, the synchronization of nonidentical dynamical units that are coupled attractively in a small-world network can be improved significantly by the introduction of just a small fraction of phase-repulsive couplings [42]. Dynamics of propagation in coupled neuronal networks with excitatory and inhibitory synapses has been investigated in detail by means of integrate-and-fire neurons [43,44]. By analyzing a canard mechanism, it has also been shown that synaptic coupling can synchronize neurons at low firing frequencies [45]. However, synchronization on scale-free neuronal networks with phase-repulsive coupling and delay has not yet been investigated.
Here, we aim to extend the scope of research by studying the dependence of synchronization transitions on the information transmission delay over scale-free neuronal networks with attractive or repulsive coupling, respectively. Since a power law distribution of the degree of neurons has been found applicable for the coherence among activated voxels using functional magnetic resonance imaging [46], and moreover, the robustness against simulated lesions of anatomic cortical networks was also found to be most similar to that of a scale-free network [47], our study addresses a relevant system setup which is still amenable to new research. We report several non-trivial effects induced by finite (non-zero) delay lengths, foremost the ability of its fine-tuning towards highly synchronized fronts of excitations. We find that the delay-induced synchronization transitions manifest as well-expressed minima in the measure for spatiotemporal synchrony. Depending on the type of coupling, however, these minima appear every integer multiple of the average oscillation period of bursting oscillations in case of attractive coupling, or they appear every odd multiple of the half of the average oscillation period for repulsive coupling. The results are robust to variations of neuronal dynamics and system size, and appear to be primarily due to the emergence of phase locking between the delay and the time scales, which are inherent to each individual neuron constituting the scale-free network.

Results
Firstly, we present in Fig. 1 space-time plots to have a look at characteristic synchronization transitions that can be induced by different information transmission delays. To do so, we employ attractive coupling as a case of example, but note that qualitatively identical space-time plots can be obtained also for repulsive coupling. We set a~2:3, for which individual neurons exhibit simple single-burst excitations. Results presented in Fig. 1(a) indicate that the spatiotemporal dynamics is synchronous if t~0, which can be attributed to sufficiently strong attractive coupling. However, if the information transmission delay is increased to t~270 the synchrony deteriorates rather drastically, as can be observed in Fig. 1(b). Interestingly, synchronization seems again fully restored at t~850, as depicted in Fig. 1(c), but then again disappears for t~1290 and reappears for t~1700, as shown in Figs. 1(d) and (e), respectively. Indeed, we find that such a succession repeats itself for higher values of t, from which we conclude that the information transmission delay can either promote or impair synchronization of neuronal activity on scalefree networks. If inspecting the values of t warranting near-perfect synchronization closely, we can observe that they equal roughly integer multiples of 850, which hints towards an underlying mechanism that can explain our observations.
In order to investigate the impact of different values of t quantitatively, and separately for attractive and repulsive coupling, we calculate the synchronization parameter s as defined by Eq.
(3). Results presented in Figs. 2(a) and (b) were obtained for attractive coupling and three different values of a. It can be observed that certain values of t significantly facilitate spatiotemporal synchronization of excitatory fronts on neuronal scale-free networks. In particular, the three minima of s appear at t&850~T, t&1700~2T and t&2550~3T if a~2:3. For a~3:0, we can observe two minima of s appearing at t&1200~T and 2400~2T. Furthermore, several more minima can be observed for a~4:1 within the considered span of information transmission delays, as depicted in Fig. 2(b). Again it is clear that they appear at integer multiples of the first minimum. This confirms the fact that delay-induced transitions to spatiotemporally synchronized neuronal activity appear intermittently, at integer multiples of a given value of t. On the other hand, values of t outside these regions impair synchronization significantly, as can be inferred from the rather sharp ascends towards larger values of s beyond the optimal delays.
Performing the same analysis for repulsive coupling reveals several similarities, but also significant differences. Results presented in Figs. 3(a) and (b) indeed have a qualitatively identical outlay with the minima of s appearing intermittently as t increases, yet the precise values warranting optimal neuronal synchrony are different if compared to the case of attractive coupling. Specifically, the three minima of s appear at t&425~T=2, t&1275~3T=2 and t&2125~5T=2 if a~2:3, while for a~3:0 and a~4:1 we can observe similar variations with odd integer multiples of half of T constituting optimal information transmission delays where s is minimal. As for attractive coupling, values of t outside these bounds impair synchronization significantly and fast. Altogether, results presented in Figs. 2 and 3 indicate that simple scaling laws account for the description of optimal information transmission delays that warrant near-perfect synchronization of neuronal activity on scale-free networks. While for attractive coupling integer multiples of a given constant period are optimal, for repulsive coupling odd integer multiples of half of the same period have the best effect. Irrespective of the coupling type, delays outside the narrow optimal span impair synchronization significantly.
It is next of interest to explore and determine the mechanisms behind these observations. We will do this by means of the duration of bursting periods of individual neurons constituting the scale-free network. The top three panels of Fig. 4 depict time courses of the membrane potential x (i) (n) for the values of a we have used in Figs. 2 and 3 above. It can be observed that, depending on a, the duration of bursts within a given trace may vary (chaotic bursting [48]), but also that the duration of bursts changes due to different a values. This is highlighted by labels T 1 , T 2 and T 3 (where applicable) in the top three panels of Fig. 4. From this it is straightforward to determine the average oscillation period of bursting T for each particular value of a, simply as the average over a large enough ensemble L as T~L {1 P i~1:::L T i . The bottom panel of Fig. 4 shows how the average period T varies with a. It can be observed that upon exceeding the Hopf bifurcation at a~2:0 the period T increases fairly linearly, but then drops rather sharply when a exceeds 4:0.
Upon connecting the values of T with the optimal information transmission delays observed in Figs. 2 and 3 for the corresponding values of a, we can establish a good understanding of the mechanism behind the observed synchronization transitions for attractive as well as for repulsive coupling. In particular, from results presented in Fig. 4 it follows that if a~2:3 then T&850, which is exactly the value of t corresponding to the first minimum of s for attractive coupling. Conversely, one half and three times one half of T&850 correspond to the first and second minima of s if a~2:3 and the coupling is repulsive. For the other two considered values of a, namely 3:0 and 4:1, an identical linkage can be established easily from the results presented in Figs. 2 and 3, depending on the type of coupling one is interested in, and the bottom panel of Fig. 4. Apparently, the average period of individual bursts determines the optimal information transmission delay that warrants the best synchrony, i.e. minimal s, of neuronal firings on the scale-free network. We therefore conclude that for attractive coupling the delay-induced transitions to spatiotemporal synchronization of neuronal activity are due to the locking between t and the average oscillation period of individual neurons constituting the scale-free network. Importantly, because the repulsive coupling can pull adjacent neurons into antiphase synchronization, the optimal delay warranting best synchronization is not equal to full integer multiples of T. Thus, it is exactly odd integer multiples of one half of the average oscillation period of an individual neuron, where the phase locking between antiphased bursts occurs.
Merging these observation into an overall insight about delayinduced synchronization transition on scale-free networks with attractive and repulsive coupling, we show in Fig. 5 contour plots of s, which depend on the two main parameters a and t for the    two types of coupling separately. The emergence of highly synchronous tongue-like regions in the two-dimensional parameter plane agrees perfectly with the reasoning we have outlined above. As the information transmission delay increases the neuronal activity enters and exits synchronous regions in an intermittent fashion. Simultaneously, as a increases, the average period of bursting increases nearly linearly according to the results presented in the bottom panel of Fig. 4, thus giving an upward momentum to the white regions. However, when aw4:0 the average oscillation period drops sharply, which terminates the white ''tongues'' of synchrony rather abruptly and shifts the optima toward much smaller t. Altogether the presented results are in agreement with those presented in Figs. 2 and 3.
In what follows, in order to test the generality of the above results, we investigate the impact of different system sizes N and different models of neuronal dynamics, including those of type I and type II. Firstly, for different system sizes, results depicted in Figs. 6(a) and (b) show clearly that the variations of N do not notably influence the outcome of our simulations. In fact, the minima of s remain located at about the same values of t irrespective of N. In order to validate our conclusions for different types of neuronal dynamics, we choose the famous Hodgkin-Huxley model (type II) and the Morris-Lecar model (type I) to describe the dynamics of individual network nodes (both models are given in the Methods section under ''Alternative models of neuronal dynamics''). Using these two models, we investigate the synchronization transition when the delay is varied. It is shown in Figs. 7 and 8 that irrespectively of the type of the governing neuronal dynamics, intermittent synchronization transitions can still be observed for both the attractive as well as repulsive coupling when the delay is increased. More importantly, the phase locking between the delay and the period of oscillators persists in a way that is identical to what we reported above for the Rulkov model. Hence, the obtained results are also deemed robust against the variations of the neuronal dynamics.
Lastly, we construct a square lattice occupying 128|128 neurons, whose nodes are modeled by the Rulkov map. Here we set the parameter a~1:99, so that every neuron operates in the excitable regime. Starting with random initial conditions, the results in Fig. 9(a) evidence that as the delay equals t~0, there is no pattern formation observable and each neuron approaches its excitable steady state value. On the other hand, however, Fig. 9(b) features coherent waves of excitation that appear as the delay equals t~50, which emerge due to the locking between the delay length and the characteristic transient time of the local neuronal dynamics. Hence, it can be concluded that appropriate information transmission delays can also evoke ordered waves of excitation in the spatial domain, thus adding to their importance for the functioning of neuronal tissue.

Discussion
We have studied delay-induced synchronization transitions on attractively and repulsively coupled scale-free neuronal networks that were locally modeled by the Rulkov map. We have shown that, irrespective of the type of couplings, information transmission delays play a pivotal role in ensuring synchronized neuronal activity. By attractive and repulsive couplings, the synchronization of bursting oscillations was found undulating intermittently as the delay was increased. However, while for attractive coupling the regions of high synchronization appeared every integer multiple of the average oscillation period, for the repulsive coupling they appeared every odd multiple of the half of the average oscillation period. Aiming to explain these observation, we have argued that by attractive coupling the intermittent outlay of synchronized regions emerges due to the locking between the delay length and the average oscillation period of bursting oscillations of individual neurons constituting the scale-free network. Conversely, by repulsive coupling the emergence of antiphase synchronization indicates locking between the delay and odd multiples of one half of the average oscillation period. Our results indicate that information transmission delays can either promote or impair synchrony among neurons and can thus effectively supplement other mechanisms of synchronization [49,50] on scale-free networks, which arguably constitutes an important ingredient of interneuronal communication. These conclusions seem to be supported by actual biological  data, stating that conduction velocities along axons connecting neurons vary from 20 to 60 m/s [51]. Real-life transmission delays are thus within the range of milliseconds, suggesting that substantially lower or higher values may be preclusive for optimal functioning of neuronal tissue. Repulsive coupling, as we have considered it in this study, is in fact an inherent ingredient of several biological systems, in particular those that contain dynamical units that are in ''competition'' with each other. Known examples are the inhibitory couplings is present in neuronal circuits associated with a synchronized behavior in central pattern generators or calcium  oscillations in epileptic human astrocyte cultures [42]. We hope that these results will foster our understanding of the observed neuronal activity.

Methods
The map proposed by Rulkov [52,53] determines the dynamics of individual nodes forming the scale-free network. It captures succinctly the main dynamical features of the more complex timecontinuous neuronal models, but simultaneously allows an efficient numerical treatment of large systems [54]. Accordingly, the spatiotemporal evolution of the studied network with information transmission delay is governed by the following iteration equations where f (x)~1 1zx 2 is a nonlinear function warranting the essential ingredients of neuronal dynamics, x (i) (n) is the membrane potential of the i-th neuron and y (i) (n) is the variation of the ion concentration, the two representing the fast and the slow variable of the map, respectively. The slow temporal evolution of y (i) (n) is due to the small values of the two parameters b and c that are here both set equal to 0:001. Moreover, n is the discrete time index, while a is the main bifurcation parameter determining the dynamics of individual neurons constituting the scale-free network. In [52] it was shown that for av2:0 all neurons are situated in excitable steady states ½x Ã~{ 1,y Ã~{ 1{(a=2), whereas if aw2:0 complex oscillatory and bursting patterns can emerge via a Hopf bifurcation. Importantly, we set the coupling strength equal to either D~0:01, corresponding to attractive coupling, or D~{0:01, corresponding to repulsive coupling. Parameter t is the information transmission delay that together with a represents the two crucial parameters that are varied in the realm of this study.
As the interaction network between neurons we use the scale-free network generated via growth and preferential attachment as proposed by Barabási and Albert [55], typically consisting of N~200 nodes or more. Each node corresponds to one neuron, whose dynamics is governed by the Rulkov map, as described above. In Eq. (1) e i,j~1 if neuron i is coupled to neuron j and e i,j~0 otherwise. Following [55], the preferential attachment is introduced via the probability P, which states that a new node will be connected to node i depending on its connectivity k i according to P(k i )~k i = P j k j . Here, k i is the degree of node i (the degree of a node is the number of links adjacent to it). This growth and preferential attachment scheme yields a network with an average degree k av~P i ki N , and a power-law degree distribution with the slope of the line equaling &{3 on a double-logarithmic graph. We will use Barabási-Albert scale-free networks having k av~4 throughout this work.
In order to study synchronization transitions quantitatively, we introduce, by means of the standard deviation, a synchronization parameter s (see e.g. [56]), which can be calculated effectively according to: In particular, s is an excellent quantity for numerically effectively measuring the spatiotemporal synchronization of excitations, hence revealing different synchronization levels and with it related transitions. From Eq. (3) it is evident that the more synchronous the neuronal network the smaller the synchronization parameter s. Accordingly, in the event of complete synchrony we have s~0. Presented results were averaged over 20 independent runs for each set of parameter values to warrant appropriate statistical accuracy with respect to the scale-free network generation and numerical simulations.

Alternative models of neuronal dynamics
The full Hodgkin-Huxley model is given by the following equations [57]: where V is the transmembrane potential of the neuron, and m, h and n are the corresponding gating variables (probabilities) characterized by a two-state, opening or closing dynamics. The voltage-dependent opening and closing rates are given explicitly by the following expressions: The membrane capacity C = 1 (mF/cm 2 ), parameters g Na , g K and g L are maximal sodium, potassium and leakage conductances, g Na~1 20 mF/cm 2 , g K~3 6 mF/cm 2 and g L~0 :3 mF/cm 2 , respectively, V Na , V K , and V L are the reversal potentials, V Na~5 0 mV, V K~{ 77 mV, V L~{ 54:4 mV and I~10 m.
The dynamics of the type I Morris-Lecar neuron is described by the following equations [58]: where V is the cell membrane potential in mV, I Ca is the depolarizing calcium current, I L is the passive leak current, respectively, v is the activation of the repolarizing potassium current I K , t is time in ms, and Iapp~14 m A/cm 2 is the applied current. The remaining parameters are V Ca = 120 mV, V K~{ 84 mV, V L~{ 60 mV, g Ca = 4 mS/cm 2 , g K = 8 mS/ cm 2 , g L = 2 mS/cm 2 . The steady state activation of the calcium current is: The potassium current activation amplitude and activation rate are: