# Frequency-Domain Order Parameters for the Burst And Spike Synchronization Transitions of Bursting Neurons

## Abstract

We are interested in characterization of synchronization transitions of bursting neurons in the frequency domain. Instantaneous population firing rate (IPFR) , which is directly obtained from the raster plot of neural spikes, is often used as a realistic collective quantity describing population activities in both the computational and the experimental neuroscience. For the case of spiking neurons, a realistic time-domain order parameter, based on , was introduced in our recent work to characterize the spike synchronization transition. Unlike the case of spiking neurons, the IPFR of bursting neurons exhibits population behaviors with both the slow bursting and the fast spiking timescales. For our aim, we decompose the IPFR into the instantaneous population bursting rate (describing the bursting behavior) and the instantaneous population spike rate (describing the spiking behavior) via frequency filtering, and extend the realistic order parameter to the case of bursting neurons. Thus, we develop the frequency-domain bursting and spiking order parameters which are just the bursting and spiking “coherence factors” and of the bursting and spiking peaks in the power spectral densities of and (i.e., “signal to noise” ratio of the spectral peak height and its relative width). Through calculation of and , we obtain the bursting and spiking thresholds beyond which the burst and spike synchronizations break up, respectively. Consequently, it is shown in explicit examples that the frequency-domain bursting and spiking order parameters may be usefully used for characterization of the bursting and the spiking transitions, respectively.

###### pacs:

87.19.lm, 87.19.lc## I Introduction

Recently, much attention has been paid to brain rhythms, observed in electrical recordings of firing activity (1). These brain rhythms emerge via synchronization between firings of individual neurons. This kind of neural synchronization may be used for efficient sensory and cognitive processing (2); (3), and it is also correlated with pathological rhythms associated with neural diseases (4); (5); (6). Here, we are interested in characterization of population synchronization of bursting neurons in terms of neural synchrony measures (7). Bursting occurs when neuronal activity alternates, on a slow timescale, between a silent phase and an active (bursting) phase of fast repetitive spikings (8); (9); (10); (11); (12). Due to the slow and fast timescales of bursting activity, bursting neurons exhibit two types of burst and spike synchronizations. Burst synchronization on the slow bursting timescale refers to a coherence between the active phase onset or offset times of bursting neurons, while spike synchronization on the fast spike timescale characterizes a coherence between intraburst spikes fired by bursting neurons (13); (14). Many recent studies on the burst and spike synchronizations have been made in several aspects (e.g., chaotic phase synchronization, transitions between different states of burst synchronization, effect of network topology, effect on information transmission, suppression of bursting synchronization, and effect of noise and coupling on the burst and spike synchronization) (15); (16); (17); (18); (19); (20); (21); (22); (23); (24); (25); (26); (27); (28); (29).

In this paper, we are interested in practical characterization of the burst and spike synchronization transitions of bursting neurons in the frequency domain. Population synchronization may be well visualized in the raster plot of neural spikes which can be obtained in experiments. Instantaneous population firing rate (IPFR), , which is directly obtained from the raster plot of spikes, is a realistic population quantity describing collective behaviors in both the computational and the experimental neuroscience (2); (30); (31); (32); (33); (34); (35). This experimentally-obtainable is in contrast to the ensemble-averaged potential which is often used as a population quantity in the computational neuroscience, because to directly get in real experiments is very difficult. To overcome this difficulty, instead of , we used as a population quantity, and developed a realistic order parameter, based on , to make practical characterization of synchronization of spiking neurons in both the computational and the experimental neuroscience (36). The mean square deviation of plays the role of the realistic order parameter used to determine the threshold value for the synchronization transition of spiking neurons. In this way, synchronization transition of spiking neurons may be well characterized in terms of the realistic order parameter , based on the IPFR .

In the field of neuroscience, power-spectral analysis of a time-series (e.g., neuronal membrane potential) is often made to examine how the variance of the data is distributed over the frequency components into which may be decomposed (30); (31); (32); (33); (34); (35). Following this conventional direction, we investigate the synchronization transitions of bursting neurons in the frequency domain. The main purpose of our works is to characterize the synchronization transitions of bursting neurons in terms of “frequency-domain” order parameters by extending the realistic “time-domain” order parameter of spiking neurons (36) to the case of bursting neurons. This extension work on the frequency-domain order parameters is in contrast to another extension work where the synchronization transitions of bursting neurons are characterized in terms of the time-domain order parameters (37). The IPFR shows the whole combined population behaviors with both the slow and fast timescales. To clearly investigate the synchronization transitions of bursting neurons, we separate the slow and fast timescales of the bursting activity via frequency filtering, and decompose the IPFR into (the instantaneous population burst rate (IPBR) describing the bursting behavior) and (the instantaneous population spike rate (IPSR) describing the intraburst spiking behavior). In presence of the burst and spike synchronizations, and exhibit regular oscillations, independently of (the number of the bursting neurons). On the other hand, in absence of the burst and spike synchronizations, and become stationary as goes to the infinity. The synchronous oscillations of and in the time domain may be well characterized by the bursting and spiking peaks of their power spectral densities. As in the case of the coherence resonance (38), each spectral “resonance” (i.e., peak) in the power spectral density may be well analyzed in terms of a “coherence factor” (i.e., a measure of spectral coherence) which is defined by a “signal to noise” ratio of the spectral peak height and its relative width (39); (40). We also note that the signal to noise ratio has a long history of being used in neuroscience as a measure of the fidelity of signal transmission and detection by neurons and synapses (41). Then, the bursting and spiking coherence factors and of the bursting and spiking peaks in the power spectral densities of and are shown to play the role of the bursting and spiking order parameters in the frequency domain which are used to determine the bursting and spiking thresholds for the bursting and spiking transitions, respectively. We also consider another raster plot of bursting onset or offset times which visualizes the bursting behaviors more directly. From this type of raster plot, we may directly obtain the IPBR, or , without frequency filtering. Then, the bursting onset and offset coherence factors, and , of the bursting onset and offset peaks in the power spectral densities of and are also shown to play the role of the frequency-domain bursting order parameters for the bursting transition. The frequency-domain order parameters and yield the same bursting threshold which is obtained through calculation of , and they are more direct ones than because they may be directly obtained without frequency filtering. Consequently, all the frequency-domain bursting and spiking order parameters may be usefully used for characterization of the burst and and spike synchronization transitions of the bursting neurons in the frequency domain.

This paper is organized as follows. In Sec. II, as an example for characterization we describe an inhibitory network of bursting Hindmarsh-Rose (HR) neurons (42); (43); (44); (45). In Sec. III, we separate the slow bursting and the fast spiking timescales via frequency filtering, and develop realistic frequency-domain bursting and spiking order parameters (i.e., the bursting and spiking coherence factors), based on the power spectral densities of the IPBR and the IPSR, which are applicable in both the computational and the experimental neuroscience. Their usefulness for characterization of the burst and spike synchronization transitions is shown in explicit examples of bursting HR neurons. Finally, a summary is given in Section IV.

## Ii A Network of Inhibitory Bursting Hindmarsh-Rose Neurons

As an example for characterization, we consider an inhibitory network of globally-coupled bursting HR neurons. The representative bursting HR neuron model was originally introduced to describe the time evolution of the membrane potential for the pond snails (42); (43); (44); (45). The population dynamics in this inhibitory network is governed by the following set of ordinary differential equations:

(1) | |||||

(2) | |||||

(3) | |||||

(4) |

where

(5) | |||||

(6) |

Here, the state of the th HR neuron at a time (measured in units of milliseconds) is described by four state variables: the fast membrane potential , the fast recovery current the slow adaptation current , and the synaptic gate variable representing the fraction of open synaptic ion channels. The parameters in the single HR neuron are taken as and (40).

Each bursting HR neuron is stimulated by using the common DC current and an independent Gaussian white noise [see the 5th and the 6th terms in Eq. (1)] satisfying and , where denotes the ensemble average. The noise is a parametric one that randomly perturbs the strength of the applied current , and its intensity is controlled by using the parameter . As passes a threshold in the absence of noise, each single HR neuron exhibits a transition from a resting state to a bursting state. Throughout this paper, we consider the suprathreshold case of where each HR neuron exhibits spontaneous bursting activity without noise. Figures 1(a)-1(b) show the time series of the fast membrane potential and the fast recovery current , while the time series of the slow adaptation current is shown in Fig. 1(c). As seen well in the time series of and , the bursting activity alternates, on a slow timescale, between a silent phase and an active (bursting) phase of fast repetitive spikings. For this case, the slow bursting timescale is ms [corresponding to the slow bursting frequency Hz)], while the fast spiking timescale is ms [corresponding to the fast spiking frequency Hz)].

The last term in Eq. (1) represents the synaptic coupling of the network. of Eq. (5) represents a synaptic current injected into the th neuron. Here the coupling strength is controlled by the parameter and is the synaptic reversal potential. Here, we use for the inhibitory synapse. The synaptic gate variable obeys the 1st order kinetics of Eq. (4) (46); (47). Here, the normalized concentration of synaptic transmitters, activating the synapse, is assumed to be an instantaneous sigmoidal function of the membrane potential with a threshold in Eq. (6), where we set and (48). The transmitter release occurs only when the neuron emits a spike (i.e., its potential is larger than ). For the inhibitory GABAergic synapse (involving the receptors), the synaptic channel opening rate, corresponding to the inverse of the synaptic rise time , is , and the synaptic closing rate , which is the inverse of the synaptic decay time , is ((49); (50)). Hence, rises fast and decays slowly.

## Iii Frequency-domain order parameters for the burst and spike synchronization transitions

In this section, we extend the realistic order parameter of spiking neurons to the case of bursting neurons for characterization of population synchronization transition in the frequency domain. For our aim, we separate the slow bursting and the fast spiking timescales through frequency filtering, and decompose the IPFR into the IPBR (describing the bursting behavior) and the IPSR (describing the intraburst spiking behavior). Then, we develop realistic frequency-domain bursting and spiking order parameters, based on the power spectral densities of the IPBR and the IPSR , and show their usefulness for characterization of the burst and spike synchronization transitions in explicit examples of bursting HR neurons.

As an example for characterization, we consider an inhibitory network of globally-coupled bursting HR neurons, and characterize the synchronization transitions of bursting HR neurons in the frequency domain by varying the noise intensity . To compare our results in the frequency domain with those in the time domain, we fix the DC current strength and the coupling strength at and , as in the time-domain work (37). In computational neuroscience, a population-averaged global potential,

(7) |

is often used for describing emergence of population synchronization. In this study, we consider the population behaviors after the transient time of ms. Although the global potential is an important ensemble-averaged quantity to describe synchronization in computational neuroscience, it is practically difficult to directly get in real experiments. To overcome this difficulty, instead of , we use the IPFR which is an experimentally-obtainable population quantity used in both the experimental and the computational neuroscience (2); (30); (31); (32); (33); (34); (35). The IPFR is obtained from the raster plot of spikes which is a collection of spike trains of individual neurons. Such raster plots of spikes, where population synchronization may be well visualized, are fundamental data in the experimental neuroscience. The raster plots of spikes in Figs. 2(a1)-2(a5) show population states for various values of noise intensity . To get a smooth IPFR from the raster plot of spikes, we employ the kernel density estimation (kernel smoother) (52). Each spike in the raster plot is convoluted (or blurred) with a kernel function to get a smooth estimate of IPFR, :

(8) |

where is the th spiking time of the th neuron, is the total number of spikes for the th neuron, and we use a Gaussian kernel function of band width :

(9) |

Figures 2(b1)-2(b5) show smooth IPFR kernel estimates of band width ms for , 0.01, 0.04, 0.06 and 0.08, respectively. For , clear “bursting bands,” each of which is composed of “stripes” of spikes, appear successively at nearly regular time intervals [see Fig. 2(a1)]; a magnified 1st intraburst band is given in Fig. 3(a1). For the case of , both the burst synchronization [synchrony on the slow bursting timescale ( ms)] and the spike synchronization [synchrony on the fast spike timescale ms)] occur in each bursting band. As a result of this complete synchronization, the IPFR kernel estimate shows a bursting activity [i.e., fast spikings appear on a slow wave in ], as shown in Fig. 2(b1). However, as is increased, loss of spike synchronization occurs in each bursting band because spiking stripes become smeared due to a destructive role of noise. As an example, see the case of where the raster plot of spikes and the IPFR kernel estimate are shown in Figs. 2(a2) and 2(b2), respectively. The magnified 1st bursting band in Fig. 3(a3) shows smearing of the spiking stripes well. Consequently, the amplitude of decreases, as shown in Fig. 2(b2). As is further increased and passes a spiking noise threshold , complete loss of spike synchronization occurs in each bursting band. Then, only the burst synchronization (without spike synchronization) occurs, as shown in the case of in Figs. 2(a3) and 2(b3). For this case, shows a slow-wave oscillation without spikes. With increase in , such “incoherent” bursting bands become more and more smeared, and hence the degree of burst synchronization decreases [e.g., see the case of in Fig. 2(a4)]. As a result, the amplitude of is further decreased, as shown in Fig. 2(b4) for . With further increasing , incoherent bursting bands begin to overlap, which eventually results in the complete loss of burst synchronization as passes another larger bursting noise threshold . Consequently, for , completely unsynchronized states with nearly stationary appear, as shown in the case of in Figs. 2(a5) and 2(b5).

The (above) IPFR kernel estimate is a population quantity describing the “whole” combined collective behaviors of bursting neurons with both the slow bursting and the fast spiking timescales. Through frequency filtering, we separate the slow and the fast timescales, and decompose the IPFR kernel estimate into the IPBR and the IPSR for more clear investigation of the burst and spike synchronizations. Through band-pass filtering of [with the lower and the higher cut-off frequencies of 3 Hz (high-pass filter) and 7 Hz (low-pass filer)], we get the regularly-oscillating IPBR (containing only the slow wave without spikes) in Figs. 2(c1)-2(c5) for , 0.01, 0.04, 0.06, and 0.08. As is increased, the amplitude of decreases gradually, and eventually becomes nearly stationary when passes the bursting noise threshold . We note that synchronous oscillations of in the time domain are characterized by the bursting peaks in the power spectral densities of , where the overline represents the time average. Figures 2(d1)-2(d5) show distinct bursting peaks in the power spectra of ; each power spectrum is made of data points and smoothed through the Daniell filters of length 3 and 5 (53). Then, each bursting peak may be analyzed well in terms of a bursting coherence factor defined by the product of the height and the factor of the peak (38); (39); (40):

(10) |

Here, and are the frequency of the bursting peak and the width of the bursting peak at the height of , respectively. For more accurate results, we repeat the process to get the bursting coherence factor for multiple realizations. Thus, we obtain (average bursting coherence factor) through average over 20 realizations. Figure 2(e) shows plots of the average bursting coherence factor versus . For ), synchronized bursting states exist because the values of become saturated to non-zero limit values in the thermodynamic limit of (i.e., bursting peaks persist, independently of ). However, as passes the bursting noise threshold , the average bursting coherence factor tends to zero as (i.e., eventually bursting peaks disappear in the thermodynamic limit), and hence a transition to unsynchronized bursting states occurs because the noise spoils the burst synchronization completely. In this way, the average bursting coherence factor describes the burst synchronization transition well in the frequency domain, and hence it plays the role of the realistic frequency-domain bursting order parameter for the bursting transition (i.e., one can determine the bursting noise threshold through calculation of ). This frequency-domain bursting order parameter is in contrast to the time-domain bursting order parameter, based on the time-averaged fluctuation of the IPBR (37). In spite of their difference, calculations of both the frequency-domain and the time-domain bursting order parameters result in the same bursting noise threshold (compare Fig. 2(e) with Fig. 3(a) in (37)). Consequently, the frequency-domain bursting order parameter may be used effectively to determine for the bursting transition, like the case of the time-domain bursting order parameter.

From now on, we investigate the intraburst spike synchronization transition of bursting HR neurons in the frequency domain by varying the noise intensity . Figures 3(a1)-3(a5) and Figures 3(b1)-3(b5) show the raster plots of intraburst spikes and the corresponding (band-pass filtered) IPSR during the 1st global bursting cycle of the IPBR , respectively for various values of : synchronized spiking states for , 0.005, 0.01, and 0.02, and unsynchronized spiking state for . Here, the IPSRs are obtained through band-pass filtering of the IPFR kernel estimate [with the lower and the higher cut-off frequencies of 30 Hz (high-pass filter) and 90 Hz (low-pass filer)]. Then, the intraburst spike synchronization may be well described in terms of . For , clear 8 spiking stripes (composed of spikes and indicating population spike synchronization) appear in the intraburst band of the 1st global bursting cycle of in Fig. 3(a1), and the band-pass filtered IPSR shows only the fast spiking oscillations (without a slow wave) with the population spiking frequency Hz) in Fig. 3(b1). However, as is increased, spiking stripes in the intraburst band become more and more smeared (e.g., see the cases of , 0.01, and 0.02). Consequently, the amplitude of decreases due to loss of spike synchronization. Eventually, when passes the spiking noise threshold , spikes become completely scattered within the intraburst band, and becomes nearly stationary. Consequently, for , complete loss of spike synchronization occurs in the intraburst band, as shown in Fig. 3(b5) for . Figures 3(c1)-3(c5) show the power spectra of in the 1st global bursting cycle of : each power spectrum is made of data points and smoothed through the Daniell filters of length 3 and 5. Spiking peaks in their power spectra are analyzed in terms of the spiking coherence factors (defined by the product of the height and the factor of the peak). For more accurate results, we repeat the process to get for multiple realizations. In each realization we follow the 20 global bursting cycles of , and get the double-averaged spiking coherence factor through average over 20 realizations. Figure 3(d) shows plots of the double-averaged spiking coherence factor versus . For ), synchronized spiking states exist because the values of become saturated to non-zero limit values as (i.e., spiking peaks persist, irrespectively of ). However, when passes the spiking noise threshold , tends to zero in the thermodynamic limit of (i.e., eventually spiking peaks disappear in the thermodynamic limit), and hence a transition to unsynchronized spiking states occurs because the noise spoils the intraburst spike synchronization completely. In this way, the double-averaged spiking coherence factor describes the intraburst spike synchronization transition well in the frequency domain, and hence it plays the role of the realistic frequency-domain spiking order parameter for the spiking transition (i.e., one can determine the spiking noise threshold through calculation of ). This frequency-domain spiking order parameter is also in contrast to the time-domain spiking order parameter, based on the time-averaged fluctuation of the IPSR (37). We also note that both the frequency-domain and the time-domain spiking order parameters yield the same spiking noise threshold (compare Fig. 3(d) with Fig. 6(d) in (37)). Consequently, the frequency-domain spiking order parameter may also be used effectively to determine for the spiking transition, as in the case of the time-domain spiking order parameter.

Finally, we consider another raster plot of bursting onset or offset times for more direct visualization of bursting behavior. [At the onset (offset) times of the th bursting HR neuron, its individual potential passes the threshold of from below (above).] Without frequency filtering, we can directly obtain the IPBR kernel estimate, [] from the raster plot of the bursting onset (offset) times. Figures 4(a1)-4(a5) show the raster plots of the bursting onset times for various values of , while the raster plots of the bursting offset times are shown in Figs. 4(c1)-4(c5). From these raster plots of the bursting onset (offset) times, we obtain smooth IPBR kernel estimates, [], of band width ms in Figs. 4(b1)[(d1)]-4(b5)[(d5)] for , 0.01, 0.04, 0.06, and 0.08. For , clear bursting “stripes” [composed of bursting onset (offset) times and indicating burst synchronization] appear successively at nearly regular time intervals; the bursting onset and offset stripes are time-shifted [see Figs. 4(a1) and 4(c1)]. The corresponding IPBR kernel estimates, and , for show regular oscillations with the same population bursting frequency Hz), although they are phase-shifted [see Figs. 4(b1) and 4(d1)]. With increasing , the bursting onset and offset stripes in the raster plots become smeared and begin to overlap, and thus the degree of the burst synchronization decreases. As a result, the amplitudes of both and decrease gradually (e.g., see the cases of , 0.04, and 0.06). Eventually, as passes the bursting noise threshold , bursting onset and offset times become completely scattered in the raster plots, and the corresponding IPBR kernel estimates, and , become nearly stationary, as shown in Figs. 4(b5) and 4(d5) for . Figures 4(e1)-4(e5) show the power spectra of , while Figures 4(f1)-4(f5) show the power spectra of ; each power spectrum is made of data points and smoothed through the Daniell filters of length 3 and 5. Bursting onset and offset peaks in these power spectra are analyzed in terms of the bursting onset and offset coherence factors and (each coherence factor of a peak is defined by the product of the height and the factor of the peak). For more accurate results, we repeat the process to obtain and for multiple realizations. Thus, we obtain (average bursting onset coherence factor) and (average bursting offset coherence factor) through average over 20 realizations. Figures 4(g1) and 4(g2) show plots of the average bursting onset and offset coherence factor and versus , respectively. As in the case of the (above) average bursting coherence factor , when passing the same bursting noise threshold , both the average bursting onset and offset coherence factors and go to zero as (i.e., eventually bursting onset and offset peaks disappear in the thermodynamic limit), and hence a transition to burst unsynchronization occurs for , because the noise breaks up the burst synchronization completely. In this way, both the average bursting onset and offset coherence factors, and , describe the burst synchronization transition well in the frequency domain, and hence they also play the role of the realistic frequency-domain bursting order parameters for the bursting transition together with . These frequency-domain bursting order parameters are also in contrast to the time-domain bursting order parameters, based on the time-averaged fluctuations of and (37). We note that both the frequency-domain and the time-domain bursting order parameters yield the same bursting noise threshold (compare Figs. 4(g1) and 4(g2) with Figs. 3(b) and 3(c) in (37)). Consequently, along with , the frequency-domain bursting order parameters, and , may also be used effectively to determine for the bursting transition, as in the case of the time-domain bursting order parameters.

## Iv Summary

We have extended the realistic time-domain order parameter of spiking neurons to the case of bursting neurons. Their usefulness for characterization of the burst and spike synchronization transitions in the frequency domain has been shown in explicit examples of bursting HR neurons by varying the noise intensity . Population synchronization may be well visualized in the raster plot of neural spikes which may be obtained in experiments. The IPFR kernel estimate , which is obtained from the raster plot of spikes, is a realistic collective quantity describing the whole combined population behaviors with the slow bursting and the fast spiking timescales. Through frequency filtering, we have decomposed the IPFR kernel estimate into the IPBR and the IPSR . We note that both and may be used to effectively characterize the burst and spike synchronizations, respectively. For synchronous cases, oscillations of and in the time domain are characterized by the bursting and spiking peaks in their power spectral densities. Similar to the case of coherence resonance, each spectral resonance (i.e., peak) may be well analyzed in terms of the coherence factor , defined by a “signal to noise” ratio of the spectral peak height and its relative width. The average bursting and spiking coherence factors and of the bursting and spiking peaks in the power spectral densities of and have been found to play the role of the frequency-domain bursting and spiking order parameters for the burst and spike synchronization transitions, respectively. Through calculation of and , we have determined the noise bursting and spiking thresholds, and , beyond which the burst and spike synchronizations break up, respectively. For more direct visualization of bursting behavior, we consider another raster plot of bursting onset or offset times, from which the IPBR, or , can be directly obtained without frequency filtering. Then, the average bursting onset and offset coherence factors, and of the bursting onset and offset peaks in the power spectral densities of and have also been shown to play the role of the frequency-domain bursting order parameters for the bursting transition. These frequency-domain order parameters and yield the same bursting noise threshold which is obtained via calculation of , and they are more direct ones than because they may be directly obtained without frequency filtering. We also note that all these bursting and spiking noise thresholds are the same as those obtained through calculations of the time-domain bursting and spiking order parameters (37). Consequently, the frequency-domain bursting and spiking order parameters may be usefully used for characterizing the burst and spike synchronization transitions of the bursting neurons, as in the case of the time-domain bursting and spiking order parameters.

###### Acknowledgements.

This research was supported by Basic Science Research Program through the National Research Foundation of Korea (NRF) funded by the Ministry of Education (Grant No. 2013057789).### References

- G. Buzski, Rhythms of the Brain (Oxford University Press, New York, 2006).
- X.-J. Wang, Physiol. Rev. 90, 1195 (2010).
- X.-J. Wang, in Encyclopedia of Cognitive Science, edited by L. Nadel (MacMillan, London, 2003), pp. 272-280.
- P. J. Uhlhaas and W. Singer, Neuron 52, 155 (2006).
- R. D. Traub and M. A. Whittington, Cortical Oscillations in Health and Diseases (Oxford University Press, New York, 2010).
- T. J. Kaper, M. A. Kramer, and H. G. Rotstein, Chaos 23, 046001 (2013).
- D. Golomb, Scholarpedia 2(1), 1347 (2007).
- J. Rinzel, in Ordinary and Partial Differential Equations, edited by B.D. Sleeman and R.J. Jarvis, Lecture Notes in Mathematics Vol. 1151 (Springer, Berlin, 1985), pp. 304-316.
- J. Rinzel, in Mathematical Topics in Population Biology, Morphogenesis, and Neurosciences, edited by E. Teramoto and M. Yamaguti, Lecture Notes in Biomathematics Vol. 71 (Springer, Berlin, 1987), pp. 267-281.
- Bursting: The Genesis of Rhythm in the Nervous System, edited by S. Coombes and P. C. Bressloff (World Scientific, Singapore, 2005).
- E. M. Izhikevich, Scholarpedia 1(3), 1300 (2006).
- E. M. Izhikevich, Dynamical Systems in Neuroscience (MIT Press, Cambridge, 2007).
- J. E. Rubin, Scholarpedia 2(10), 1666 (2007).
- I. Omelchenko, M. Rosenblum, and A. Pikovsky, Eur. Phys. J. 191, 3 (2010).
- X. Sun, J. Lei, M. Perc, J. Kurths, and G. Chen, Chaos 21, 016110 (2011).
- C. van Vreeswijk and D. Hansel, Neural Comput. 13, 959 (2001).
- M. Dhamala, V. Jirsa, and M. Ding, Phys. Rev. Lett. 92, 028101 (2004).
- M.V. Ivanchenko, G. Osipov, V. Shalfeev, and J. Kurths, Phys. Rev. Lett. 93, 134101 (2004).
- T. Pereira, M. Baptista, and J. Kurths, Eur. Phys. J. Spec. Top. 146, 155 (2007).
- H. Yu, J. Wang, B. Deng, X. Wei, Y.K. Wong, W.L. Chan, K.M. Tsang, and Z. Yu, Chaos 21, 013127 (2011).
- G. Tanaka, B. Ibarz, M.A. Sanjuan, and K. Aihara, Chaos 16, 013113 (2006).
- X. Shi and Q. Lu, Physica A 388, 2410 (2009).
- X. Shi and Q. Lu, Chin. Phys. 14, 77 (2005).
- C.A.S. Batista, A.M. Batista, J.A.C. de Pontes, R.L. Viana, and S.R. Lopes, Phys. Rev. E 76, 016218 (2007).
- C. A. Batista, E.L. Lameu, A.M. Batista, S.R. Lopes, T. Pereira, G.Zamora-Lopez, J. Kurths, and R.L. Viana, Phys. Rev. E 86, 016211 (2012).
- E. L. Lameu, C. A. S. Batista, A. M. Batista, K. Larosz, R. L. Viana, S. R. Lopes, and J. Kurths, Chaos 22, 043149 (2012).
- S.-Y. Kim, Y. Kim, D.-G. Hong, J. Kim, and W. Lim, J. Korean Phys. Soc. 60, 1441 (2012).
- S.-Y. Kim and W. Lim, Cognitive Neurodynamics 7, 495 (2013).
- I. Belykh, E. de Lange, and M. Hasler, Phys. Rev. Lett. 94 188101 (2005).
- N. Brunel and V. Hakim, Chaos 18, 015113 (2008).
- N. Brunel and V. Hakim, Neural Comput. 11, 1621 (1999).
- N. Brunel, J. Comput. Neurosci. 8, 183 (2000).
- N. Brunel and X.-J. Wang, J. Neurophysiol. 90, 415 (2003).
- C. Geisler, N. Brunel, and X.-J. Wang, J. Neurophysiol. 94, 4344 (2005).
- N. Brunel and D. Hansel, Neural Comp. 18, 1066 (2006).
- S.-Y. Kim and W. Lim, J. Neurosci. Methods 226, 161 (2014).
- S.-Y. Kim and W. Lim, e-print: arXiv:1403.3994 [q-bio.NC].
- A. Neiman, Scholarpedia 2(11), 1442 (2007).
- H. Gang, T. Ditzinger, C.Z. Ning, and H. Haken, Phys. Rev. Lett. 71, 807 (1993).
- A. Longtin, Phys. Rev. E 55, 868 (1997).
- S. R. Schultz, Scholarpedia 2(6), 2046 (2007).
- J.L. Hindmarsh and R.M. Rose, Nature 296, 162 (1982).
- J.L. Hindmarsh and R.M. Rose, Proc. R. Soc. London, Ser. B 221, 87 (1984).
- R. M. Rose and J.L. Hindmarsh, Proc. R. Soc. London, Ser. B 225, 161 (1985).
- A. Shilnikov and M. Kolomiets, Int. J. Bifur. Chaos 18, 2141 (2008).
- D. Golomb and J. Rinzel, Physica D 72, 259 (1994).
- X.-J. Wang and G. Buzski, J. Neurosci. 16, 6402 (1996).
- X. Liang, M. Tang, M. Dhamala and Z. Liu, Phys. Rev. E 80, 066202 (2009).
- C. Brgers and N. Kopell, Neural Comput. 15, 509 (2003).
- C. Brgers and N. Kopell, Neural Comput. 17, 557 (2005).
- M. San Miguel and R. Toral, in Instabilities and Nonequilibrium Structures VI, edited by J. Martinez, R. Tiemann, and E. Tirapegui (Kluwer Academic Publisher, Dordrecht, 2000), p. 35.
- H. Shimazaki and S. Shinomoto, J. Comput. Neurosci. 29, 171 (2010).
- P. Bloomfield, Fourier analysis of time series: an introduction, 2nd edition (New York, John Wiley Sons Inc., 2000), p. 261.